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