Warnings corrected.
[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  = 0.;
944       Float_t eta = 0.;
945       Float_t phi = 0.;
946       for (i=0; i< np; i++) {
947         TParticle* iparticle = (TParticle *) fParticles.At(i);
948         if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg && 
949            iparticle->Pt() > fElectronMinPt){
950           pt = iparticle->Pt();
951           phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
952           eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax    
953           if(IsInEMCAL(phi,eta))
954             ok =kTRUE;
955         }
956       }
957       if(!ok)
958         return 0;
959       AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
960     }
961     
962     // Check for minimum multiplicity
963     if (fTriggerMultiplicity > 0) {
964       Int_t multiplicity = 0;
965       for (i = 0; i < np; i++) {
966         TParticle *  iparticle = (TParticle *) fParticles.At(i);
967         
968         Int_t statusCode = iparticle->GetStatusCode();
969         
970         // Initial state particle
971         if (statusCode > 20)
972           continue;
973         
974         // skip quarks and gluons
975         Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
976         if (pdgCode <= 10 || pdgCode == 21)
977           continue;
978         
979         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
980           continue;
981         
982         TParticlePDG* pdgPart = iparticle->GetPDG();
983         if (pdgPart && pdgPart->Charge() == 0)
984           continue;
985         
986         ++multiplicity;
987       }
988
989       if (multiplicity < fTriggerMultiplicity) {
990         delete [] pParent;
991         return 0;
992       }
993       
994       Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
995     }    
996     
997      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
998     if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
999
1000       Bool_t okd = kFALSE;
1001
1002       Int_t pdg  = 22; 
1003       Int_t iphcand = -1;
1004       for (i=0; i< np; i++) {
1005          TParticle* iparticle = (TParticle *) fParticles.At(i);
1006          Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1007          Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax 
1008          
1009          if(iparticle->GetStatusCode() == 1 
1010             && iparticle->GetPdgCode() == pdg   
1011             && iparticle->Pt() > fPhotonMinPt    
1012             && eta < fPHOSEta){                 
1013             
1014             // first check if the photon is in PHOS phi
1015             if(IsInPHOS(phi,eta)){ 
1016                 okd = kTRUE;
1017                 break;
1018             } 
1019             if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1020              
1021          }
1022       }
1023       
1024       if(!okd && iphcand != -1) // execute rotation in phi 
1025           RotatePhi(iphcand,okd);
1026       
1027       if(!okd)
1028         return 0;
1029     }
1030     
1031     if (fTriggerParticle) {
1032         Bool_t triggered = kFALSE;
1033         for (i = 0; i < np; i++) {
1034             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1035             kf = CheckPDGCode(iparticle->GetPdgCode());
1036             if (kf != fTriggerParticle) continue;
1037             if (iparticle->Pt() == 0.) continue;
1038             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1039             triggered = kTRUE;
1040             break;
1041         }
1042         if (!triggered) {
1043           delete [] pParent;
1044           return 0;
1045         }
1046     }
1047         
1048
1049     // Check if there is a ccbar or bbbar pair with at least one of the two
1050     // in fYMin < y < fYMax
1051     //
1052     // TO BE CHECKED:  Should we require this for beauty jets?
1053     //
1054     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1055       TParticle *partCheck;
1056       TParticle *mother;
1057       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1058       Bool_t  theChild=kFALSE;
1059       Float_t y;  
1060       Int_t   pdg,mpdg,mpdgUpperFamily;
1061       for(i=0; i<np; i++) {
1062         partCheck = (TParticle*)fParticles.At(i);
1063         pdg = partCheck->GetPdgCode();  
1064         if(TMath::Abs(pdg) == fFlavorSelect) { // quark  
1065           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1066           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1067                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
1068           if(y>fYMin && y<fYMax) inYcut=kTRUE;
1069         }
1070         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1071           Int_t mi = partCheck->GetFirstMother() - 1;
1072           if(mi<0) continue;
1073           mother = (TParticle*)fParticles.At(mi);
1074           mpdg=TMath::Abs(mother->GetPdgCode());
1075           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1076           if ( ParentSelected(mpdg) || 
1077               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1078             if (KinematicSelection(partCheck,1)) {
1079               theChild=kTRUE;
1080             }
1081           }
1082         }
1083       }
1084       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1085         delete[] pParent;
1086         return 0;
1087       }
1088       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1089         delete[] pParent;
1090         return 0;       
1091       }
1092
1093     }
1094
1095     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1096     if ( (fProcess == kPyW ||
1097           fProcess == kPyZ ||
1098           fProcess == kPyMbDefault ||
1099           fProcess == kPyMb ||
1100           fProcess == kPyMbWithDirectPhoton ||
1101           fProcess == kPyMbNonDiffr)  
1102          && (fCutOnChild == 1) ) {
1103       if ( !CheckKinematicsOnChild() ) {
1104         delete[] pParent;
1105         return 0;
1106       }
1107     }
1108   
1109
1110     for (i = 0; i < np; i++) {
1111         Int_t trackIt = 0;
1112         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1113         kf = CheckPDGCode(iparticle->GetPdgCode());
1114         Int_t ks = iparticle->GetStatusCode();
1115         Int_t km = iparticle->GetFirstMother();
1116         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
1117             (ks != 1) ||
1118             //
1119             //TO BE CHECKED:  Should we require this for beauty jets?
1120             //
1121             ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
1122             nc++;
1123             if (ks == 1) trackIt = 1;
1124             Int_t ipa = iparticle->GetFirstMother()-1;
1125             
1126             iparent = (ipa > -1) ? pParent[ipa] : -1;
1127             
1128 //
1129 // store track information
1130             p[0] = iparticle->Px();
1131             p[1] = iparticle->Py();
1132             p[2] = iparticle->Pz();
1133             p[3] = iparticle->Energy();
1134
1135             
1136             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1137             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1138             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1139             
1140             Float_t tof = fEventTime + kconv * iparticle->T();
1141
1142             PushTrack(fTrackIt*trackIt, iparent, kf, 
1143                       p[0], p[1], p[2], p[3], 
1144                       origin[0], origin[1], origin[2], tof, 
1145                       polar[0], polar[1], polar[2],
1146                       kPPrimary, nt, 1., ks);
1147             fNprimaries++;
1148             KeepTrack(nt);
1149             pParent[i] = nt;
1150             SetHighWaterMark(nt);
1151             
1152         } // select particle
1153     } // particle loop 
1154
1155     delete[] pParent;
1156     
1157     return 1;
1158 }
1159
1160
1161 void AliGenPythia::FinishRun()
1162 {
1163 // Print x-section summary
1164     fPythia->Pystat(1);
1165
1166     if (fNev > 0.) {
1167         fQ  /= fNev;
1168         fX1 /= fNev;
1169         fX2 /= fNev;    
1170     }
1171     
1172     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1173     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1174 }
1175
1176 void AliGenPythia::AdjustWeights() const
1177 {
1178 // Adjust the weights after generation of all events
1179 //
1180     if (gAlice) {
1181         TParticle *part;
1182         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1183         for (Int_t i=0; i<ntrack; i++) {
1184             part= gAlice->GetMCApp()->Particle(i);
1185             part->SetWeight(part->GetWeight()*fKineBias);
1186         }
1187     }
1188 }
1189     
1190 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1191 {
1192 // Treat protons as inside nuclei with mass numbers a1 and a2  
1193
1194     fAProjectile = a1;
1195     fATarget     = a2;
1196     fNucPdf      = pdfset;  // 0 EKS98 1 EPS08
1197     fSetNuclei   = kTRUE;
1198 }
1199
1200
1201 void AliGenPythia::MakeHeader()
1202 {
1203 //
1204 // Make header for the simulated event
1205 // 
1206   if (gAlice) {
1207     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1208         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1209   }
1210
1211 // Builds the event header, to be called after each event
1212     if (fHeader) delete fHeader;
1213     fHeader = new AliGenPythiaEventHeader("Pythia");
1214 //
1215 // Event type  
1216     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1217 //
1218 // Number of trials
1219     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1220 //
1221 // Event Vertex 
1222     fHeader->SetPrimaryVertex(fVertex);
1223     
1224 //
1225 // Number of primaries
1226     fHeader->SetNProduced(fNprimaries);
1227 //
1228 // Jets that have triggered
1229
1230     //Need to store jets for b-jet studies too!
1231     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets)
1232     {
1233         Int_t ntrig, njet;
1234         Float_t jets[4][10];
1235         GetJets(njet, ntrig, jets);
1236
1237         
1238         for (Int_t i = 0; i < ntrig; i++) {
1239             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
1240                                                         jets[3][i]);
1241         }
1242     }
1243 //
1244 // Copy relevant information from external header, if present.
1245 //
1246     Float_t uqJet[4];
1247     
1248     if (fRL) {
1249         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1250         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1251         {
1252             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1253             
1254             
1255             exHeader->TriggerJet(i, uqJet);
1256             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1257         }
1258     }
1259 //
1260 // Store quenching parameters
1261 //
1262     if (fQuench){
1263         Double_t z[4];
1264         Double_t xp, yp;
1265         if (fQuench == 1) {
1266             // Pythia::Quench()
1267             fPythia->GetQuenchingParameters(xp, yp, z);
1268         } else if (fQuench == 2){
1269             // Pyquen
1270             Double_t r1 = PARIMP.rb1;
1271             Double_t r2 = PARIMP.rb2;
1272             Double_t b  = PARIMP.b1;
1273             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1274             Double_t phi = PARIMP.psib1;
1275             xp = r * TMath::Cos(phi);
1276             yp = r * TMath::Sin(phi);
1277             
1278         } else if (fQuench == 4) {
1279             // QPythia
1280             Double_t xy[2];
1281             Double_t i0i1[2];
1282             AliFastGlauber::Instance()->GetSavedXY(xy);
1283             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1284             xp = xy[0];
1285             yp = xy[1];
1286             ((AliGenPythiaEventHeader*) fHeader)->SetInMediumLength(2. * i0i1[1] / i0i1[0]);
1287         }
1288         
1289             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1290             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1291     }
1292 //
1293 // Store pt^hard 
1294     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1295 //
1296 //  Pass header
1297 //
1298     AddHeader(fHeader);
1299     fHeader = 0x0;
1300 }
1301
1302 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1303 {
1304 // Check the kinematic trigger condition
1305 //
1306     Double_t eta[2];
1307     eta[0] = jet1->Eta();
1308     eta[1] = jet2->Eta();
1309     Double_t phi[2];
1310     phi[0] = jet1->Phi();
1311     phi[1] = jet2->Phi();
1312     Int_t    pdg[2]; 
1313     pdg[0] = jet1->GetPdgCode();
1314     pdg[1] = jet2->GetPdgCode();    
1315     Bool_t   triggered = kFALSE;
1316
1317     //
1318     //TO BE CHECKED: If we call this method for kPyBeautyJets, we need it here
1319     //
1320     if (fProcess == kPyJets || fProcess == kPyBeautyJets) {
1321         Int_t njets = 0;
1322         Int_t ntrig = 0;
1323         Float_t jets[4][10];
1324 //
1325 // Use Pythia clustering on parton level to determine jet axis
1326 //
1327         GetJets(njets, ntrig, jets);
1328         
1329         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1330 //
1331     } else {
1332         Int_t ij = 0;
1333         Int_t ig = 1;
1334         if (pdg[0] == kGamma) {
1335             ij = 1;
1336             ig = 0;
1337         }
1338         //Check eta range first...
1339         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1340             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1341         {
1342             //Eta is okay, now check phi range
1343             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1344                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1345             {
1346                 triggered = kTRUE;
1347             }
1348         }
1349     }
1350     return triggered;
1351 }
1352
1353
1354
1355 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1356 //
1357 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1358 //
1359     Bool_t checking = kFALSE;
1360     Int_t j, kcode, ks, km;
1361     Int_t nPartAcc = 0; //number of particles in the acceptance range
1362     Int_t numberOfAcceptedParticles = 1;
1363     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1364     Int_t npart = fParticles.GetEntriesFast();
1365     
1366     for (j = 0; j<npart; j++) {
1367         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1368         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1369         ks = jparticle->GetStatusCode();
1370         km = jparticle->GetFirstMother(); 
1371         
1372         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1373             nPartAcc++;
1374         }
1375         if( numberOfAcceptedParticles <= nPartAcc){
1376           checking = kTRUE;
1377           break;
1378         }
1379     }
1380
1381     return checking;
1382 }
1383
1384 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1385 {
1386   //
1387   // Load event into Pythia Common Block
1388   //
1389   
1390   Int_t npart = stack -> GetNprimary();
1391   Int_t n0 = 0;
1392   
1393   if (!flag) {
1394     (fPythia->GetPyjets())->N = npart;
1395   } else {
1396     n0 = (fPythia->GetPyjets())->N;
1397     (fPythia->GetPyjets())->N = n0 + npart;
1398   }
1399   
1400   
1401   for (Int_t part = 0; part < npart; part++) {
1402     TParticle *mPart = stack->Particle(part);
1403     
1404     Int_t kf     =  mPart->GetPdgCode();
1405     Int_t ks     =  mPart->GetStatusCode();
1406     Int_t idf    =  mPart->GetFirstDaughter();
1407     Int_t idl    =  mPart->GetLastDaughter();
1408     
1409     if (reHadr) {
1410             if (ks == 11 || ks == 12) {
1411         ks  -= 10;
1412         idf  = -1;
1413         idl  = -1;
1414             }
1415     }
1416     
1417     Float_t px = mPart->Px();
1418     Float_t py = mPart->Py();
1419     Float_t pz = mPart->Pz();
1420     Float_t e  = mPart->Energy();
1421     Float_t m  = mPart->GetCalcMass();
1422     
1423     
1424     (fPythia->GetPyjets())->P[0][part+n0] = px;
1425     (fPythia->GetPyjets())->P[1][part+n0] = py;
1426     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1427     (fPythia->GetPyjets())->P[3][part+n0] = e;
1428     (fPythia->GetPyjets())->P[4][part+n0] = m;
1429     
1430     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1431     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1432     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1433     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1434     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1435   }
1436 }
1437
1438 void  AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
1439 {
1440   //
1441   // Load event into Pythia Common Block
1442   //
1443   
1444   Int_t npart = stack -> GetEntries();
1445   Int_t n0 = 0;
1446   
1447   if (!flag) {
1448     (fPythia->GetPyjets())->N = npart;
1449   } else {
1450     n0 = (fPythia->GetPyjets())->N;
1451     (fPythia->GetPyjets())->N = n0 + npart;
1452   }
1453   
1454   
1455   for (Int_t part = 0; part < npart; part++) {
1456     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1457     Int_t kf     =  mPart->GetPdgCode();
1458     Int_t ks     =  mPart->GetStatusCode();
1459     Int_t idf    =  mPart->GetFirstDaughter();
1460     Int_t idl    =  mPart->GetLastDaughter();
1461     
1462     if (reHadr) {
1463             if (ks == 11 || ks == 12) {
1464         ks  -= 10;
1465         idf  = -1;
1466         idl  = -1;
1467             }
1468     }
1469     
1470     Float_t px = mPart->Px();
1471     Float_t py = mPart->Py();
1472     Float_t pz = mPart->Pz();
1473     Float_t e  = mPart->Energy();
1474     Float_t m  = mPart->GetCalcMass();
1475     
1476     
1477     (fPythia->GetPyjets())->P[0][part+n0] = px;
1478     (fPythia->GetPyjets())->P[1][part+n0] = py;
1479     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1480     (fPythia->GetPyjets())->P[3][part+n0] = e;
1481     (fPythia->GetPyjets())->P[4][part+n0] = m;
1482     
1483     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1484     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1485     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1486     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1487     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1488   }
1489 }
1490
1491
1492 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1493 {
1494 //
1495 //  Calls the Pythia jet finding algorithm to find jets in the current event
1496 //
1497 //
1498 //
1499 //  Save jets
1500     Int_t n     = fPythia->GetN();
1501
1502 //
1503 //  Run Jet Finder
1504     fPythia->Pycell(njets);
1505     Int_t i;
1506     for (i = 0; i < njets; i++) {
1507         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1508         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1509         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1510         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1511
1512         jets[0][i] = px;
1513         jets[1][i] = py;
1514         jets[2][i] = pz;
1515         jets[3][i] = e;
1516     }
1517 }
1518
1519
1520
1521 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1522 {
1523 //
1524 //  Calls the Pythia clustering algorithm to find jets in the current event
1525 //
1526     Int_t n     = fPythia->GetN();
1527     nJets       = 0;
1528     nJetsTrig   = 0;
1529     if (fJetReconstruction == kCluster) {
1530 //
1531 //  Configure cluster algorithm
1532 //    
1533         fPythia->SetPARU(43, 2.);
1534         fPythia->SetMSTU(41, 1);
1535 //
1536 //  Call cluster algorithm
1537 //    
1538         fPythia->Pyclus(nJets);
1539 //
1540 //  Loading jets from common block
1541 //
1542     } else {
1543
1544 //
1545 //  Run Jet Finder
1546         fPythia->Pycell(nJets);
1547     }
1548
1549     Int_t i;
1550     for (i = 0; i < nJets; i++) {
1551         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1552         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1553         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1554         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1555         Float_t pt    = TMath::Sqrt(px * px + py * py);
1556         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);  
1557         Float_t theta = TMath::ATan2(pt,pz);
1558         Float_t et    = e * TMath::Sin(theta);
1559         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1560         if (
1561             eta > fEtaMinJet && eta < fEtaMaxJet && 
1562             phi > fPhiMinJet && phi < fPhiMaxJet &&
1563             et  > fEtMinJet  && et  < fEtMaxJet     
1564             ) 
1565         {
1566             jets[0][nJetsTrig] = px;
1567             jets[1][nJetsTrig] = py;
1568             jets[2][nJetsTrig] = pz;
1569             jets[3][nJetsTrig] = e;
1570             nJetsTrig++;
1571 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1572         } else {
1573 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1574         }
1575     }
1576 }
1577
1578 void AliGenPythia::GetSubEventTime()
1579 {
1580   // Calculates time of the next subevent
1581   fEventTime = 0.;
1582   if (fEventsTime) {
1583     TArrayF &array = *fEventsTime;
1584     fEventTime = array[fCurSubEvent++];
1585   }
1586   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1587   return;
1588 }
1589
1590 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1591 {
1592   // Is particle in EMCAL acceptance? 
1593   // phi in degrees, etamin=-etamax
1594   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi && 
1595      eta < fEMCALEta  ) 
1596     return kTRUE;
1597   else 
1598     return kFALSE;
1599 }
1600
1601 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1602 {
1603   // Is particle in PHOS acceptance? 
1604   // Acceptance slightly larger considered.
1605   // phi in degrees, etamin=-etamax
1606   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi && 
1607      eta < fPHOSEta  ) 
1608     return kTRUE;
1609   else 
1610     return kFALSE;
1611 }
1612
1613 void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1614 {
1615   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi 
1616   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1617   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1618   Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1619   
1620   //calculate deltaphi
1621   TParticle* ph = (TParticle *) fParticles.At(iphcand);
1622   Double_t phphi = ph->Phi();
1623   Double_t deltaphi = phiPHOS - phphi;
1624
1625   
1626   
1627   //loop for all particles and produce the phi rotation
1628   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1629   Double_t oldphi, newphi;
1630   Double_t newVx, newVy, R, Vz, time; 
1631   Double_t newPx, newPy, pt, Pz, e;
1632   for(Int_t i=0; i< np; i++) {
1633       TParticle* iparticle = (TParticle *) fParticles.At(i);
1634       oldphi = iparticle->Phi();
1635       newphi = oldphi + deltaphi;
1636       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle 
1637       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1638       
1639       R = iparticle->R();
1640       newVx = R*TMath::Cos(newphi);
1641       newVy = R*TMath::Sin(newphi);
1642       Vz = iparticle->Vz(); // don't transform
1643       time = iparticle->T(); // don't transform
1644       
1645       pt = iparticle->Pt();
1646       newPx = pt*TMath::Cos(newphi);
1647       newPy = pt*TMath::Sin(newphi);
1648       Pz = iparticle->Pz(); // don't transform
1649       e = iparticle->Energy(); // don't transform
1650       
1651       // apply rotation 
1652       iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1653       iparticle->SetMomentum(newPx, newPy, Pz, e);
1654       
1655   } //end particle loop 
1656   
1657    // now let's check that we put correctly the candidate photon in PHOS
1658    Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1659    Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax 
1660    if(IsInPHOS(phi,eta)) 
1661       okdd = kTRUE;
1662 }
1663
1664
1665 #ifdef never
1666 void AliGenPythia::Streamer(TBuffer &R__b)
1667 {
1668    // Stream an object of class AliGenPythia.
1669
1670    if (R__b.IsReading()) {
1671       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1672       AliGenerator::Streamer(R__b);
1673       R__b >> (Int_t&)fProcess;
1674       R__b >> (Int_t&)fStrucFunc;
1675       R__b >> (Int_t&)fForceDecay;
1676       R__b >> fEnergyCMS;
1677       R__b >> fKineBias;
1678       R__b >> fTrials;
1679       fParentSelect.Streamer(R__b);
1680       fChildSelect.Streamer(R__b);
1681       R__b >> fXsection;
1682 //      (AliPythia::Instance())->Streamer(R__b);
1683       R__b >> fPtHardMin;
1684       R__b >> fPtHardMax;
1685 //      if (fDecayer) fDecayer->Streamer(R__b);
1686    } else {
1687       R__b.WriteVersion(AliGenPythia::IsA());
1688       AliGenerator::Streamer(R__b);
1689       R__b << (Int_t)fProcess;
1690       R__b << (Int_t)fStrucFunc;
1691       R__b << (Int_t)fForceDecay;
1692       R__b << fEnergyCMS;
1693       R__b << fKineBias;
1694       R__b << fTrials;
1695       fParentSelect.Streamer(R__b);
1696       fChildSelect.Streamer(R__b);
1697       R__b << fXsection;
1698 //      R__b << fPythia;
1699       R__b << fPtHardMin;
1700       R__b << fPtHardMax;
1701       //     fDecayer->Streamer(R__b);
1702    }
1703 }
1704 #endif
1705
1706
1707