68d59d5b5429026e636cfc5c0b62cc27cbee266d
[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     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 }