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      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      fNoGammas(0),
  85      fLHC(0),
  86      fPi0Decay(0),
  87      fGenImpPar(0.),
  88      fProcess(kDpmMb)
  89 {
  90 // Default PbPb collisions at 5. 5 TeV
  91 //
  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      fEnergyCMS(5500.),
 104      fMinImpactParam(0.),
 105      fMaxImpactParam(5.),
 106      fICentr(0),
 107      fSelectAll(0),
 108      fFlavor(0),
 109      fTrials(0),
 110      fSpectators(1),
 111      fSpecn(0),
 112      fSpecp(0),
 113      fDPMjet(0),
 114      fNoGammas(0),
 115      fLHC(0),
 116      fPi0Decay(0),
 117      fGenImpPar(0.),
 118      fProcess(kDpmMb)
 119 {
 120     // Dummy copy constructor
 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 && 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) 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) ? newPos[imo] : -1;
 296               } // if has mother
 297
 298               Bool_t tFlag = (fTrackIt && (ks == 1));
 299
 300               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
 301               KeepTrack(nt);
 302               newPos[i] = nt;
 303           } // if selected
 304       } // particle loop
 305       delete[] newPos;
 306       delete[] pSelected;
 307       if (nc>0) {
 308           jev += nc;
 309           if (jev >= fNpart || fNpart == -1) {
 310               break;
 311           }
 312       }
 313   } // event loop
 314   MakeHeader();
 315   SetHighWaterMark(nt);
 316 }
 317
 318 //______________________________________________________________________________
 319 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 320 {
 321 //
 322 // Looks recursively if one of the daughters has been selected
 323 //
 324 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 325     Int_t imin = -1;
 326     Int_t imax = -1;
 327     Int_t i;
 328     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 329     Bool_t selected = kFALSE;
 330     if (hasDaughters) {
 331         imin = iparticle->GetFirstDaughter();
 332         imax = iparticle->GetLastDaughter();
 333         for (i = imin; i <= imax; i++){
 334             TParticle *  jparticle = (TParticle *) fParticles.At(i);
 335             Int_t ip = jparticle->GetPdgCode();
 336             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 337                 selected=kTRUE; break;
 338             }
 339             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 340         }
 341     } else {
 342         return kFALSE;
 343     }
 344     return selected;
 345 }
 346
 347
 348
 349 //______________________________________________________________________________
 350 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
 351 {
 352 // Select flavor of particle
 353 // 0: all
 354 // 4: charm and beauty
 355 // 5: beauty
 356     Bool_t res = 0;
 357
 358     if (fFlavor == 0) {
 359         res = kTRUE;
 360     } else {
 361         Int_t ifl = TMath::Abs(pid/100);
 362         if (ifl > 10) ifl/=10;
 363         res = (fFlavor == ifl);
 364     }
 365 //
 366 //  This part if gamma writing is inhibited
 367     if (fNoGammas)
 368         res = res && (pid != kGamma && pid != kPi0);
 369 //
 370     return res;
 371 }
 372
 373 //______________________________________________________________________________
 374 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 375 {
 376 // Return true for a stable particle
 377 //
 378
 379 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
 380     if (particle->GetStatusCode() == 1) return kTRUE;
 381     else return kFALSE;
 382
 383 }
 384
 385 //______________________________________________________________________________
 386 void AliGenDPMjet::MakeHeader()
 387 {
 388 // Builds the event header, to be called after each event
 389     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
 390     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
 391     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
 392     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
 393     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
 394                                                          fDPMjet->GetfIt());
 395  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 396 // Bookkeeping for kinematic bias
 397     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 398 // Event Vertex
 399     header->SetPrimaryVertex(fVertex);
 400     gAlice->SetGenEventHeader(header);
 401     AddHeader(header);
 402 }
 403
 404 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 405 {
 406     // Add header to container or runloader
 407     if (fContainer) {
 408         fContainer->AddHeader(header);
 409     } else {
 410         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
 411     }
 412 }
 413
 414
 415 //______________________________________________________________________________
 416 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 417 {
 418 // Assignment operator
 419     return *this;
 420 }
 421
 422
 423 void AliGenDPMjet::FinishRun()
 424 {
 425     // Print run statistics
 426     fDPMjet->Dt_Dtuout();
 427 }
 428
 429
 430
 431 //______________________________________________________________________________