TDPMjet/AliGenDPMjet.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 // Generator using DPMJET as an external generator
  18 // The main DPMJET options are accessable for the user through this interface.
  19 // Uses the TDPMjet implementation of TGenerator.
  20
  21 #include <TDPMjet.h>
  22 #include <TRandom.h>
  23 #include <TArrayI.h>
  24 #include <TParticle.h>
  25 #include <TGraph.h>
  26 #include <TDatabasePDG.h>
  27 #include <TParticlePDG.h>
  28 #include <TParticleClassPDG.h>
  29 #include <TPDGCode.h>
  30 #include <TLorentzVector.h>
  31 #include <TClonesArray.h>
  32 #include "AliRunLoader.h"
  33 #include "AliGenDPMjet.h"
  34 #include "AliGenDPMjetEventHeader.h"
  35 #include "AliRun.h"
  36 #include "AliDpmJetRndm.h"
  37 #include "AliHeader.h"
  38 #include "AliStack.h"
  39 #include "AliMC.h"
  40
  41 ClassImp(AliGenDPMjet)
  42
  43 //______________________________________________________________________________
  44 AliGenDPMjet::AliGenDPMjet()
  45     :AliGenMC(),
  46      fBeamEn(2750.),
  47      fMinImpactParam(0.),
  48      fMaxImpactParam(5.),
  49      fICentr(0),
  50      fSelectAll(0),
  51      fFlavor(0),
  52      fTrials(0),
  53      fSpectators(1),
  54      fSpecn(0),
  55      fSpecp(0),
  56      fDPMjet(0),
  57      fNoGammas(0),
  58      fLHC(0),
  59      fPi0Decay(1),
  60      fDecayAll(0),
  61      fGenImpPar(0.),
  62      fProcess(kDpmMb)
  63 {
  64 // Constructor
  65     fEnergyCMS = 5500.;
  66     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
  67 }
  68
  69
  70 //______________________________________________________________________________
  71 AliGenDPMjet::AliGenDPMjet(Int_t npart)
  72     :AliGenMC(npart),
  73      fBeamEn(2750.),
  74      fMinImpactParam(0.),
  75      fMaxImpactParam(5.),
  76      fICentr(0),
  77      fSelectAll(0),
  78      fFlavor(0),
  79      fTrials(0),
  80      fSpectators(1),
  81      fSpecn(0),
  82      fSpecp(0),
  83      fDPMjet(0),
  84      fNoGammas(0),
  85      fLHC(0),
  86      fPi0Decay(1),
  87      fDecayAll(0),
  88      fGenImpPar(0.),
  89      fProcess(kDpmMb)
  90 {
  91 // Default PbPb collisions at 5. 5 TeV
  92 //
  93     fEnergyCMS = 5500.;
  94     fName = "DPMJET";
  95     fTitle= "Particle Generator using DPMJET";
  96     SetTarget();
  97     SetProjectile();
  98     fVertex.Set(3);
  99     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 100 }
 101
 102 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
 103     :AliGenMC(),
 104      fBeamEn(2750.),
 105      fMinImpactParam(0.),
 106      fMaxImpactParam(5.),
 107      fICentr(0),
 108      fSelectAll(0),
 109      fFlavor(0),
 110      fTrials(0),
 111      fSpectators(1),
 112      fSpecn(0),
 113      fSpecp(0),
 114      fDPMjet(0),
 115      fNoGammas(0),
 116      fLHC(0),
 117      fPi0Decay(1),
 118      fDecayAll(0),
 119      fGenImpPar(0.),
 120      fProcess(kDpmMb)
 121 {
 122     // Dummy copy constructor
 123     fEnergyCMS = 5500.;
 124 }
 125
 126 //______________________________________________________________________________
 127 AliGenDPMjet::~AliGenDPMjet()
 128 {
 129 // Destructor
 130 }
 131 //______________________________________________________________________________
 132 void AliGenDPMjet::Init()
 133 {
 134 // Initialization
 135
 136     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
 137                       fBeamEn,fEnergyCMS));
 138
 139     fDPMjet=(TDPMjet*) fMCEvGen;
 140     //
 141     // **** Flag to force central production
 142     // fICentr=1. central production forced
 143     // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
 144     // fICentr<-99 -> fraction of x-sec. = XSFRAC
 145     // fICentr=-1. -> evaporation/fzc suppressed
 146     // fICentr<-1. -> evaporation/fzc suppressed
 147     if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
 148
 149     fDPMjet->SetfFCentr(fICentr);
 150     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
 151     fDPMjet->SetPi0Decay(fPi0Decay);
 152     fDPMjet->SetDecayAll(fDecayAll);
 153 //
 154 //  Initialize DPMjet
 155 //
 156     fDPMjet->Initialize();
 157 }
 158
 159
 160 //______________________________________________________________________________
 161 void AliGenDPMjet::Generate()
 162 {
 163 // Generate one event
 164
 165   Float_t polar[3]    =   {0,0,0};
 166   Float_t origin[3]   =   {0,0,0};
 167   Float_t p[3];
 168   Float_t tof;
 169
 170 //  converts from mm/c to s
 171   const Float_t kconv = 0.001/2.999792458e8;
 172   Int_t nt  = 0;
 173   Int_t jev = 0;
 174   Int_t kf, ks, imo;
 175   kf = 0;
 176   fTrials = 0;
 177   //  Set collision vertex position
 178   if (fVertexSmear == kPerEvent) Vertex();
 179
 180   while(1)
 181   {
 182 //    Generate one event
 183 // --------------------------------------------------------------------------
 184       fSpecn = 0;
 185       fSpecp = 0;
 186 // --------------------------------------------------------------------------
 187       fDPMjet->GenerateEvent();
 188
 189       fTrials++;
 190
 191       fDPMjet->ImportParticles(&fParticles,"All");
 192       if (fLHC) Boost();
 193
 194       // Temporaneo
 195       fGenImpPar = fDPMjet->GetBImpac();
 196
 197       Int_t np = fParticles.GetEntriesFast();
 198       printf("\n **************************************************%d\n",np);
 199       Int_t nc=0;
 200       if (np==0) continue;
 201       Int_t i;
 202       Int_t* newPos     = new Int_t[np];
 203       Int_t* pSelected  = new Int_t[np];
 204
 205       for (i = 0; i<np; i++) {
 206           newPos[i]    = i;
 207           pSelected[i] = 0;
 208       }
 209
 210 //      First select parent particles
 211
 212       for (i = 0; i<np; i++) {
 213           TParticle *iparticle = (TParticle *) fParticles.At(i);
 214
 215 // Is this a parent particle ?
 216
 217           if (Stable(iparticle)) continue;
 218
 219           Bool_t  selected             =  kTRUE;
 220           Bool_t  hasSelectedDaughters =  kFALSE;
 221
 222           kf = iparticle->GetPdgCode();
 223           if (kf == 92 || kf == 99999) continue;
 224           ks = iparticle->GetStatusCode();
 225 // No initial state partons
 226           if (ks==21) continue;
 227
 228           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
 229                                SelectFlavor(kf);
 230
 231
 232           hasSelectedDaughters = DaughtersSelection(iparticle);
 233
 234
 235 // Put particle on the stack if it is either selected or
 236 // it is the mother of at least one seleted particle
 237
 238           if (selected || hasSelectedDaughters) {
 239               nc++;
 240               pSelected[i] = 1;
 241           } // selected
 242       } // particle loop parents
 243
 244 // Now select the final state particles
 245
 246
 247       for (i=0; i<np; i++) {
 248           TParticle *iparticle = (TParticle *) fParticles.At(i);
 249
 250 // Is this a final state particle ?
 251
 252           if (!Stable(iparticle)) continue;
 253
 254           Bool_t  selected =  kTRUE;
 255           kf = iparticle->GetPdgCode();
 256           ks = iparticle->GetStatusCode();
 257
 258 // --------------------------------------------------------------------------
 259 // Count spectator neutrons and protons (ks == 13, 14)
 260           if(ks == 13 || ks == 14){
 261               if(kf == kNeutron) fSpecn += 1;
 262               if(kf == kProton)  fSpecp += 1;
 263           }
 264 // --------------------------------------------------------------------------
 265
 266           if (!fSelectAll) {
 267               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 268               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
 269           }
 270
 271 // Put particle on the stack if selected
 272
 273           if (selected) {
 274               nc++;
 275               pSelected[i] = 1;
 276           } // selected
 277       } // particle loop final state
 278
 279 // Write particles to stack
 280
 281       for (i = 0; i<np; i++) {
 282           TParticle *  iparticle = (TParticle *) fParticles.At(i);
 283           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
 284           if (pSelected[i]) {
 285
 286               kf   = iparticle->GetPdgCode();
 287               ks   = iparticle->GetStatusCode();
 288
 289               p[0] = iparticle->Px();
 290               p[1] = iparticle->Py();
 291               p[2] = iparticle->Pz();
 292               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 293               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 294               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 295
 296               tof = kconv*iparticle->T();
 297
 298               imo = -1;
 299               TParticle* mother = 0;
 300               if (hasMother) {
 301                   imo = iparticle->GetFirstMother();
 302                   mother = (TParticle *) fParticles.At(imo);
 303                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
 304               } // if has mother
 305
 306
 307
 308               Bool_t tFlag = (fTrackIt && (ks == 1));
 309               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 310               KeepTrack(nt);
 311               newPos[i] = nt;
 312           } // if selected
 313       } // particle loop
 314       delete[] newPos;
 315       delete[] pSelected;
 316       if (nc>0) {
 317           jev += nc;
 318           if (jev >= fNpart || fNpart == -1) {
 319               break;
 320           }
 321       }
 322   } // event loop
 323   MakeHeader();
 324   SetHighWaterMark(nt);
 325 }
 326
 327 //______________________________________________________________________________
 328 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 329 {
 330 //
 331 // Looks recursively if one of the daughters has been selected
 332 //
 333 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 334     Int_t imin = -1;
 335     Int_t imax = -1;
 336     Int_t i;
 337     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 338     Bool_t selected = kFALSE;
 339     if (hasDaughters) {
 340         imin = iparticle->GetFirstDaughter();
 341         imax = iparticle->GetLastDaughter();
 342         for (i = imin; i <= imax; i++){
 343             TParticle *  jparticle = (TParticle *) fParticles.At(i);
 344             Int_t ip = jparticle->GetPdgCode();
 345             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 346                 selected=kTRUE; break;
 347             }
 348             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 349         }
 350     } else {
 351         return kFALSE;
 352     }
 353     return selected;
 354 }
 355
 356
 357
 358 //______________________________________________________________________________
 359 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
 360 {
 361 // Select flavor of particle
 362 // 0: all
 363 // 4: charm and beauty
 364 // 5: beauty
 365     Bool_t res = 0;
 366
 367     if (fFlavor == 0) {
 368         res = kTRUE;
 369     } else {
 370         Int_t ifl = TMath::Abs(pid/100);
 371         if (ifl > 10) ifl/=10;
 372         res = (fFlavor == ifl);
 373     }
 374 //
 375 //  This part if gamma writing is inhibited
 376     if (fNoGammas)
 377         res = res && (pid != kGamma && pid != kPi0);
 378 //
 379     return res;
 380 }
 381
 382 //______________________________________________________________________________
 383 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 384 {
 385 // Return true for a stable particle
 386 //
 387
 388 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
 389     if (particle->GetStatusCode() == 1) return kTRUE;
 390     else return kFALSE;
 391
 392 }
 393
 394 //______________________________________________________________________________
 395 void AliGenDPMjet::MakeHeader()
 396 {
 397 // Builds the event header, to be called after each event
 398     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
 399     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
 400     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
 401     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
 402     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
 403                                                          fDPMjet->GetfIt());
 404  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 405 // Bookkeeping for kinematic bias
 406     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 407 // Event Vertex
 408     header->SetPrimaryVertex(fVertex);
 409     gAlice->SetGenEventHeader(header);
 410     AddHeader(header);
 411 }
 412
 413 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 414 {
 415     // Add header to container or runloader
 416     if (fContainer) {
 417         fContainer->AddHeader(header);
 418     } else {
 419         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
 420     }
 421 }
 422
 423
 424 //______________________________________________________________________________
 425 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 426 {
 427 // Assignment operator
 428     return *this;
 429 }
 430
 431
 432 void AliGenDPMjet::FinishRun()
 433 {
 434     // Print run statistics
 435     fDPMjet->Dt_Dtuout();
 436 }
 437
 438
 439
 440 //______________________________________________________________________________