32624e68615d9b2c7996123fc376eda906eb275d
[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(0),
60      fGenImpPar(0.),
61      fProcess(kDpmMb)
62 {
63 // Constructor
64     fEnergyCMS = 5500.;
65     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
66 }
67
68
69 //______________________________________________________________________________
70 AliGenDPMjet::AliGenDPMjet(Int_t npart)
71     :AliGenMC(npart),
72      fBeamEn(2750.),
73      fMinImpactParam(0.),
74      fMaxImpactParam(5.),
75      fICentr(0),
76      fSelectAll(0),
77      fFlavor(0),
78      fTrials(0),
79      fSpectators(1),
80      fSpecn(0),
81      fSpecp(0),
82      fDPMjet(0),
83      fNoGammas(0),
84      fLHC(0),
85      fPi0Decay(0),
86      fGenImpPar(0.),
87      fProcess(kDpmMb)
88 {
89 // Default PbPb collisions at 5. 5 TeV
90 //
91     fEnergyCMS = 5500.;
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      fMinImpactParam(0.),
104      fMaxImpactParam(5.),
105      fICentr(0),
106      fSelectAll(0),
107      fFlavor(0),
108      fTrials(0),
109      fSpectators(1),
110      fSpecn(0),
111      fSpecp(0),
112      fDPMjet(0),
113      fNoGammas(0),
114      fLHC(0),
115      fPi0Decay(0),
116      fGenImpPar(0.),
117      fProcess(kDpmMb)
118 {
119     // Dummy copy constructor
120     fEnergyCMS = 5500.;
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 && TMath::Abs(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 || kf == 99999) 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 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
296               } // if has mother   
297
298               Bool_t tFlag = (fTrackIt && (ks == 1));
299               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
300               KeepTrack(nt);
301               newPos[i] = nt;
302           } // if selected
303       } // particle loop
304       delete[] newPos;
305       delete[] pSelected;
306       if (nc>0) {
307           jev += nc;
308           if (jev >= fNpart || fNpart == -1) {
309               break;
310           }
311       }
312   } // event loop
313   MakeHeader();
314   SetHighWaterMark(nt);
315 }
316
317 //______________________________________________________________________________
318 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
319 {
320 //
321 // Looks recursively if one of the daughters has been selected
322 //
323 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
324     Int_t imin = -1;
325     Int_t imax = -1;
326     Int_t i;
327     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
328     Bool_t selected = kFALSE;
329     if (hasDaughters) {
330         imin = iparticle->GetFirstDaughter();
331         imax = iparticle->GetLastDaughter();       
332         for (i = imin; i <= imax; i++){
333             TParticle *  jparticle = (TParticle *) fParticles.At(i);    
334             Int_t ip = jparticle->GetPdgCode();
335             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
336                 selected=kTRUE; break;
337             }
338             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
339         }
340     } else {
341         return kFALSE;
342     }
343     return selected;
344 }
345
346
347
348 //______________________________________________________________________________
349 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
350 {
351 // Select flavor of particle
352 // 0: all
353 // 4: charm and beauty
354 // 5: beauty
355     Bool_t res = 0;
356     
357     if (fFlavor == 0) {
358         res = kTRUE;
359     } else {
360         Int_t ifl = TMath::Abs(pid/100);
361         if (ifl > 10) ifl/=10;
362         res = (fFlavor == ifl);
363     }
364 //
365 //  This part if gamma writing is inhibited
366     if (fNoGammas) 
367         res = res && (pid != kGamma && pid != kPi0);
368 //
369     return res;
370 }
371
372 //______________________________________________________________________________
373 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
374 {
375 // Return true for a stable particle
376 //
377     
378 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
379     if (particle->GetStatusCode() == 1) return kTRUE;
380     else return kFALSE;
381
382 }
383
384 //______________________________________________________________________________
385 void AliGenDPMjet::MakeHeader()
386 {
387 // Builds the event header, to be called after each event
388     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
389     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
390     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
391     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
392     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
393                                                          fDPMjet->GetfIt());
394  ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
395 // Bookkeeping for kinematic bias
396     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
397 // Event Vertex
398     header->SetPrimaryVertex(fVertex);
399     gAlice->SetGenEventHeader(header);    
400     AddHeader(header);
401 }
402
403 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
404 {
405     // Add header to container or runloader
406     if (fContainer) {
407         fContainer->AddHeader(header);
408     } else {
409         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
410     }
411 }
412
413
414 //______________________________________________________________________________
415 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
416 {
417 // Assignment operator
418     return *this;
419 }
420
421
422 void AliGenDPMjet::FinishRun()
423 {
424     // Print run statistics
425     fDPMjet->Dt_Dtuout();
426 }
427
428     
429
430 //______________________________________________________________________________