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