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