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