Implement decaying via external decayer for long-lived resonances that are not know...
[u/mrichter/AliRoot.git] / TAmpt / AliGenAmpt.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 /* $Id$ */
17
18 // Generator using AMPT as an external generator
19
20 #include "AliGenAmpt.h"
21
22 #include <TClonesArray.h>
23 #include <TGraph.h>
24 #include <TAmpt.h>
25 #include <TLorentzVector.h>
26 #include <TPDGCode.h>
27 #include <TParticle.h>
28 #include <TVirtualMC.h>
29 #include <TParticlePDG.h>
30 #include "AliGenHijingEventHeader.h"
31 #define AliGenAmptEventHeader AliGenHijingEventHeader
32 #include "AliAmptRndm.h"
33 #include "AliLog.h"
34 #include "AliRun.h"
35 #include "AliDecayer.h"
36
37 ClassImp(AliGenAmpt)
38
39 AliGenAmpt::AliGenAmpt() 
40   : AliGenMC(),
41     fFrame("CMS"),
42     fMinImpactParam(0.),
43     fMaxImpactParam(5.),
44     fKeep(0),
45     fQuench(0),
46     fShadowing(1),
47     fDecaysOff(1),
48     fTrigger(0),     
49     fEvaluate(0),
50     fSelectAll(0),
51     fFlavor(0),
52     fKineBias(0.),
53     fTrials(0),
54     fXsection(0.),
55     fAmpt(0),
56     fPtHardMin(2.0),
57     fPtHardMax(-1),
58     fSpectators(1),
59     fDsigmaDb(0),
60     fDnDb(0),
61     fPtMinJet(-2.5),
62     fEtaMinJet(-20.),
63     fEtaMaxJet(+20.),
64     fPhiMinJet(0.),
65     fPhiMaxJet(TMath::TwoPi()),
66     fRadiation(3),
67     fSimpleJet(kFALSE),
68     fNoGammas(kFALSE),
69     fProjectileSpecn(0),
70     fProjectileSpecp(0),
71     fTargetSpecn(0),
72     fTargetSpecp(0),
73     fLHC(kFALSE),
74     fRandomPz(kFALSE),
75     fNoHeavyQuarks(kFALSE),
76     fEventTime(0.),
77     fIsoft(1),
78     fNtMax(150),
79     fIpop(1),
80     fXmu(3.2264),
81     fAlpha(1./3),
82     fStringA(0.5),
83     fStringB(0.9),
84     fHeader(new AliGenAmptEventHeader("Ampt")),
85     fDecay(kTRUE)
86 {
87   // Constructor
88   fEnergyCMS = 2760.;
89   AliAmptRndm::SetAmptRandom(GetRandom());
90 }
91
92 AliGenAmpt::AliGenAmpt(Int_t npart)
93   : AliGenMC(npart),
94     fFrame("CMS"),
95     fMinImpactParam(0.),
96     fMaxImpactParam(5.),
97     fKeep(0),
98     fQuench(0),
99     fShadowing(1),
100     fDecaysOff(1),
101     fTrigger(0),     
102     fEvaluate(0),
103     fSelectAll(0),
104     fFlavor(0),
105     fKineBias(0.),
106     fTrials(0),
107     fXsection(0.),
108     fAmpt(0),
109     fPtHardMin(2.0),
110     fPtHardMax(-1),
111     fSpectators(1),
112     fDsigmaDb(0),
113     fDnDb(0),
114     fPtMinJet(-2.5),
115     fEtaMinJet(-20.),
116     fEtaMaxJet(+20.),
117     fPhiMinJet(0.),
118     fPhiMaxJet(2. * TMath::Pi()),
119     fRadiation(3),
120     fSimpleJet(kFALSE),
121     fNoGammas(kFALSE),
122     fProjectileSpecn(0),
123     fProjectileSpecp(0),
124     fTargetSpecn(0),
125     fTargetSpecp(0),
126     fLHC(kFALSE),
127     fRandomPz(kFALSE),
128     fNoHeavyQuarks(kFALSE),
129     fEventTime(0.),
130     fIsoft(1),
131     fNtMax(150),
132     fIpop(1),
133     fXmu(3.2264),
134     fAlpha(1./3),
135     fStringA(0.5),
136     fStringB(0.9),
137     fHeader(new AliGenAmptEventHeader("Ampt")),
138     fDecay(kTRUE)
139 {
140   // Default PbPb collisions at 2.76 TeV
141
142   fEnergyCMS = 2760.;
143   fName = "Ampt";
144   fTitle= "Particle Generator using AMPT";
145   AliAmptRndm::SetAmptRandom(GetRandom());
146 }
147
148 AliGenAmpt::~AliGenAmpt()
149 {
150   // Destructor
151   if ( fDsigmaDb) delete fDsigmaDb;  
152   if ( fDnDb)     delete fDnDb;
153   if ( fHeader)   delete fHeader;
154 }
155
156 void AliGenAmpt::Init()
157 {
158   // Initialisation
159
160   fFrame.Resize(8);
161   fTarget.Resize(8);
162   fProjectile.Resize(8);
163
164   fAmpt = new TAmpt(fEnergyCMS, fFrame, fProjectile, fTarget, 
165                     fAProjectile, fZProjectile, fATarget, fZTarget, 
166                     fMinImpactParam, fMaxImpactParam);
167   SetMC(fAmpt);
168
169   fAmpt->SetIHPR2(2,  fRadiation);
170   fAmpt->SetIHPR2(3,  fTrigger);
171   fAmpt->SetIHPR2(6,  fShadowing);
172   fAmpt->SetIHPR2(12, fDecaysOff);    
173   fAmpt->SetIHPR2(21, fKeep);
174   fAmpt->SetHIPR1(8,  fPtHardMin);      
175   fAmpt->SetHIPR1(9,  fPtHardMax);      
176   fAmpt->SetHIPR1(10, fPtMinJet);       
177   fAmpt->SetHIPR1(50, fSimpleJet);
178
179   //  Quenching
180   //  fQuench = 0:  no quenching
181   //  fQuench = 1:  Hijing default
182   //  fQuench = 2:  new LHC  parameters for HIPR1(11) and HIPR1(14)
183   //  fQuench = 3:  new RHIC parameters for HIPR1(11) and HIPR1(14)
184   //  fQuench = 4:  new LHC  parameters with log(e) dependence
185   //  fQuench = 5:  new RHIC parameters with log(e) dependence
186   fAmpt->SetIHPR2(50, 0);
187   if (fQuench > 0) 
188     fAmpt->SetIHPR2(4,  1);
189   else
190     fAmpt->SetIHPR2(4,  0);
191
192   if (fQuench == 2) {
193     fAmpt->SetHIPR1(14, 1.1);
194     fAmpt->SetHIPR1(11, 3.7);
195   } else if (fQuench == 3) {
196     fAmpt->SetHIPR1(14, 0.20);
197     fAmpt->SetHIPR1(11, 2.5);
198   } else if (fQuench == 4) {
199     fAmpt->SetIHPR2(50, 1);
200     fAmpt->SetHIPR1(14, 4.*0.34);
201     fAmpt->SetHIPR1(11, 3.7);
202   } else if (fQuench == 5) {
203     fAmpt->SetIHPR2(50, 1);
204     fAmpt->SetHIPR1(14, 0.34);
205     fAmpt->SetHIPR1(11, 2.5);
206   }
207     
208   // Heavy quarks
209   if (fNoHeavyQuarks) {
210     fAmpt->SetIHPR2(49, 1);
211   } else {
212     fAmpt->SetIHPR2(49, 0);
213   }
214
215   // Ampt specific
216   fAmpt->SetIsoft(fIsoft);
217   fAmpt->SetNtMax(fNtMax);
218   fAmpt->SetIpop(fIpop);
219   fAmpt->SetXmu(fXmu);
220   fAmpt->SetAlpha(fAlpha);
221   fAmpt->SetStringFrag(fStringA, fStringB);
222
223   AliGenMC::Init();
224     
225   // Initialize Ampt  
226   fAmpt->Initialize();
227   if (fEvaluate) 
228     EvaluateCrossSections();
229 }
230
231 void AliGenAmpt::Generate()
232 {
233   // Generate one event
234
235   Float_t polar[3]    =   {0,0,0};
236   Float_t origin[3]   =   {0,0,0};
237   Float_t origin0[3]  =   {0,0,0};
238   Float_t p[3];
239   Float_t tof;
240
241   //  converts from mm/c to s
242   const Float_t kconv = 0.001/2.99792458e8;
243
244   Int_t nt  = 0;
245   Int_t jev = 0;
246   Int_t j, kf, ks, ksp, imo;
247   kf = 0;
248     
249   fTrials = 0;
250   for (j = 0;j < 3; j++) 
251     origin0[j] = fOrigin[j];
252
253   if(fVertexSmear == kPerEvent) {
254     Vertex();
255     for (j=0; j < 3; j++) 
256       origin0[j] = fVertex[j];
257   } 
258
259   Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
260
261   while(1) {
262     // Generate one event
263     Int_t fpemask = gSystem->GetFPEMask();
264     gSystem->SetFPEMask(0);
265     fAmpt->GenerateEvent();
266     gSystem->SetFPEMask(fpemask);
267     fTrials++;
268     fNprimaries = 0;
269     fAmpt->ImportParticles(&fParticles,"All");
270     Int_t np = fParticles.GetEntriesFast();
271     if (np == 0 ) 
272       continue;
273
274     if (fTrigger != kNoTrigger) {
275       if (!CheckTrigger()) 
276         continue;
277     }
278
279     AliDecayer *decayer = 0;
280     if (gMC)
281       decayer = gMC->GetDecayer();
282
283     if (decayer&&fDecay) {
284       TClonesArray arr("TParticle",100);
285       Int_t np2 = np;
286       for (Int_t i = 0; i < np; i++) {
287         TParticle *iparticle = (TParticle *)fParticles.At(i);
288         if (!Stable(iparticle)) 
289           continue;
290         kf        = iparticle->GetPdgCode();
291         if (kf==92)
292           continue;
293         if (0) { // this turned out to be too cumbersome!
294           if (kf!=331&&kf!=3114&&kf!=3114&&kf!=411&&kf!=-4122&&kf!=-3324&&kf!=-3312&&kf!=-3114&&
295               kf!=-311&&kf!=3214&&kf!=-3214&&kf!=-433&&kf!=413&&kf!=3122&&kf!=-3122&&kf!=-413&&
296               kf!=-421&&kf!=-423&&kf!=3324&&kf!=-313&&kf!=213&&kf!=-213&&kf!=3314&&kf!=3222&&
297               kf!=-3222&&kf!=3224&&kf!=-3224&&kf!=-4212&&kf!=4212&&kf!=433&&kf!=423&&kf!=-3322&&
298               kf!=3322&&kf!=-3314)
299             continue; //decay eta',Sigma*+,Sigma*-,D+,Lambda_c-,Xi*0_bar,Xi-_bar,Sigma*-,
300                       //      K0_bar,Sigma*0,Sigma*0_bar,D*_s-,D*+,Lambda0,Lambda0_bar,D*-
301                       //      D0_bar,D*0_bar,Xi*0,K*0_bar,rho+,rho-,Xi*-,Sigma-,
302                       //      Sigma+,Sigma*+,Sigma*-,Sigma_c-,Sigma_c+,D*_s+,D*0,Xi0_bar
303                       //      Xi0,Xi*+
304         } else { // really only decay particles if there are not known to Geant3
305           if (gMC->IdFromPDG(kf)>0)
306             continue;
307         }
308         if (0) { // defining the particle for Geant3 leads to a floating point exception.
309           TParticlePDG *pdg = iparticle->GetPDG(1);
310           //pdg->Print(); printf("%s\n",pdg->ParticleClass());
311           TString ptype(pdg->ParticleClass());
312           TMCParticleType mctype(kPTUndefined);
313           if (ptype=="Baryon" || ptype=="Meson")
314             mctype = kPTHadron;
315           gMC->DefineParticle(pdg->PdgCode(), pdg->GetName(), mctype, pdg->Mass(), pdg->Charge(), pdg->Lifetime(),
316                               ptype,pdg->Width(), (Int_t)pdg->Spin(), (Int_t)pdg->Parity(), 0, 
317                               (Int_t)pdg->Isospin(), 0, 0, 0, 0, pdg->Stable());
318           gMC->SetUserDecay(pdg->PdgCode());
319           continue;
320         }
321
322         TLorentzVector pmom(iparticle->Px(),iparticle->Py(),iparticle->Pz(),iparticle->Energy());
323         decayer->Decay(kf,&pmom);
324         decayer->ImportParticles(&arr);
325         Int_t ndecayed = arr.GetEntries();
326         if (ndecayed>1) {
327           if (np2+ndecayed>fParticles.GetSize())
328             fParticles.Expand(2*fParticles.GetSize());
329           //arr.Print();
330           iparticle->SetStatusCode(2);
331           iparticle->SetFirstDaughter(np2);
332           for (Int_t jj = 1; jj < ndecayed; jj++) {
333             TParticle *jp = (TParticle *)arr.At(jj);
334             if (jp->GetFirstMother()!=1)
335               continue;
336             TParticle *newp = new(fParticles[np2]) TParticle(jp->GetPdgCode(),
337                                                              1,
338                                                              i,
339                                                              -1,
340                                                              -1,
341                                                              -1,
342                                                              jp->Px(),jp->Py(),jp->Pz(),jp->Energy(),
343                                                              jp->Vx(),jp->Vy(),jp->Vz(),jp->T());
344             newp->SetUniqueID(999);
345             np2++;
346           }
347           iparticle->SetLastDaughter(np2-1);
348         }
349       }
350       np = fParticles.GetEntries();
351       if (np!=np2) {
352         AliError(Form("Something is fishy: %d %d\n", np,np2));
353       }
354     } else {
355       if (fDecay)
356         AliError("No decayer found, but fDecay==kTRUE!");
357     }
358
359     if (fLHC) 
360       Boost();
361       
362     Int_t nc = 0;
363     Int_t* newPos     = new Int_t[np];
364     Int_t* pSelected  = new Int_t[np];
365
366     for (Int_t i = 0; i < np; i++) {
367       newPos[i]    = i;
368       pSelected[i] = 0;
369     }
370       
371     // Get event vertex
372     //TParticle *  iparticle = (TParticle *) fParticles.At(0);
373     fVertex[0] = origin0[0];
374     fVertex[1] = origin0[1];    
375     fVertex[2] = origin0[2];
376       
377     // First select parent particles
378     for (Int_t i = 0; i < np; i++) {
379       TParticle *iparticle = (TParticle *) fParticles.At(i);
380
381     // Is this a parent particle ?
382       if (Stable(iparticle)) continue;
383       Bool_t  selected             =  kTRUE;
384       Bool_t  hasSelectedDaughters =  kFALSE;
385       kf        = iparticle->GetPdgCode();
386       ks        = iparticle->GetStatusCode();
387       if (kf == 92) 
388         continue;
389             
390       if (!fSelectAll) 
391         selected = KinematicSelection(iparticle, 0) && SelectFlavor(kf);
392       hasSelectedDaughters = DaughtersSelection(iparticle);
393
394       // Put particle on the stack if it is either selected or 
395       // it is the mother of at least one seleted particle
396       if (selected || hasSelectedDaughters) {
397         nc++;
398         pSelected[i] = 1;
399       } // selected
400     } // particle loop parents
401
402     // Now select the final state particles
403     fProjectileSpecn    = 0;  
404     fProjectileSpecp    = 0;
405     fTargetSpecn        = 0;  
406     fTargetSpecp        = 0;
407     for (Int_t i = 0; i<np; i++) {
408       TParticle *iparticle = (TParticle *) fParticles.At(i);
409       // Is this a final state particle ?
410       if (!Stable(iparticle)) 
411         continue;
412       Bool_t  selected             =  kTRUE;
413       kf        = iparticle->GetPdgCode();
414       if (kf == 92) 
415         continue;
416       ks        = iparticle->GetStatusCode();
417       ksp       = iparticle->GetUniqueID();
418           
419       // --------------------------------------------------------------------------
420       // Count spectator neutrons and protons
421       if(ksp == 0 || ksp == 1) {
422         if(kf == kNeutron) fProjectileSpecn += 1;
423         if(kf == kProton)  fProjectileSpecp += 1;
424       } else if(ksp == 10 || ksp == 11) {
425         if(kf == kNeutron) fTargetSpecn += 1;
426         if(kf == kProton)  fTargetSpecp += 1;
427       }
428       // --------------------------------------------------------------------------
429       if (!fSelectAll) {
430         selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf);
431         if (!fSpectators && selected) 
432           selected = (ksp != 0 && ksp != 1 && ksp != 10 && ksp != 11);
433       }
434
435      // Put particle on the stack if selected
436       if (selected) {
437         nc++;
438         pSelected[i] = 1;
439         if (1) printf("---> %d %d %d %s\n",i,nc,kf,iparticle->GetName());
440       } // selected
441     } // particle loop final state
442
443     //Time of the interactions
444     Float_t tInt = 0.;
445     if (fPileUpTimeWindow > 0.) 
446       tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
447
448     // Write particles to stack
449     for (Int_t i = 0; i<np; i++) {
450       TParticle *iparticle = (TParticle *) fParticles.At(i);
451       Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
452       Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
453       if (pSelected[i]) {
454         kf   = iparticle->GetPdgCode();
455         ks   = iparticle->GetStatusCode();
456         p[0] = iparticle->Px();
457         p[1] = iparticle->Py();
458         p[2] = iparticle->Pz() * sign;
459         origin[0] = origin0[0]+iparticle->Vx()/10;
460         origin[1] = origin0[1]+iparticle->Vy()/10;
461         origin[2] = origin0[2]+iparticle->Vz()/10;
462         fEventTime = 0.;
463               
464         if (TestBit(kVertexRange)) {
465           fEventTime = sign * origin0[2] / 2.99792458e10;
466           tof = kconv * iparticle->T() + fEventTime;
467         } else {
468           tof = kconv * iparticle->T();
469           fEventTime = tInt;
470           if (fPileUpTimeWindow > 0.) tof += tInt;
471         }
472         imo = -1;
473         TParticle* mother = 0;
474         if (hasMother) {
475           imo = iparticle->GetFirstMother();
476           mother = (TParticle *) fParticles.At(imo);
477           imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
478         } // if has mother   
479         Bool_t tFlag = (fTrackIt && !hasDaughter);
480         PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
481         fNprimaries++;
482         KeepTrack(nt);
483         newPos[i] = nt;
484       } // if selected
485     } // particle loop
486     delete[] newPos;
487     delete[] pSelected;
488       
489     AliInfo(Form("\n I've put %i particles on the stack \n",nc));
490     if (nc > 0) {
491       jev += nc;
492       if (jev >= fNpart || fNpart == -1) {
493         fKineBias = Float_t(fNpart)/Float_t(fTrials);
494         AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
495         break;
496       }
497     }
498   } // event loop
499   MakeHeader();
500   SetHighWaterMark(nt);
501 }
502
503 void AliGenAmpt::EvaluateCrossSections()
504 {
505   // Glauber Calculation of geometrical x-section
506
507   Float_t xTot       = 0.;          // barn
508   Float_t xTotHard   = 0.;          // barn 
509   Float_t xPart      = 0.;          // barn
510   Float_t xPartHard  = 0.;          // barn 
511   Float_t sigmaHard  = 0.1;         // mbarn
512   Float_t bMin       = 0.;
513   Float_t bMax       = fAmpt->GetHIPR1(34)+fAmpt->GetHIPR1(35);
514   const Float_t kdib = 0.2;
515   Int_t   kMax       = Int_t((bMax-bMin)/kdib)+1;
516
517   printf("\n Projectile Radius (fm): %f \n",fAmpt->GetHIPR1(34));
518   printf("\n Target     Radius (fm): %f \n",fAmpt->GetHIPR1(35));    
519
520   Int_t i;
521   Float_t oldvalue= 0.;
522   Float_t* b   = new Float_t[kMax]; memset(b,0,kMax*sizeof(Float_t));
523   Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
524   Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
525   for (i = 0; i < kMax; i++) {
526     Float_t xb  = bMin+i*kdib;
527     Float_t ov=fAmpt->Profile(xb);
528     Float_t gb  =  2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fAmpt->GetHINT1(12)*ov));
529     Float_t gbh =  2.*0.01*fAmpt->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
530     xTot+=gb;
531     xTotHard += gbh;
532     printf("profile %f %f %f\n", xb, ov, fAmpt->GetHINT1(12));
533         
534     if (xb > fMinImpactParam && xb < fMaxImpactParam) {
535       xPart += gb;
536       xPartHard += gbh;
537     }
538         
539     if ((oldvalue) && ((xTot-oldvalue)/oldvalue<0.0001)) 
540       break;
541     oldvalue = xTot;
542     printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
543     printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
544     if (i>0) {
545       si1[i] = gb/kdib;
546       si2[i] = gbh/gb;
547       b[i]  = xb;
548     }
549   }
550
551   printf("\n Total cross section (barn): %f \n",xTot);
552   printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
553   printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
554   printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
555
556   //  Store result as a graph
557   b[0] = 0;
558   si1[0] = 0;
559   si2[0]=si2[1];
560   delete fDsigmaDb;
561   fDsigmaDb  = new TGraph(i, b, si1);
562   delete fDnDb;
563   fDnDb      = new TGraph(i, b, si2);
564 }
565
566 Bool_t AliGenAmpt::DaughtersSelection(TParticle* iparticle)
567 {
568   // Looks recursively if one of the daughters has been selected
569   //printf("\n Consider daughters %d:",iparticle->GetPdgCode());
570   Int_t imin = -1;
571   Int_t imax = -1;
572   Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0);
573   Bool_t selected = kFALSE;
574   if (hasDaughters) {
575     imin = iparticle->GetFirstDaughter();
576     imax = iparticle->GetLastDaughter();       
577     for (Int_t i = imin; i <= imax; i++){
578       TParticle *  jparticle = (TParticle *) fParticles.At(i);  
579       Int_t ip = jparticle->GetPdgCode();
580       if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
581         selected=kTRUE; break;
582       }
583       if (DaughtersSelection(jparticle)) {selected=kTRUE; break; }
584     }
585   } else {
586     return kFALSE;
587   }
588   return selected;
589 }
590
591 Bool_t AliGenAmpt::SelectFlavor(Int_t pid)
592 {
593   // Select flavor of particle
594   // 0: all
595   // 4: charm and beauty
596   // 5: beauty
597   Bool_t res = 0;
598     
599   if (fFlavor == 0) {
600     res = kTRUE;
601   } else {
602     Int_t ifl = TMath::Abs(pid/100);
603     if (ifl > 10) ifl/=10;
604     res = (fFlavor == ifl);
605   }
606
607   //  This part if gamma writing is inhibited
608   if (fNoGammas) 
609     res = res && (pid != kGamma && pid != kPi0);
610
611   return res;
612 }
613
614 Bool_t AliGenAmpt::Stable(TParticle* particle) const
615 {
616   // Return true for a stable particle
617
618   if (!particle)
619     return kFALSE;
620   if (particle->GetFirstDaughter() < 0 )
621     return kTRUE;
622   return kFALSE;
623 }
624
625 void AliGenAmpt::MakeHeader()
626 {
627   // Fills the event header, to be called after each event
628
629   fHeader->SetNProduced(fNprimaries);
630   fHeader->SetImpactParameter(fAmpt->GetHINT1(19));
631   fHeader->SetTotalEnergy(fAmpt->GetEATT());
632   fHeader->SetHardScatters(fAmpt->GetJATT());
633   fHeader->SetParticipants(fAmpt->GetNP(), fAmpt->GetNT());
634   fHeader->SetCollisions(fAmpt->GetN0(),
635                         fAmpt->GetN01(),
636                         fAmpt->GetN10(),
637                         fAmpt->GetN11());
638   fHeader->SetSpectators(fProjectileSpecn, fProjectileSpecp,
639                         fTargetSpecn,fTargetSpecp);
640   fHeader->SetReactionPlaneAngle(fAmpt->GetHINT1(20));
641   //printf("Impact Parameter %13.3f \n", fAmpt->GetHINT1(19));
642
643   // 4-momentum vectors of the triggered jets.
644   // Before final state gluon radiation.
645   TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(21), 
646                                             fAmpt->GetHINT1(22),
647                                             fAmpt->GetHINT1(23),
648                                             fAmpt->GetHINT1(24));
649
650   TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(31), 
651                                             fAmpt->GetHINT1(32),
652                                             fAmpt->GetHINT1(33),
653                                             fAmpt->GetHINT1(34));
654   // After final state gluon radiation.
655   TLorentzVector* jet3 = new TLorentzVector(fAmpt->GetHINT1(26), 
656                                             fAmpt->GetHINT1(27),
657                                             fAmpt->GetHINT1(28),
658                                             fAmpt->GetHINT1(29));
659
660   TLorentzVector* jet4 = new TLorentzVector(fAmpt->GetHINT1(36), 
661                                             fAmpt->GetHINT1(37),
662                                             fAmpt->GetHINT1(38),
663                                             fAmpt->GetHINT1(39));
664   fHeader->SetJets(jet1, jet2, jet3, jet4);
665   // Bookkeeping for kinematic bias
666   fHeader->SetTrials(fTrials);
667   // Event Vertex
668   fHeader->SetPrimaryVertex(fVertex);
669   fHeader->SetInteractionTime(fEventTime);
670
671   fCollisionGeometry = fHeader;
672   AddHeader(fHeader);
673 }
674
675
676 Bool_t AliGenAmpt::CheckTrigger()
677 {
678   // Check the kinematic trigger condition
679
680   Bool_t   triggered = kFALSE;
681  
682   if (fTrigger == 1) {
683     // jet-jet Trigger  
684     TLorentzVector* jet1 = new TLorentzVector(fAmpt->GetHINT1(26), 
685                                               fAmpt->GetHINT1(27),
686                                               fAmpt->GetHINT1(28),
687                                               fAmpt->GetHINT1(29));
688         
689     TLorentzVector* jet2 = new TLorentzVector(fAmpt->GetHINT1(36), 
690                                               fAmpt->GetHINT1(37),
691                                               fAmpt->GetHINT1(38),
692                                               fAmpt->GetHINT1(39));
693     Double_t eta1      = jet1->Eta();
694     Double_t eta2      = jet2->Eta();
695     Double_t phi1      = jet1->Phi();
696     Double_t phi2      = jet2->Phi();
697     //printf("\n Trigger: %f %f %f %f", fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet);
698     if ( (eta1 < fEtaMaxJet && eta1 > fEtaMinJet &&  
699           phi1 < fPhiMaxJet && phi1 > fPhiMinJet) 
700          ||
701          (eta2 < fEtaMaxJet && eta2 > fEtaMinJet &&  
702           phi2 < fPhiMaxJet && phi2 > fPhiMinJet)
703       ) 
704       triggered = kTRUE;
705   } else if (fTrigger == 2) {
706   // Gamma Jet
707     Int_t np = fParticles.GetEntriesFast();
708     for (Int_t i = 0; i < np; i++) {
709       TParticle* part = (TParticle*) fParticles.At(i);
710       Int_t kf = part->GetPdgCode();
711       Int_t ksp = part->GetUniqueID();
712       if (kf == 22 && ksp == 40) {
713         Float_t phi = part->Phi();
714         Float_t eta = part->Eta();
715         if  (eta < fEtaMaxJet && 
716              eta > fEtaMinJet &&
717              phi < fPhiMaxJet && 
718              phi > fPhiMinJet) {
719           triggered = 1;
720           break;
721         } // check phi,eta within limits
722       } // direct gamma ? 
723     } // particle loop
724   } // fTrigger == 2
725   return triggered;
726 }