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