d7c462dc67a80d100c05803e408b6ebdbeb62c62
[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 "PyquenCommon.h"
48
49 ClassImp(AliGenPythia)
50
51
52 AliGenPythia::AliGenPythia():
53     AliGenMC(),
54     fProcess(kPyCharm),          
55     fStrucFunc(kCTEQ5L), 
56     fKineBias(0.),
57     fTrials(0),
58     fTrialsRun(0),
59     fQ(0.),
60     fX1(0.),
61     fX2(0.),
62     fEventTime(0.),
63     fInteractionRate(0.),
64     fTimeWindow(0.),
65     fCurSubEvent(0),
66     fEventsTime(0),
67     fNev(0),
68     fFlavorSelect(0),
69     fXsection(0.),
70     fPythia(0),
71     fPtHardMin(0.),
72     fPtHardMax(1.e4),
73     fYHardMin(-1.e10),
74     fYHardMax(1.e10),
75     fGinit(1),
76     fGfinal(1),
77     fHadronisation(1),
78     fNpartons(0),
79     fReadFromFile(0),
80     fQuench(0),
81     fQhat(0.),
82     fLength(0.),
83     fPtKick(1.),
84     fFullEvent(kTRUE),
85     fDecayer(new AliDecayerPythia()),
86     fDebugEventFirst(-1),
87     fDebugEventLast(-1),
88     fEtMinJet(0.),      
89     fEtMaxJet(1.e4),      
90     fEtaMinJet(-20.),     
91     fEtaMaxJet(20.),     
92     fPhiMinJet(0.),     
93     fPhiMaxJet(2.* TMath::Pi()),     
94     fJetReconstruction(kCell),
95     fEtaMinGamma(-20.),      
96     fEtaMaxGamma(20.),      
97     fPhiMinGamma(0.),      
98     fPhiMaxGamma(2. * TMath::Pi()),      
99     fPycellEtaMax(2.),     
100     fPycellNEta(274),       
101     fPycellNPhi(432),       
102     fPycellThreshold(0.),  
103     fPycellEtSeed(4.),     
104     fPycellMinEtJet(10.),  
105     fPycellMaxRadius(1.), 
106     fStackFillOpt(kFlavorSelection),   
107     fFeedDownOpt(kTRUE),    
108     fFragmentation(kTRUE),
109     fSetNuclei(kFALSE),
110     fNewMIS(kFALSE),   
111     fHFoff(kFALSE),    
112     fNucPdf(0),
113     fTriggerParticle(0),
114     fTriggerEta(0.9),     
115     fTriggerMultiplicity(0),
116     fTriggerMultiplicityEta(0),
117     fCountMode(kCountAll),      
118     fHeader(0),  
119     fRL(0),      
120     fFileName(0),
121     fFragPhotonInCalo(kFALSE),
122     fPi0InCalo(kFALSE) ,
123     fPhotonInCalo(kFALSE),
124     fEleInEMCAL(kFALSE),
125     fCheckEMCAL(kFALSE),
126     fCheckPHOS(kFALSE),
127     fCheckPHOSeta(kFALSE),
128     fFragPhotonOrPi0MinPt(0), 
129     fPhotonMinPt(0), 
130     fElectronMinPt(0), 
131     fPHOSMinPhi(219.),
132     fPHOSMaxPhi(321.),
133     fPHOSEta(0.13),
134     fEMCALMinPhi(79.),
135     fEMCALMaxPhi(191.),
136     fEMCALEta(0.71)
137
138 {
139 // Default Constructor
140   fEnergyCMS = 5500.;
141   SetNuclei(0,0);
142   if (!AliPythiaRndm::GetPythiaRandom()) 
143       AliPythiaRndm::SetPythiaRandom(GetRandom());
144 }
145
146 AliGenPythia::AliGenPythia(Int_t npart)
147     :AliGenMC(npart),
148      fProcess(kPyCharm),          
149      fStrucFunc(kCTEQ5L), 
150      fKineBias(0.),
151      fTrials(0),
152      fTrialsRun(0),
153      fQ(0.),
154      fX1(0.),
155      fX2(0.),
156      fEventTime(0.),
157      fInteractionRate(0.),
158      fTimeWindow(0.),
159      fCurSubEvent(0),
160      fEventsTime(0),
161      fNev(0),
162      fFlavorSelect(0),
163      fXsection(0.),
164      fPythia(0),
165      fPtHardMin(0.),
166      fPtHardMax(1.e4),
167      fYHardMin(-1.e10),
168      fYHardMax(1.e10),
169      fGinit(kTRUE),
170      fGfinal(kTRUE),
171      fHadronisation(kTRUE),
172      fNpartons(0),
173      fReadFromFile(kFALSE),
174      fQuench(kFALSE),
175      fQhat(0.),
176      fLength(0.),
177      fPtKick(1.),
178      fFullEvent(kTRUE),
179      fDecayer(new AliDecayerPythia()),
180      fDebugEventFirst(-1),
181      fDebugEventLast(-1),
182      fEtMinJet(0.),      
183      fEtMaxJet(1.e4),      
184      fEtaMinJet(-20.),     
185      fEtaMaxJet(20.),     
186      fPhiMinJet(0.),     
187      fPhiMaxJet(2.* TMath::Pi()),     
188      fJetReconstruction(kCell),
189      fEtaMinGamma(-20.),      
190      fEtaMaxGamma(20.),      
191      fPhiMinGamma(0.),      
192      fPhiMaxGamma(2. * TMath::Pi()),      
193      fPycellEtaMax(2.),     
194      fPycellNEta(274),       
195      fPycellNPhi(432),       
196      fPycellThreshold(0.),  
197      fPycellEtSeed(4.),     
198      fPycellMinEtJet(10.),  
199      fPycellMaxRadius(1.), 
200      fStackFillOpt(kFlavorSelection),   
201      fFeedDownOpt(kTRUE),    
202      fFragmentation(kTRUE),
203      fSetNuclei(kFALSE),
204      fNewMIS(kFALSE),   
205      fHFoff(kFALSE),    
206      fNucPdf(0),
207      fTriggerParticle(0),
208      fTriggerEta(0.9),     
209      fTriggerMultiplicity(0),
210      fTriggerMultiplicityEta(0),
211      fCountMode(kCountAll),      
212      fHeader(0),  
213      fRL(0),      
214      fFileName(0),
215      fFragPhotonInCalo(kFALSE),
216      fPi0InCalo(kFALSE) ,
217      fPhotonInCalo(kFALSE),
218      fEleInEMCAL(kFALSE),
219      fCheckEMCAL(kFALSE),
220      fCheckPHOS(kFALSE),
221      fCheckPHOSeta(kFALSE),
222      fFragPhotonOrPi0MinPt(0),
223      fPhotonMinPt(0),
224      fElectronMinPt(0),
225      fPHOSMinPhi(219.),
226      fPHOSMaxPhi(321.),
227      fPHOSEta(0.13),
228      fEMCALMinPhi(79.),
229      fEMCALMaxPhi(191.),
230      fEMCALEta(0.71)
231 {
232 // default charm production at 5. 5 TeV
233 // semimuonic decay
234 // structure function GRVHO
235 //
236     fEnergyCMS = 5500.;
237     fName = "Pythia";
238     fTitle= "Particle Generator using PYTHIA";
239     SetForceDecay();
240     // Set random number generator 
241     if (!AliPythiaRndm::GetPythiaRandom()) 
242       AliPythiaRndm::SetPythiaRandom(GetRandom());
243     SetNuclei(0,0);
244  }
245
246 AliGenPythia::~AliGenPythia()
247 {
248 // Destructor
249   if(fEventsTime) delete fEventsTime;
250 }
251
252 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
253 {
254 // Generate pileup using user specified rate
255     fInteractionRate = rate;
256     fTimeWindow = timewindow;
257     GeneratePileup();
258 }
259
260 void AliGenPythia::GeneratePileup()
261 {
262 // Generate sub events time for pileup
263     fEventsTime = 0;
264     if(fInteractionRate == 0.) {
265       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
266       return;
267     }
268
269     Int_t npart = NumberParticles();
270     if(npart < 0) {
271       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
272       return;
273     }
274
275     if(fEventsTime) delete fEventsTime;
276     fEventsTime = new TArrayF(npart);
277     TArrayF &array = *fEventsTime;
278     for(Int_t ipart = 0; ipart < npart; ipart++)
279       array[ipart] = 0.;
280
281     Float_t eventtime = 0.;
282     while(1)
283       {
284         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
285         if(eventtime > fTimeWindow) break;
286         array.Set(array.GetSize()+1);
287         array[array.GetSize()-1] = eventtime;
288       }
289
290     eventtime = 0.;
291     while(1)
292       {
293         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
294         if(TMath::Abs(eventtime) > fTimeWindow) break;
295         array.Set(array.GetSize()+1);
296         array[array.GetSize()-1] = eventtime;
297       }
298
299     SetNumberParticles(fEventsTime->GetSize());
300 }
301
302 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
303                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
304 {
305 // Set pycell parameters
306     fPycellEtaMax    =  etamax;
307     fPycellNEta      =  neta;
308     fPycellNPhi      =  nphi;
309     fPycellThreshold =  thresh;
310     fPycellEtSeed    =  etseed;
311     fPycellMinEtJet  =  minet;
312     fPycellMaxRadius =  r;
313 }
314
315
316
317 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
318 {
319   // Set a range of event numbers, for which a table
320   // of generated particle will be printed
321   fDebugEventFirst = eventFirst;
322   fDebugEventLast  = eventLast;
323   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
324 }
325
326 void AliGenPythia::Init()
327 {
328 // Initialisation
329     
330     SetMC(AliPythia::Instance());
331     fPythia=(AliPythia*) fMCEvGen;
332     
333 //
334     fParentWeight=1./Float_t(fNpart);
335 //
336
337
338     fPythia->SetCKIN(3,fPtHardMin);
339     fPythia->SetCKIN(4,fPtHardMax);
340     fPythia->SetCKIN(7,fYHardMin);
341     fPythia->SetCKIN(8,fYHardMax);
342     
343     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);  
344     // Fragmentation?
345     if (fFragmentation) {
346       fPythia->SetMSTP(111,1);
347     } else {
348       fPythia->SetMSTP(111,0);
349     }
350
351
352 //  initial state radiation   
353     fPythia->SetMSTP(61,fGinit);
354 //  final state radiation
355     fPythia->SetMSTP(71,fGfinal);
356 //  pt - kick
357     if (fPtKick > 0.) {
358         fPythia->SetMSTP(91,1);
359         fPythia->SetPARP(91,fPtKick);   
360         fPythia->SetPARP(93, 4. * fPtKick);
361     } else {
362         fPythia->SetMSTP(91,0);
363     }
364
365
366     if (fReadFromFile) {
367         fRL  =  AliRunLoader::Open(fFileName, "Partons");
368         fRL->LoadKinematics();
369         fRL->LoadHeader();
370     } else {
371         fRL = 0x0;
372     }
373  //
374     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
375     //  Forward Paramters to the AliPythia object
376     fDecayer->SetForceDecay(fForceDecay);    
377 // Switch off Heavy Flavors on request  
378     if (fHFoff) {
379         // Maximum number of quark flavours used in pdf 
380         fPythia->SetMSTP(58, 3);
381         // Maximum number of flavors that can be used in showers
382         fPythia->SetMSTJ(45, 3);        
383         // Switch off g->QQbar splitting in decay table
384         ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
385     }
386
387     fDecayer->Init();
388
389
390 //  Parent and Children Selection
391     switch (fProcess) 
392     {
393     case kPyOldUEQ2ordered:
394     case kPyOldUEQ2ordered2:
395     case kPyOldPopcorn:
396       break;
397     case kPyCharm:
398     case kPyCharmUnforced:
399     case kPyCharmPbPbMNR:
400     case kPyCharmpPbMNR:
401     case kPyCharmppMNR:
402     case kPyCharmppMNRwmi:
403         fParentSelect[0] =   411;
404         fParentSelect[1] =   421;
405         fParentSelect[2] =   431;
406         fParentSelect[3] =  4122;
407         fParentSelect[4] =  4232;
408         fParentSelect[5] =  4132;
409         fParentSelect[6] =  4332;
410         fFlavorSelect    =  4;  
411         break;
412     case kPyD0PbPbMNR:
413     case kPyD0pPbMNR:
414     case kPyD0ppMNR:
415         fParentSelect[0] =   421;
416         fFlavorSelect    =   4; 
417         break;
418     case kPyDPlusPbPbMNR:
419     case kPyDPluspPbMNR:
420     case kPyDPlusppMNR:
421         fParentSelect[0] =   411;
422         fFlavorSelect    =   4; 
423         break;
424     case kPyDPlusStrangePbPbMNR:
425     case kPyDPlusStrangepPbMNR:
426     case kPyDPlusStrangeppMNR:
427         fParentSelect[0] =   431;
428         fFlavorSelect    =   4; 
429         break;
430     case kPyBeauty:
431     case kPyBeautyJets:
432     case kPyBeautyPbPbMNR:
433     case kPyBeautypPbMNR:
434     case kPyBeautyppMNR:
435     case kPyBeautyppMNRwmi:
436         fParentSelect[0]=  511;
437         fParentSelect[1]=  521;
438         fParentSelect[2]=  531;
439         fParentSelect[3]= 5122;
440         fParentSelect[4]= 5132;
441         fParentSelect[5]= 5232;
442         fParentSelect[6]= 5332;
443         fFlavorSelect   = 5;    
444         break;
445     case kPyBeautyUnforced:
446         fParentSelect[0] =  511;
447         fParentSelect[1] =  521;
448         fParentSelect[2] =  531;
449         fParentSelect[3] = 5122;
450         fParentSelect[4] = 5132;
451         fParentSelect[5] = 5232;
452         fParentSelect[6] = 5332;
453         fFlavorSelect    = 5;   
454         break;
455     case kPyJpsiChi:
456     case kPyJpsi:
457         fParentSelect[0] = 443;
458         break;
459     case kPyMbDefault:
460     case kPyMb:
461     case kPyMbWithDirectPhoton:
462     case kPyMbNonDiffr:
463     case kPyMbMSEL1:
464     case kPyJets:
465     case kPyDirectGamma:
466     case kPyLhwgMb:     
467         break;
468     case kPyW:
469     case kPyZ:
470         break;
471     }
472 //
473 //
474 //  JetFinder for Trigger
475 //
476 //  Configure detector (EMCAL like)
477 //
478     fPythia->SetPARU(51, fPycellEtaMax);
479     fPythia->SetMSTU(51, fPycellNEta);
480     fPythia->SetMSTU(52, fPycellNPhi);
481 //
482 //  Configure Jet Finder
483 //  
484     fPythia->SetPARU(58,  fPycellThreshold);
485     fPythia->SetPARU(52,  fPycellEtSeed);
486     fPythia->SetPARU(53,  fPycellMinEtJet);
487     fPythia->SetPARU(54,  fPycellMaxRadius);
488     fPythia->SetMSTU(54,  2);
489 //
490 //  This counts the total number of calls to Pyevnt() per run.
491     fTrialsRun = 0;
492     fQ         = 0.;
493     fX1        = 0.;
494     fX2        = 0.;    
495     fNev       = 0 ;
496 //    
497 //
498 //
499     AliGenMC::Init();
500 //
501 //
502 //  
503     if (fSetNuclei) {
504         fDyBoost = 0;
505         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
506     }
507     
508     fPythia->SetPARJ(200, 0.0);
509     fPythia->SetPARJ(199, 0.0);
510     fPythia->SetPARJ(198, 0.0);
511     fPythia->SetPARJ(197, 0.0);
512
513     if (fQuench == 1) {
514         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
515     }
516
517     if (fQuench == 3) {
518         // Nestor's change of the splittings
519         fPythia->SetPARJ(200, 0.8);
520         fPythia->SetMSTJ(41, 1);  // QCD radiation only
521         fPythia->SetMSTJ(42, 2);  // angular ordering
522         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
523         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
524         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
525         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
526     } else if (fQuench == 4) {
527         // Armesto-Cunqueiro-Salgado change of the splittings.
528         AliFastGlauber* glauber = AliFastGlauber::Instance();
529         glauber->Init(2);
530         //read and store transverse almonds corresponding to differnt
531         //impact parameters.
532         glauber->SetCentralityClass(0.,0.1);
533         fPythia->SetPARJ(200, 1.);
534         fPythia->SetPARJ(198, fQhat);
535         fPythia->SetPARJ(199, fLength);
536         
537         fPythia->SetMSTJ(41, 1);  // QCD radiation only
538         fPythia->SetMSTJ(42, 2);  // angular ordering
539         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
540         //fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
541         //fPythia->SetMSTJ(50, 0);  // No coherence in first branching
542         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
543         //    MSTJ(41) must NOT be 11 or 12, as then FSR may go through PYPTFS
544         //   (kt-ordered cascade) in which medium effects have not been introduced.
545     }
546 }
547
548 void AliGenPythia::Generate()
549 {
550 // Generate one event
551     if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
552     fDecayer->ForceDecay();
553
554     Float_t polar[3]   =   {0,0,0};
555     Float_t origin[3]  =   {0,0,0};
556     Float_t p[4];
557 //  converts from mm/c to s
558     const Float_t kconv=0.001/2.999792458e8;
559 //
560     Int_t nt=0;
561     Int_t jev=0;
562     Int_t j, kf;
563     fTrials=0;
564     fEventTime = 0.;
565     
566     
567
568     //  Set collision vertex position 
569     if (fVertexSmear == kPerEvent) Vertex();
570     
571 //  event loop    
572     while(1)
573     {
574 //
575 // Produce event
576 //
577 //
578 // Switch hadronisation off
579 //
580         fPythia->SetMSTJ(1, 0);
581
582         if (fQuench ==4){
583             Double_t bimp;
584             // Quenching comes through medium-modified splitting functions.
585             AliFastGlauber::Instance()->GetRandomBHard(bimp);
586             fPythia->SetPARJ(197,bimp);
587         } 
588 //
589 // Either produce new event or read partons from file
590 //      
591         if (!fReadFromFile) {
592             if (!fNewMIS) {
593                 fPythia->Pyevnt();
594             } else {
595                 fPythia->Pyevnw();
596             }
597             fNpartons = fPythia->GetN();
598         } else {
599             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
600             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
601             fPythia->SetN(0);
602             LoadEvent(fRL->Stack(), 0 , 1);
603             fPythia->Pyedit(21);
604         }
605         
606 //
607 //  Run quenching routine 
608 //
609         if (fQuench == 1) {
610             fPythia->Quench();
611         } else if (fQuench == 2){
612             fPythia->Pyquen(208., 0, 0.);
613         } else if (fQuench == 3) {
614             // Quenching is via multiplicative correction of the splittings
615         }
616         
617 //
618 // Switch hadronisation on
619 //
620         if (fHadronisation) {
621             fPythia->SetMSTJ(1, 1);
622 //
623 // .. and perform hadronisation
624 //      printf("Calling hadronisation %d\n", fPythia->GetN());
625             fPythia->Pyexec();  
626         }
627         fTrials++;
628         fPythia->ImportParticles(&fParticles,"All");
629         Boost();
630 //
631 //
632 //
633         Int_t i;
634         
635         fNprimaries = 0;
636         Int_t np = fParticles.GetEntriesFast();
637         
638         if (np == 0) continue;
639 //
640         
641 //
642         Int_t* pParent   = new Int_t[np];
643         Int_t* pSelected = new Int_t[np];
644         Int_t* trackIt   = new Int_t[np];
645         for (i = 0; i < np; i++) {
646             pParent[i]   = -1;
647             pSelected[i] =  0;
648             trackIt[i]   =  0;
649         }
650
651         Int_t nc = 0;        // Total n. of selected particles
652         Int_t nParents = 0;  // Selected parents
653         Int_t nTkbles = 0;   // Trackable particles
654         if (fProcess != kPyMbDefault && 
655             fProcess != kPyMb && 
656             fProcess != kPyMbWithDirectPhoton && 
657             fProcess != kPyJets && 
658             fProcess != kPyDirectGamma &&
659             fProcess != kPyMbNonDiffr  &&
660             fProcess != kPyMbMSEL1     &&
661             fProcess != kPyW && 
662             fProcess != kPyZ &&
663             fProcess != kPyCharmppMNRwmi && 
664             fProcess != kPyBeautyppMNRwmi &&
665             fProcess != kPyBeautyJets) {
666             
667             for (i = 0; i < np; i++) {
668                 TParticle* iparticle = (TParticle *) fParticles.At(i);
669                 Int_t ks = iparticle->GetStatusCode();
670                 kf = CheckPDGCode(iparticle->GetPdgCode());
671 // No initial state partons
672                 if (ks==21) continue;
673 //
674 // Heavy Flavor Selection
675 //
676                 // quark ?
677                 kf = TMath::Abs(kf);
678                 Int_t kfl = kf;
679                 // Resonance
680
681                 if (kfl > 100000) kfl %= 100000;
682                 if (kfl > 10000)  kfl %= 10000;
683                 // meson ?
684                 if  (kfl > 10) kfl/=100;
685                 // baryon
686                 if (kfl > 10) kfl/=10;
687                 Int_t ipa = iparticle->GetFirstMother()-1;
688                 Int_t kfMo = 0;
689 //
690 // Establish mother daughter relation between heavy quarks and mesons
691 //
692                 if (kf >= fFlavorSelect && kf <= 6) {
693                     Int_t idau = iparticle->GetFirstDaughter() - 1;
694                     if (idau > -1) {
695                         TParticle* daughter = (TParticle *) fParticles.At(idau);
696                         Int_t pdgD = daughter->GetPdgCode();
697                         if (pdgD == 91 || pdgD == 92) {
698                             Int_t jmin = daughter->GetFirstDaughter() - 1;
699                             Int_t jmax = daughter->GetLastDaughter()  - 1;                          
700                             for (Int_t jp = jmin; jp <= jmax; jp++)
701                                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
702                         } // is string or cluster
703                     } // has daughter
704                 } // heavy quark
705                 
706
707                 if (ipa > -1) {
708                     TParticle *  mother = (TParticle *) fParticles.At(ipa);
709                     kfMo = TMath::Abs(mother->GetPdgCode());
710                 }
711                 
712                 // What to keep in Stack?
713                 Bool_t flavorOK = kFALSE;
714                 Bool_t selectOK = kFALSE;
715                 if (fFeedDownOpt) {
716                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
717                 } else {
718                     if (kfl > fFlavorSelect) {
719                         nc = -1;
720                         break;
721                     }
722                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
723                 }
724                 switch (fStackFillOpt) {
725                 case kFlavorSelection:
726                     selectOK = kTRUE;
727                     break;
728                 case kParentSelection:
729                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
730                     break;
731                 }
732                 if (flavorOK && selectOK) { 
733 //
734 // Heavy flavor hadron or quark
735 //
736 // Kinematic seletion on final state heavy flavor mesons
737                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
738                     {
739                         continue;
740                     }
741                     pSelected[i] = 1;
742                     if (ParentSelected(kf)) ++nParents; // Update parent count
743 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
744                 } else {
745 // Kinematic seletion on decay products
746                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
747                         && !KinematicSelection(iparticle, 1)) 
748                     {
749                         continue;
750                     }
751 //
752 // Decay products 
753 // Select if mother was selected and is not tracked
754
755                     if (pSelected[ipa] && 
756                         !trackIt[ipa]  &&     // mother will be  tracked ?
757                         kfMo !=  5 &&         // mother is b-quark, don't store fragments          
758                         kfMo !=  4 &&         // mother is c-quark, don't store fragments 
759                         kf   != 92)           // don't store string
760                     {
761 //
762 // Semi-stable or de-selected: diselect decay products:
763 // 
764 //
765                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
766                         {
767                             Int_t ipF = iparticle->GetFirstDaughter();
768                             Int_t ipL = iparticle->GetLastDaughter();   
769                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
770                         }
771 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
772                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
773                     }
774                 }
775                 if (pSelected[i] == -1) pSelected[i] = 0;
776                 if (!pSelected[i]) continue;
777                 // Count quarks only if you did not include fragmentation
778                 if (fFragmentation && kf <= 10) continue;
779
780                 nc++;
781 // Decision on tracking
782                 trackIt[i] = 0;
783 //
784 // Track final state particle
785                 if (ks == 1) trackIt[i] = 1;
786 // Track semi-stable particles
787                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
788 // Track particles selected by process if undecayed. 
789                 if (fForceDecay == kNoDecay) {
790                     if (ParentSelected(kf)) trackIt[i] = 1;
791                 } else {
792                     if (ParentSelected(kf)) trackIt[i] = 0;
793                 }
794                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
795 //
796 //
797
798             } // particle selection loop
799             if (nc > 0) {
800                 for (i = 0; i<np; i++) {
801                     if (!pSelected[i]) continue;
802                     TParticle *  iparticle = (TParticle *) fParticles.At(i);
803                     kf = CheckPDGCode(iparticle->GetPdgCode());
804                     Int_t ks = iparticle->GetStatusCode();  
805                     p[0] = iparticle->Px();
806                     p[1] = iparticle->Py();
807                     p[2] = iparticle->Pz();
808                     p[3] = iparticle->Energy();
809                     
810                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
811                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
812                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
813                     
814                     Float_t tof   = kconv*iparticle->T();
815                     Int_t ipa     = iparticle->GetFirstMother()-1;
816                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
817  
818                     PushTrack(fTrackIt*trackIt[i], iparent, kf, 
819                               p[0], p[1], p[2], p[3], 
820                               origin[0], origin[1], origin[2], tof, 
821                               polar[0], polar[1], polar[2],
822                               kPPrimary, nt, 1., ks);
823                     pParent[i] = nt;
824                     KeepTrack(nt);
825                     fNprimaries++;
826                 } //  PushTrack loop
827             }
828         } else {
829             nc = GenerateMB();
830         } // mb ?
831         
832         GetSubEventTime();
833
834         delete[] pParent;
835         delete[] pSelected;
836         delete[] trackIt;
837
838         if (nc > 0) {
839           switch (fCountMode) {
840           case kCountAll:
841             // printf(" Count all \n");
842             jev += nc;
843             break;
844           case kCountParents:
845             // printf(" Count parents \n");
846             jev += nParents;
847             break;
848           case kCountTrackables:
849             // printf(" Count trackable \n");
850             jev += nTkbles;
851             break;
852           }
853             if (jev >= fNpart || fNpart == -1) {
854                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
855                 
856                 fQ  += fPythia->GetVINT(51);
857                 fX1 += fPythia->GetVINT(41);
858                 fX2 += fPythia->GetVINT(42);
859                 fTrialsRun += fTrials;
860                 fNev++;
861                 MakeHeader();
862                 break;
863             }
864         }
865     } // event loop
866     SetHighWaterMark(nt);
867 //  adjust weight due to kinematic selection
868 //    AdjustWeights();
869 //  get cross-section
870     fXsection=fPythia->GetPARI(1);
871 }
872
873 Int_t  AliGenPythia::GenerateMB()
874 {
875 //
876 // Min Bias selection and other global selections
877 //
878     Int_t i, kf, nt, iparent;
879     Int_t nc = 0;
880     Float_t p[4];
881     Float_t polar[3]   =   {0,0,0};
882     Float_t origin[3]  =   {0,0,0};
883 //  converts from mm/c to s
884     const Float_t kconv=0.001/2.999792458e8;
885     
886
887     
888     Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
889
890
891
892     Int_t* pParent = new Int_t[np];
893     for (i=0; i< np; i++) pParent[i] = -1;
894     //
895     //TO BE CHECKED:  Should we require this for Beauty Jets?
896     //
897     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets) {
898         TParticle* jet1 = (TParticle *) fParticles.At(6);
899         TParticle* jet2 = (TParticle *) fParticles.At(7);
900         if (!CheckTrigger(jet1, jet2)) {
901           delete [] pParent;
902           return 0;
903         }
904     }
905
906     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
907     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
908
909       Bool_t ok = kFALSE;
910
911       Int_t pdg  = 0; 
912       if (fFragPhotonInCalo) pdg = 22   ; // Photon
913       else if (fPi0InCalo) pdg = 111 ;    // Pi0
914
915       for (i=0; i< np; i++) {
916         TParticle* iparticle = (TParticle *) fParticles.At(i);
917         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg && 
918            iparticle->Pt() > fFragPhotonOrPi0MinPt){
919           Int_t imother = iparticle->GetFirstMother() - 1;
920           TParticle* pmother = (TParticle *) fParticles.At(imother);
921           if(pdg == 111 || 
922              (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
923             {
924               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
925               Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax        
926               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
927                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
928                 ok =kTRUE;
929             }
930         }
931       }
932       if(!ok)
933         return 0;
934     }
935
936     // Select beauty jets with electron in EMCAL
937     if (fProcess == kPyBeautyJets && fEleInEMCAL) {
938
939       Bool_t ok = kFALSE;
940
941       Int_t pdg  = 11; //electron
942
943       Float_t pt, eta, phi;
944       for (i=0; i< np; i++) {
945         TParticle* iparticle = (TParticle *) fParticles.At(i);
946         if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg && 
947            iparticle->Pt() > fElectronMinPt){
948           pt = iparticle->Pt();
949           phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
950           eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax    
951           if(IsInEMCAL(phi,eta))
952             ok =kTRUE;
953         }
954       }
955       if(!ok)
956         return 0;
957       AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
958     }
959     
960     // Check for minimum multiplicity
961     if (fTriggerMultiplicity > 0) {
962       Int_t multiplicity = 0;
963       for (i = 0; i < np; i++) {
964         TParticle *  iparticle = (TParticle *) fParticles.At(i);
965         
966         Int_t statusCode = iparticle->GetStatusCode();
967         
968         // Initial state particle
969         if (statusCode > 20)
970           continue;
971         
972         // skip quarks and gluons
973         Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
974         if (pdgCode <= 10 || pdgCode == 21)
975           continue;
976         
977         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
978           continue;
979         
980         TParticlePDG* pdgPart = iparticle->GetPDG();
981         if (pdgPart && pdgPart->Charge() == 0)
982           continue;
983         
984         ++multiplicity;
985       }
986
987       if (multiplicity < fTriggerMultiplicity) {
988         delete [] pParent;
989         return 0;
990       }
991       
992       Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
993     }    
994     
995      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
996     if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
997
998       Bool_t okd = kFALSE;
999
1000       Int_t pdg  = 22; 
1001       Int_t iphcand = -1;
1002       for (i=0; i< np; i++) {
1003          TParticle* iparticle = (TParticle *) fParticles.At(i);
1004          Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1005          Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax 
1006          
1007          if(iparticle->GetStatusCode() == 1 
1008             && iparticle->GetPdgCode() == pdg   
1009             && iparticle->Pt() > fPhotonMinPt    
1010             && eta < fPHOSEta){                 
1011             
1012             // first check if the photon is in PHOS phi
1013             if(IsInPHOS(phi,eta)){ 
1014                 okd = kTRUE;
1015                 break;
1016             } 
1017             if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1018              
1019          }
1020       }
1021       
1022       if(!okd && iphcand != -1) // execute rotation in phi 
1023           RotatePhi(iphcand,okd);
1024       
1025       if(!okd)
1026         return 0;
1027     }
1028     
1029     if (fTriggerParticle) {
1030         Bool_t triggered = kFALSE;
1031         for (i = 0; i < np; i++) {
1032             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1033             kf = CheckPDGCode(iparticle->GetPdgCode());
1034             if (kf != fTriggerParticle) continue;
1035             if (iparticle->Pt() == 0.) continue;
1036             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1037             triggered = kTRUE;
1038             break;
1039         }
1040         if (!triggered) {
1041           delete [] pParent;
1042           return 0;
1043         }
1044     }
1045         
1046
1047     // Check if there is a ccbar or bbbar pair with at least one of the two
1048     // in fYMin < y < fYMax
1049     //
1050     // TO BE CHECKED:  Should we require this for beauty jets?
1051     //
1052     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1053       TParticle *partCheck;
1054       TParticle *mother;
1055       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1056       Bool_t  theChild=kFALSE;
1057       Float_t y;  
1058       Int_t   pdg,mpdg,mpdgUpperFamily;
1059       for(i=0; i<np; i++) {
1060         partCheck = (TParticle*)fParticles.At(i);
1061         pdg = partCheck->GetPdgCode();  
1062         if(TMath::Abs(pdg) == fFlavorSelect) { // quark  
1063           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1064           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1065                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
1066           if(y>fYMin && y<fYMax) inYcut=kTRUE;
1067         }
1068         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1069           Int_t mi = partCheck->GetFirstMother() - 1;
1070           if(mi<0) continue;
1071           mother = (TParticle*)fParticles.At(mi);
1072           mpdg=TMath::Abs(mother->GetPdgCode());
1073           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1074           if ( ParentSelected(mpdg) || 
1075               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1076             if (KinematicSelection(partCheck,1)) {
1077               theChild=kTRUE;
1078             }
1079           }
1080         }
1081       }
1082       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1083         delete[] pParent;
1084         return 0;
1085       }
1086       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1087         delete[] pParent;
1088         return 0;       
1089       }
1090
1091     }
1092
1093     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1094     if ( (fProcess == kPyW ||
1095           fProcess == kPyZ ||
1096           fProcess == kPyMbDefault ||
1097           fProcess == kPyMb ||
1098           fProcess == kPyMbWithDirectPhoton ||
1099           fProcess == kPyMbNonDiffr)  
1100          && (fCutOnChild == 1) ) {
1101       if ( !CheckKinematicsOnChild() ) {
1102         delete[] pParent;
1103         return 0;
1104       }
1105     }
1106   
1107
1108     for (i = 0; i < np; i++) {
1109         Int_t trackIt = 0;
1110         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1111         kf = CheckPDGCode(iparticle->GetPdgCode());
1112         Int_t ks = iparticle->GetStatusCode();
1113         Int_t km = iparticle->GetFirstMother();
1114         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
1115             (ks != 1) ||
1116             //
1117             //TO BE CHECKED:  Should we require this for beauty jets?
1118             //
1119             ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
1120             nc++;
1121             if (ks == 1) trackIt = 1;
1122             Int_t ipa = iparticle->GetFirstMother()-1;
1123             
1124             iparent = (ipa > -1) ? pParent[ipa] : -1;
1125             
1126 //
1127 // store track information
1128             p[0] = iparticle->Px();
1129             p[1] = iparticle->Py();
1130             p[2] = iparticle->Pz();
1131             p[3] = iparticle->Energy();
1132
1133             
1134             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1135             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1136             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1137             
1138             Float_t tof = fEventTime + kconv * iparticle->T();
1139
1140             PushTrack(fTrackIt*trackIt, iparent, kf, 
1141                       p[0], p[1], p[2], p[3], 
1142                       origin[0], origin[1], origin[2], tof, 
1143                       polar[0], polar[1], polar[2],
1144                       kPPrimary, nt, 1., ks);
1145             fNprimaries++;
1146             KeepTrack(nt);
1147             pParent[i] = nt;
1148             SetHighWaterMark(nt);
1149             
1150         } // select particle
1151     } // particle loop 
1152
1153     delete[] pParent;
1154     
1155     return 1;
1156 }
1157
1158
1159 void AliGenPythia::FinishRun()
1160 {
1161 // Print x-section summary
1162     fPythia->Pystat(1);
1163
1164     if (fNev > 0.) {
1165         fQ  /= fNev;
1166         fX1 /= fNev;
1167         fX2 /= fNev;    
1168     }
1169     
1170     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1171     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1172 }
1173
1174 void AliGenPythia::AdjustWeights() const
1175 {
1176 // Adjust the weights after generation of all events
1177 //
1178     if (gAlice) {
1179         TParticle *part;
1180         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1181         for (Int_t i=0; i<ntrack; i++) {
1182             part= gAlice->GetMCApp()->Particle(i);
1183             part->SetWeight(part->GetWeight()*fKineBias);
1184         }
1185     }
1186 }
1187     
1188 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1189 {
1190 // Treat protons as inside nuclei with mass numbers a1 and a2  
1191
1192     fAProjectile = a1;
1193     fATarget     = a2;
1194     fNucPdf      = pdfset;  // 0 EKS98 1 EPS08
1195     fSetNuclei   = kTRUE;
1196 }
1197
1198
1199 void AliGenPythia::MakeHeader()
1200 {
1201 //
1202 // Make header for the simulated event
1203 // 
1204   if (gAlice) {
1205     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1206         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1207   }
1208
1209 // Builds the event header, to be called after each event
1210     if (fHeader) delete fHeader;
1211     fHeader = new AliGenPythiaEventHeader("Pythia");
1212 //
1213 // Event type  
1214     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1215 //
1216 // Number of trials
1217     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1218 //
1219 // Event Vertex 
1220     fHeader->SetPrimaryVertex(fVertex);
1221     
1222 //
1223 // Number of primaries
1224     fHeader->SetNProduced(fNprimaries);
1225 //
1226 // Jets that have triggered
1227
1228     //Need to store jets for b-jet studies too!
1229     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets)
1230     {
1231         Int_t ntrig, njet;
1232         Float_t jets[4][10];
1233         GetJets(njet, ntrig, jets);
1234
1235         
1236         for (Int_t i = 0; i < ntrig; i++) {
1237             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
1238                                                         jets[3][i]);
1239         }
1240     }
1241 //
1242 // Copy relevant information from external header, if present.
1243 //
1244     Float_t uqJet[4];
1245     
1246     if (fRL) {
1247         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1248         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1249         {
1250             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1251             
1252             
1253             exHeader->TriggerJet(i, uqJet);
1254             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1255         }
1256     }
1257 //
1258 // Store quenching parameters
1259 //
1260     if (fQuench){
1261         Double_t z[4];
1262         Double_t xp, yp;
1263         if (fQuench == 1) {
1264             // Pythia::Quench()
1265             fPythia->GetQuenchingParameters(xp, yp, z);
1266         } else if (fQuench == 2){
1267             // Pyquen
1268             Double_t r1 = PARIMP.rb1;
1269             Double_t r2 = PARIMP.rb2;
1270             Double_t b  = PARIMP.b1;
1271             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1272             Double_t phi = PARIMP.psib1;
1273             xp = r * TMath::Cos(phi);
1274             yp = r * TMath::Sin(phi);
1275             
1276         } else if (fQuench == 4) {
1277             // QPythia
1278             Double_t xy[2];
1279             Double_t i0i1[2];
1280             AliFastGlauber::Instance()->GetSavedXY(xy);
1281             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1282             xp = xy[0];
1283             yp = xy[1];
1284             ((AliGenPythiaEventHeader*) fHeader)->SetInMediumLength(2. * i0i1[1] / i0i1[0]);
1285         }
1286         
1287             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1288             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1289     }
1290 //
1291 // Store pt^hard 
1292     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1293 //
1294 //  Pass header
1295 //
1296     AddHeader(fHeader);
1297     fHeader = 0x0;
1298 }
1299
1300 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1301 {
1302 // Check the kinematic trigger condition
1303 //
1304     Double_t eta[2];
1305     eta[0] = jet1->Eta();
1306     eta[1] = jet2->Eta();
1307     Double_t phi[2];
1308     phi[0] = jet1->Phi();
1309     phi[1] = jet2->Phi();
1310     Int_t    pdg[2]; 
1311     pdg[0] = jet1->GetPdgCode();
1312     pdg[1] = jet2->GetPdgCode();    
1313     Bool_t   triggered = kFALSE;
1314
1315     //
1316     //TO BE CHECKED: If we call this method for kPyBeautyJets, we need it here
1317     //
1318     if (fProcess == kPyJets || fProcess == kPyBeautyJets) {
1319         Int_t njets = 0;
1320         Int_t ntrig = 0;
1321         Float_t jets[4][10];
1322 //
1323 // Use Pythia clustering on parton level to determine jet axis
1324 //
1325         GetJets(njets, ntrig, jets);
1326         
1327         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1328 //
1329     } else {
1330         Int_t ij = 0;
1331         Int_t ig = 1;
1332         if (pdg[0] == kGamma) {
1333             ij = 1;
1334             ig = 0;
1335         }
1336         //Check eta range first...
1337         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1338             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1339         {
1340             //Eta is okay, now check phi range
1341             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1342                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1343             {
1344                 triggered = kTRUE;
1345             }
1346         }
1347     }
1348     return triggered;
1349 }
1350
1351
1352
1353 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1354 //
1355 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1356 //
1357     Bool_t checking = kFALSE;
1358     Int_t j, kcode, ks, km;
1359     Int_t nPartAcc = 0; //number of particles in the acceptance range
1360     Int_t numberOfAcceptedParticles = 1;
1361     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1362     Int_t npart = fParticles.GetEntriesFast();
1363     
1364     for (j = 0; j<npart; j++) {
1365         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1366         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1367         ks = jparticle->GetStatusCode();
1368         km = jparticle->GetFirstMother(); 
1369         
1370         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1371             nPartAcc++;
1372         }
1373         if( numberOfAcceptedParticles <= nPartAcc){
1374           checking = kTRUE;
1375           break;
1376         }
1377     }
1378
1379     return checking;
1380 }
1381
1382 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1383 {
1384   //
1385   // Load event into Pythia Common Block
1386   //
1387   
1388   Int_t npart = stack -> GetNprimary();
1389   Int_t n0 = 0;
1390   
1391   if (!flag) {
1392     (fPythia->GetPyjets())->N = npart;
1393   } else {
1394     n0 = (fPythia->GetPyjets())->N;
1395     (fPythia->GetPyjets())->N = n0 + npart;
1396   }
1397   
1398   
1399   for (Int_t part = 0; part < npart; part++) {
1400     TParticle *mPart = stack->Particle(part);
1401     
1402     Int_t kf     =  mPart->GetPdgCode();
1403     Int_t ks     =  mPart->GetStatusCode();
1404     Int_t idf    =  mPart->GetFirstDaughter();
1405     Int_t idl    =  mPart->GetLastDaughter();
1406     
1407     if (reHadr) {
1408             if (ks == 11 || ks == 12) {
1409         ks  -= 10;
1410         idf  = -1;
1411         idl  = -1;
1412             }
1413     }
1414     
1415     Float_t px = mPart->Px();
1416     Float_t py = mPart->Py();
1417     Float_t pz = mPart->Pz();
1418     Float_t e  = mPart->Energy();
1419     Float_t m  = mPart->GetCalcMass();
1420     
1421     
1422     (fPythia->GetPyjets())->P[0][part+n0] = px;
1423     (fPythia->GetPyjets())->P[1][part+n0] = py;
1424     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1425     (fPythia->GetPyjets())->P[3][part+n0] = e;
1426     (fPythia->GetPyjets())->P[4][part+n0] = m;
1427     
1428     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1429     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1430     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1431     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1432     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1433   }
1434 }
1435
1436 void  AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
1437 {
1438   //
1439   // Load event into Pythia Common Block
1440   //
1441   
1442   Int_t npart = stack -> GetEntries();
1443   Int_t n0 = 0;
1444   
1445   if (!flag) {
1446     (fPythia->GetPyjets())->N = npart;
1447   } else {
1448     n0 = (fPythia->GetPyjets())->N;
1449     (fPythia->GetPyjets())->N = n0 + npart;
1450   }
1451   
1452   
1453   for (Int_t part = 0; part < npart; part++) {
1454     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
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::IsInEMCAL(Float_t phi, Float_t eta)
1589 {
1590   // Is particle in EMCAL acceptance? 
1591   // phi in degrees, etamin=-etamax
1592   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi && 
1593      eta < fEMCALEta  ) 
1594     return kTRUE;
1595   else 
1596     return kFALSE;
1597 }
1598
1599 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1600 {
1601   // Is particle in PHOS acceptance? 
1602   // Acceptance slightly larger considered.
1603   // phi in degrees, etamin=-etamax
1604   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi && 
1605      eta < fPHOSEta  ) 
1606     return kTRUE;
1607   else 
1608     return kFALSE;
1609 }
1610
1611 void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1612 {
1613   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi 
1614   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1615   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1616   Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1617   
1618   //calculate deltaphi
1619   TParticle* ph = (TParticle *) fParticles.At(iphcand);
1620   Double_t phphi = ph->Phi();
1621   Double_t deltaphi = phiPHOS - phphi;
1622
1623   
1624   
1625   //loop for all particles and produce the phi rotation
1626   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1627   Double_t oldphi, newphi;
1628   Double_t newVx, newVy, R, Vz, time; 
1629   Double_t newPx, newPy, pt, Pz, e;
1630   for(Int_t i=0; i< np; i++) {
1631       TParticle* iparticle = (TParticle *) fParticles.At(i);
1632       oldphi = iparticle->Phi();
1633       newphi = oldphi + deltaphi;
1634       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle 
1635       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1636       
1637       R = iparticle->R();
1638       newVx = R*TMath::Cos(newphi);
1639       newVy = R*TMath::Sin(newphi);
1640       Vz = iparticle->Vz(); // don't transform
1641       time = iparticle->T(); // don't transform
1642       
1643       pt = iparticle->Pt();
1644       newPx = pt*TMath::Cos(newphi);
1645       newPy = pt*TMath::Sin(newphi);
1646       Pz = iparticle->Pz(); // don't transform
1647       e = iparticle->Energy(); // don't transform
1648       
1649       // apply rotation 
1650       iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1651       iparticle->SetMomentum(newPx, newPy, Pz, e);
1652       
1653   } //end particle loop 
1654   
1655    // now let's check that we put correctly the candidate photon in PHOS
1656    Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1657    Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax 
1658    if(IsInPHOS(phi,eta)) 
1659       okdd = kTRUE;
1660 }
1661
1662
1663 #ifdef never
1664 void AliGenPythia::Streamer(TBuffer &R__b)
1665 {
1666    // Stream an object of class AliGenPythia.
1667
1668    if (R__b.IsReading()) {
1669       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1670       AliGenerator::Streamer(R__b);
1671       R__b >> (Int_t&)fProcess;
1672       R__b >> (Int_t&)fStrucFunc;
1673       R__b >> (Int_t&)fForceDecay;
1674       R__b >> fEnergyCMS;
1675       R__b >> fKineBias;
1676       R__b >> fTrials;
1677       fParentSelect.Streamer(R__b);
1678       fChildSelect.Streamer(R__b);
1679       R__b >> fXsection;
1680 //      (AliPythia::Instance())->Streamer(R__b);
1681       R__b >> fPtHardMin;
1682       R__b >> fPtHardMax;
1683 //      if (fDecayer) fDecayer->Streamer(R__b);
1684    } else {
1685       R__b.WriteVersion(AliGenPythia::IsA());
1686       AliGenerator::Streamer(R__b);
1687       R__b << (Int_t)fProcess;
1688       R__b << (Int_t)fStrucFunc;
1689       R__b << (Int_t)fForceDecay;
1690       R__b << fEnergyCMS;
1691       R__b << fKineBias;
1692       R__b << fTrials;
1693       fParentSelect.Streamer(R__b);
1694       fChildSelect.Streamer(R__b);
1695       R__b << fXsection;
1696 //      R__b << fPythia;
1697       R__b << fPtHardMin;
1698       R__b << fPtHardMax;
1699       //     fDecayer->Streamer(R__b);
1700    }
1701 }
1702 #endif
1703
1704
1705