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