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