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