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