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 }