- Clean-up
[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
32 #include "AliGenDPMjet.h"
33 #include "AliGenDPMjetEventHeader.h"
34 #include "AliPythia.h"
35 #include "AliRun.h"
36 #include "AliDpmJetRndm.h"
37
38 ClassImp(AliGenDPMjet)
39
40 //______________________________________________________________________________
41 AliGenDPMjet::AliGenDPMjet()
42     :AliGenMC(), 
43      fBeamEn(2750.),
44      fEnergyCMS(5500.),
45      fMinImpactParam(0.),
46      fMaxImpactParam(5.),
47      fICentr(0),
48      fSelectAll(0),
49      fFlavor(0),
50      fTrials(0),
51      fSpectators(1),
52      fSpecn(0),
53      fSpecp(0),
54      fDPMjet(0),
55      fParticles(0),
56      fNoGammas(0),
57      fLHC(0),
58      fPi0Decay(0),
59      fGenImpPar(0.),
60      fProcess(kDpmMb)
61 {
62 // Constructor
63     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
64 }
65
66
67 //______________________________________________________________________________
68 AliGenDPMjet::AliGenDPMjet(Int_t npart)
69     :AliGenMC(npart),
70      fBeamEn(2750.),
71      fEnergyCMS(5500.),
72      fMinImpactParam(0.),
73      fMaxImpactParam(5.),
74      fICentr(0),
75      fSelectAll(0),
76      fFlavor(0),
77      fTrials(0),
78      fSpectators(1),
79      fSpecn(0),
80      fSpecp(0),
81      fDPMjet(0),
82      fParticles(new TClonesArray("TParticle",10000)),
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     fName = "DPMJET";
92     fTitle= "Particle Generator using DPMJET";
93     SetTarget();
94     SetProjectile();
95     fVertex.Set(3);
96     // Instance AliPythia
97     AliPythia::Instance(); 
98     AliDpmJetRndm::SetDpmJetRandom(GetRandom());
99 }
100
101 AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/)
102     :AliGenMC(),
103      fBeamEn(2750.),
104      fEnergyCMS(5500.),
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      fParticles(0),
116      fNoGammas(0),
117      fLHC(0),
118      fPi0Decay(0),
119      fGenImpPar(0.),
120      fProcess(kDpmMb)
121 {
122     // Dummy copy constructor
123 }
124
125 //______________________________________________________________________________
126 AliGenDPMjet::~AliGenDPMjet()
127 {
128 // Destructor
129     delete fParticles;
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 && fZProjectile == 1) fDPMjet->SetfIdp(1);
148     
149     fDPMjet->SetfFCentr(fICentr);  
150     fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); 
151     fDPMjet->SetPi0Decay(fPi0Decay);
152 //
153 //  Initialize DPMjet  
154 //    
155     fDPMjet->Initialize();
156 }
157
158
159 //______________________________________________________________________________
160 void AliGenDPMjet::Generate()
161 {
162 // Generate one event
163
164   Float_t polar[3]    =   {0,0,0};
165   Float_t origin[3]   =   {0,0,0};
166   Float_t origin0[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 //      Get event vertex
210       
211       fVertex[0] = origin0[0];
212       fVertex[1] = origin0[1];  
213       fVertex[2] = origin0[2];
214       
215 //      First select parent particles
216
217       for (i = 0; i<np; i++) {
218           TParticle *iparticle = (TParticle *) fParticles->At(i);
219
220 // Is this a parent particle ?
221
222           if (Stable(iparticle)) continue;
223
224           Bool_t  selected             =  kTRUE;
225           Bool_t  hasSelectedDaughters =  kFALSE;
226           
227           kf = iparticle->GetPdgCode();
228           if (kf == 92) continue;
229           ks = iparticle->GetStatusCode();
230 // No initial state partons
231           if (ks==21) continue;
232             
233           if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && 
234                                SelectFlavor(kf);
235           hasSelectedDaughters = DaughtersSelection(iparticle);
236
237 // Put particle on the stack if it is either selected or 
238 // it is the mother of at least one seleted particle
239
240           if (selected || hasSelectedDaughters) {
241               nc++;
242               pSelected[i] = 1;
243           } // selected
244       } // particle loop parents
245
246 // Now select the final state particles
247
248
249       for (i=0; i<np; i++) {
250           TParticle *iparticle = (TParticle *) fParticles->At(i);
251
252 // Is this a final state particle ?
253
254           if (!Stable(iparticle)) continue;
255       
256           Bool_t  selected =  kTRUE;
257           kf = iparticle->GetPdgCode();
258           ks = iparticle->GetStatusCode();
259
260 // --------------------------------------------------------------------------
261 // Count spectator neutrons and protons (ks == 13, 14)
262           if(ks == 13 || ks == 14){
263               if(kf == kNeutron) fSpecn += 1;
264               if(kf == kProton)  fSpecp += 1;
265           }
266 // --------------------------------------------------------------------------
267
268           if (!fSelectAll) {
269               selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
270               if (!fSpectators && selected) selected = (ks == 13 || ks == 14);
271           }
272
273 // Put particle on the stack if selected
274
275           if (selected) {
276               nc++;
277               pSelected[i] = 1;
278           } // selected
279       } // particle loop final state
280
281 // Write particles to stack
282
283       for (i = 0; i<np; i++) {
284           TParticle *  iparticle = (TParticle *) fParticles->At(i);
285           Bool_t  hasMother   = (iparticle->GetFirstMother()>=0);
286           if (pSelected[i]) {
287               
288               kf   = iparticle->GetPdgCode();         
289               ks   = iparticle->GetStatusCode();              
290               
291               p[0] = iparticle->Px();
292               p[1] = iparticle->Py();
293               p[2] = iparticle->Pz();
294               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
295               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
296               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
297                     
298               tof = kconv*iparticle->T();
299               
300               imo = -1;
301               TParticle* mother = 0;
302               if (hasMother) {
303                   imo = iparticle->GetFirstMother();
304                   mother = (TParticle *) fParticles->At(imo);
305                   imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1;
306               } // if has mother   
307
308               Bool_t tFlag = (fTrackIt && (ks == 1));
309               
310               PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
311               KeepTrack(nt);
312               newPos[i] = nt;
313           } // if selected
314       } // particle loop
315       delete[] newPos;
316       delete[] pSelected;
317       if (nc>0) {
318           jev += nc;
319           if (jev >= fNpart || fNpart == -1) {
320               break;
321           }
322       }
323   } // event loop
324   MakeHeader();
325   SetHighWaterMark(nt);
326 }
327
328 //______________________________________________________________________________
329 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
330 {
331 //
332 // Looks recursively if one of the daughters has been selected
333 //
334 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
335     Int_t imin = -1;
336     Int_t imax = -1;
337     Int_t i;
338     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
339     Bool_t selected = kFALSE;
340     if (hasDaughters) {
341         imin = iparticle->GetFirstDaughter();
342         imax = iparticle->GetLastDaughter();       
343         for (i = imin; i <= imax; i++){
344             TParticle *  jparticle = (TParticle *) fParticles->At(i);   
345             Int_t ip = jparticle->GetPdgCode();
346             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
347                 selected=kTRUE; break;
348             }
349             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
350         }
351     } else {
352         return kFALSE;
353     }
354     return selected;
355 }
356
357
358
359 //______________________________________________________________________________
360 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
361 {
362 // Select flavor of particle
363 // 0: all
364 // 4: charm and beauty
365 // 5: beauty
366     Bool_t res = 0;
367     
368     if (fFlavor == 0) {
369         res = kTRUE;
370     } else {
371         Int_t ifl = TMath::Abs(pid/100);
372         if (ifl > 10) ifl/=10;
373         res = (fFlavor == ifl);
374     }
375 //
376 //  This part if gamma writing is inhibited
377     if (fNoGammas) 
378         res = res && (pid != kGamma && pid != kPi0);
379 //
380     return res;
381 }
382
383 //______________________________________________________________________________
384 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
385 {
386 // Return true for a stable particle
387 //
388     
389 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
390     if (particle->GetStatusCode() == 1) return kTRUE;
391     else return kFALSE;
392
393 }
394
395 //______________________________________________________________________________
396 void AliGenDPMjet::MakeHeader()
397 {
398 // Builds the event header, to be called after each event
399     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
400     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
401     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
402     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
403     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), 
404                                                          fDPMjet->GetfIt());
405
406 // Bookkeeping for kinematic bias
407     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
408 // Event Vertex
409     header->SetPrimaryVertex(fVertex);
410     gAlice->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
424 //______________________________________________________________________________