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()
66 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
70 //______________________________________________________________________________
71 AliGenDPMjet::AliGenDPMjet(Int_t npart)
85 fParticles(new TClonesArray("TParticle",10000)),
92 // Default PbPb collisions at 5. 5 TeV
95 fTitle= "Particle Generator using DPMJET";
99 AliDpmJetRndm::SetDpmJetRandom(GetRandom());
102 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
123 // Dummy copy constructor
126 //______________________________________________________________________________
127 AliGenDPMjet::~AliGenDPMjet()
132 //______________________________________________________________________________
133 void AliGenDPMjet::Init()
137 SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget,
138 fBeamEn,fEnergyCMS));
140 fDPMjet=(TDPMjet*) fMCEvGen;
142 // **** Flag to force central production
143 // fICentr=1. central production forced
144 // fICentr<0 && fICentr>-100 -> bmin = fMinImpactParam, bmax = fMaxImpactParam
145 // fICentr<-99 -> fraction of x-sec. = XSFRAC
146 // fICentr=-1. -> evaporation/fzc suppressed
147 // fICentr<-1. -> evaporation/fzc suppressed
148 if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1);
150 fDPMjet->SetfFCentr(fICentr);
151 fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam);
152 fDPMjet->SetPi0Decay(fPi0Decay);
156 fDPMjet->Initialize();
160 //______________________________________________________________________________
161 void AliGenDPMjet::Generate()
163 // Generate one event
165 Float_t polar[3] = {0,0,0};
166 Float_t origin[3] = {0,0,0};
167 Float_t origin0[3] = {0,0,0};
171 // converts from mm/c to s
172 const Float_t kconv = 0.001/2.999792458e8;
178 // Set collision vertex position
179 if (fVertexSmear == kPerEvent) Vertex();
183 // Generate one event
184 // --------------------------------------------------------------------------
187 // --------------------------------------------------------------------------
188 fDPMjet->GenerateEvent();
191 fDPMjet->ImportParticles(fParticles,"All");
195 fGenImpPar = fDPMjet->GetBImpac();
197 Int_t np = fParticles->GetEntriesFast();
198 printf("\n **************************************************%d\n",np);
202 Int_t* newPos = new Int_t[np];
203 Int_t* pSelected = new Int_t[np];
205 for (i = 0; i<np; i++) {
212 fVertex[0] = origin0[0];
213 fVertex[1] = origin0[1];
214 fVertex[2] = origin0[2];
216 // First select parent particles
218 for (i = 0; i<np; i++) {
219 TParticle *iparticle = (TParticle *) fParticles->At(i);
221 // Is this a parent particle ?
223 if (Stable(iparticle)) continue;
225 Bool_t selected = kTRUE;
226 Bool_t hasSelectedDaughters = kFALSE;
228 kf = iparticle->GetPdgCode();
229 if (kf == 92) continue;
230 ks = iparticle->GetStatusCode();
231 // No initial state partons
232 if (ks==21) continue;
234 if (!fSelectAll) selected = KinematicSelection(iparticle, 0) &&
236 hasSelectedDaughters = DaughtersSelection(iparticle);
238 // Put particle on the stack if it is either selected or
239 // it is the mother of at least one seleted particle
241 if (selected || hasSelectedDaughters) {
245 } // particle loop parents
247 // Now select the final state particles
250 for (i=0; i<np; i++) {
251 TParticle *iparticle = (TParticle *) fParticles->At(i);
253 // Is this a final state particle ?
255 if (!Stable(iparticle)) continue;
257 Bool_t selected = kTRUE;
258 kf = iparticle->GetPdgCode();
259 ks = iparticle->GetStatusCode();
261 // --------------------------------------------------------------------------
262 // Count spectator neutrons and protons (ks == 13, 14)
263 if(ks == 13 || ks == 14){
264 if(kf == kNeutron) fSpecn += 1;
265 if(kf == kProton) fSpecp += 1;
267 // --------------------------------------------------------------------------
270 selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
271 if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
274 // Put particle on the stack if selected
280 } // particle loop final state
282 // Write particles to stack
284 for (i = 0; i<np; i++) {
285 TParticle * iparticle = (TParticle *) fParticles->At(i);
286 Bool_t hasMother = (iparticle->GetFirstMother()>=0);
289 kf = iparticle->GetPdgCode();
290 ks = iparticle->GetStatusCode();
292 p[0] = iparticle->Px();
293 p[1] = iparticle->Py();
294 p[2] = iparticle->Pz();
295 origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
296 origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
297 origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
299 tof = kconv*iparticle->T();
302 TParticle* mother = 0;
304 imo = iparticle->GetFirstMother();
305 mother = (TParticle *) fParticles->At(imo);
306 imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
309 Bool_t tFlag = (fTrackIt && (ks == 1));
311 PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
320 if (jev >= fNpart || fNpart == -1) {
326 SetHighWaterMark(nt);
329 //______________________________________________________________________________
330 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
333 // Looks recursively if one of the daughters has been selected
335 // printf("\n Consider daughters %d:",iparticle->GetPdgCode());
339 Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
340 Bool_t selected = kFALSE;
342 imin = iparticle->GetFirstDaughter();
343 imax = iparticle->GetLastDaughter();
344 for (i = imin; i <= imax; i++){
345 TParticle * jparticle = (TParticle *) fParticles->At(i);
346 Int_t ip = jparticle->GetPdgCode();
347 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
348 selected=kTRUE; break;
350 if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
360 //______________________________________________________________________________
361 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
363 // Select flavor of particle
365 // 4: charm and beauty
372 Int_t ifl = TMath::Abs(pid/100);
373 if (ifl > 10) ifl/=10;
374 res = (fFlavor == ifl);
377 // This part if gamma writing is inhibited
379 res = res && (pid != kGamma && pid != kPi0);
384 //______________________________________________________________________________
385 Bool_t AliGenDPMjet::Stable(TParticle* particle)
387 // Return true for a stable particle
390 // if (particle->GetFirstDaughter() < 0 ) return kTRUE;
391 if (particle->GetStatusCode() == 1) return kTRUE;
396 //______________________________________________________________________________
397 void AliGenDPMjet::MakeHeader()
399 // Builds the event header, to be called after each event
400 AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
401 ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
402 ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
403 ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
404 ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(),
406 ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
407 // Bookkeeping for kinematic bias
408 ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
410 header->SetPrimaryVertex(fVertex);
411 gAlice->SetGenEventHeader(header);
415 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
417 // Add header to container or runloader
419 fContainer->AddHeader(header);
421 AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
426 //______________________________________________________________________________
427 AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/)
429 // Assignment operator
435 //______________________________________________________________________________