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