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