TAmpt/AliGenAmpt.cxx

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