A fix for being able to compile the par file of ANALYSISalice if requested
[u/mrichter/AliRoot.git] / 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       Int_t np = fParticles.GetEntriesFast();
207       //
208       // Multiplicity Trigger
209       if (fTriggerMultiplicity > 0) {
210         Int_t multiplicity = 0;
211         for (Int_t i = 0; i < np; i++) {
212           TParticle *  iparticle = (TParticle *) fParticles.At(i);
213         
214           Int_t statusCode = iparticle->GetStatusCode();
215         
216           // Initial state particle
217           if (statusCode != 1)
218             continue;
219           // eta cut
220           if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
221             continue;
222           // pt cut
223           if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
224             continue;
225           
226           TParticlePDG* pdgPart = iparticle->GetPDG();
227           if (pdgPart && pdgPart->Charge() == 0)
228             continue;
229           ++multiplicity;
230         }
231         //
232         //
233         if (multiplicity < fTriggerMultiplicity) continue;
234         Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
235       }    
236
237       Int_t nc = 0;
238       if (np == 0) continue;
239
240       Int_t i;
241       Int_t* newPos     = new Int_t[np];
242       Int_t* pSelected  = new Int_t[np];
243
244       for (i = 0; i<np; i++) {
245           newPos[i]    = i;
246           pSelected[i] = 0;
247       }
248       
249 //      First select parent particles
250
251       for (i = 0; i<np; i++) {
252           TParticle *iparticle = (TParticle *) fParticles.At(i);
253
254 // Is this a parent particle ?
255
256           if (Stable(iparticle)) continue;
257
258           Bool_t  selected             =  kTRUE;
259           Bool_t  hasSelectedDaughters =  kFALSE;
260           
261           kf = iparticle->GetPdgCode();
262           if (kf == 92 || kf == 99999) continue;
263           ks = iparticle->GetStatusCode();
264 // No initial state partons
265           if (ks==21) continue;
266             
267           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
268                                SelectFlavor(kf);
269
270           
271           hasSelectedDaughters = DaughtersSelection(iparticle);
272
273
274 // Put particle on the stack if it is either selected or 
275 // it is the mother of at least one seleted particle
276
277           if (selected || hasSelectedDaughters) {
278               nc++;
279               pSelected[i] = 1;
280           } // selected
281       } // particle loop parents
282
283 // Now select the final state particles
284
285
286       for (i=0; i<np; i++) {
287           TParticle *iparticle = (TParticle *) fParticles.At(i);
288
289 // Is this a final state particle ?
290
291           if (!Stable(iparticle)) continue;
292       
293           Bool_t  selected =  kTRUE;
294           kf = iparticle->GetPdgCode();
295           ks = iparticle->GetStatusCode();
296
297 // --------------------------------------------------------------------------
298 // Count spectator neutrons and protons (ks == 13, 14)
299           if(ks == 13 || ks == 14){
300               if(kf == kNeutron) fSpecn += 1;
301               if(kf == kProton)  fSpecp += 1;
302           }
303 // --------------------------------------------------------------------------
304
305           if (!fSelectAll) {
306               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
307               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
308           }
309
310 // Put particle on the stack if selected
311
312           if (selected) {
313               nc++;
314               pSelected[i] = 1;
315           } // selected
316       } // particle loop final state
317
318 // Write particles to stack
319
320       for (i = 0; i<np; i++) {
321           TParticle *  iparticle = (TParticle *) fParticles.At(i);
322           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
323           if (pSelected[i]) {
324               
325               kf   = iparticle->GetPdgCode();         
326               ks   = iparticle->GetStatusCode();              
327               
328               p[0] = iparticle->Px();
329               p[1] = iparticle->Py();
330               p[2] = iparticle->Pz();
331               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
332               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
333               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
334                     
335               tof = kconv*iparticle->T();
336               
337               imo = -1;
338               TParticle* mother = 0;
339               if (hasMother) {
340                   imo = iparticle->GetFirstMother();
341                   mother = (TParticle *) fParticles.At(imo);
342                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
343               } // if has mother   
344
345
346               
347               Bool_t tFlag = (fTrackIt && (ks == 1));
348               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
349               KeepTrack(nt);
350               newPos[i] = nt;
351           } // if selected
352       } // particle loop
353       delete[] newPos;
354       delete[] pSelected;
355       if (nc>0) {
356           jev += nc;
357           if (jev >= fNpart || fNpart == -1) {
358               break;
359           }
360       }
361   } // event loop
362   MakeHeader();
363   SetHighWaterMark(nt);
364 }
365
366 //______________________________________________________________________________
367 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
368 {
369 //
370 // Looks recursively if one of the daughters has been selected
371 //
372 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
373     Int_t imin = -1;
374     Int_t imax = -1;
375     Int_t i;
376     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
377     Bool_t selected = kFALSE;
378     if (hasDaughters) {
379         imin = iparticle->GetFirstDaughter();
380         imax = iparticle->GetLastDaughter();       
381         for (i = imin; i <= imax; i++){
382             TParticle *  jparticle = (TParticle *) fParticles.At(i);    
383             Int_t ip = jparticle->GetPdgCode();
384             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
385                 selected=kTRUE; break;
386             }
387             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
388         }
389     } else {
390         return kFALSE;
391     }
392     return selected;
393 }
394
395
396
397 //______________________________________________________________________________
398 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
399 {
400 // Select flavor of particle
401 // 0: all
402 // 4: charm and beauty
403 // 5: beauty
404     Bool_t res = 0;
405     
406     if (fFlavor == 0) {
407         res = kTRUE;
408     } else {
409         Int_t ifl = TMath::Abs(pid/100);
410         if (ifl > 10) ifl/=10;
411         res = (fFlavor == ifl);
412     }
413 //
414 //  This part if gamma writing is inhibited
415     if (fNoGammas) 
416         res = res && (pid != kGamma && pid != kPi0);
417 //
418     return res;
419 }
420
421 //______________________________________________________________________________
422 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
423 {
424 // Return true for a stable particle
425 //
426     
427 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
428     if (particle->GetStatusCode() == 1) return kTRUE;
429     else return kFALSE;
430
431 }
432
433 //______________________________________________________________________________
434 void AliGenDPMjet::MakeHeader()
435 {
436 // Builds the event header, to be called after each event
437     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
438     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
439     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
440     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
441     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
442                                                          fDPMjet->GetfIt());
443  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
444 // Bookkeeping for kinematic bias
445     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
446 // Event Vertex
447     header->SetPrimaryVertex(fVertex);
448     gAlice->SetGenEventHeader(header);    
449     AddHeader(header);
450 }
451
452 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
453 {
454     // Add header to container or runloader
455     if (fContainer) {
456         fContainer->AddHeader(header);
457     } else {
458         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
459     }
460 }
461
462
463 //______________________________________________________________________________
464 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
465 {
466 // Assignment operator
467     return *this;
468 }
469
470
471 void AliGenDPMjet::FinishRun()
472 {
473     // Print run statistics
474     fDPMjet->Dt_Dtuout();
475 }
476
477     
478
479 //______________________________________________________________________________