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