Finish run with print-out of run statistics added.
[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 p[3];
164   Float_t tof;
165
166 //  converts from mm/c to s
167   const Float_t kconv = 0.001/2.999792458e8;
168   Int_t nt  = 0;
169   Int_t jev = 0;
170   Int_t kf, ks, imo;
171   kf = 0;
172   fTrials = 0;
173   //  Set collision vertex position 
174   if (fVertexSmear == kPerEvent) Vertex();
175   
176   while(1)
177   {
178 //    Generate one event
179 // --------------------------------------------------------------------------
180       fSpecn = 0;  
181       fSpecp = 0;
182 // --------------------------------------------------------------------------
183       fDPMjet->GenerateEvent();
184       fTrials++;
185
186       fDPMjet->ImportParticles(&fParticles,"All");      
187       if (fLHC) Boost();
188
189       // Temporaneo
190       fGenImpPar = fDPMjet->GetBImpac();
191       
192       Int_t np = fParticles.GetEntriesFast();
193       printf("\n **************************************************%d\n",np);
194       Int_t nc=0;
195       if (np==0) continue;
196       Int_t i;
197       Int_t* newPos     = new Int_t[np];
198       Int_t* pSelected  = new Int_t[np];
199
200       for (i = 0; i<np; i++) {
201           newPos[i]    = i;
202           pSelected[i] = 0;
203       }
204       
205 //      First select parent particles
206
207       for (i = 0; i<np; i++) {
208           TParticle *iparticle = (TParticle *) fParticles.At(i);
209
210 // Is this a parent particle ?
211
212           if (Stable(iparticle)) continue;
213
214           Bool_t  selected             =  kTRUE;
215           Bool_t  hasSelectedDaughters =  kFALSE;
216           
217           kf = iparticle->GetPdgCode();
218           if (kf == 92) continue;
219           ks = iparticle->GetStatusCode();
220 // No initial state partons
221           if (ks==21) continue;
222             
223           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
224                                SelectFlavor(kf);
225           hasSelectedDaughters = DaughtersSelection(iparticle);
226
227 // Put particle on the stack if it is either selected or 
228 // it is the mother of at least one seleted particle
229
230           if (selected || hasSelectedDaughters) {
231               nc++;
232               pSelected[i] = 1;
233           } // selected
234       } // particle loop parents
235
236 // Now select the final state particles
237
238
239       for (i=0; i<np; i++) {
240           TParticle *iparticle = (TParticle *) fParticles.At(i);
241
242 // Is this a final state particle ?
243
244           if (!Stable(iparticle)) continue;
245       
246           Bool_t  selected =  kTRUE;
247           kf = iparticle->GetPdgCode();
248           ks = iparticle->GetStatusCode();
249
250 // --------------------------------------------------------------------------
251 // Count spectator neutrons and protons (ks == 13, 14)
252           if(ks == 13 || ks == 14){
253               if(kf == kNeutron) fSpecn += 1;
254               if(kf == kProton)  fSpecp += 1;
255           }
256 // --------------------------------------------------------------------------
257
258           if (!fSelectAll) {
259               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
260               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
261           }
262
263 // Put particle on the stack if selected
264
265           if (selected) {
266               nc++;
267               pSelected[i] = 1;
268           } // selected
269       } // particle loop final state
270
271 // Write particles to stack
272
273       for (i = 0; i<np; i++) {
274           TParticle *  iparticle = (TParticle *) fParticles.At(i);
275           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
276           if (pSelected[i]) {
277               
278               kf   = iparticle->GetPdgCode();         
279               ks   = iparticle->GetStatusCode();              
280               
281               p[0] = iparticle->Px();
282               p[1] = iparticle->Py();
283               p[2] = iparticle->Pz();
284               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
285               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
286               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
287                     
288               tof = kconv*iparticle->T();
289               
290               imo = -1;
291               TParticle* mother = 0;
292               if (hasMother) {
293                   imo = iparticle->GetFirstMother();
294                   mother = (TParticle *) fParticles.At(imo);
295                   imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
296               } // if has mother   
297
298               Bool_t tFlag = (fTrackIt && (ks == 1));
299               
300               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
301               KeepTrack(nt);
302               newPos[i] = nt;
303           } // if selected
304       } // particle loop
305       delete[] newPos;
306       delete[] pSelected;
307       if (nc>0) {
308           jev += nc;
309           if (jev >= fNpart || fNpart == -1) {
310               break;
311           }
312       }
313   } // event loop
314   MakeHeader();
315   SetHighWaterMark(nt);
316 }
317
318 //______________________________________________________________________________
319 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
320 {
321 //
322 // Looks recursively if one of the daughters has been selected
323 //
324 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
325     Int_t imin = -1;
326     Int_t imax = -1;
327     Int_t i;
328     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
329     Bool_t selected = kFALSE;
330     if (hasDaughters) {
331         imin = iparticle->GetFirstDaughter();
332         imax = iparticle->GetLastDaughter();       
333         for (i = imin; i <= imax; i++){
334             TParticle *  jparticle = (TParticle *) fParticles.At(i);    
335             Int_t ip = jparticle->GetPdgCode();
336             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
337                 selected=kTRUE; break;
338             }
339             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
340         }
341     } else {
342         return kFALSE;
343     }
344     return selected;
345 }
346
347
348
349 //______________________________________________________________________________
350 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
351 {
352 // Select flavor of particle
353 // 0: all
354 // 4: charm and beauty
355 // 5: beauty
356     Bool_t res = 0;
357     
358     if (fFlavor == 0) {
359         res = kTRUE;
360     } else {
361         Int_t ifl = TMath::Abs(pid/100);
362         if (ifl > 10) ifl/=10;
363         res = (fFlavor == ifl);
364     }
365 //
366 //  This part if gamma writing is inhibited
367     if (fNoGammas) 
368         res = res && (pid != kGamma && pid != kPi0);
369 //
370     return res;
371 }
372
373 //______________________________________________________________________________
374 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
375 {
376 // Return true for a stable particle
377 //
378     
379 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
380     if (particle->GetStatusCode() == 1) return kTRUE;
381     else return kFALSE;
382
383 }
384
385 //______________________________________________________________________________
386 void AliGenDPMjet::MakeHeader()
387 {
388 // Builds the event header, to be called after each event
389     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
390     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
391     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
392     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
393     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
394                                                          fDPMjet->GetfIt());
395  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
396 // Bookkeeping for kinematic bias
397     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
398 // Event Vertex
399     header->SetPrimaryVertex(fVertex);
400     gAlice->SetGenEventHeader(header);    
401     AddHeader(header);
402 }
403
404 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
405 {
406     // Add header to container or runloader
407     if (fContainer) {
408         fContainer->AddHeader(header);
409     } else {
410         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
411     }
412 }
413
414
415 //______________________________________________________________________________
416 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
417 {
418 // Assignment operator
419     return *this;
420 }
421
422
423 void AliGenDPMjet::FinishRun()
424 {
425     // Print run statistics
426     fDPMjet->Dt_Dtuout();
427 }
428
429     
430
431 //______________________________________________________________________________