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