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