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