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