Introducing the interaction time into the aliroot generators. In case of gaussian...
[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   Double_t polar[3]    =   {0,0,0};
175   Double_t origin[3]   =   {0,0,0};
176   Double_t p[4]        =   {0};
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               p[3] = iparticle->Energy();
334               origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
335               origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
336               origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
337                     
338               tof = fTime + kconv*iparticle->T();
339               
340               imo = -1;
341               TParticle* mother = 0;
342               if (hasMother) {
343                   imo = iparticle->GetFirstMother();
344                   mother = (TParticle *) fParticles.At(imo);
345                   imo = (mother->GetPdgCode() != 92 && mother->GetPdgCode() != 99999) ? newPos[imo] : -1;
346               } // if has mother   
347
348
349               
350               Bool_t tFlag = (fTrackIt && (ks == 1));
351               PushTrack(tFlag, imo, kf, 
352                         p[0], p[1], p[2], p[3], 
353                         origin[0], origin[1], origin[2], tof,
354                         polar[0], polar[1], polar[2],
355                         kPNoProcess, nt, 1., ks);
356               KeepTrack(nt);
357               newPos[i] = nt;
358           } // if selected
359       } // particle loop
360       delete[] newPos;
361       delete[] pSelected;
362       if (nc>0) {
363           jev += nc;
364           if (jev >= fNpart || fNpart == -1) {
365               break;
366           }
367       }
368   } // event loop
369   MakeHeader();
370   SetHighWaterMark(nt);
371 }
372
373 //______________________________________________________________________________
374 Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle)
375 {
376 //
377 // Looks recursively if one of the daughters has been selected
378 //
379 //    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
380     Int_t imin = -1;
381     Int_t imax = -1;
382     Int_t i;
383     Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
384     Bool_t selected = kFALSE;
385     if (hasDaughters) {
386         imin = iparticle->GetFirstDaughter();
387         imax = iparticle->GetLastDaughter();       
388         for (i = imin; i <= imax; i++){
389             TParticle *  jparticle = (TParticle *) fParticles.At(i);    
390             Int_t ip = jparticle->GetPdgCode();
391             if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
392                 selected=kTRUE; break;
393             }
394             if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
395         }
396     } else {
397         return kFALSE;
398     }
399     return selected;
400 }
401
402
403
404 //______________________________________________________________________________
405 Bool_t AliGenDPMjet::SelectFlavor(Int_t pid)
406 {
407 // Select flavor of particle
408 // 0: all
409 // 4: charm and beauty
410 // 5: beauty
411     Bool_t res = 0;
412     
413     if (fFlavor == 0) {
414         res = kTRUE;
415     } else {
416         Int_t ifl = TMath::Abs(pid/100);
417         if (ifl > 10) ifl/=10;
418         res = (fFlavor == ifl);
419     }
420 //
421 //  This part if gamma writing is inhibited
422     if (fNoGammas) 
423         res = res && (pid != kGamma && pid != kPi0);
424 //
425     return res;
426 }
427
428 //______________________________________________________________________________
429 Bool_t AliGenDPMjet::Stable(TParticle*  particle)
430 {
431 // Return true for a stable particle
432 //
433     
434 //    if (particle->GetFirstDaughter() < 0 ) return kTRUE;
435     if (particle->GetStatusCode() == 1) return kTRUE;
436     else return kFALSE;
437
438 }
439
440 //______________________________________________________________________________
441 void AliGenDPMjet::MakeHeader()
442 {
443   printf("MakeHeader %13.3f \n", fDPMjet->GetBImpac());
444 // Builds the event header, to be called after each event
445     AliGenEventHeader* header = new AliGenDPMjetEventHeader("DPMJET");
446     ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc());
447     ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac());
448     ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy());
449     ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetProjParticipants(), 
450                                                          fDPMjet->GetTargParticipants());
451      ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode());
452     // Bookkeeping for kinematic bias
453     ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials);
454     // Event Vertex
455     header->SetPrimaryVertex(fVertex);
456     header->SetInteractionTime(fTime);
457     gAlice->SetGenEventHeader(header);    
458     AddHeader(header);
459 }
460
461 void AliGenDPMjet::AddHeader(AliGenEventHeader* header)
462 {
463     // Add header to container or runloader
464     if (fContainer) {
465         fContainer->AddHeader(header);
466     } else {
467         AliRunLoader::Instance()->GetHeader()->SetGenEventHeader(header);
468     }
469 }
470
471
472 //______________________________________________________________________________
473 AliGenDPMjet& AliGenDPMjet::operator=(const  AliGenDPMjet& /*rhs*/)
474 {
475 // Assignment operator
476     return *this;
477 }
478
479
480 void AliGenDPMjet::FinishRun()
481 {
482     // Print run statistics
483     fDPMjet->Dt_Dtuout();
484 }
485
486     
487
488 //______________________________________________________________________________