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