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      fTriggerMultiplicity(0),
  64      fTriggerMultiplicityEta(0),
  65      fTriggerMultiplicityPtMin(0)
  66 {
  67 // Constructor
  68     fEnergyCMS = 5500.;
  69     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
  70 }
  71
  72
  73 //______________________________________________________________________________
  74 AliGenDPMjet::AliGenDPMjet(Int_t npart)
  75     :AliGenMC(npart),
  76      fBeamEn(2750.),
  77      fMinImpactParam(0.),
  78      fMaxImpactParam(5.),
  79      fICentr(0),
  80      fSelectAll(0),
  81      fFlavor(0),
  82      fTrials(0),
  83      fSpectators(1),
  84      fSpecn(0),
  85      fSpecp(0),
  86      fDPMjet(0),
  87      fNoGammas(0),
  88      fLHC(0),
  89      fPi0Decay(1),
  90      fDecayAll(0),
  91      fGenImpPar(0.),
  92      fProcess(kDpmMb),
  93      fTriggerMultiplicity(0),
  94      fTriggerMultiplicityEta(0),
  95      fTriggerMultiplicityPtMin(0)
  96 {
  97 // Default PbPb collisions at 5. 5 TeV
  98 //
  99     fEnergyCMS = 5500.;
 100     fName = "DPMJET";
 101     fTitle= "Particle Generator using DPMJET";
 102     SetTarget();
 103     SetProjectile();
 104     fVertex.Set(3);
 105     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
 106 }
 107
 108 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
 109     :AliGenMC(),
 110      fBeamEn(2750.),
 111      fMinImpactParam(0.),
 112      fMaxImpactParam(5.),
 113      fICentr(0),
 114      fSelectAll(0),
 115      fFlavor(0),
 116      fTrials(0),
 117      fSpectators(1),
 118      fSpecn(0),
 119      fSpecp(0),
 120      fDPMjet(0),
 121      fNoGammas(0),
 122      fLHC(0),
 123      fPi0Decay(1),
 124      fDecayAll(0),
 125      fGenImpPar(0.),
 126      fProcess(kDpmMb),
 127      fTriggerMultiplicity(0),
 128      fTriggerMultiplicityEta(0),
 129      fTriggerMultiplicityPtMin(0)
 130 {
 131     // Dummy copy constructor
 132     fEnergyCMS = 5500.;
 133 }
 134
 135 //______________________________________________________________________________
 136 AliGenDPMjet::~AliGenDPMjet()
 137 {
 138 // Destructor
 139 }
 140 //______________________________________________________________________________
 141 void AliGenDPMjet::Init()
 142 {
 143 // Initialization
 144
 145     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
 146                       fBeamEn,fEnergyCMS));
 147
 148     fDPMjet=(TDPMjet*) fMCEvGen;
 149     //
 150     // **** Flag to force central production
 151     // fICentr=1. central production forced
 152     // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
 153     // fICentr<-99 -> fraction of x-sec. = XSFRAC
 154     // fICentr=-1. -> evaporation/fzc suppressed
 155     // fICentr<-1. -> evaporation/fzc suppressed
 156     if (fAProjectile == 1 && TMath::Abs(fZProjectile == 1)) fDPMjet->SetfIdp(1);
 157
 158     fDPMjet->SetfFCentr(fICentr);
 159     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
 160     fDPMjet->SetPi0Decay(fPi0Decay);
 161     fDPMjet->SetDecayAll(fDecayAll);
 162 //
 163 //  Initialize DPMjet
 164 //
 165     fDPMjet->Initialize();
 166 }
 167
 168
 169 //______________________________________________________________________________
 170 void AliGenDPMjet::Generate()
 171 {
 172 // Generate one event
 173
 174   Float_t polar[3]    =   {0,0,0};
 175   Float_t origin[3]   =   {0,0,0};
 176   Float_t p[3];
 177   Float_t tof;
 178
 179 //  converts from mm/c to s
 180   const Float_t kconv = 0.001/2.999792458e8;
 181   Int_t nt  = 0;
 182   Int_t jev = 0;
 183   Int_t kf, ks, imo;
 184   kf = 0;
 185   fTrials = 0;
 186   //  Set collision vertex position
 187   if (fVertexSmear == kPerEvent) Vertex();
 188
 189   while(1)
 190   {
 191 //    Generate one event
 192 // --------------------------------------------------------------------------
 193       fSpecn = 0;
 194       fSpecp = 0;
 195 // --------------------------------------------------------------------------
 196       fDPMjet->GenerateEvent();
 197
 198       fTrials++;
 199
 200       fDPMjet->ImportParticles(&fParticles,"All");
 201       if (fLHC) Boost();
 202
 203       // Temporaneo
 204       fGenImpPar = fDPMjet->GetBImpac();
 205
 206       Int_t np = fParticles.GetEntriesFast();
 207       //
 208       // Multiplicity Trigger
 209       if (fTriggerMultiplicity > 0) {
 210         Int_t multiplicity = 0;
 211         for (Int_t i = 0; i < np; i++) {
 212           TParticle *  iparticle = (TParticle *) fParticles.At(i);
 213
 214           Int_t statusCode = iparticle->GetStatusCode();
 215
 216           // Initial state particle
 217           if (statusCode != 1)
 218             continue;
 219           // eta cut
 220           if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
 221             continue;
 222           // pt cut
 223           if (iparticle->Pt() < fTriggerMultiplicityPtMin)
 224             continue;
 225
 226           TParticlePDG* pdgPart = iparticle->GetPDG();
 227           if (pdgPart && pdgPart->Charge() == 0)
 228             continue;
 229           ++multiplicity;
 230         }
 231         //
 232         //
 233         if (multiplicity < fTriggerMultiplicity) continue;
 234         Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
 235       }
 236
 237       Int_t nc = 0;
 238       if (np == 0) continue;
 239
 240       Int_t i;
 241       Int_t* newPos     = new Int_t[np];
 242       Int_t* pSelected  = new Int_t[np];
 243
 244       for (i = 0; i<np; i++) {
 245           newPos[i]    = i;
 246           pSelected[i] = 0;
 247       }
 248
 249 //      First select parent particles
 250
 251       for (i = 0; i<np; i++) {
 252           TParticle *iparticle = (TParticle *) fParticles.At(i);
 253
 254 // Is this a parent particle ?
 255
 256           if (Stable(iparticle)) continue;
 257
 258           Bool_t  selected             =  kTRUE;
 259           Bool_t  hasSelectedDaughters =  kFALSE;
 260
 261           kf = iparticle->GetPdgCode();
 262           if (kf == 92 || kf == 99999) continue;
 263           ks = iparticle->GetStatusCode();
 264 // No initial state partons
 265           if (ks==21) continue;
 266
 267           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
 268                                SelectFlavor(kf);
 269
 270
 271           hasSelectedDaughters = DaughtersSelection(iparticle);
 272
 273
 274 // Put particle on the stack if it is either selected or
 275 // it is the mother of at least one seleted particle
 276
 277           if (selected || hasSelectedDaughters) {
 278               nc++;
 279               pSelected[i] = 1;
 280           } // selected
 281       } // particle loop parents
 282
 283 // Now select the final state particles
 284
 285
 286       for (i=0; i<np; i++) {
 287           TParticle *iparticle = (TParticle *) fParticles.At(i);
 288
 289 // Is this a final state particle ?
 290
 291           if (!Stable(iparticle)) continue;
 292
 293           Bool_t  selected =  kTRUE;
 294           kf = iparticle->GetPdgCode();
 295           ks = iparticle->GetStatusCode();
 296
 297 // --------------------------------------------------------------------------
 298 // Count spectator neutrons and protons (ks == 13, 14)
 299           if(ks == 13 || ks == 14){
 300               if(kf == kNeutron) fSpecn += 1;
 301               if(kf == kProton)  fSpecp += 1;
 302           }
 303 // --------------------------------------------------------------------------
 304
 305           if (!fSelectAll) {
 306               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 307               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
 308           }
 309
 310 // Put particle on the stack if selected
 311
 312           if (selected) {
 313               nc++;
 314               pSelected[i] = 1;
 315           } // selected
 316       } // particle loop final state
 317
 318 // Write particles to stack
 319
 320       for (i = 0; i<np; i++) {
 321           TParticle *  iparticle = (TParticle *) fParticles.At(i);
 322           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
 323           if (pSelected[i]) {
 324
 325               kf   = iparticle->GetPdgCode();
 326               ks   = iparticle->GetStatusCode();
 327
 328               p[0] = iparticle->Px();
 329               p[1] = iparticle->Py();
 330               p[2] = iparticle->Pz();
 331               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 332               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 333               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 334
 335               tof = kconv*iparticle->T();
 336
 337               imo = -1;
 338               TParticle* mother = 0;
 339               if (hasMother) {
 340                   imo = iparticle->GetFirstMother();
 341                   mother = (TParticle *) fParticles.At(imo);
 342                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
 343               } // if has mother
 344
 345
 346
 347               Bool_t tFlag = (fTrackIt && (ks == 1));
 348               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 349               KeepTrack(nt);
 350               newPos[i] = nt;
 351           } // if selected
 352       } // particle loop
 353       delete[] newPos;
 354       delete[] pSelected;
 355       if (nc>0) {
 356           jev += nc;
 357           if (jev >= fNpart || fNpart == -1) {
 358               break;
 359           }
 360       }
 361   } // event loop
 362   MakeHeader();
 363   SetHighWaterMark(nt);
 364 }
 365
 366 //______________________________________________________________________________
 367 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 368 {
 369 //
 370 // Looks recursively if one of the daughters has been selected
 371 //
 372 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 373     Int_t imin = -1;
 374     Int_t imax = -1;
 375     Int_t i;
 376     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 377     Bool_t selected = kFALSE;
 378     if (hasDaughters) {
 379         imin = iparticle->GetFirstDaughter();
 380         imax = iparticle->GetLastDaughter();
 381         for (i = imin; i <= imax; i++){
 382             TParticle *  jparticle = (TParticle *) fParticles.At(i);
 383             Int_t ip = jparticle->GetPdgCode();
 384             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 385                 selected=kTRUE; break;
 386             }
 387             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 388         }
 389     } else {
 390         return kFALSE;
 391     }
 392     return selected;
 393 }
 394
 395
 396
 397 //______________________________________________________________________________
 398 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
 399 {
 400 // Select flavor of particle
 401 // 0: all
 402 // 4: charm and beauty
 403 // 5: beauty
 404     Bool_t res = 0;
 405
 406     if (fFlavor == 0) {
 407         res = kTRUE;
 408     } else {
 409         Int_t ifl = TMath::Abs(pid/100);
 410         if (ifl > 10) ifl/=10;
 411         res = (fFlavor == ifl);
 412     }
 413 //
 414 //  This part if gamma writing is inhibited
 415     if (fNoGammas)
 416         res = res && (pid != kGamma && pid != kPi0);
 417 //
 418     return res;
 419 }
 420
 421 //______________________________________________________________________________
 422 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 423 {
 424 // Return true for a stable particle
 425 //
 426
 427 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
 428     if (particle->GetStatusCode() == 1) return kTRUE;
 429     else return kFALSE;
 430
 431 }
 432
 433 //______________________________________________________________________________
 434 void AliGenDPMjet::MakeHeader()
 435 {
 436 // Builds the event header, to be called after each event
 437     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
 438     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
 439     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
 440     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
 441     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
 442                                                          fDPMjet->GetfIt());
 443  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 444 // Bookkeeping for kinematic bias
 445     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 446 // Event Vertex
 447     header->SetPrimaryVertex(fVertex);
 448     gAlice->SetGenEventHeader(header);
 449     AddHeader(header);
 450 }
 451
 452 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 453 {
 454     // Add header to container or runloader
 455     if (fContainer) {
 456         fContainer->AddHeader(header);
 457     } else {
 458         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
 459     }
 460 }
 461
 462
 463 //______________________________________________________________________________
 464 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 465 {
 466 // Assignment operator
 467     return *this;
 468 }
 469
 470
 471 void AliGenDPMjet::FinishRun()
 472 {
 473     // Print run statistics
 474     fDPMjet->Dt_Dtuout();
 475 }
 476
 477
 478
 479 //______________________________________________________________________________