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