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      fEnergyCMS(5500.),
  48      fMinImpactParam(0.),
  49      fMaxImpactParam(5.),
  50      fICentr(0),
  51      fSelectAll(0),
  52      fFlavor(0),
  53      fTrials(0),
  54      fSpectators(1),
  55      fSpecn(0),
  56      fSpecp(0),
  57      fDPMjet(0),
  58      fParticles(0),
  59      fNoGammas(0),
  60      fLHC(0),
  61      fPi0Decay(0),
  62      fGenImpPar(0.),
  63      fProcess(kDpmMb)
  64 {
  65 // Constructor
  66     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
  67 }
  68
  69
  70 //______________________________________________________________________________
  71 AliGenDPMjet::AliGenDPMjet(Int_t npart)
  72     :AliGenMC(npart),
  73      fBeamEn(2750.),
  74      fEnergyCMS(5500.),
  75      fMinImpactParam(0.),
  76      fMaxImpactParam(5.),
  77      fICentr(0),
  78      fSelectAll(0),
  79      fFlavor(0),
  80      fTrials(0),
  81      fSpectators(1),
  82      fSpecn(0),
  83      fSpecp(0),
  84      fDPMjet(0),
  85      fParticles(new TClonesArray("TParticle",10000)),
  86      fNoGammas(0),
  87      fLHC(0),
  88      fPi0Decay(0),
  89      fGenImpPar(0.),
  90      fProcess(kDpmMb)
  91 {
  92 // Default PbPb collisions at 5. 5 TeV
  93 //
  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      fEnergyCMS(5500.),
 106      fMinImpactParam(0.),
 107      fMaxImpactParam(5.),
 108      fICentr(0),
 109      fSelectAll(0),
 110      fFlavor(0),
 111      fTrials(0),
 112      fSpectators(1),
 113      fSpecn(0),
 114      fSpecp(0),
 115      fDPMjet(0),
 116      fParticles(0),
 117      fNoGammas(0),
 118      fLHC(0),
 119      fPi0Decay(0),
 120      fGenImpPar(0.),
 121      fProcess(kDpmMb)
 122 {
 123     // Dummy copy constructor
 124 }
 125
 126 //______________________________________________________________________________
 127 AliGenDPMjet::~AliGenDPMjet()
 128 {
 129 // Destructor
 130     delete fParticles;
 131 }
 132 //______________________________________________________________________________
 133 void AliGenDPMjet::Init()
 134 {
 135 // Initialization
 136
 137     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
 138                       fBeamEn,fEnergyCMS));
 139
 140     fDPMjet=(TDPMjet*) fMCEvGen;
 141     //
 142     // **** Flag to force central production
 143     // fICentr=1. central production forced
 144     // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
 145     // fICentr<-99 -> fraction of x-sec. = XSFRAC
 146     // fICentr=-1. -> evaporation/fzc suppressed
 147     // fICentr<-1. -> evaporation/fzc suppressed
 148     if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
 149
 150     fDPMjet->SetfFCentr(fICentr);
 151     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
 152     fDPMjet->SetPi0Decay(fPi0Decay);
 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 origin0[3]  =   {0,0,0};
 168   Float_t p[3];
 169   Float_t tof;
 170
 171 //  converts from mm/c to s
 172   const Float_t kconv = 0.001/2.999792458e8;
 173   Int_t nt  = 0;
 174   Int_t jev = 0;
 175   Int_t kf, ks, imo;
 176   kf = 0;
 177   fTrials = 0;
 178   //  Set collision vertex position
 179   if (fVertexSmear == kPerEvent) Vertex();
 180
 181   while(1)
 182   {
 183 //    Generate one event
 184 // --------------------------------------------------------------------------
 185       fSpecn = 0;
 186       fSpecp = 0;
 187 // --------------------------------------------------------------------------
 188       fDPMjet->GenerateEvent();
 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 //      Get event vertex
 211
 212       fVertex[0] = origin0[0];
 213       fVertex[1] = origin0[1];
 214       fVertex[2] = origin0[2];
 215
 216 //      First select parent particles
 217
 218       for (i = 0; i<np; i++) {
 219           TParticle *iparticle = (TParticle *) fParticles->At(i);
 220
 221 // Is this a parent particle ?
 222
 223           if (Stable(iparticle)) continue;
 224
 225           Bool_t  selected             =  kTRUE;
 226           Bool_t  hasSelectedDaughters =  kFALSE;
 227
 228           kf = iparticle->GetPdgCode();
 229           if (kf == 92) continue;
 230           ks = iparticle->GetStatusCode();
 231 // No initial state partons
 232           if (ks==21) continue;
 233
 234           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
 235                                SelectFlavor(kf);
 236           hasSelectedDaughters = DaughtersSelection(iparticle);
 237
 238 // Put particle on the stack if it is either selected or
 239 // it is the mother of at least one seleted particle
 240
 241           if (selected || hasSelectedDaughters) {
 242               nc++;
 243               pSelected[i] = 1;
 244           } // selected
 245       } // particle loop parents
 246
 247 // Now select the final state particles
 248
 249
 250       for (i=0; i<np; i++) {
 251           TParticle *iparticle = (TParticle *) fParticles->At(i);
 252
 253 // Is this a final state particle ?
 254
 255           if (!Stable(iparticle)) continue;
 256
 257           Bool_t  selected =  kTRUE;
 258           kf = iparticle->GetPdgCode();
 259           ks = iparticle->GetStatusCode();
 260
 261 // --------------------------------------------------------------------------
 262 // Count spectator neutrons and protons (ks == 13, 14)
 263           if(ks == 13 || ks == 14){
 264               if(kf == kNeutron) fSpecn += 1;
 265               if(kf == kProton)  fSpecp += 1;
 266           }
 267 // --------------------------------------------------------------------------
 268
 269           if (!fSelectAll) {
 270               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 271               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
 272           }
 273
 274 // Put particle on the stack if selected
 275
 276           if (selected) {
 277               nc++;
 278               pSelected[i] = 1;
 279           } // selected
 280       } // particle loop final state
 281
 282 // Write particles to stack
 283
 284       for (i = 0; i<np; i++) {
 285           TParticle *  iparticle = (TParticle *) fParticles->At(i);
 286           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
 287           if (pSelected[i]) {
 288
 289               kf   = iparticle->GetPdgCode();
 290               ks   = iparticle->GetStatusCode();
 291
 292               p[0] = iparticle->Px();
 293               p[1] = iparticle->Py();
 294               p[2] = iparticle->Pz();
 295               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 296               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 297               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 298
 299               tof = kconv*iparticle->T();
 300
 301               imo = -1;
 302               TParticle* mother = 0;
 303               if (hasMother) {
 304                   imo = iparticle->GetFirstMother();
 305                   mother = (TParticle *) fParticles->At(imo);
 306                   imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
 307               } // if has mother
 308
 309               Bool_t tFlag = (fTrackIt && (ks == 1));
 310
 311               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 312               KeepTrack(nt);
 313               newPos[i] = nt;
 314           } // if selected
 315       } // particle loop
 316       delete[] newPos;
 317       delete[] pSelected;
 318       if (nc>0) {
 319           jev += nc;
 320           if (jev >= fNpart || fNpart == -1) {
 321               break;
 322           }
 323       }
 324   } // event loop
 325   MakeHeader();
 326   SetHighWaterMark(nt);
 327 }
 328
 329 //______________________________________________________________________________
 330 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 331 {
 332 //
 333 // Looks recursively if one of the daughters has been selected
 334 //
 335 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 336     Int_t imin = -1;
 337     Int_t imax = -1;
 338     Int_t i;
 339     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 340     Bool_t selected = kFALSE;
 341     if (hasDaughters) {
 342         imin = iparticle->GetFirstDaughter();
 343         imax = iparticle->GetLastDaughter();
 344         for (i = imin; i <= imax; i++){
 345             TParticle *  jparticle = (TParticle *) fParticles->At(i);
 346             Int_t ip = jparticle->GetPdgCode();
 347             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 348                 selected=kTRUE; break;
 349             }
 350             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 351         }
 352     } else {
 353         return kFALSE;
 354     }
 355     return selected;
 356 }
 357
 358
 359
 360 //______________________________________________________________________________
 361 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
 362 {
 363 // Select flavor of particle
 364 // 0: all
 365 // 4: charm and beauty
 366 // 5: beauty
 367     Bool_t res = 0;
 368
 369     if (fFlavor == 0) {
 370         res = kTRUE;
 371     } else {
 372         Int_t ifl = TMath::Abs(pid/100);
 373         if (ifl > 10) ifl/=10;
 374         res = (fFlavor == ifl);
 375     }
 376 //
 377 //  This part if gamma writing is inhibited
 378     if (fNoGammas)
 379         res = res && (pid != kGamma && pid != kPi0);
 380 //
 381     return res;
 382 }
 383
 384 //______________________________________________________________________________
 385 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 386 {
 387 // Return true for a stable particle
 388 //
 389
 390 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
 391     if (particle->GetStatusCode() == 1) return kTRUE;
 392     else return kFALSE;
 393
 394 }
 395
 396 //______________________________________________________________________________
 397 void AliGenDPMjet::MakeHeader()
 398 {
 399 // Builds the event header, to be called after each event
 400     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
 401     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
 402     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
 403     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
 404     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
 405                                                          fDPMjet->GetfIt());
 406  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 407 // Bookkeeping for kinematic bias
 408     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 409 // Event Vertex
 410     header->SetPrimaryVertex(fVertex);
 411     gAlice->SetGenEventHeader(header);
 412  AddHeader(header);
 413 }
 414
 415 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 416 {
 417     // Add header to container or runloader
 418     if (fContainer) {
 419         fContainer->AddHeader(header);
 420     } else {
 421         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
 422     }
 423 }
 424
 425
 426 //______________________________________________________________________________
 427 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 428 {
 429 // Assignment operator
 430     return *this;
 431 }
 432
 433
 434
 435 //______________________________________________________________________________