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