Empirical cut to increase robustness
[u/mrichter/AliRoot.git] / TAmpt / AliGenAmpt.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 // Generator using AMPT as an external generator
19
20 #include "AliGenAmpt.h"
21
22 #include <TClonesArray.h>
23 #include <TGraph.h>
24 #include <TAmpt.h>
25 #include <TLorentzVector.h>
26 #include <TPDGCode.h>
27 #include <TParticle.h>
28
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]; memset(b,0,kMax*sizeof(Float_t));
435   Float_t* si1 = new Float_t[kMax]; memset(si1,0,kMax*sizeof(Float_t));
436   Float_t* si2 = new Float_t[kMax]; memset(si2,0,kMax*sizeof(Float_t));
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 }