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