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