1844a3d3a8ac8e7dd003af5112eec5898bbe4456
[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           if(pdg == 111 || 
844              (pdg == 22 && iparticle->GetFirstMother() != 11))//No photon from hadron decay
845             {
846               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
847               Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax         
848               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
849                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
850                 ok =kTRUE;
851             }
852         }
853       }
854       if(!ok)
855         return 0;
856     }
857     
858     if (fTriggerParticle) {
859         Bool_t triggered = kFALSE;
860         for (i = 0; i < np; i++) {
861             TParticle *  iparticle = (TParticle *) fParticles->At(i);
862             kf = CheckPDGCode(iparticle->GetPdgCode());
863             if (kf != fTriggerParticle) continue;
864             if (iparticle->Pt() == 0.) continue;
865             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
866             triggered = kTRUE;
867             break;
868         }
869         if (!triggered) {
870           delete [] pParent;
871           return 0;
872         }
873     }
874         
875
876     // Check if there is a ccbar or bbbar pair with at least one of the two
877     // in fYMin < y < fYMax
878     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
879       TParticle *hvq;
880       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
881       Float_t yQ;  
882       Int_t   pdgQ;
883       for(i=0; i<np; i++) {
884         hvq = (TParticle*)fParticles->At(i);
885         pdgQ = hvq->GetPdgCode();  
886         if(TMath::Abs(pdgQ) != fFlavorSelect) continue; 
887         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
888         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
889                             (hvq->Energy()-hvq->Pz()+1.e-13));
890         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
891       }
892       if (!theQ || !theQbar || !inYcut) {
893         delete[] pParent;
894         return 0;
895       }
896     }
897
898     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
899     if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)  
900          && (fCutOnChild == 1) ) {
901       if ( !CheckKinematicsOnChild() ) {
902         delete[] pParent;
903         return 0;
904       }
905     }
906   
907
908     for (i = 0; i < np; i++) {
909         Int_t trackIt = 0;
910         TParticle *  iparticle = (TParticle *) fParticles->At(i);
911         kf = CheckPDGCode(iparticle->GetPdgCode());
912         Int_t ks = iparticle->GetStatusCode();
913         Int_t km = iparticle->GetFirstMother();
914         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
915             (ks != 1) ||
916             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
917             nc++;
918             if (ks == 1) trackIt = 1;
919             Int_t ipa = iparticle->GetFirstMother()-1;
920             
921             iparent = (ipa > -1) ? pParent[ipa] : -1;
922             
923 //
924 // store track information
925             p[0] = iparticle->Px();
926             p[1] = iparticle->Py();
927             p[2] = iparticle->Pz();
928             p[3] = iparticle->Energy();
929
930             
931             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
932             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
933             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
934             
935             Float_t tof = fEventTime + kconv * iparticle->T();
936
937             PushTrack(fTrackIt*trackIt, iparent, kf, 
938                       p[0], p[1], p[2], p[3], 
939                       origin[0], origin[1], origin[2], tof, 
940                       polar[0], polar[1], polar[2],
941                       kPPrimary, nt, 1., ks);
942             //
943             // Special Treatment to store color-flow
944             //
945             if (ks == 3 || ks == 13 || ks == 14) {
946                 TParticle* particle = 0;
947                 if (fStack) {
948                     particle = fStack->Particle(nt);
949                 } else {
950                     particle = gAlice->Stack()->Particle(nt);
951                 }
952                 particle->SetFirstDaughter(fPythia->GetK(2, i));
953                 particle->SetLastDaughter(fPythia->GetK(3, i));         
954             }
955             
956             KeepTrack(nt);
957             pParent[i] = nt;
958             SetHighWaterMark(nt);
959             
960         } // select particle
961     } // particle loop 
962
963     delete[] pParent;
964     
965     return 1;
966 }
967
968
969 void AliGenPythia::FinishRun()
970 {
971 // Print x-section summary
972     fPythia->Pystat(1);
973
974     if (fNev > 0.) {
975         fQ  /= fNev;
976         fX1 /= fNev;
977         fX2 /= fNev;    
978     }
979     
980     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
981     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
982 }
983
984 void AliGenPythia::AdjustWeights() const
985 {
986 // Adjust the weights after generation of all events
987 //
988     if (gAlice) {
989         TParticle *part;
990         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
991         for (Int_t i=0; i<ntrack; i++) {
992             part= gAlice->GetMCApp()->Particle(i);
993             part->SetWeight(part->GetWeight()*fKineBias);
994         }
995     }
996 }
997     
998 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
999 {
1000 // Treat protons as inside nuclei with mass numbers a1 and a2  
1001
1002     fAProjectile = a1;
1003     fATarget     = a2;
1004     fSetNuclei   = kTRUE;
1005 }
1006
1007
1008 void AliGenPythia::MakeHeader()
1009 {
1010 //
1011 // Make header for the simulated event
1012 // 
1013   if (gAlice) {
1014     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1015         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1016   }
1017
1018 // Builds the event header, to be called after each event
1019     if (fHeader) delete fHeader;
1020     fHeader = new AliGenPythiaEventHeader("Pythia");
1021 //
1022 // Event type  
1023     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1024 //
1025 // Number of trials
1026     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1027 //
1028 // Event Vertex 
1029     fHeader->SetPrimaryVertex(fVertex);
1030 //
1031 // Jets that have triggered
1032
1033     if (fProcess == kPyJets)
1034     {
1035         Int_t ntrig, njet;
1036         Float_t jets[4][10];
1037         GetJets(njet, ntrig, jets);
1038
1039         
1040         for (Int_t i = 0; i < ntrig; i++) {
1041             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
1042                                                         jets[3][i]);
1043         }
1044     }
1045 //
1046 // Copy relevant information from external header, if present.
1047 //
1048     Float_t uqJet[4];
1049     
1050     if (fRL) {
1051         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1052         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1053         {
1054             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1055             
1056             
1057             exHeader->TriggerJet(i, uqJet);
1058             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1059         }
1060     }
1061 //
1062 // Store quenching parameters
1063 //
1064     if (fQuench){
1065         Double_t z[4];
1066         Double_t xp, yp;
1067         if (fQuench == 1) {
1068             // Pythia::Quench()
1069             fPythia->GetQuenchingParameters(xp, yp, z);
1070         } else {
1071             // Pyquen
1072             Double_t r1 = PARIMP.rb1;
1073             Double_t r2 = PARIMP.rb2;
1074             Double_t b  = PARIMP.b1;
1075             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1076             Double_t phi = PARIMP.psib1;
1077             xp = r * TMath::Cos(phi);
1078             yp = r * TMath::Sin(phi);
1079             
1080         }
1081             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1082             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1083         }
1084 //
1085 // Store pt^hard 
1086     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1087 //
1088 //  Pass header
1089 //
1090     AddHeader(fHeader);
1091 }
1092
1093 void AliGenPythia::AddHeader(AliGenEventHeader* header)
1094 {
1095     // Add header to container or runloader
1096     if (fContainer) {
1097         fContainer->AddHeader(header);
1098     } else {
1099         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);   
1100     }
1101 }
1102
1103
1104 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1105 {
1106 // Check the kinematic trigger condition
1107 //
1108     Double_t eta[2];
1109     eta[0] = jet1->Eta();
1110     eta[1] = jet2->Eta();
1111     Double_t phi[2];
1112     phi[0] = jet1->Phi();
1113     phi[1] = jet2->Phi();
1114     Int_t    pdg[2]; 
1115     pdg[0] = jet1->GetPdgCode();
1116     pdg[1] = jet2->GetPdgCode();    
1117     Bool_t   triggered = kFALSE;
1118
1119     if (fProcess == kPyJets) {
1120         Int_t njets = 0;
1121         Int_t ntrig = 0;
1122         Float_t jets[4][10];
1123 //
1124 // Use Pythia clustering on parton level to determine jet axis
1125 //
1126         GetJets(njets, ntrig, jets);
1127         
1128         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1129 //
1130     } else {
1131         Int_t ij = 0;
1132         Int_t ig = 1;
1133         if (pdg[0] == kGamma) {
1134             ij = 1;
1135             ig = 0;
1136         }
1137         //Check eta range first...
1138         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1139             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1140         {
1141             //Eta is okay, now check phi range
1142             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1143                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1144             {
1145                 triggered = kTRUE;
1146             }
1147         }
1148     }
1149     return triggered;
1150 }
1151
1152
1153
1154 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1155 //
1156 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1157 //
1158     Bool_t checking = kFALSE;
1159     Int_t j, kcode, ks, km;
1160     Int_t nPartAcc = 0; //number of particles in the acceptance range
1161     Int_t numberOfAcceptedParticles = 1;
1162     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1163     Int_t npart = fParticles->GetEntriesFast();
1164     
1165     for (j = 0; j<npart; j++) {
1166         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1167         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1168         ks = jparticle->GetStatusCode();
1169         km = jparticle->GetFirstMother(); 
1170         
1171         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1172             nPartAcc++;
1173         }
1174         if( numberOfAcceptedParticles <= nPartAcc){
1175           checking = kTRUE;
1176           break;
1177         }
1178     }
1179
1180     return checking;
1181 }
1182
1183 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1184 {
1185 //
1186 // Load event into Pythia Common Block
1187 //
1188
1189     Int_t npart = stack -> GetNprimary();
1190     Int_t n0 = 0;
1191     
1192     if (!flag) {
1193         (fPythia->GetPyjets())->N = npart;
1194     } else {
1195         n0 = (fPythia->GetPyjets())->N;
1196         (fPythia->GetPyjets())->N = n0 + npart;
1197     }
1198     
1199     
1200     for (Int_t part = 0; part < npart; part++) {
1201         TParticle *mPart = stack->Particle(part);
1202         
1203         Int_t kf     =  mPart->GetPdgCode();
1204         Int_t ks     =  mPart->GetStatusCode();
1205         Int_t idf    =  mPart->GetFirstDaughter();
1206         Int_t idl    =  mPart->GetLastDaughter();
1207         
1208         if (reHadr) {
1209             if (ks == 11 || ks == 12) {
1210                 ks  -= 10;
1211                 idf  = -1;
1212                 idl  = -1;
1213             }
1214         }
1215         
1216         Float_t px = mPart->Px();
1217         Float_t py = mPart->Py();
1218         Float_t pz = mPart->Pz();
1219         Float_t e  = mPart->Energy();
1220         Float_t m  = mPart->GetCalcMass();
1221         
1222         
1223         (fPythia->GetPyjets())->P[0][part+n0] = px;
1224         (fPythia->GetPyjets())->P[1][part+n0] = py;
1225         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1226         (fPythia->GetPyjets())->P[3][part+n0] = e;
1227         (fPythia->GetPyjets())->P[4][part+n0] = m;
1228         
1229         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1230         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1231         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1232         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1233         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1234     }
1235 }
1236
1237
1238 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1239 {
1240 //
1241 //  Calls the Pythia jet finding algorithm to find jets in the current event
1242 //
1243 //
1244 //
1245 //  Save jets
1246     Int_t n     = fPythia->GetN();
1247
1248 //
1249 //  Run Jet Finder
1250     fPythia->Pycell(njets);
1251     Int_t i;
1252     for (i = 0; i < njets; i++) {
1253         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1254         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1255         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1256         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1257
1258         jets[0][i] = px;
1259         jets[1][i] = py;
1260         jets[2][i] = pz;
1261         jets[3][i] = e;
1262     }
1263 }
1264
1265
1266
1267 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1268 {
1269 //
1270 //  Calls the Pythia clustering algorithm to find jets in the current event
1271 //
1272     Int_t n     = fPythia->GetN();
1273     nJets       = 0;
1274     nJetsTrig   = 0;
1275     if (fJetReconstruction == kCluster) {
1276 //
1277 //  Configure cluster algorithm
1278 //    
1279         fPythia->SetPARU(43, 2.);
1280         fPythia->SetMSTU(41, 1);
1281 //
1282 //  Call cluster algorithm
1283 //    
1284         fPythia->Pyclus(nJets);
1285 //
1286 //  Loading jets from common block
1287 //
1288     } else {
1289
1290 //
1291 //  Run Jet Finder
1292         fPythia->Pycell(nJets);
1293     }
1294
1295     Int_t i;
1296     for (i = 0; i < nJets; i++) {
1297         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1298         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1299         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1300         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1301         Float_t pt    = TMath::Sqrt(px * px + py * py);
1302         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);  
1303         Float_t theta = TMath::ATan2(pt,pz);
1304         Float_t et    = e * TMath::Sin(theta);
1305         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1306         if (
1307             eta > fEtaMinJet && eta < fEtaMaxJet && 
1308             phi > fPhiMinJet && phi < fPhiMaxJet &&
1309             et  > fEtMinJet  && et  < fEtMaxJet     
1310             ) 
1311         {
1312             jets[0][nJetsTrig] = px;
1313             jets[1][nJetsTrig] = py;
1314             jets[2][nJetsTrig] = pz;
1315             jets[3][nJetsTrig] = e;
1316             nJetsTrig++;
1317 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1318         } else {
1319 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1320         }
1321     }
1322 }
1323
1324 void AliGenPythia::GetSubEventTime()
1325 {
1326   // Calculates time of the next subevent
1327   fEventTime = 0.;
1328   if (fEventsTime) {
1329     TArrayF &array = *fEventsTime;
1330     fEventTime = array[fCurSubEvent++];
1331   }
1332   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1333   return;
1334 }
1335
1336 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1337 {
1338   // Is particle in EMCAL acceptance? 
1339   // phi in degrees, etamin=-etamax
1340   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi && 
1341      eta < fEMCALEta  ) 
1342     return kTRUE;
1343   else 
1344     return kFALSE;
1345 }
1346
1347 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1348 {
1349   // Is particle in PHOS acceptance? 
1350   // Acceptance slightly larger considered.
1351   // phi in degrees, etamin=-etamax
1352   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi && 
1353      eta < fPHOSEta  ) 
1354     return kTRUE;
1355   else 
1356     return kFALSE;
1357 }
1358
1359
1360
1361 #ifdef never
1362 void AliGenPythia::Streamer(TBuffer &R__b)
1363 {
1364    // Stream an object of class AliGenPythia.
1365
1366    if (R__b.IsReading()) {
1367       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1368       AliGenerator::Streamer(R__b);
1369       R__b >> (Int_t&)fProcess;
1370       R__b >> (Int_t&)fStrucFunc;
1371       R__b >> (Int_t&)fForceDecay;
1372       R__b >> fEnergyCMS;
1373       R__b >> fKineBias;
1374       R__b >> fTrials;
1375       fParentSelect.Streamer(R__b);
1376       fChildSelect.Streamer(R__b);
1377       R__b >> fXsection;
1378 //      (AliPythia::Instance())->Streamer(R__b);
1379       R__b >> fPtHardMin;
1380       R__b >> fPtHardMax;
1381 //      if (fDecayer) fDecayer->Streamer(R__b);
1382    } else {
1383       R__b.WriteVersion(AliGenPythia::IsA());
1384       AliGenerator::Streamer(R__b);
1385       R__b << (Int_t)fProcess;
1386       R__b << (Int_t)fStrucFunc;
1387       R__b << (Int_t)fForceDecay;
1388       R__b << fEnergyCMS;
1389       R__b << fKineBias;
1390       R__b << fTrials;
1391       fParentSelect.Streamer(R__b);
1392       fChildSelect.Streamer(R__b);
1393       R__b << fXsection;
1394 //      R__b << fPythia;
1395       R__b << fPtHardMin;
1396       R__b << fPtHardMax;
1397       //     fDecayer->Streamer(R__b);
1398    }
1399 }
1400 #endif
1401
1402
1403