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   Double_t polar[3]    =   {0,0,0};
 175   Double_t origin[3]   =   {0,0,0};
 176   Double_t p[4]        =   {0};
 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               p[3] = iparticle->Energy();
 334               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 335               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 336               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 337
 338               tof = kconv*iparticle->T();
 339
 340               imo = -1;
 341               TParticle* mother = 0;
 342               if (hasMother) {
 343                   imo = iparticle->GetFirstMother();
 344                   mother = (TParticle *) fParticles.At(imo);
 345                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
 346               } // if has mother
 347
 348
 349
 350               Bool_t tFlag = (fTrackIt && (ks == 1));
 351               PushTrack(tFlag, imo, kf,
 352                         p[0], p[1], p[2], p[3],
 353                         origin[0], origin[1], origin[2], tof,
 354                         polar[0], polar[1], polar[2],
 355                         kPNoProcess, nt, 1., ks);
 356               KeepTrack(nt);
 357               newPos[i] = nt;
 358           } // if selected
 359       } // particle loop
 360       delete[] newPos;
 361       delete[] pSelected;
 362       if (nc>0) {
 363           jev += nc;
 364           if (jev >= fNpart || fNpart == -1) {
 365               break;
 366           }
 367       }
 368   } // event loop
 369   MakeHeader();
 370   SetHighWaterMark(nt);
 371 }
 372
 373 //______________________________________________________________________________
 374 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
 375 {
 376 //
 377 // Looks recursively if one of the daughters has been selected
 378 //
 379 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
 380     Int_t imin = -1;
 381     Int_t imax = -1;
 382     Int_t i;
 383     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
 384     Bool_t selected = kFALSE;
 385     if (hasDaughters) {
 386         imin = iparticle->GetFirstDaughter();
 387         imax = iparticle->GetLastDaughter();
 388         for (i = imin; i <= imax; i++){
 389             TParticle *  jparticle = (TParticle *) fParticles.At(i);
 390             Int_t ip = jparticle->GetPdgCode();
 391             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
 392                 selected=kTRUE; break;
 393             }
 394             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
 395         }
 396     } else {
 397         return kFALSE;
 398     }
 399     return selected;
 400 }
 401
 402
 403
 404 //______________________________________________________________________________
 405 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
 406 {
 407 // Select flavor of particle
 408 // 0: all
 409 // 4: charm and beauty
 410 // 5: beauty
 411     Bool_t res = 0;
 412
 413     if (fFlavor == 0) {
 414         res = kTRUE;
 415     } else {
 416         Int_t ifl = TMath::Abs(pid/100);
 417         if (ifl > 10) ifl/=10;
 418         res = (fFlavor == ifl);
 419     }
 420 //
 421 //  This part if gamma writing is inhibited
 422     if (fNoGammas)
 423         res = res && (pid != kGamma && pid != kPi0);
 424 //
 425     return res;
 426 }
 427
 428 //______________________________________________________________________________
 429 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
 430 {
 431 // Return true for a stable particle
 432 //
 433
 434 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
 435     if (particle->GetStatusCode() == 1) return kTRUE;
 436     else return kFALSE;
 437
 438 }
 439
 440 //______________________________________________________________________________
 441 void AliGenDPMjet::MakeHeader()
 442 {
 443   printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
 444 // Builds the event header, to be called after each event
 445     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
 446     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
 447     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
 448     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
 449     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(),
 450                                                          fDPMjet->GetTargParticipants());
 451      ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
 452     // Bookkeeping for kinematic bias
 453     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
 454     // Event Vertex
 455     header->SetPrimaryVertex(fVertex);
 456     gAlice->SetGenEventHeader(header);
 457     AddHeader(header);
 458 }
 459
 460 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
 461 {
 462     // Add header to container or runloader
 463     if (fContainer) {
 464         fContainer->AddHeader(header);
 465     } else {
 466         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
 467     }
 468 }
 469
 470
 471 //______________________________________________________________________________
 472 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
 473 {
 474 // Assignment operator
 475     return *this;
 476 }
 477
 478
 479 void AliGenDPMjet::FinishRun()
 480 {
 481     // Print run statistics
 482     fDPMjet->Dt_Dtuout();
 483 }
 484
 485
 486
 487 //______________________________________________________________________________