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