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