1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
24 #include <TParticle.h>
26 #include <TDatabasePDG.h>
27 #include <TParticlePDG.h>
28 #include <TParticleClassPDG.h>
30 #include <TLorentzVector.h>
31 #include <TClonesArray.h>
32 #include "AliRunLoader.h"
33 #include "AliGenDPMjet.h"
34 #include "AliGenDPMjetEventHeader.h"
36 #include "AliDpmJetRndm.h"
37 #include "AliHeader.h"
41 ClassImp(AliGenDPMjet)
43 //______________________________________________________________________________
44 AliGenDPMjet::AliGenDPMjet()
63 fTriggerMultiplicity(0),
64 fTriggerMultiplicityEta(0),
65 fTriggerMultiplicityPtMin(0)
69 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
73 //______________________________________________________________________________
74 AliGenDPMjet::AliGenDPMjet(Int_t npart)
93 fTriggerMultiplicity(0),
94 fTriggerMultiplicityEta(0),
95 fTriggerMultiplicityPtMin(0)
97 // Default PbPb collisions at 5. 5 TeV
101 fTitle= "Particle Generator using DPMJET";
105 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
108 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
127 fTriggerMultiplicity(0),
128 fTriggerMultiplicityEta(0),
129 fTriggerMultiplicityPtMin(0)
131 // Dummy copy constructor
135 //______________________________________________________________________________
136 AliGenDPMjet::~AliGenDPMjet()
140 //______________________________________________________________________________
141 void AliGenDPMjet::Init()
145 SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
146 fBeamEn,fEnergyCMS));
148 fDPMjet=(TDPMjet*) fMCEvGen;
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);
158 fDPMjet->SetfFCentr(fICentr);
159 fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
160 fDPMjet->SetPi0Decay(fPi0Decay);
161 fDPMjet->SetDecayAll(fDecayAll);
165 fDPMjet->Initialize();
169 //______________________________________________________________________________
170 void AliGenDPMjet::Generate()
172 // Generate one event
174 Double_t polar[3] = {0,0,0};
175 Double_t origin[3] = {0,0,0};
179 // converts from mm/c to s
180 const Float_t kconv = 0.001/2.999792458e8;
186 // Set collision vertex position
187 if (fVertexSmear == kPerEvent) Vertex();
191 // Generate one event
192 // --------------------------------------------------------------------------
195 // --------------------------------------------------------------------------
196 fDPMjet->GenerateEvent();
200 fDPMjet->ImportParticles(&fParticles,"All");
204 fGenImpPar = fDPMjet->GetBImpac();
206 if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
208 Int_t np = fParticles.GetEntriesFast();
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);
216 Int_t statusCode = iparticle->GetStatusCode();
218 // Initial state particle
222 if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
225 if (iparticle->Pt() < fTriggerMultiplicityPtMin)
228 TParticlePDG* pdgPart = iparticle->GetPDG();
229 if (pdgPart && pdgPart->Charge() == 0)
235 if (multiplicity < fTriggerMultiplicity) continue;
236 Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
240 if (np == 0) continue;
243 Int_t* newPos = new Int_t[np];
244 Int_t* pSelected = new Int_t[np];
246 for (i = 0; i<np; i++) {
251 // First select parent particles
253 for (i = 0; i<np; i++) {
254 TParticle *iparticle = (TParticle *) fParticles.At(i);
256 // Is this a parent particle ?
258 if (Stable(iparticle)) continue;
260 Bool_t selected = kTRUE;
261 Bool_t hasSelectedDaughters = kFALSE;
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;
269 if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
273 hasSelectedDaughters = DaughtersSelection(iparticle);
276 // Put particle on the stack if it is either selected or
277 // it is the mother of at least one seleted particle
279 if (selected || hasSelectedDaughters) {
283 } // particle loop parents
285 // Now select the final state particles
288 for (i=0; i<np; i++) {
289 TParticle *iparticle = (TParticle *) fParticles.At(i);
291 // Is this a final state particle ?
293 if (!Stable(iparticle)) continue;
295 Bool_t selected = kTRUE;
296 kf = iparticle->GetPdgCode();
297 ks = iparticle->GetStatusCode();
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;
305 // --------------------------------------------------------------------------
308 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
309 if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
312 // Put particle on the stack if selected
318 } // particle loop final state
320 // Write particles to stack
322 for (i = 0; i<np; i++) {
323 TParticle * iparticle = (TParticle *) fParticles.At(i);
324 Bool_t hasMother = (iparticle->GetFirstMother()>=0);
327 kf = iparticle->GetPdgCode();
328 ks = iparticle->GetStatusCode();
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]
338 tof = kconv*iparticle->T();
341 TParticle* mother = 0;
343 imo = iparticle->GetFirstMother();
344 mother = (TParticle *) fParticles.At(imo);
345 imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
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);
364 if (jev >= fNpart || fNpart == -1) {
370 SetHighWaterMark(nt);
373 //______________________________________________________________________________
374 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
377 // Looks recursively if one of the daughters has been selected
379 // printf("\n Consider daughters %d:",iparticle->GetPdgCode());
383 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
384 Bool_t selected = kFALSE;
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;
394 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
404 //______________________________________________________________________________
405 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
407 // Select flavor of particle
409 // 4: charm and beauty
416 Int_t ifl = TMath::Abs(pid/100);
417 if (ifl > 10) ifl/=10;
418 res = (fFlavor == ifl);
421 // This part if gamma writing is inhibited
423 res = res && (pid != kGamma && pid != kPi0);
428 //______________________________________________________________________________
429 Bool_t AliGenDPMjet::Stable(TParticle* particle)
431 // Return true for a stable particle
434 // if (particle->GetFirstDaughter() < 0 ) return kTRUE;
435 if (particle->GetStatusCode() == 1) return kTRUE;
440 //______________________________________________________________________________
441 void AliGenDPMjet::MakeHeader()
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);
455 header->SetPrimaryVertex(fVertex);
456 gAlice->SetGenEventHeader(header);
460 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
462 // Add header to container or runloader
464 fContainer->AddHeader(header);
466 AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
471 //______________________________________________________________________________
472 AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
474 // Assignment operator
479 void AliGenDPMjet::FinishRun()
481 // Print run statistics
482 fDPMjet->Dt_Dtuout();
487 //______________________________________________________________________________