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      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(0),
  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(0),
 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       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 //      First select parent particles
 210
 211       for (i = 0; i<np; i++) {
 212           TParticle *iparticle = (TParticle *) fParticles.At(i);
 213
 214 // Is this a parent particle ?
 215
 216           if (Stable(iparticle)) continue;
 217
 218           Bool_t  selected             =  kTRUE;
 219           Bool_t  hasSelectedDaughters =  kFALSE;
 220
 221           kf = iparticle->GetPdgCode();
 222           if (kf == 92 || kf == 99999) continue;
 223           ks = iparticle->GetStatusCode();
 224 // No initial state partons
 225           if (ks==21) continue;
 226
 227           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
 228                                SelectFlavor(kf);
 229           hasSelectedDaughters = DaughtersSelection(iparticle);
 230
 231 // Put particle on the stack if it is either selected or
 232 // it is the mother of at least one seleted particle
 233
 234           if (selected || hasSelectedDaughters) {
 235               nc++;
 236               pSelected[i] = 1;
 237           } // selected
 238       } // particle loop parents
 239
 240 // Now select the final state particles
 241
 242
 243       for (i=0; i<np; i++) {
 244           TParticle *iparticle = (TParticle *) fParticles.At(i);
 245
 246 // Is this a final state particle ?
 247
 248           if (!Stable(iparticle)) continue;
 249
 250           Bool_t  selected =  kTRUE;
 251           kf = iparticle->GetPdgCode();
 252           ks = iparticle->GetStatusCode();
 253
 254 // --------------------------------------------------------------------------
 255 // Count spectator neutrons and protons (ks == 13, 14)
 256           if(ks == 13 || ks == 14){
 257               if(kf == kNeutron) fSpecn += 1;
 258               if(kf == kProton)  fSpecp += 1;
 259           }
 260 // --------------------------------------------------------------------------
 261
 262           if (!fSelectAll) {
 263               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
 264               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
 265           }
 266
 267 // Put particle on the stack if selected
 268
 269           if (selected) {
 270               nc++;
 271               pSelected[i] = 1;
 272           } // selected
 273       } // particle loop final state
 274
 275 // Write particles to stack
 276
 277       for (i = 0; i<np; i++) {
 278           TParticle *  iparticle = (TParticle *) fParticles.At(i);
 279           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
 280           if (pSelected[i]) {
 281
 282               kf   = iparticle->GetPdgCode();
 283               ks   = iparticle->GetStatusCode();
 284
 285               p[0] = iparticle->Px();
 286               p[1] = iparticle->Py();
 287               p[2] = iparticle->Pz();
 288               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 289               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 290               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 291
 292               tof = kconv*iparticle->T();
 293
 294               imo = -1;
 295               TParticle* mother = 0;
 296               if (hasMother) {
 297                   imo = iparticle->GetFirstMother();
 298                   mother = (TParticle *) fParticles.At(imo);
 299                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
 300               } // if has mother
 301
 302               Bool_t tFlag = (fTrackIt && (ks == 1));
 303               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 304               KeepTrack(nt);
 305               newPos[i] = nt;
 306           } // if selected
 307       } // particle loop
 308       delete[] newPos;
 309       delete[] pSelected;
 310       if (nc>0) {
 311           jev += nc;
 312           if (jev >= fNpart || fNpart == -1) {
 313               break;
 314           }
 315       }
 316   } // event loop
 317   MakeHeader();
 318   SetHighWaterMark(nt);
 319 }
 320
 321 //______________________________________________________________________________
 322 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 323 {
 324 //
 325 // Looks recursively if one of the daughters has been selected
 326 //
 327 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 328     Int_t imin = -1;
 329     Int_t imax = -1;
 330     Int_t i;
 331     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 332     Bool_t selected = kFALSE;
 333     if (hasDaughters) {
 334         imin = iparticle->GetFirstDaughter();
 335         imax = iparticle->GetLastDaughter();
 336         for (i = imin; i <= imax; i++){
 337             TParticle *  jparticle = (TParticle *) fParticles.At(i);
 338             Int_t ip = jparticle->GetPdgCode();
 339             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 340                 selected=kTRUE; break;
 341             }
 342             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 343         }
 344     } else {
 345         return kFALSE;
 346     }
 347     return selected;
 348 }
 349
 350
 351
 352 //______________________________________________________________________________
 353 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
 354 {
 355 // Select flavor of particle
 356 // 0: all
 357 // 4: charm and beauty
 358 // 5: beauty
 359     Bool_t res = 0;
 360
 361     if (fFlavor == 0) {
 362         res = kTRUE;
 363     } else {
 364         Int_t ifl = TMath::Abs(pid/100);
 365         if (ifl > 10) ifl/=10;
 366         res = (fFlavor == ifl);
 367     }
 368 //
 369 //  This part if gamma writing is inhibited
 370     if (fNoGammas)
 371         res = res && (pid != kGamma && pid != kPi0);
 372 //
 373     return res;
 374 }
 375
 376 //______________________________________________________________________________
 377 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 378 {
 379 // Return true for a stable particle
 380 //
 381
 382 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
 383     if (particle->GetStatusCode() == 1) return kTRUE;
 384     else return kFALSE;
 385
 386 }
 387
 388 //______________________________________________________________________________
 389 void AliGenDPMjet::MakeHeader()
 390 {
 391 // Builds the event header, to be called after each event
 392     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
 393     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
 394     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
 395     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
 396     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
 397                                                          fDPMjet->GetfIt());
 398  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 399 // Bookkeeping for kinematic bias
 400     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 401 // Event Vertex
 402     header->SetPrimaryVertex(fVertex);
 403     gAlice->SetGenEventHeader(header);
 404     AddHeader(header);
 405 }
 406
 407 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 408 {
 409     // Add header to container or runloader
 410     if (fContainer) {
 411         fContainer->AddHeader(header);
 412     } else {
 413         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
 414     }
 415 }
 416
 417
 418 //______________________________________________________________________________
 419 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 420 {
 421 // Assignment operator
 422     return *this;
 423 }
 424
 425
 426 void AliGenDPMjet::FinishRun()
 427 {
 428     // Print run statistics
 429     fDPMjet->Dt_Dtuout();
 430 }
 431
 432
 433
 434 //______________________________________________________________________________