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