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