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