756273c4befe7d804cf58a25eb7a47aba28c75e6
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.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 //
19 // Generator using the TPythia interface (via AliPythia)
20 // to generate pp collisions.
21 // Using SetNuclei() also nuclear modifications to the structure functions
22 // can be taken into account. This makes, of course, only sense for the
23 // generation of the products of hard processes (heavy flavor, jets ...)
24 //
25 // andreas.morsch@cern.ch
26 //
27
28 #include <TClonesArray.h>
29 #include <TDatabasePDG.h>
30 #include <TParticle.h>
31 #include <TPDGCode.h>
32 #include <TObjArray.h>
33 #include <TSystem.h>
34 #include <TTree.h>
35 #include "AliConst.h"
36 #include "AliDecayerPythia.h"
37 #include "AliGenPythia.h"
38 #include "AliFastGlauber.h"
39 #include "AliHeader.h"
40 #include "AliGenPythiaEventHeader.h"
41 #include "AliPythia.h"
42 #include "AliPythiaRndm.h"
43 #include "AliRun.h"
44 #include "AliStack.h"
45 #include "AliRunLoader.h"
46 #include "AliMC.h"
47 #include "AliLog.h"
48 #include "PyquenCommon.h"
49 #include "AliLog.h"
50
51 ClassImp(AliGenPythia)
52
53
54 AliGenPythia::AliGenPythia():
55     AliGenMC(),
56     fProcess(kPyCharm),          
57     fItune(-1),
58     fStrucFunc(kCTEQ5L), 
59     fKineBias(0.),
60     fTrials(0),
61     fTrialsRun(0),
62     fQ(0.),
63     fX1(0.),
64     fX2(0.),
65     fEventTime(0.),
66     fInteractionRate(0.),
67     fTimeWindow(0.),
68     fCurSubEvent(0),
69     fEventsTime(0),
70     fNev(0),
71     fFlavorSelect(0),
72     fXsection(0.),
73     fPythia(0),
74     fPtHardMin(0.),
75     fPtHardMax(1.e4),
76     fYHardMin(-1.e10),
77     fYHardMax(1.e10),
78     fGinit(1),
79     fGfinal(1),
80     fHadronisation(1),
81     fPatchOmegaDalitz(0), 
82     fNpartons(0),
83     fReadFromFile(0),
84     fQuench(0),
85     fQhat(0.),
86     fLength(0.),
87     fpyquenT(1.),
88     fpyquenTau(0.1),
89     fpyquenNf(0),
90     fpyquenEloss(0),
91     fpyquenAngle(0), 
92     fImpact(0.),
93     fPtKick(1.),
94     fFullEvent(kTRUE),
95     fDecayer(new AliDecayerPythia()),
96     fDebugEventFirst(-1),
97     fDebugEventLast(-1),
98     fEtMinJet(0.),      
99     fEtMaxJet(1.e4),      
100     fEtaMinJet(-20.),     
101     fEtaMaxJet(20.),     
102     fPhiMinJet(0.),     
103     fPhiMaxJet(2.* TMath::Pi()),     
104     fJetReconstruction(kCell),
105     fEtaMinGamma(-20.),      
106     fEtaMaxGamma(20.),      
107     fPhiMinGamma(0.),      
108     fPhiMaxGamma(2. * TMath::Pi()),      
109     fPycellEtaMax(2.),     
110     fPycellNEta(274),       
111     fPycellNPhi(432),       
112     fPycellThreshold(0.),  
113     fPycellEtSeed(4.),     
114     fPycellMinEtJet(10.),  
115     fPycellMaxRadius(1.), 
116     fStackFillOpt(kFlavorSelection),   
117     fFeedDownOpt(kTRUE),    
118     fFragmentation(kTRUE),
119     fSetNuclei(kFALSE),
120     fNewMIS(kFALSE),   
121     fHFoff(kFALSE),    
122     fNucPdf(0),
123     fTriggerParticle(0),
124     fTriggerEta(0.9),     
125     fTriggerMinPt(-1),  
126     fTriggerMaxPt(1000),  
127     fTriggerMultiplicity(0),
128     fTriggerMultiplicityEta(0),
129     fTriggerMultiplicityPtMin(0),
130     fCountMode(kCountAll),      
131     fHeader(0),  
132     fRL(0),      
133     fkFileName(0),
134     fFragPhotonInCalo(kFALSE),
135     fHadronInCalo(kFALSE) ,
136     fPi0InCalo(kFALSE) ,
137     fEtaInCalo(kFALSE) ,
138     fPhotonInCalo(kFALSE), // not in use
139     fDecayPhotonInCalo(kFALSE),
140     fForceNeutralMeson2PhotonDecay(kFALSE),
141     fEleInCalo(kFALSE),
142     fEleInEMCAL(kFALSE), // not in use
143     fCheckBarrel(kFALSE),
144     fCheckEMCAL(kFALSE),
145     fCheckPHOS(kFALSE),
146     fCheckPHOSeta(kFALSE),
147     fPHOSRotateCandidate(-1),
148     fTriggerParticleMinPt(0), 
149     fPhotonMinPt(0), // not in use
150     fElectronMinPt(0), // not in use
151     fPHOSMinPhi(219.),
152     fPHOSMaxPhi(321.),
153     fPHOSEta(0.13),
154     fEMCALMinPhi(79.),
155     fEMCALMaxPhi(191.),
156     fEMCALEta(0.71),
157     fkTuneForDiff(0),
158     fProcDiff(0)
159 {
160 // Default Constructor
161   fEnergyCMS = 5500.;
162   if (!AliPythiaRndm::GetPythiaRandom()) 
163       AliPythiaRndm::SetPythiaRandom(GetRandom());
164 }
165
166 AliGenPythia::AliGenPythia(Int_t npart)
167     :AliGenMC(npart),
168      fProcess(kPyCharm),          
169      fItune(-1),
170      fStrucFunc(kCTEQ5L), 
171      fKineBias(0.),
172      fTrials(0),
173      fTrialsRun(0),
174      fQ(0.),
175      fX1(0.),
176      fX2(0.),
177      fEventTime(0.),
178      fInteractionRate(0.),
179      fTimeWindow(0.),
180      fCurSubEvent(0),
181      fEventsTime(0),
182      fNev(0),
183      fFlavorSelect(0),
184      fXsection(0.),
185      fPythia(0),
186      fPtHardMin(0.),
187      fPtHardMax(1.e4),
188      fYHardMin(-1.e10),
189      fYHardMax(1.e10),
190      fGinit(kTRUE),
191      fGfinal(kTRUE),
192      fHadronisation(kTRUE),
193      fPatchOmegaDalitz(0), 
194      fNpartons(0),
195      fReadFromFile(kFALSE),
196      fQuench(kFALSE),
197      fQhat(0.),
198      fLength(0.),
199      fpyquenT(1.),
200      fpyquenTau(0.1),
201      fpyquenNf(0),
202      fpyquenEloss(0),
203      fpyquenAngle(0), 
204      fImpact(0.),
205      fPtKick(1.),
206      fFullEvent(kTRUE),
207      fDecayer(new AliDecayerPythia()),
208      fDebugEventFirst(-1),
209      fDebugEventLast(-1),
210      fEtMinJet(0.),      
211      fEtMaxJet(1.e4),      
212      fEtaMinJet(-20.),     
213      fEtaMaxJet(20.),     
214      fPhiMinJet(0.),     
215      fPhiMaxJet(2.* TMath::Pi()),     
216      fJetReconstruction(kCell),
217      fEtaMinGamma(-20.),      
218      fEtaMaxGamma(20.),      
219      fPhiMinGamma(0.),      
220      fPhiMaxGamma(2. * TMath::Pi()),      
221      fPycellEtaMax(2.),     
222      fPycellNEta(274),       
223      fPycellNPhi(432),       
224      fPycellThreshold(0.),  
225      fPycellEtSeed(4.),     
226      fPycellMinEtJet(10.),  
227      fPycellMaxRadius(1.), 
228      fStackFillOpt(kFlavorSelection),   
229      fFeedDownOpt(kTRUE),    
230      fFragmentation(kTRUE),
231      fSetNuclei(kFALSE),
232      fNewMIS(kFALSE),   
233      fHFoff(kFALSE),    
234      fNucPdf(0),
235      fTriggerParticle(0),
236      fTriggerEta(0.9), 
237      fTriggerMinPt(-1),  
238      fTriggerMaxPt(1000),      
239      fTriggerMultiplicity(0),
240      fTriggerMultiplicityEta(0),
241      fTriggerMultiplicityPtMin(0),
242      fCountMode(kCountAll),      
243      fHeader(0),  
244      fRL(0),      
245      fkFileName(0),
246      fFragPhotonInCalo(kFALSE),
247      fHadronInCalo(kFALSE) ,
248      fPi0InCalo(kFALSE) ,
249      fEtaInCalo(kFALSE) ,
250      fPhotonInCalo(kFALSE), // not in use
251      fDecayPhotonInCalo(kFALSE),
252      fForceNeutralMeson2PhotonDecay(kFALSE),
253      fEleInCalo(kFALSE),
254      fEleInEMCAL(kFALSE), // not in use
255      fCheckBarrel(kFALSE),
256      fCheckEMCAL(kFALSE),
257      fCheckPHOS(kFALSE),
258      fCheckPHOSeta(kFALSE),
259      fPHOSRotateCandidate(-1),
260      fTriggerParticleMinPt(0),
261      fPhotonMinPt(0), // not in use
262      fElectronMinPt(0), // not in use
263      fPHOSMinPhi(219.),
264      fPHOSMaxPhi(321.),
265      fPHOSEta(0.13),
266      fEMCALMinPhi(79.),
267      fEMCALMaxPhi(191.),
268      fEMCALEta(0.71),
269      fkTuneForDiff(0),
270      fProcDiff(0)
271 {
272 // default charm production at 5. 5 TeV
273 // semimuonic decay
274 // structure function GRVHO
275 //
276     fEnergyCMS = 5500.;
277     fName = "Pythia";
278     fTitle= "Particle Generator using PYTHIA";
279     SetForceDecay();
280     // Set random number generator 
281     if (!AliPythiaRndm::GetPythiaRandom()) 
282       AliPythiaRndm::SetPythiaRandom(GetRandom());
283  }
284
285 AliGenPythia::~AliGenPythia()
286 {
287 // Destructor
288   if(fEventsTime) delete fEventsTime;
289 }
290
291 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
292 {
293 // Generate pileup using user specified rate
294     fInteractionRate = rate;
295     fTimeWindow = timewindow;
296     GeneratePileup();
297 }
298
299 void AliGenPythia::GeneratePileup()
300 {
301 // Generate sub events time for pileup
302     fEventsTime = 0;
303     if(fInteractionRate == 0.) {
304       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
305       return;
306     }
307
308     Int_t npart = NumberParticles();
309     if(npart < 0) {
310       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
311       return;
312     }
313
314     if(fEventsTime) delete fEventsTime;
315     fEventsTime = new TArrayF(npart);
316     TArrayF &array = *fEventsTime;
317     for(Int_t ipart = 0; ipart < npart; ipart++)
318       array[ipart] = 0.;
319
320     Float_t eventtime = 0.;
321     while(1)
322       {
323         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
324         if(eventtime > fTimeWindow) break;
325         array.Set(array.GetSize()+1);
326         array[array.GetSize()-1] = eventtime;
327       }
328
329     eventtime = 0.;
330     while(1)
331       {
332         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
333         if(TMath::Abs(eventtime) > fTimeWindow) break;
334         array.Set(array.GetSize()+1);
335         array[array.GetSize()-1] = eventtime;
336       }
337
338     SetNumberParticles(fEventsTime->GetSize());
339 }
340
341 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
342                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
343 {
344 // Set pycell parameters
345     fPycellEtaMax    =  etamax;
346     fPycellNEta      =  neta;
347     fPycellNPhi      =  nphi;
348     fPycellThreshold =  thresh;
349     fPycellEtSeed    =  etseed;
350     fPycellMinEtJet  =  minet;
351     fPycellMaxRadius =  r;
352 }
353
354
355
356 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
357 {
358   // Set a range of event numbers, for which a table
359   // of generated particle will be printed
360   fDebugEventFirst = eventFirst;
361   fDebugEventLast  = eventLast;
362   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
363 }
364
365 void AliGenPythia::Init()
366 {
367 // Initialisation
368     
369     SetMC(AliPythia::Instance());
370     fPythia=(AliPythia*) fMCEvGen;
371     
372 //
373     fParentWeight=1./Float_t(fNpart);
374 //
375
376
377     fPythia->SetCKIN(3,fPtHardMin);
378     fPythia->SetCKIN(4,fPtHardMax);
379     fPythia->SetCKIN(7,fYHardMin);
380     fPythia->SetCKIN(8,fYHardMax);
381     
382     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
383     if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
384   
385     // Fragmentation?
386     if (fFragmentation) {
387       fPythia->SetMSTP(111,1);
388     } else {
389       fPythia->SetMSTP(111,0);
390     }
391
392
393 //  initial state radiation   
394     fPythia->SetMSTP(61,fGinit);
395 //  final state radiation
396     fPythia->SetMSTP(71,fGfinal);
397 //  pt - kick
398     if (fPtKick > 0.) {
399         fPythia->SetMSTP(91,1);
400         fPythia->SetPARP(91,fPtKick);   
401         fPythia->SetPARP(93, 4. * fPtKick);
402     } else {
403         fPythia->SetMSTP(91,0);
404     }
405
406
407     if (fReadFromFile) {
408         fRL  =  AliRunLoader::Open(fkFileName, "Partons");
409         fRL->LoadKinematics();
410         fRL->LoadHeader();
411     } else {
412         fRL = 0x0;
413     }
414  //
415     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
416     //  Forward Paramters to the AliPythia object
417     fDecayer->SetForceDecay(fForceDecay);    
418 // Switch off Heavy Flavors on request  
419     if (fHFoff) {
420         // Maximum number of quark flavours used in pdf 
421         fPythia->SetMSTP(58, 3);
422         // Maximum number of flavors that can be used in showers
423         fPythia->SetMSTJ(45, 3);        
424         // Switch off g->QQbar splitting in decay table
425         ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
426     }
427
428     fDecayer->Init();
429
430
431 //  Parent and Children Selection
432     switch (fProcess) 
433     {
434     case kPyOldUEQ2ordered:
435     case kPyOldUEQ2ordered2:
436     case kPyOldPopcorn:
437       break;
438     case kPyCharm:
439     case kPyCharmUnforced:
440     case kPyCharmPbPbMNR:
441     case kPyCharmpPbMNR:
442     case kPyCharmppMNR:
443     case kPyCharmppMNRwmi:
444         fParentSelect[0] =   411;
445         fParentSelect[1] =   421;
446         fParentSelect[2] =   431;
447         fParentSelect[3] =  4122;
448         fParentSelect[4] =  4232;
449         fParentSelect[5] =  4132;
450         fParentSelect[6] =  4332;
451         fFlavorSelect    =  4;  
452         break;
453     case kPyD0PbPbMNR:
454     case kPyD0pPbMNR:
455     case kPyD0ppMNR:
456         fParentSelect[0] =   421;
457         fFlavorSelect    =   4; 
458         break;
459     case kPyDPlusPbPbMNR:
460     case kPyDPluspPbMNR:
461     case kPyDPlusppMNR:
462         fParentSelect[0] =   411;
463         fFlavorSelect    =   4; 
464         break;
465     case kPyDPlusStrangePbPbMNR:
466     case kPyDPlusStrangepPbMNR:
467     case kPyDPlusStrangeppMNR:
468         fParentSelect[0] =   431;
469         fFlavorSelect    =   4; 
470         break;
471     case kPyLambdacppMNR:
472         fParentSelect[0] =  4122;
473         fFlavorSelect    =   4; 
474         break;      
475     case kPyBeauty:
476     case kPyBeautyJets:
477     case kPyBeautyPbPbMNR:
478     case kPyBeautypPbMNR:
479     case kPyBeautyppMNR:
480     case kPyBeautyppMNRwmi:
481         fParentSelect[0]=  511;
482         fParentSelect[1]=  521;
483         fParentSelect[2]=  531;
484         fParentSelect[3]= 5122;
485         fParentSelect[4]= 5132;
486         fParentSelect[5]= 5232;
487         fParentSelect[6]= 5332;
488         fFlavorSelect   = 5;    
489         break;
490     case kPyBeautyUnforced:
491         fParentSelect[0] =  511;
492         fParentSelect[1] =  521;
493         fParentSelect[2] =  531;
494         fParentSelect[3] = 5122;
495         fParentSelect[4] = 5132;
496         fParentSelect[5] = 5232;
497         fParentSelect[6] = 5332;
498         fFlavorSelect    = 5;   
499         break;
500     case kPyJpsiChi:
501     case kPyJpsi:
502         fParentSelect[0] = 443;
503         break;
504     case kPyMbDefault:
505     case kPyMbAtlasTuneMC09:
506     case kPyMb:
507     case kPyMbWithDirectPhoton:
508     case kPyMbNonDiffr:
509     case kPyMbMSEL1:
510     case kPyJets:
511     case kPyDirectGamma:
512     case kPyLhwgMb:     
513         break;
514     case kPyW:
515     case kPyZ:
516     case kPyMBRSingleDiffraction:
517     case kPyMBRDoubleDiffraction:
518     case kPyMBRCentralDiffraction:
519         break;
520     }
521 //
522 //
523 //  JetFinder for Trigger
524 //
525 //  Configure detector (EMCAL like)
526 //
527     fPythia->SetPARU(51, fPycellEtaMax);
528     fPythia->SetMSTU(51, fPycellNEta);
529     fPythia->SetMSTU(52, fPycellNPhi);
530 //
531 //  Configure Jet Finder
532 //  
533     fPythia->SetPARU(58,  fPycellThreshold);
534     fPythia->SetPARU(52,  fPycellEtSeed);
535     fPythia->SetPARU(53,  fPycellMinEtJet);
536     fPythia->SetPARU(54,  fPycellMaxRadius);
537     fPythia->SetMSTU(54,  2);
538 //
539 //  This counts the total number of calls to Pyevnt() per run.
540     fTrialsRun = 0;
541     fQ         = 0.;
542     fX1        = 0.;
543     fX2        = 0.;    
544     fNev       = 0 ;
545 //    
546 //
547 //
548     AliGenMC::Init();
549 //
550 //
551 //  
552     if (fSetNuclei) {
553         fDyBoost = 0;
554         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
555     }
556     
557     fPythia->SetPARJ(200, 0.0);
558     fPythia->SetPARJ(199, 0.0);
559     fPythia->SetPARJ(198, 0.0);
560     fPythia->SetPARJ(197, 0.0);
561
562     if (fQuench == 1) {
563         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
564     }
565
566     if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
567   
568     if (fQuench == 3) {
569         // Nestor's change of the splittings
570         fPythia->SetPARJ(200, 0.8);
571         fPythia->SetMSTJ(41, 1);  // QCD radiation only
572         fPythia->SetMSTJ(42, 2);  // angular ordering
573         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
574         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
575         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
576         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
577     } else if (fQuench == 4) {
578         // Armesto-Cunqueiro-Salgado change of the splittings.
579         AliFastGlauber* glauber = AliFastGlauber::Instance();
580         glauber->Init(2);
581         //read and store transverse almonds corresponding to differnt
582         //impact parameters.
583         glauber->SetCentralityClass(0.,0.1);
584         fPythia->SetPARJ(200, 1.);
585         fPythia->SetPARJ(198, fQhat);
586         fPythia->SetPARJ(199, fLength);
587         fPythia->SetMSTJ(42, 2);  // angular ordering
588         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
589         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
590     }
591 }
592
593 void AliGenPythia::Generate()
594 {
595 // Generate one event
596     if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
597     fDecayer->ForceDecay();
598
599     Double_t polar[3]   =   {0,0,0};
600     Double_t origin[3]  =   {0,0,0};
601     Double_t p[4];
602 //  converts from mm/c to s
603     const Float_t kconv=0.001/2.999792458e8;
604 //
605     Int_t nt=0;
606     Int_t jev=0;
607     Int_t j, kf;
608     fTrials=0;
609     fEventTime = 0.;
610     
611     
612
613     //  Set collision vertex position 
614     if (fVertexSmear == kPerEvent) Vertex();
615     
616 //  event loop    
617     while(1)
618     {
619 //
620 // Produce event
621 //
622 //
623 // Switch hadronisation off
624 //
625         fPythia->SetMSTJ(1, 0);
626
627         if (fQuench ==4){
628             Double_t bimp;
629             // Quenching comes through medium-modified splitting functions.
630             AliFastGlauber::Instance()->GetRandomBHard(bimp);
631             fPythia->SetPARJ(197, bimp);
632             fImpact = bimp;
633             fPythia->Qpygin0();
634         } 
635 //
636 // Either produce new event or read partons from file
637 //      
638         if (!fReadFromFile) {
639             if (!fNewMIS) {
640                 fPythia->Pyevnt();
641             } else {
642                 fPythia->Pyevnw();
643             }
644             fNpartons = fPythia->GetN();
645         } else {
646             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
647             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
648             fPythia->SetN(0);
649             LoadEvent(fRL->Stack(), 0 , 1);
650             fPythia->Pyedit(21);
651         }
652         
653 //
654 //  Run quenching routine 
655 //
656         if (fQuench == 1) {
657             fPythia->Quench();
658         } else if (fQuench == 2){
659             fPythia->Pyquen(208., 0, 0.);
660         } else if (fQuench == 3) {
661             // Quenching is via multiplicative correction of the splittings
662         }
663         
664 //
665 // Switch hadronisation on
666 //
667         if (fHadronisation) {
668             fPythia->SetMSTJ(1, 1);
669 //
670 // .. and perform hadronisation
671 //      printf("Calling hadronisation %d\n", fPythia->GetN());
672
673             if (fPatchOmegaDalitz) {
674               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
675               fPythia->Pyexec();
676               fPythia->DalitzDecays();
677               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
678             } 
679             fPythia->Pyexec();
680         }
681         fTrials++;
682         fPythia->ImportParticles(&fParticles,"All");
683         
684         if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
685         if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
686         
687 //
688 //
689 //
690         Int_t i;
691         
692         fNprimaries = 0;
693         Int_t np = fParticles.GetEntriesFast();
694         
695         if (np == 0) continue;
696 //
697         
698 //
699         Int_t* pParent   = new Int_t[np];
700         Int_t* pSelected = new Int_t[np];
701         Int_t* trackIt   = new Int_t[np];
702         for (i = 0; i < np; i++) {
703             pParent[i]   = -1;
704             pSelected[i] =  0;
705             trackIt[i]   =  0;
706         }
707
708         Int_t nc = 0;        // Total n. of selected particles
709         Int_t nParents = 0;  // Selected parents
710         Int_t nTkbles = 0;   // Trackable particles
711         if (fProcess != kPyMbDefault && 
712             fProcess != kPyMb && 
713             fProcess != kPyMbAtlasTuneMC09 && 
714             fProcess != kPyMbWithDirectPhoton && 
715             fProcess != kPyJets && 
716             fProcess != kPyDirectGamma &&
717             fProcess != kPyMbNonDiffr  &&
718             fProcess != kPyMbMSEL1     &&
719             fProcess != kPyW && 
720             fProcess != kPyZ &&
721             fProcess != kPyCharmppMNRwmi && 
722             fProcess != kPyBeautyppMNRwmi &&
723             fProcess != kPyBeautyJets) {
724             
725             for (i = 0; i < np; i++) {
726                 TParticle* iparticle = (TParticle *) fParticles.At(i);
727                 Int_t ks = iparticle->GetStatusCode();
728                 kf = CheckPDGCode(iparticle->GetPdgCode());
729 // No initial state partons
730                 if (ks==21) continue;
731 //
732 // Heavy Flavor Selection
733 //
734                 // quark ?
735                 kf = TMath::Abs(kf);
736                 Int_t kfl = kf;
737                 // Resonance
738
739                 if (kfl > 100000) kfl %= 100000;
740                 if (kfl > 10000)  kfl %= 10000;
741                 // meson ?
742                 if  (kfl > 10) kfl/=100;
743                 // baryon
744                 if (kfl > 10) kfl/=10;
745                 Int_t ipa = iparticle->GetFirstMother()-1;
746                 Int_t kfMo = 0;
747 //
748 // Establish mother daughter relation between heavy quarks and mesons
749 //
750                 if (kf >= fFlavorSelect && kf <= 6) {
751                     Int_t idau = iparticle->GetFirstDaughter() - 1;
752                     if (idau > -1) {
753                         TParticle* daughter = (TParticle *) fParticles.At(idau);
754                         Int_t pdgD = daughter->GetPdgCode();
755                         if (pdgD == 91 || pdgD == 92) {
756                             Int_t jmin = daughter->GetFirstDaughter() - 1;
757                             Int_t jmax = daughter->GetLastDaughter()  - 1;                          
758                             for (Int_t jp = jmin; jp <= jmax; jp++)
759                                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
760                         } // is string or cluster
761                     } // has daughter
762                 } // heavy quark
763                 
764
765                 if (ipa > -1) {
766                     TParticle *  mother = (TParticle *) fParticles.At(ipa);
767                     kfMo = TMath::Abs(mother->GetPdgCode());
768                 }
769                 
770                 // What to keep in Stack?
771                 Bool_t flavorOK = kFALSE;
772                 Bool_t selectOK = kFALSE;
773                 if (fFeedDownOpt) {
774                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
775                 } else {
776                     if (kfl > fFlavorSelect) {
777                         nc = -1;
778                         break;
779                     }
780                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
781                 }
782                 switch (fStackFillOpt) {
783                 case kHeavyFlavor:
784                 case kFlavorSelection:
785                     selectOK = kTRUE;
786                     break;
787                 case kParentSelection:
788                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
789                     break;
790                 }
791                 if (flavorOK && selectOK) { 
792 //
793 // Heavy flavor hadron or quark
794 //
795 // Kinematic seletion on final state heavy flavor mesons
796                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
797                     {
798                         continue;
799                     }
800                     pSelected[i] = 1;
801                     if (ParentSelected(kf)) ++nParents; // Update parent count
802 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
803                 } else {
804 // Kinematic seletion on decay products
805                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
806                         && !KinematicSelection(iparticle, 1)) 
807                     {
808                         continue;
809                     }
810 //
811 // Decay products 
812 // Select if mother was selected and is not tracked
813
814                     if (pSelected[ipa] && 
815                         !trackIt[ipa]  &&     // mother will be  tracked ?
816                         kfMo !=  5 &&         // mother is b-quark, don't store fragments          
817                         kfMo !=  4 &&         // mother is c-quark, don't store fragments 
818                         kf   != 92)           // don't store string
819                     {
820 //
821 // Semi-stable or de-selected: diselect decay products:
822 // 
823 //
824                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
825                         {
826                             Int_t ipF = iparticle->GetFirstDaughter();
827                             Int_t ipL = iparticle->GetLastDaughter();   
828                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
829                         }
830 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
831                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
832                     }
833                 }
834                 if (pSelected[i] == -1) pSelected[i] = 0;
835                 if (!pSelected[i]) continue;
836                 // Count quarks only if you did not include fragmentation
837                 if (fFragmentation && kf <= 10) continue;
838
839                 nc++;
840 // Decision on tracking
841                 trackIt[i] = 0;
842 //
843 // Track final state particle
844                 if (ks == 1) trackIt[i] = 1;
845 // Track semi-stable particles
846                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
847 // Track particles selected by process if undecayed. 
848                 if (fForceDecay == kNoDecay) {
849                     if (ParentSelected(kf)) trackIt[i] = 1;
850                 } else {
851                     if (ParentSelected(kf)) trackIt[i] = 0;
852                 }
853                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
854 //
855 //
856
857             } // particle selection loop
858             if (nc > 0) {
859                 for (i = 0; i<np; i++) {
860                     if (!pSelected[i]) continue;
861                     TParticle *  iparticle = (TParticle *) fParticles.At(i);
862                     kf = CheckPDGCode(iparticle->GetPdgCode());
863                     Int_t ks = iparticle->GetStatusCode();  
864                     p[0] = iparticle->Px();
865                     p[1] = iparticle->Py();
866                     p[2] = iparticle->Pz();
867                     p[3] = iparticle->Energy();
868                     
869                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
870                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
871                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
872                     
873                     Float_t tof   = fTime + kconv*iparticle->T();
874                     Int_t ipa     = iparticle->GetFirstMother()-1;
875                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
876  
877                     PushTrack(fTrackIt*trackIt[i], iparent, kf, 
878                               p[0], p[1], p[2], p[3], 
879                               origin[0], origin[1], origin[2], tof, 
880                               polar[0], polar[1], polar[2],
881                               kPPrimary, nt, 1., ks);
882                     pParent[i] = nt;
883                     KeepTrack(nt);
884                     fNprimaries++;
885                 } //  PushTrack loop
886             }
887         } else {
888             nc = GenerateMB();
889         } // mb ?
890         
891         GetSubEventTime();
892
893         delete[] pParent;
894         delete[] pSelected;
895         delete[] trackIt;
896
897         if (nc > 0) {
898           switch (fCountMode) {
899           case kCountAll:
900             // printf(" Count all \n");
901             jev += nc;
902             break;
903           case kCountParents:
904             // printf(" Count parents \n");
905             jev += nParents;
906             break;
907           case kCountTrackables:
908             // printf(" Count trackable \n");
909             jev += nTkbles;
910             break;
911           }
912             if (jev >= fNpart || fNpart == -1) {
913                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
914                 
915                 fQ  += fPythia->GetVINT(51);
916                 fX1 += fPythia->GetVINT(41);
917                 fX2 += fPythia->GetVINT(42);
918                 fTrialsRun += fTrials;
919                 fNev++;
920                 MakeHeader();
921                 break;
922             }
923         }
924     } // event loop
925     SetHighWaterMark(nt);
926 //  adjust weight due to kinematic selection
927 //    AdjustWeights();
928 //  get cross-section
929     fXsection=fPythia->GetPARI(1);
930 }
931
932 Int_t  AliGenPythia::GenerateMB()
933 {
934 //
935 // Min Bias selection and other global selections
936 //
937     Int_t i, kf, nt, iparent;
938     Int_t nc = 0;
939     Double_t p[4];
940     Double_t polar[3]   =   {0,0,0};
941     Double_t origin[3]  =   {0,0,0};
942 //  converts from mm/c to s
943     const Float_t kconv=0.001/2.999792458e8;
944     
945     
946     Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
947
948
949
950     Int_t* pParent = new Int_t[np];
951     for (i=0; i< np; i++) pParent[i] = -1;
952     if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) 
953         && fEtMinJet > 0.) {
954         TParticle* jet1 = (TParticle *) fParticles.At(6);
955         TParticle* jet2 = (TParticle *) fParticles.At(7);
956         
957         if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
958           delete [] pParent;
959           return 0;
960         }
961     }
962
963   
964   // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
965   // implemented primaryly for kPyJets, but extended to any kind of process.
966   if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) && 
967       (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
968     Bool_t ok = TriggerOnSelectedParticles(np);
969     
970     if(!ok) {
971       delete[] pParent;
972       return 0;
973     }
974   }
975
976     // Check for diffraction
977     if(fkTuneForDiff) {
978       if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
979         if(!CheckDiffraction()) {
980           delete [] pParent;
981           return 0;
982         }
983       }
984     }
985
986     // Check for minimum multiplicity
987     if (fTriggerMultiplicity > 0) {
988       Int_t multiplicity = 0;
989       for (i = 0; i < np; i++) {
990         TParticle *  iparticle = (TParticle *) fParticles.At(i);
991         
992         Int_t statusCode = iparticle->GetStatusCode();
993         
994         // Initial state particle
995         if (statusCode != 1)
996           continue;
997         // eta cut
998         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
999           continue;
1000         // pt cut
1001         if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
1002             continue;
1003
1004         TParticlePDG* pdgPart = iparticle->GetPDG();
1005         if (pdgPart && pdgPart->Charge() == 0)
1006           continue;
1007         
1008         ++multiplicity;
1009       }
1010
1011       if (multiplicity < fTriggerMultiplicity) {
1012         delete [] pParent;
1013         return 0;
1014       }
1015       Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1016     }    
1017     
1018     
1019     if (fTriggerParticle) {
1020         Bool_t triggered = kFALSE;
1021         for (i = 0; i < np; i++) {
1022             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1023             kf = CheckPDGCode(iparticle->GetPdgCode());
1024             if (kf != fTriggerParticle) continue;
1025             if (iparticle->Pt() == 0.) continue;
1026             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1027             if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
1028             triggered = kTRUE;
1029             break;
1030         }
1031         if (!triggered) {
1032           delete [] pParent;
1033           return 0;
1034         }
1035     }
1036         
1037
1038     // Check if there is a ccbar or bbbar pair with at least one of the two
1039     // in fYMin < y < fYMax
1040
1041     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1042       TParticle *partCheck;
1043       TParticle *mother;
1044       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1045       Bool_t  theChild=kFALSE;
1046       Bool_t  triggered=kFALSE;
1047       Float_t y;  
1048       Int_t   pdg,mpdg,mpdgUpperFamily;
1049       for(i=0; i<np; i++) {
1050         partCheck = (TParticle*)fParticles.At(i);
1051         pdg = partCheck->GetPdgCode();  
1052         if(TMath::Abs(pdg) == fFlavorSelect) { // quark  
1053           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1054           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1055                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
1056           if(y>fYMin && y<fYMax) inYcut=kTRUE;
1057         }
1058         if(fTriggerParticle) {
1059           if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1060         }
1061         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1062           Int_t mi = partCheck->GetFirstMother() - 1;
1063           if(mi<0) continue;
1064           mother = (TParticle*)fParticles.At(mi);
1065           mpdg=TMath::Abs(mother->GetPdgCode());
1066           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1067           if ( ParentSelected(mpdg) || 
1068               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1069             if (KinematicSelection(partCheck,1)) {
1070               theChild=kTRUE;
1071             }
1072           }
1073         }
1074       }
1075       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1076         delete[] pParent;
1077         return 0;
1078       }
1079       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1080         delete[] pParent;
1081         return 0;       
1082       }
1083       if(fTriggerParticle && !triggered) { // particle requested is not produced
1084         delete[] pParent;
1085         return 0;       
1086       }
1087
1088     }
1089
1090     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1091     if ( (fProcess == kPyW ||
1092           fProcess == kPyZ ||
1093           fProcess == kPyMbDefault ||
1094           fProcess == kPyMb ||
1095           fProcess == kPyMbAtlasTuneMC09 ||
1096           fProcess == kPyMbWithDirectPhoton ||
1097           fProcess == kPyMbNonDiffr)  
1098          && (fCutOnChild == 1) ) {
1099       if ( !CheckKinematicsOnChild() ) {
1100         delete[] pParent;
1101         return 0;
1102       }
1103     }
1104   
1105
1106     for (i = 0; i < np; i++) {
1107         Int_t trackIt = 0;
1108         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1109         kf = CheckPDGCode(iparticle->GetPdgCode());
1110         Int_t ks = iparticle->GetStatusCode();
1111         Int_t km = iparticle->GetFirstMother();
1112         if (
1113             (((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1114             ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
1115             ) 
1116           {
1117             nc++;
1118             if (ks == 1) trackIt = 1;
1119             Int_t ipa = iparticle->GetFirstMother()-1;
1120             
1121             iparent = (ipa > -1) ? pParent[ipa] : -1;
1122             
1123 //
1124 // store track information
1125             p[0] = iparticle->Px();
1126             p[1] = iparticle->Py();
1127             p[2] = iparticle->Pz();
1128             p[3] = iparticle->Energy();
1129
1130             
1131             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1132             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1133             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1134             
1135             Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1136
1137             PushTrack(fTrackIt*trackIt, iparent, kf, 
1138                       p[0], p[1], p[2], p[3], 
1139                       origin[0], origin[1], origin[2], tof, 
1140                       polar[0], polar[1], polar[2],
1141                       kPPrimary, nt, 1., ks);
1142             fNprimaries++;
1143             KeepTrack(nt);
1144             pParent[i] = nt;
1145             SetHighWaterMark(nt);
1146             
1147         } // select particle
1148     } // particle loop 
1149
1150     delete[] pParent;
1151     
1152     return 1;
1153 }
1154
1155
1156 void AliGenPythia::FinishRun()
1157 {
1158 // Print x-section summary
1159     fPythia->Pystat(1);
1160
1161     if (fNev > 0.) {
1162         fQ  /= fNev;
1163         fX1 /= fNev;
1164         fX2 /= fNev;    
1165     }
1166     
1167     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1168     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1169 }
1170
1171 void AliGenPythia::AdjustWeights() const
1172 {
1173 // Adjust the weights after generation of all events
1174 //
1175     if (gAlice) {
1176         TParticle *part;
1177         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1178         for (Int_t i=0; i<ntrack; i++) {
1179             part= gAlice->GetMCApp()->Particle(i);
1180             part->SetWeight(part->GetWeight()*fKineBias);
1181         }
1182     }
1183 }
1184     
1185 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1186 {
1187 // Treat protons as inside nuclei with mass numbers a1 and a2  
1188
1189     fAProjectile = a1;
1190     fATarget     = a2;
1191     fNucPdf      = pdfset;  // 0 EKS98 9 EPS09LO 19 EPS09NLO
1192     fSetNuclei   = kTRUE;
1193 }
1194
1195
1196 void AliGenPythia::MakeHeader()
1197 {
1198 //
1199 // Make header for the simulated event
1200 // 
1201   if (gAlice) {
1202     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1203         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1204   }
1205
1206 // Builds the event header, to be called after each event
1207     if (fHeader) delete fHeader;
1208     fHeader = new AliGenPythiaEventHeader("Pythia");
1209 //
1210 // Event type  
1211     if(fProcDiff>0){
1212       //      if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1213       //      printf("fPythia->GetMSTI(1) = %d   fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1214     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1215     }
1216     else 
1217     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1218 //
1219 // Number of trials
1220     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1221 //
1222 // Event Vertex 
1223     fHeader->SetPrimaryVertex(fVertex);
1224     fHeader->SetInteractionTime(fTime+fEventTime);
1225 //
1226 // Number of primaries
1227     fHeader->SetNProduced(fNprimaries);
1228 //
1229 // Jets that have triggered
1230
1231     //Need to store jets for b-jet studies too!
1232     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1233     {
1234         Int_t ntrig, njet;
1235         Float_t jets[4][10];
1236         GetJets(njet, ntrig, jets);
1237
1238         
1239         for (Int_t i = 0; i < ntrig; i++) {
1240             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
1241                                                         jets[3][i]);
1242         }
1243     }
1244 //
1245 // Copy relevant information from external header, if present.
1246 //
1247     Float_t uqJet[4];
1248     
1249     if (fRL) {
1250         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1251         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1252         {
1253             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1254             
1255             
1256             exHeader->TriggerJet(i, uqJet);
1257             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1258         }
1259     }
1260 //
1261 // Store quenching parameters
1262 //
1263     if (fQuench){
1264         Double_t z[4] = {0.};
1265         Double_t xp = 0.;
1266         Double_t yp = 0.;
1267         
1268         if (fQuench == 1) {
1269             // Pythia::Quench()
1270             fPythia->GetQuenchingParameters(xp, yp, z);
1271         } else if (fQuench == 2){
1272             // Pyquen
1273             Double_t r1 = PARIMP.rb1;
1274             Double_t r2 = PARIMP.rb2;
1275             Double_t b  = PARIMP.b1;
1276             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1277             Double_t phi = PARIMP.psib1;
1278             xp = r * TMath::Cos(phi);
1279             yp = r * TMath::Sin(phi);
1280             
1281         } else if (fQuench == 4) {
1282             // QPythia
1283             Double_t xy[2];
1284             Double_t i0i1[2];
1285             AliFastGlauber::Instance()->GetSavedXY(xy);
1286             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1287             xp = xy[0];
1288             yp = xy[1];
1289             ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1290         }
1291         
1292             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1293             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1294     }
1295 //
1296 // Store pt^hard 
1297     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1298 //
1299 //  Pass header
1300 //
1301     AddHeader(fHeader);
1302     fHeader = 0x0;
1303 }
1304
1305 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1306 {
1307 // Check the kinematic trigger condition
1308 //
1309     Double_t eta[2];
1310     eta[0] = jet1->Eta();
1311     eta[1] = jet2->Eta();
1312     Double_t phi[2];
1313     phi[0] = jet1->Phi();
1314     phi[1] = jet2->Phi();
1315     Int_t    pdg[2]; 
1316     pdg[0] = jet1->GetPdgCode();
1317     pdg[1] = jet2->GetPdgCode();    
1318     Bool_t   triggered = kFALSE;
1319
1320     if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi) {
1321         Int_t njets = 0;
1322         Int_t ntrig = 0;
1323         Float_t jets[4][10];
1324 //
1325 // Use Pythia clustering on parton level to determine jet axis
1326 //
1327         GetJets(njets, ntrig, jets);
1328         
1329         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1330 //
1331     } else {
1332         Int_t ij = 0;
1333         Int_t ig = 1;
1334         if (pdg[0] == kGamma) {
1335             ij = 1;
1336             ig = 0;
1337         }
1338         //Check eta range first...
1339         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1340             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1341         {
1342             //Eta is okay, now check phi range
1343             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1344                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1345             {
1346                 triggered = kTRUE;
1347             }
1348         }
1349     }
1350     return triggered;
1351 }
1352
1353
1354
1355 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1356 //
1357 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1358 //
1359     Bool_t checking = kFALSE;
1360     Int_t j, kcode, ks, km;
1361     Int_t nPartAcc = 0; //number of particles in the acceptance range
1362     Int_t numberOfAcceptedParticles = 1;
1363     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1364     Int_t npart = fParticles.GetEntriesFast();
1365     
1366     for (j = 0; j<npart; j++) {
1367         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1368         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1369         ks = jparticle->GetStatusCode();
1370         km = jparticle->GetFirstMother(); 
1371         
1372         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1373             nPartAcc++;
1374         }
1375         if( numberOfAcceptedParticles <= nPartAcc){
1376           checking = kTRUE;
1377           break;
1378         }
1379     }
1380
1381     return checking;
1382 }
1383
1384 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1385 {
1386   //
1387   // Load event into Pythia Common Block
1388   //
1389   
1390   Int_t npart = stack -> GetNprimary();
1391   Int_t n0 = 0;
1392   
1393   if (!flag) {
1394     (fPythia->GetPyjets())->N = npart;
1395   } else {
1396     n0 = (fPythia->GetPyjets())->N;
1397     (fPythia->GetPyjets())->N = n0 + npart;
1398   }
1399   
1400   
1401   for (Int_t part = 0; part < npart; part++) {
1402     TParticle *mPart = stack->Particle(part);
1403     
1404     Int_t kf     =  mPart->GetPdgCode();
1405     Int_t ks     =  mPart->GetStatusCode();
1406     Int_t idf    =  mPart->GetFirstDaughter();
1407     Int_t idl    =  mPart->GetLastDaughter();
1408     
1409     if (reHadr) {
1410             if (ks == 11 || ks == 12) {
1411         ks  -= 10;
1412         idf  = -1;
1413         idl  = -1;
1414             }
1415     }
1416     
1417     Float_t px = mPart->Px();
1418     Float_t py = mPart->Py();
1419     Float_t pz = mPart->Pz();
1420     Float_t e  = mPart->Energy();
1421     Float_t m  = mPart->GetCalcMass();
1422     
1423     
1424     (fPythia->GetPyjets())->P[0][part+n0] = px;
1425     (fPythia->GetPyjets())->P[1][part+n0] = py;
1426     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1427     (fPythia->GetPyjets())->P[3][part+n0] = e;
1428     (fPythia->GetPyjets())->P[4][part+n0] = m;
1429     
1430     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1431     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1432     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1433     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1434     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1435   }
1436 }
1437
1438 void  AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1439 {
1440   //
1441   // Load event into Pythia Common Block
1442   //
1443   
1444   Int_t npart = stack -> GetEntries();
1445   Int_t n0 = 0;
1446   
1447   if (!flag) {
1448     (fPythia->GetPyjets())->N = npart;
1449   } else {
1450     n0 = (fPythia->GetPyjets())->N;
1451     (fPythia->GetPyjets())->N = n0 + npart;
1452   }
1453   
1454   
1455   for (Int_t part = 0; part < npart; part++) {
1456     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1457     if (!mPart) continue;
1458     
1459     Int_t kf     =  mPart->GetPdgCode();
1460     Int_t ks     =  mPart->GetStatusCode();
1461     Int_t idf    =  mPart->GetFirstDaughter();
1462     Int_t idl    =  mPart->GetLastDaughter();
1463     
1464     if (reHadr) {
1465         if (ks == 11 || ks == 12) {
1466             ks  -= 10;
1467             idf  = -1;
1468             idl  = -1;
1469         }
1470     }
1471     
1472     Float_t px = mPart->Px();
1473     Float_t py = mPart->Py();
1474     Float_t pz = mPart->Pz();
1475     Float_t e  = mPart->Energy();
1476     Float_t m  = mPart->GetCalcMass();
1477     
1478     
1479     (fPythia->GetPyjets())->P[0][part+n0] = px;
1480     (fPythia->GetPyjets())->P[1][part+n0] = py;
1481     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1482     (fPythia->GetPyjets())->P[3][part+n0] = e;
1483     (fPythia->GetPyjets())->P[4][part+n0] = m;
1484     
1485     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1486     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1487     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1488     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1489     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1490   }
1491 }
1492
1493
1494 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1495 {
1496 //
1497 //  Calls the Pythia jet finding algorithm to find jets in the current event
1498 //
1499 //
1500 //
1501 //  Save jets
1502     Int_t n     = fPythia->GetN();
1503
1504 //
1505 //  Run Jet Finder
1506     fPythia->Pycell(njets);
1507     Int_t i;
1508     for (i = 0; i < njets; i++) {
1509         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1510         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1511         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1512         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1513
1514         jets[0][i] = px;
1515         jets[1][i] = py;
1516         jets[2][i] = pz;
1517         jets[3][i] = e;
1518     }
1519 }
1520
1521
1522
1523 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1524 {
1525 //
1526 //  Calls the Pythia clustering algorithm to find jets in the current event
1527 //
1528     Int_t n     = fPythia->GetN();
1529     nJets       = 0;
1530     nJetsTrig   = 0;
1531     if (fJetReconstruction == kCluster) {
1532 //
1533 //  Configure cluster algorithm
1534 //    
1535         fPythia->SetPARU(43, 2.);
1536         fPythia->SetMSTU(41, 1);
1537 //
1538 //  Call cluster algorithm
1539 //    
1540         fPythia->Pyclus(nJets);
1541 //
1542 //  Loading jets from common block
1543 //
1544     } else {
1545
1546 //
1547 //  Run Jet Finder
1548         fPythia->Pycell(nJets);
1549     }
1550
1551     Int_t i;
1552     for (i = 0; i < nJets; i++) {
1553         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1554         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1555         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1556         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1557         Float_t pt    = TMath::Sqrt(px * px + py * py);
1558         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);  
1559         Float_t theta = TMath::ATan2(pt,pz);
1560         Float_t et    = e * TMath::Sin(theta);
1561         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1562         if (
1563             eta > fEtaMinJet && eta < fEtaMaxJet && 
1564             phi > fPhiMinJet && phi < fPhiMaxJet &&
1565             et  > fEtMinJet  && et  < fEtMaxJet     
1566             ) 
1567         {
1568             jets[0][nJetsTrig] = px;
1569             jets[1][nJetsTrig] = py;
1570             jets[2][nJetsTrig] = pz;
1571             jets[3][nJetsTrig] = e;
1572             nJetsTrig++;
1573 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1574         } else {
1575 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1576         }
1577     }
1578 }
1579
1580 void AliGenPythia::GetSubEventTime()
1581 {
1582   // Calculates time of the next subevent
1583   fEventTime = 0.;
1584   if (fEventsTime) {
1585     TArrayF &array = *fEventsTime;
1586     fEventTime = array[fCurSubEvent++];
1587   }
1588   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1589   return;
1590 }
1591
1592 Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
1593 {
1594   // Is particle in Central Barrel acceptance? 
1595   // etamin=-etamax
1596   if( eta < fTriggerEta  ) 
1597     return kTRUE;
1598   else 
1599     return kFALSE;
1600 }
1601
1602 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
1603 {
1604   // Is particle in EMCAL acceptance? 
1605   // phi in degrees, etamin=-etamax
1606   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi && 
1607      eta < fEMCALEta  ) 
1608     return kTRUE;
1609   else 
1610     return kFALSE;
1611 }
1612
1613 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
1614 {
1615   // Is particle in PHOS acceptance? 
1616   // Acceptance slightly larger considered.
1617   // phi in degrees, etamin=-etamax
1618   // iparticle is the index of the particle to be checked, for PHOS rotation case
1619
1620   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi && 
1621      eta < fPHOSEta  ) 
1622     return kTRUE;
1623   else 
1624   {
1625     if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1626
1627     return kFALSE;
1628   }
1629 }
1630
1631 void AliGenPythia::RotatePhi(Bool_t& okdd)
1632 {
1633   //Rotate event in phi to enhance events in PHOS acceptance
1634   
1635   if(fPHOSRotateCandidate < 0) return ; 
1636   
1637   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi 
1638   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1639   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1640   Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1641   
1642   //calculate deltaphi
1643   TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
1644   Double_t phphi = ph->Phi();
1645   Double_t deltaphi = phiPHOS - phphi;
1646   
1647   
1648   
1649   //loop for all particles and produce the phi rotation
1650   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1651   Double_t oldphi, newphi;
1652   Double_t newVx, newVy, r, vZ, time; 
1653   Double_t newPx, newPy, pt, pz, e;
1654   for(Int_t i=0; i< np; i++) {
1655     TParticle* iparticle = (TParticle *) fParticles.At(i);
1656     oldphi = iparticle->Phi();
1657     newphi = oldphi + deltaphi;
1658     if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle 
1659     if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1660     
1661     r = iparticle->R();
1662     newVx = r * TMath::Cos(newphi);
1663     newVy = r * TMath::Sin(newphi);
1664     vZ   = iparticle->Vz(); // don't transform
1665     time = iparticle->T(); // don't transform
1666     
1667     pt = iparticle->Pt();
1668     newPx = pt * TMath::Cos(newphi);
1669     newPy = pt * TMath::Sin(newphi);
1670     pz = iparticle->Pz(); // don't transform
1671     e  = iparticle->Energy(); // don't transform
1672     
1673     // apply rotation 
1674     iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1675     iparticle->SetMomentum(newPx, newPy, pz, e);
1676     
1677   } //end particle loop 
1678   
1679   // now let's check that we put correctly the candidate photon in PHOS
1680   Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1681   Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax 
1682   if(IsInPHOS(phi,eta,-1)) 
1683     okdd = kTRUE;
1684   
1685   // reset the value for next event
1686   fPHOSRotateCandidate = -1;
1687   
1688 }
1689
1690
1691 Bool_t AliGenPythia::CheckDiffraction()
1692 {
1693   // use this method only with Perugia-0 tune!
1694
1695   //  printf("AAA\n");
1696
1697    Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1698
1699    Int_t iPart1=-1;
1700    Int_t iPart2=-1;
1701
1702    Double_t y1 = 1e10;
1703    Double_t y2 = -1e10;
1704
1705   const Int_t kNstable=20;
1706   const Int_t pdgStable[20] = {
1707     22,             // Photon
1708     11,             // Electron
1709     12,             // Electron Neutrino 
1710     13,             // Muon 
1711     14,             // Muon Neutrino
1712     15,             // Tau 
1713     16,             // Tau Neutrino
1714     211,            // Pion
1715     321,            // Kaon
1716     311,            // K0
1717     130,            // K0s
1718     310,            // K0l
1719     2212,           // Proton 
1720     2112,           // Neutron
1721     3122,           // Lambda_0
1722     3112,           // Sigma Minus
1723     3222,           // Sigma Plus
1724     3312,           // Xsi Minus 
1725     3322,           // Xsi0
1726     3334            // Omega
1727   };
1728     
1729      for (Int_t i = 0; i < np; i++) {
1730         TParticle *  part = (TParticle *) fParticles.At(i);
1731         
1732         Int_t statusCode = part->GetStatusCode();
1733         
1734         // Initial state particle
1735         if (statusCode != 1)
1736           continue;
1737
1738         Int_t pdg = TMath::Abs(part->GetPdgCode());
1739         Bool_t isStable = kFALSE;
1740         for (Int_t i1 = 0; i1 < kNstable; i1++) {
1741           if (pdg == pdgStable[i1]) {
1742             isStable = kTRUE;
1743             break;
1744           }
1745         }
1746         if(!isStable) 
1747           continue;
1748
1749         Double_t y = part->Y();
1750
1751         if (y < y1)
1752           {
1753             y1 = y;
1754             iPart1 = i;
1755           }
1756         if (y > y2)
1757         {
1758           y2 = y;
1759           iPart2 = i;
1760         }
1761      }
1762
1763      if(iPart1<0 || iPart2<0) return kFALSE;
1764
1765      y1=TMath::Abs(y1);
1766      y2=TMath::Abs(y2);
1767
1768      TParticle *  part1 = (TParticle *) fParticles.At(iPart1);
1769      TParticle *  part2 = (TParticle *) fParticles.At(iPart2);
1770
1771      Int_t pdg1 = part1->GetPdgCode();
1772      Int_t pdg2 = part2->GetPdgCode();
1773
1774
1775      Int_t iPart = -1;
1776      if (pdg1 == 2212 && pdg2 == 2212)
1777        {
1778          if(y1 > y2) 
1779            iPart = iPart1;
1780          else if(y1 < y2) 
1781            iPart = iPart2;
1782          else {
1783            iPart = iPart1;
1784            if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1785          }
1786        }
1787      else if (pdg1 == 2212)
1788        iPart = iPart1;
1789      else if (pdg2 == 2212)
1790        iPart = iPart2;
1791
1792
1793
1794
1795
1796      Double_t M=-1.;
1797      if(iPart>0) {
1798        TParticle *  part = (TParticle *) fParticles.At(iPart);
1799        Double_t E= part->Energy();
1800        Double_t P= part->P();
1801        M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1802      }
1803
1804      Double_t Mmin, Mmax, wSD, wDD, wND;
1805      if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1806
1807      if(M>-1 && M<Mmin) return kFALSE;
1808      if(M>Mmax) M=-1;
1809
1810      Int_t procType=fPythia->GetMSTI(1);
1811      Int_t proc0=2;
1812      if(procType== 94) proc0=1;
1813      if(procType== 92 || procType== 93) proc0=0;
1814
1815      Int_t proc=2;
1816      if(M>0) proc=0;
1817      else if(proc0==1) proc=1;
1818
1819      if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1820      if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1821
1822
1823     //     if(proc==1 || proc==2) return kFALSE;
1824
1825      if(proc!=0) {
1826        if(proc0!=0) fProcDiff = procType;
1827        else       fProcDiff = 95;
1828        return kTRUE;
1829      }
1830
1831     if(wSD<0)  AliError("wSD<0 ! \n");
1832
1833     if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1834
1835     //    printf("iPart = %d\n", iPart);
1836
1837     if(iPart==iPart1) fProcDiff=93;
1838     else if(iPart==iPart2) fProcDiff=92;
1839     else {
1840       printf("EROOR:  iPart!=iPart1 && iPart!=iPart2\n");
1841
1842     }
1843
1844     return kTRUE;
1845 }
1846
1847
1848
1849 Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax, 
1850                                                        Double_t &wSD, Double_t &wDD, Double_t &wND)
1851 {
1852
1853   // 900 GeV
1854   if(TMath::Abs(fEnergyCMS-900)<1 ){
1855
1856 const Int_t nbin=400;
1857 Double_t bin[]={
1858 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500, 
1859 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300, 
1860 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100, 
1861 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900, 
1862 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700, 
1863 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500, 
1864 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300, 
1865 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100, 
1866 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900, 
1867 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700, 
1868 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500, 
1869 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300, 
1870 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100, 
1871 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900, 
1872 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700, 
1873 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500, 
1874 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300, 
1875 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100, 
1876 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900, 
1877 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700, 
1878 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500, 
1879 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300, 
1880 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100, 
1881 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900, 
1882 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700, 
1883 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500, 
1884 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300, 
1885 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100, 
1886 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900, 
1887 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700, 
1888 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500, 
1889 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300, 
1890 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100, 
1891 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900, 
1892 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700, 
1893 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500, 
1894 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300, 
1895 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100, 
1896 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900, 
1897 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700, 
1898 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500, 
1899 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300, 
1900 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100, 
1901 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900, 
1902 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700, 
1903 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500, 
1904 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300, 
1905 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100, 
1906 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900, 
1907 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700, 
1908 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500, 
1909 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300, 
1910 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100, 
1911 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900, 
1912 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700, 
1913 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500, 
1914 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300, 
1915 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100, 
1916 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900, 
1917 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700, 
1918 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500, 
1919 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300, 
1920 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100, 
1921 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900, 
1922 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700, 
1923 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500, 
1924 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
1925 Double_t w[]={
1926 1.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002, 
1927 0.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285, 
1928 0.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754, 
1929 0.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686, 
1930 0.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483, 
1931 0.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964, 
1932 0.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759, 
1933 0.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836, 
1934 0.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836, 
1935 0.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688, 
1936 0.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331, 
1937 0.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728, 
1938 0.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921, 
1939 0.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763, 
1940 0.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208, 
1941 0.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992, 
1942 0.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822, 
1943 0.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597, 
1944 0.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908, 
1945 0.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008, 
1946 0.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658, 
1947 0.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051, 
1948 0.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882, 
1949 0.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734, 
1950 0.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476, 
1951 0.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245, 
1952 0.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417, 
1953 0.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507, 
1954 0.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188, 
1955 0.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408, 
1956 0.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228, 
1957 0.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830, 
1958 0.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991, 
1959 0.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279, 
1960 0.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818, 
1961 0.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686, 
1962 0.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242, 
1963 0.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969, 
1964 0.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382, 
1965 0.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880, 
1966 0.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040, 
1967 0.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072, 
1968 0.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978, 
1969 0.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000, 
1970 0.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155, 
1971 0.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093, 
1972 0.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390, 
1973 0.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977, 
1974 0.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739, 
1975 0.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220, 
1976 0.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314, 
1977 0.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331, 
1978 0.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391, 
1979 0.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989, 
1980 0.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348, 
1981 0.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153, 
1982 0.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453, 
1983 0.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558, 
1984 0.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521, 
1985 0.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577, 
1986 0.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585, 
1987 0.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980, 
1988 0.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873, 
1989 0.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928, 
1990 0.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462, 
1991 0.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252, 
1992 0.386112, 0.364314, 0.434375, 0.334629};
1993 wDD = 0.379611;
1994 wND = 0.496961;
1995 wSD = -1;
1996
1997     Mmin = bin[0];
1998     Mmax = bin[nbin];
1999     if(M<Mmin || M>Mmax) return kTRUE;
2000
2001     Int_t ibin=nbin-1;
2002     for(Int_t i=1; i<=nbin; i++) 
2003       if(M<=bin[i]) {
2004         ibin=i-1;
2005         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2006         break;
2007       }
2008     wSD=w[ibin];
2009     return kTRUE;
2010   }
2011  else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2012
2013 const Int_t nbin=400;
2014 Double_t bin[]={
2015 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500, 
2016 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300, 
2017 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100, 
2018 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900, 
2019 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700, 
2020 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500, 
2021 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300, 
2022 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100, 
2023 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900, 
2024 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700, 
2025 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500, 
2026 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300, 
2027 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100, 
2028 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900, 
2029 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700, 
2030 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500, 
2031 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300, 
2032 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100, 
2033 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900, 
2034 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700, 
2035 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500, 
2036 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300, 
2037 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100, 
2038 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900, 
2039 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700, 
2040 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500, 
2041 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300, 
2042 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100, 
2043 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900, 
2044 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700, 
2045 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500, 
2046 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300, 
2047 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100, 
2048 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900, 
2049 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700, 
2050 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500, 
2051 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300, 
2052 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100, 
2053 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900, 
2054 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700, 
2055 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500, 
2056 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300, 
2057 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100, 
2058 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900, 
2059 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700, 
2060 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500, 
2061 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300, 
2062 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100, 
2063 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900, 
2064 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700, 
2065 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500, 
2066 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300, 
2067 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100, 
2068 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900, 
2069 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700, 
2070 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500, 
2071 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300, 
2072 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100, 
2073 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900, 
2074 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700, 
2075 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500, 
2076 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300, 
2077 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100, 
2078 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900, 
2079 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700, 
2080 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500, 
2081 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2082 Double_t w[]={
2083 1.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723, 
2084 0.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705, 
2085 0.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186, 
2086 0.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638, 
2087 0.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431, 
2088 0.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138, 
2089 0.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077, 
2090 0.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291, 
2091 0.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975, 
2092 0.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411, 
2093 0.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911, 
2094 0.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867, 
2095 0.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675, 
2096 0.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140, 
2097 0.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678, 
2098 0.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345, 
2099 0.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044, 
2100 0.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527, 
2101 0.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399, 
2102 0.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695, 
2103 0.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323, 
2104 0.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478, 
2105 0.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321, 
2106 0.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223, 
2107 0.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186, 
2108 0.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548, 
2109 0.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368, 
2110 0.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152, 
2111 0.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471, 
2112 0.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872, 
2113 0.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062, 
2114 0.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696, 
2115 0.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455, 
2116 0.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111, 
2117 0.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490, 
2118 0.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975, 
2119 0.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739, 
2120 0.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472, 
2121 0.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139, 
2122 0.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843, 
2123 0.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598, 
2124 0.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948, 
2125 0.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487, 
2126 0.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762, 
2127 0.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245, 
2128 0.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926, 
2129 0.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985, 
2130 0.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870, 
2131 0.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446, 
2132 0.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229, 
2133 0.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778, 
2134 0.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398, 
2135 0.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936, 
2136 0.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324, 
2137 0.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554, 
2138 0.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403, 
2139 0.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683, 
2140 0.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038, 
2141 0.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310, 
2142 0.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109, 
2143 0.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974, 
2144 0.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506, 
2145 0.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351, 
2146 0.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870, 
2147 0.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336, 
2148 0.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611, 
2149 0.201712, 0.242225, 0.255565, 0.258738};
2150 wDD = 0.512813;
2151 wND = 0.518820;
2152 wSD = -1;
2153
2154     Mmin = bin[0];
2155     Mmax = bin[nbin];
2156     if(M<Mmin || M>Mmax) return kTRUE;
2157
2158     Int_t ibin=nbin-1;
2159     for(Int_t i=1; i<=nbin; i++) 
2160       if(M<=bin[i]) {
2161         ibin=i-1;
2162         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2163         break;
2164       }
2165     wSD=w[ibin];
2166     return kTRUE;
2167   }
2168
2169
2170   else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2171 const Int_t nbin=400;
2172 Double_t bin[]={
2173 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500, 
2174 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300, 
2175 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100, 
2176 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900, 
2177 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700, 
2178 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500, 
2179 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300, 
2180 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100, 
2181 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900, 
2182 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700, 
2183 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500, 
2184 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300, 
2185 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100, 
2186 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900, 
2187 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700, 
2188 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500, 
2189 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300, 
2190 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100, 
2191 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900, 
2192 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700, 
2193 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500, 
2194 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300, 
2195 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100, 
2196 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900, 
2197 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700, 
2198 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500, 
2199 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300, 
2200 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100, 
2201 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900, 
2202 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700, 
2203 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500, 
2204 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300, 
2205 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100, 
2206 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900, 
2207 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700, 
2208 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500, 
2209 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300, 
2210 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100, 
2211 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900, 
2212 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700, 
2213 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500, 
2214 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300, 
2215 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100, 
2216 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900, 
2217 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700, 
2218 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500, 
2219 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300, 
2220 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100, 
2221 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900, 
2222 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700, 
2223 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500, 
2224 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300, 
2225 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100, 
2226 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900, 
2227 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700, 
2228 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500, 
2229 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300, 
2230 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100, 
2231 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900, 
2232 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700, 
2233 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500, 
2234 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300, 
2235 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100, 
2236 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900, 
2237 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700, 
2238 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500, 
2239 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2240 Double_t w[]={
2241 1.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093, 
2242 0.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548, 
2243 0.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117, 
2244 0.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274, 
2245 0.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740, 
2246 0.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602, 
2247 0.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363, 
2248 0.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014, 
2249 0.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298, 
2250 0.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501, 
2251 0.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820, 
2252 0.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242, 
2253 0.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060, 
2254 0.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033, 
2255 0.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969, 
2256 0.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421, 
2257 0.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165, 
2258 0.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295, 
2259 0.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590, 
2260 0.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245, 
2261 0.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191, 
2262 0.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485, 
2263 0.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379, 
2264 0.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470, 
2265 0.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143, 
2266 0.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866, 
2267 0.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686, 
2268 0.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773, 
2269 0.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810, 
2270 0.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225, 
2271 0.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591, 
2272 0.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213, 
2273 0.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267, 
2274 0.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150, 
2275 0.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691, 
2276 0.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534, 
2277 0.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353, 
2278 0.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136, 
2279 0.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707, 
2280 0.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390, 
2281 0.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804, 
2282 0.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828, 
2283 0.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258, 
2284 0.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801, 
2285 0.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926, 
2286 0.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133, 
2287 0.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729, 
2288 0.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415, 
2289 0.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328, 
2290 0.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690, 
2291 0.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651, 
2292 0.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718, 
2293 0.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734, 
2294 0.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402, 
2295 0.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809, 
2296 0.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570, 
2297 0.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492, 
2298 0.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162, 
2299 0.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066, 
2300 0.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092, 
2301 0.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577, 
2302 0.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618, 
2303 0.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440, 
2304 0.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519, 
2305 0.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812, 
2306 0.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903, 
2307 0.175006, 0.223482, 0.202706, 0.218108};
2308 wDD = 0.207705;
2309 wND = 0.289628;
2310 wSD = -1;
2311
2312     Mmin = bin[0];
2313     Mmax = bin[nbin];
2314
2315     if(M<Mmin || M>Mmax) return kTRUE;
2316
2317     Int_t ibin=nbin-1;
2318     for(Int_t i=1; i<=nbin; i++) 
2319       if(M<=bin[i]) {
2320         ibin=i-1;
2321         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2322         break;
2323       }
2324     wSD=w[ibin];
2325     return kTRUE;
2326   }
2327
2328   return kFALSE;
2329 }
2330
2331 Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2332 {
2333 // Check if this is a heavy flavor decay product
2334   TParticle *  part = (TParticle *) fParticles.At(ipart);
2335   Int_t mpdg = TMath::Abs(part->GetPdgCode());
2336   Int_t mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2337   //
2338   // Light hadron
2339   if (mfl >= 4 && mfl < 6) return kTRUE;
2340   
2341   Int_t imo = part->GetFirstMother()-1;
2342   TParticle* pm = part;
2343   //
2344   // Heavy flavor hadron produced by generator
2345   while (imo >  -1) {
2346     pm  =  (TParticle*)fParticles.At(imo);
2347     mpdg = TMath::Abs(pm->GetPdgCode());
2348     mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2349     if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2350     imo = pm->GetFirstMother()-1;
2351   }
2352   return kFALSE;
2353 }
2354
2355 Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
2356 {
2357   // check the eta/phi correspond to the detectors acceptance
2358   // iparticle is the index of the particle to be checked, for PHOS rotation case
2359   if     (fCheckPHOS   && IsInPHOS  (phi,eta,iparticle)) return kTRUE;
2360   else if(fCheckEMCAL  && IsInEMCAL (phi,eta)) return kTRUE;
2361   else if(fCheckBarrel && IsInBarrel(    eta)) return kTRUE;
2362   else                                         return kFALSE;
2363 }
2364
2365 Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
2366 {
2367   // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2368   // implemented primaryly for kPyJets, but extended to any kind of process.
2369   
2370   //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2371   //       fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel); 
2372   
2373   Bool_t ok = kFALSE;
2374   for (Int_t i=0; i< np; i++) {
2375     
2376     TParticle* iparticle = (TParticle *) fParticles.At(i);
2377     
2378     Int_t pdg          = iparticle->GetPdgCode();
2379     Int_t status       = iparticle->GetStatusCode();
2380     Int_t imother      = iparticle->GetFirstMother() - 1;
2381     
2382     TParticle* pmother = 0x0;
2383     Int_t momStatus    = -1;
2384     Int_t momPdg       = -1;
2385     if(imother > 0 ){  
2386       pmother = (TParticle *) fParticles.At(imother);
2387       momStatus    = pmother->GetStatusCode();
2388       momPdg       = pmother->GetPdgCode();
2389     }
2390     
2391     ok = kFALSE;
2392     
2393     //
2394     // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2395     //
2396     // Hadron
2397     if (fHadronInCalo && status == 1)
2398     {
2399       if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0 
2400         // (in case neutral mesons were declared stable)
2401         ok = kTRUE;
2402     }
2403     //Electron
2404     else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2405     {
2406         ok = kTRUE;
2407     }
2408     //Fragmentation photon
2409     else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2410     {        
2411       if(momStatus != 11) ok = kTRUE ;  // No photon from hadron decay
2412     }
2413     // Decay photon
2414     else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2415     {              
2416       if( momStatus == 11)
2417       {
2418         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2419         //                                                   pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2420         ok = kTRUE ;  // photon from hadron decay
2421         
2422         //In case only decays from pi0 or eta requested
2423         if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2424         if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2425       }
2426
2427     }
2428     // Pi0 or Eta particle
2429     else if ((fPi0InCalo || fEtaInCalo))
2430     {
2431       if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2432       
2433       if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2434       {
2435         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2436         ok = kTRUE;
2437       }      
2438       else if (fEtaInCalo && pdg == 221) 
2439       {
2440         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2441         ok = kTRUE;
2442       }
2443       
2444     }// pi0 or eta
2445     
2446     //
2447     // Check that the selected particle is in the calorimeter acceptance
2448     //
2449     if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2450     {
2451       //Just check if the selected particle falls in the acceptance
2452       if(!fForceNeutralMeson2PhotonDecay )
2453       {
2454         //printf("\t Check acceptance! \n");
2455         Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2456         Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax      
2457         
2458         if(CheckDetectorAcceptance(phi,eta,i))
2459         {
2460           ok =kTRUE;
2461           AliDebug(1,Form("Selected trigger pdg %d, status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2462           //printf("\t Accept \n");
2463           break;
2464         }
2465         else ok = kFALSE;
2466       }
2467       //Mesons have several decay modes, select only those decaying into 2 photons
2468       else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2469       {
2470         // In case we want the pi0/eta trigger, 
2471         // check the decay mode (2 photons)
2472         
2473         //printf("\t Force decay 2 gamma\n");          
2474         
2475         Int_t ndaughters = iparticle->GetNDaughters();
2476         if(ndaughters != 2){
2477           ok=kFALSE;
2478           continue;
2479         }
2480         
2481         TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2482         TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2483         if(!d1 || !d2) {
2484           ok=kFALSE;
2485           continue;
2486         }
2487         
2488         //iparticle->Print();
2489         //d1->Print();
2490         //d2->Print();
2491         
2492         Int_t pdgD1 = d1->GetPdgCode();
2493         Int_t pdgD2 = d2->GetPdgCode();
2494         //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2495         //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2496         
2497         if(pdgD1 != 22  || pdgD2 != 22){ 
2498           ok = kFALSE;
2499           continue;
2500         }
2501         
2502         //printf("\t accept decay\n");
2503         
2504         //Trigger on the meson, not on the daughter
2505         if(!fDecayPhotonInCalo){    
2506           
2507           Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2508           Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax  
2509           
2510           if(CheckDetectorAcceptance(phi,eta,i))
2511           {
2512             //printf("\t Accept meson pdg %d\n",pdg);
2513             ok =kTRUE;
2514             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2515             break;
2516           } else {
2517             ok=kFALSE;
2518             continue;
2519           }
2520         }
2521         
2522         //printf("Check daughters acceptance\n");
2523         
2524         //Trigger on the meson daughters
2525         //Photon 1
2526         Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2527         Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax     
2528         if(d1->Pt() > fTriggerParticleMinPt)
2529         {
2530           //printf("\t Check acceptance photon 1! \n");
2531           if(CheckDetectorAcceptance(phi,eta,i))
2532           {
2533             //printf("\t Accept Photon 1\n");
2534             ok =kTRUE;
2535             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2536             break;
2537           }
2538           else ok = kFALSE;
2539         } // pt cut
2540         else  ok = kFALSE;
2541         
2542         //Photon 2
2543         phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2544         eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax   
2545         
2546         if(d2->Pt() > fTriggerParticleMinPt)
2547         {
2548           //printf("\t Check acceptance photon 2! \n");
2549           if(CheckDetectorAcceptance(phi,eta,i))
2550           {
2551             //printf("\t Accept Photon 2\n");
2552             ok =kTRUE;
2553             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2554             break;
2555           } 
2556           else ok = kFALSE;         
2557         } // pt cut
2558         else ok = kFALSE;
2559       } // force 2 photon daughters in pi0/eta decays
2560       else ok = kFALSE;
2561     } else ok = kFALSE; // check acceptance
2562   } // primary loop
2563     
2564   //
2565   // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2566   // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2567   //
2568   if(fCheckPHOSeta)
2569   {
2570     RotatePhi(ok);
2571   }
2572   
2573   return ok;
2574 }
2575