Correcting warnings
[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      fEnergyCMS(5500.),
48      fMinImpactParam(0.),
49      fMaxImpactParam(5.),
50      fICentr(0),
51      fSelectAll(0),
52      fFlavor(0),
53      fTrials(0),
54      fSpectators(1),
55      fSpecn(0),
56      fSpecp(0),
57      fDPMjet(0),
58      fParticles(0),
59      fNoGammas(0),
60      fLHC(0),
61      fPi0Decay(0),
62      fGenImpPar(0.),
63      fProcess(kDpmMb)
64 {
65 // Constructor
66     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
67 }
68
69
70 //______________________________________________________________________________
71 AliGenDPMjet::AliGenDPMjet(Int_t npart)
72     :AliGenMC(npart),
73      fBeamEn(2750.),
74      fEnergyCMS(5500.),
75      fMinImpactParam(0.),
76      fMaxImpactParam(5.),
77      fICentr(0),
78      fSelectAll(0),
79      fFlavor(0),
80      fTrials(0),
81      fSpectators(1),
82      fSpecn(0),
83      fSpecp(0),
84      fDPMjet(0),
85      fParticles(new TClonesArray("TParticle",10000)),
86      fNoGammas(0),
87      fLHC(0),
88      fPi0Decay(0),
89      fGenImpPar(0.),
90      fProcess(kDpmMb)
91 {
92 // Default PbPb collisions at 5. 5 TeV
93 //
94     fName = "DPMJET";
95     fTitle= "Particle Generator using DPMJET";
96     SetTarget();
97     SetProjectile();
98     fVertex.Set(3);
99     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
100 }
101
102 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
103     :AliGenMC(),
104      fBeamEn(2750.),
105      fEnergyCMS(5500.),
106      fMinImpactParam(0.),
107      fMaxImpactParam(5.),
108      fICentr(0),
109      fSelectAll(0),
110      fFlavor(0),
111      fTrials(0),
112      fSpectators(1),
113      fSpecn(0),
114      fSpecp(0),
115      fDPMjet(0),
116      fParticles(0),
117      fNoGammas(0),
118      fLHC(0),
119      fPi0Decay(0),
120      fGenImpPar(0.),
121      fProcess(kDpmMb)
122 {
123     // Dummy copy constructor
124 }
125
126 //______________________________________________________________________________
127 AliGenDPMjet::~AliGenDPMjet()
128 {
129 // Destructor
130     delete fParticles;
131 }
132 //______________________________________________________________________________
133 void AliGenDPMjet::Init()
134 {
135 // Initialization
136     
137     SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, 
138                       fBeamEn,fEnergyCMS));
139
140     fDPMjet=(TDPMjet*) fMCEvGen;
141     //
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);
149     
150     fDPMjet->SetfFCentr(fICentr);  
151     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
152     fDPMjet->SetPi0Decay(fPi0Decay);
153 //
154 //  Initialize DPMjet  
155 //    
156     fDPMjet->Initialize();
157 }
158
159
160 //______________________________________________________________________________
161 void AliGenDPMjet::Generate()
162 {
163 // Generate one event
164
165   Float_t polar[3]    =   {0,0,0};
166   Float_t origin[3]   =   {0,0,0};
167   Float_t origin0[3]  =   {0,0,0};
168   Float_t p[3];
169   Float_t tof;
170
171 //  converts from mm/c to s
172   const Float_t kconv = 0.001/2.999792458e8;
173   Int_t nt  = 0;
174   Int_t jev = 0;
175   Int_t kf, ks, imo;
176   kf = 0;
177   fTrials = 0;
178   //  Set collision vertex position 
179   if (fVertexSmear == kPerEvent) Vertex();
180   
181   while(1)
182   {
183 //    Generate one event
184 // --------------------------------------------------------------------------
185       fSpecn = 0;  
186       fSpecp = 0;
187 // --------------------------------------------------------------------------
188       fDPMjet->GenerateEvent();
189       fTrials++;
190
191       fDPMjet->ImportParticles(fParticles,"All");      
192       if (fLHC) Boost();
193
194       // Temporaneo
195       fGenImpPar = fDPMjet->GetBImpac();
196       
197       Int_t np = fParticles->GetEntriesFast();
198       printf("\n **************************************************%d\n",np);
199       Int_t nc=0;
200       if (np==0) continue;
201       Int_t i;
202       Int_t* newPos     = new Int_t[np];
203       Int_t* pSelected  = new Int_t[np];
204
205       for (i = 0; i<np; i++) {
206           newPos[i]    = i;
207           pSelected[i] = 0;
208       }
209       
210 //      Get event vertex
211       
212       fVertex[0] = origin0[0];
213       fVertex[1] = origin0[1];  
214       fVertex[2] = origin0[2];
215       
216 //      First select parent particles
217
218       for (i = 0; i<np; i++) {
219           TParticle *iparticle = (TParticle *) fParticles->At(i);
220
221 // Is this a parent particle ?
222
223           if (Stable(iparticle)) continue;
224
225           Bool_t  selected             =  kTRUE;
226           Bool_t  hasSelectedDaughters =  kFALSE;
227           
228           kf = iparticle->GetPdgCode();
229           if (kf == 92) continue;
230           ks = iparticle->GetStatusCode();
231 // No initial state partons
232           if (ks==21) continue;
233             
234           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
235                                SelectFlavor(kf);
236           hasSelectedDaughters = DaughtersSelection(iparticle);
237
238 // Put particle on the stack if it is either selected or 
239 // it is the mother of at least one seleted particle
240
241           if (selected || hasSelectedDaughters) {
242               nc++;
243               pSelected[i] = 1;
244           } // selected
245       } // particle loop parents
246
247 // Now select the final state particles
248
249
250       for (i=0; i<np; i++) {
251           TParticle *iparticle = (TParticle *) fParticles->At(i);
252
253 // Is this a final state particle ?
254
255           if (!Stable(iparticle)) continue;
256       
257           Bool_t  selected =  kTRUE;
258           kf = iparticle->GetPdgCode();
259           ks = iparticle->GetStatusCode();
260
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;
266           }
267 // --------------------------------------------------------------------------
268
269           if (!fSelectAll) {
270               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
271               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
272           }
273
274 // Put particle on the stack if selected
275
276           if (selected) {
277               nc++;
278               pSelected[i] = 1;
279           } // selected
280       } // particle loop final state
281
282 // Write particles to stack
283
284       for (i = 0; i<np; i++) {
285           TParticle *  iparticle = (TParticle *) fParticles->At(i);
286           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
287           if (pSelected[i]) {
288               
289               kf   = iparticle->GetPdgCode();         
290               ks   = iparticle->GetStatusCode();              
291               
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]
298                     
299               tof = kconv*iparticle->T();
300               
301               imo = -1;
302               TParticle* mother = 0;
303               if (hasMother) {
304                   imo = iparticle->GetFirstMother();
305                   mother = (TParticle *) fParticles->At(imo);
306                   imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
307               } // if has mother   
308
309               Bool_t tFlag = (fTrackIt && (ks == 1));
310               
311               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
312               KeepTrack(nt);
313               newPos[i] = nt;
314           } // if selected
315       } // particle loop
316       delete[] newPos;
317       delete[] pSelected;
318       if (nc>0) {
319           jev += nc;
320           if (jev >= fNpart || fNpart == -1) {
321               break;
322           }
323       }
324   } // event loop
325   MakeHeader();
326   SetHighWaterMark(nt);
327 }
328
329 //______________________________________________________________________________
330 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
331 {
332 //
333 // Looks recursively if one of the daughters has been selected
334 //
335 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
336     Int_t imin = -1;
337     Int_t imax = -1;
338     Int_t i;
339     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
340     Bool_t selected = kFALSE;
341     if (hasDaughters) {
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;
349             }
350             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
351         }
352     } else {
353         return kFALSE;
354     }
355     return selected;
356 }
357
358
359
360 //______________________________________________________________________________
361 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
362 {
363 // Select flavor of particle
364 // 0: all
365 // 4: charm and beauty
366 // 5: beauty
367     Bool_t res = 0;
368     
369     if (fFlavor == 0) {
370         res = kTRUE;
371     } else {
372         Int_t ifl = TMath::Abs(pid/100);
373         if (ifl > 10) ifl/=10;
374         res = (fFlavor == ifl);
375     }
376 //
377 //  This part if gamma writing is inhibited
378     if (fNoGammas) 
379         res = res && (pid != kGamma && pid != kPi0);
380 //
381     return res;
382 }
383
384 //______________________________________________________________________________
385 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
386 {
387 // Return true for a stable particle
388 //
389     
390 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
391     if (particle->GetStatusCode() == 1) return kTRUE;
392     else return kFALSE;
393
394 }
395
396 //______________________________________________________________________________
397 void AliGenDPMjet::MakeHeader()
398 {
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(), 
405                                                          fDPMjet->GetfIt());
406  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
407 // Bookkeeping for kinematic bias
408     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
409 // Event Vertex
410     header->SetPrimaryVertex(fVertex);
411     gAlice->SetGenEventHeader(header);    
412  AddHeader(header);
413 }
414
415 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
416 {
417     // Add header to container or runloader
418     if (fContainer) {
419         fContainer->AddHeader(header);
420     } else {
421         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
422     }
423 }
424
425
426 //______________________________________________________________________________
427 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
428 {
429 // Assignment operator
430     return *this;
431 }
432
433
434
435 //______________________________________________________________________________