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