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