Possibility to simulate pile-up event with correct time distribution.
[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 <TDatabasePDG.h>
29 #include <TParticle.h>
30 #include <TPDGCode.h>
31 #include <TSystem.h>
32 #include <TTree.h>
33 #include "AliConst.h"
34 #include "AliDecayerPythia.h"
35 #include "AliGenPythia.h"
36 #include "AliHeader.h"
37 #include "AliGenPythiaEventHeader.h"
38 #include "AliPythia.h"
39 #include "AliPythiaRndm.h"
40 #include "AliRun.h"
41 #include "AliStack.h"
42 #include "AliRunLoader.h"
43 #include "AliMC.h"
44 #include "pyquenCommon.h"
45
46 ClassImp(AliGenPythia)
47
48 AliGenPythia::AliGenPythia()
49                  :AliGenMC()
50 {
51 // Default Constructor
52   fParticles = 0;
53   fPythia    = 0;
54   fHeader = 0;
55   fReadFromFile = 0;
56   fEventTime = 0.;
57   fDecayer   = new AliDecayerPythia();
58   SetEventListRange();
59   SetJetPhiRange();
60   SetJetEtaRange();
61   SetJetEtRange();
62   SetGammaPhiRange();
63   SetGammaEtaRange();
64   SetPtKick();
65   SetQuench();
66   SetHadronisation();  
67   fSetNuclei = kFALSE;
68   if (!AliPythiaRndm::GetPythiaRandom()) 
69     AliPythiaRndm::SetPythiaRandom(GetRandom());
70 }
71
72 AliGenPythia::AliGenPythia(Int_t npart)
73                  :AliGenMC(npart)
74 {
75 // default charm production at 5. 5 TeV
76 // semimuonic decay
77 // structure function GRVHO
78 //
79     fName = "Pythia";
80     fTitle= "Particle Generator using PYTHIA";
81     fXsection  = 0.;
82     fReadFromFile = 0;
83     fEventTime = 0.;
84     SetProcess();
85     SetStrucFunc();
86     SetForceDecay();
87     SetPtHard();
88     SetYHard();
89     SetEnergyCMS();
90     fDecayer = new AliDecayerPythia();
91     // Set random number generator 
92     if (!AliPythiaRndm::GetPythiaRandom()) 
93       AliPythiaRndm::SetPythiaRandom(GetRandom());
94     fFlavorSelect   = 0;
95     // Produced particles  
96     fParticles = new TClonesArray("TParticle",1000);
97     fHeader = 0;
98     SetEventListRange();
99     SetJetPhiRange();
100     SetJetEtaRange();
101     SetJetEtRange();
102     SetGammaPhiRange();
103     SetGammaEtaRange();
104     SetJetReconstructionMode();
105     SetQuench();
106     SetHadronisation();
107     SetPtKick();
108     // Options determining what to keep in the stack (Heavy flavour generation)
109     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
110     fFeedDownOpt = kTRUE;             // allow feed down from higher family
111     // Fragmentation on/off
112     fFragmentation = kTRUE;
113     // Default counting mode
114     fCountMode = kCountAll;
115     // Pycel
116     SetPycellParameters();
117     fSetNuclei = kFALSE;
118 }
119
120 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
121     :AliGenMC(Pythia)
122 {
123 // copy constructor
124     Pythia.Copy(*this);
125 }
126
127 AliGenPythia::~AliGenPythia()
128 {
129 // Destructor
130 }
131
132 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
133                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
134 {
135 // Set pycell parameters
136     fPycellEtaMax    =  etamax;
137     fPycellNEta      =  neta;
138     fPycellNPhi      =  nphi;
139     fPycellThreshold =  thresh;
140     fPycellEtSeed    =  etseed;
141     fPycellMinEtJet  =  minet;
142     fPycellMaxRadius =  r;
143 }
144
145
146
147 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
148 {
149   // Set a range of event numbers, for which a table
150   // of generated particle will be printed
151   fDebugEventFirst = eventFirst;
152   fDebugEventLast  = eventLast;
153   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
154 }
155
156 void AliGenPythia::Init()
157 {
158 // Initialisation
159     
160     SetMC(AliPythia::Instance());
161     fPythia=(AliPythia*) fMCEvGen;
162     
163 //
164     fParentWeight=1./Float_t(fNpart);
165 //
166 //  Forward Paramters to the AliPythia object
167     fDecayer->SetForceDecay(fForceDecay);    
168     fDecayer->Init();
169
170
171     fPythia->SetCKIN(3,fPtHardMin);
172     fPythia->SetCKIN(4,fPtHardMax);
173     fPythia->SetCKIN(7,fYHardMin);
174     fPythia->SetCKIN(8,fYHardMax);
175     
176     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);  
177     // Fragmentation?
178     if (fFragmentation) {
179       fPythia->SetMSTP(111,1);
180     } else {
181       fPythia->SetMSTP(111,0);
182     }
183
184
185 //  initial state radiation   
186     fPythia->SetMSTP(61,fGinit);
187 //  final state radiation
188     fPythia->SetMSTP(71,fGfinal);
189 //  pt - kick
190     if (fPtKick > 0.) {
191         fPythia->SetMSTP(91,1);
192         fPythia->SetPARP(91,fPtKick);
193     } else {
194         fPythia->SetMSTP(91,0);
195     }
196
197
198     if (fReadFromFile) {
199         fRL  =  AliRunLoader::Open(fFileName, "Partons");
200         fRL->LoadKinematics();
201         fRL->LoadHeader();
202     } else {
203         fRL = 0x0;
204     }
205     
206     
207  //
208     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
209
210 //  Parent and Children Selection
211     switch (fProcess) 
212     {
213     case kPyCharm:
214     case kPyCharmUnforced:
215     case kPyCharmPbPbMNR:
216     case kPyCharmppMNR:
217     case kPyCharmpPbMNR:
218         fParentSelect[0] =   411;
219         fParentSelect[1] =   421;
220         fParentSelect[2] =   431;
221         fParentSelect[3] =  4122;
222         fFlavorSelect    =  4;  
223         break;
224     case kPyD0PbPbMNR:
225     case kPyD0pPbMNR:
226     case kPyD0ppMNR:
227         fParentSelect[0] =   421;
228         fFlavorSelect    =   4; 
229         break;
230     case kPyDPlusPbPbMNR:
231     case kPyDPluspPbMNR:
232     case kPyDPlusppMNR:
233         fParentSelect[0] =   411;
234         fFlavorSelect    =   4; 
235         break;
236     case kPyBeauty:
237     case kPyBeautyPbPbMNR:
238     case kPyBeautypPbMNR:
239     case kPyBeautyppMNR:
240         fParentSelect[0]=  511;
241         fParentSelect[1]=  521;
242         fParentSelect[2]=  531;
243         fParentSelect[3]= 5122;
244         fParentSelect[4]= 5132;
245         fParentSelect[5]= 5232;
246         fParentSelect[6]= 5332;
247         fFlavorSelect   = 5;    
248         break;
249     case kPyBeautyUnforced:
250         fParentSelect[0] =  511;
251         fParentSelect[1] =  521;
252         fParentSelect[2] =  531;
253         fParentSelect[3] = 5122;
254         fParentSelect[4] = 5132;
255         fParentSelect[5] = 5232;
256         fParentSelect[6] = 5332;
257         fFlavorSelect    = 5;   
258         break;
259     case kPyJpsiChi:
260     case kPyJpsi:
261         fParentSelect[0] = 443;
262         break;
263     case kPyMb:
264     case kPyMbNonDiffr:
265     case kPyJets:
266     case kPyDirectGamma:
267         break;
268     case kPyW:
269         break;
270     }
271 //
272 //
273 //  JetFinder for Trigger
274 //
275 //  Configure detector (EMCAL like)
276 //
277         fPythia->SetPARU(51, fPycellEtaMax);
278         fPythia->SetMSTU(51, fPycellNEta);
279         fPythia->SetMSTU(52, fPycellNPhi);
280 //
281 //  Configure Jet Finder
282 //  
283         fPythia->SetPARU(58,  fPycellThreshold);
284         fPythia->SetPARU(52,  fPycellEtSeed);
285         fPythia->SetPARU(53,  fPycellMinEtJet);
286         fPythia->SetPARU(54,  fPycellMaxRadius);
287         fPythia->SetMSTU(54,  2);
288 //
289 //  This counts the total number of calls to Pyevnt() per run.
290     fTrialsRun = 0;
291     fQ         = 0.;
292     fX1        = 0.;
293     fX2        = 0.;    
294     fNev       = 0 ;
295 //    
296 //
297 //
298     AliGenMC::Init();
299 //
300 //
301 //  
302     if (fSetNuclei) {
303         fDyBoost = 0;
304         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
305     }
306
307     if (fQuench) {
308         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
309     }
310     
311 }
312
313 void AliGenPythia::Generate()
314 {
315 // Generate one event
316     
317     fDecayer->ForceDecay();
318
319     Float_t polar[3]   =   {0,0,0};
320     Float_t origin[3]  =   {0,0,0};
321     Float_t p[4];
322 //  converts from mm/c to s
323     const Float_t kconv=0.001/2.999792458e8;
324 //
325     Int_t nt=0;
326     Int_t jev=0;
327     Int_t j, kf;
328     fTrials=0;
329     fEventTime = 0.;
330     
331     
332
333     //  Set collision vertex position 
334     if (fVertexSmear == kPerEvent) Vertex();
335     
336 //  event loop    
337     while(1)
338     {
339 //
340 // Produce event
341 //
342 //
343 // Switch hadronisation off
344 //
345         fPythia->SetMSTJ(1, 0);
346 //
347 // Either produce new event or read partons from file
348 //      
349         if (!fReadFromFile) {
350             fPythia->Pyevnt();
351             fNpartons = fPythia->GetN();
352         } else {
353             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
354             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
355             fPythia->SetN(0);
356             LoadEvent(fRL->Stack(), 0 , 1);
357             fPythia->Pyedit(21);
358         }
359         
360 //
361 //  Run quenching routine 
362 //
363         if (fQuench == 1) {
364             fPythia->Quench();
365         } else if (fQuench == 2){
366             fPythia->Pyquen(208., 0, 0.);
367         }
368 //
369 // Switch hadronisation on
370 //
371         fPythia->SetMSTJ(1, 1);
372 //
373 // .. and perform hadronisation
374 //      printf("Calling hadronisation %d\n", fPythia->GetN());
375         fPythia->Pyexec();      
376         fTrials++;
377         fPythia->ImportParticles(fParticles,"All");
378         Boost();
379 //
380 //
381 //
382         Int_t i;
383         
384
385         Int_t np = fParticles->GetEntriesFast();
386         
387         if (np == 0) continue;
388 //
389         
390 //
391         Int_t* pParent   = new Int_t[np];
392         Int_t* pSelected = new Int_t[np];
393         Int_t* trackIt   = new Int_t[np];
394         for (i = 0; i < np; i++) {
395             pParent[i]   = -1;
396             pSelected[i] =  0;
397             trackIt[i]   =  0;
398         }
399
400         Int_t nc = 0;        // Total n. of selected particles
401         Int_t nParents = 0;  // Selected parents
402         Int_t nTkbles = 0;   // Trackable particles
403         if (fProcess != kPyMb && fProcess != kPyJets && 
404             fProcess != kPyDirectGamma &&
405             fProcess != kPyMbNonDiffr  &&
406             fProcess != kPyW) {
407             
408             for (i = 0; i < np; i++) {
409                 TParticle* iparticle = (TParticle *) fParticles->At(i);
410                 Int_t ks = iparticle->GetStatusCode();
411                 kf = CheckPDGCode(iparticle->GetPdgCode());
412 // No initial state partons
413                 if (ks==21) continue;
414 //
415 // Heavy Flavor Selection
416 //
417                 // quark ?
418                 kf = TMath::Abs(kf);
419                 Int_t kfl = kf;
420                 // Resonance
421
422                 if (kfl > 100000) kfl %= 100000;
423                 if (kfl > 10000)  kfl %= 10000;
424                 // meson ?
425                 if  (kfl > 10) kfl/=100;
426                 // baryon
427                 if (kfl > 10) kfl/=10;
428                 Int_t ipa = iparticle->GetFirstMother()-1;
429                 Int_t kfMo = 0;
430 //
431 // Establish mother daughter relation between heavy quarks and mesons
432 //
433                 if (kf >= fFlavorSelect && kf <= 6) {
434                     Int_t idau = iparticle->GetFirstDaughter() - 1;
435                     if (idau > -1) {
436                         TParticle* daughter = (TParticle *) fParticles->At(idau);
437                         Int_t pdgD = daughter->GetPdgCode();
438                         if (pdgD == 91 || pdgD == 92) {
439                             Int_t jmin = daughter->GetFirstDaughter() - 1;
440                             Int_t jmax = daughter->GetLastDaughter()  - 1;                          
441                             for (Int_t j = jmin; j <= jmax; j++)
442                                 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
443                         } // is string or cluster
444                     } // has daughter
445                 } // heavy quark
446                 
447
448                 if (ipa > -1) {
449                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
450                     kfMo = TMath::Abs(mother->GetPdgCode());
451                 }
452                 
453                 // What to keep in Stack?
454                 Bool_t flavorOK = kFALSE;
455                 Bool_t selectOK = kFALSE;
456                 if (fFeedDownOpt) {
457                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
458                 } else {
459                     if (kfl > fFlavorSelect) {
460                         nc = -1;
461                         break;
462                     }
463                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
464                 }
465                 switch (fStackFillOpt) {
466                 case kFlavorSelection:
467                     selectOK = kTRUE;
468                     break;
469                 case kParentSelection:
470                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
471                     break;
472                 }
473                 if (flavorOK && selectOK) { 
474 //
475 // Heavy flavor hadron or quark
476 //
477 // Kinematic seletion on final state heavy flavor mesons
478                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
479                     {
480                         continue;
481                     }
482                     pSelected[i] = 1;
483                     if (ParentSelected(kf)) ++nParents; // Update parent count
484 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
485                 } else {
486 // Kinematic seletion on decay products
487                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
488                         && !KinematicSelection(iparticle, 1)) 
489                     {
490                         continue;
491                     }
492 //
493 // Decay products 
494 // Select if mother was selected and is not tracked
495
496                     if (pSelected[ipa] && 
497                         !trackIt[ipa]  &&     // mother will be  tracked ?
498                         kfMo !=  5 &&         // mother is b-quark, don't store fragments          
499                         kfMo !=  4 &&         // mother is c-quark, don't store fragments 
500                         kf   != 92)           // don't store string
501                     {
502 //
503 // Semi-stable or de-selected: diselect decay products:
504 // 
505 //
506                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
507                         {
508                             Int_t ipF = iparticle->GetFirstDaughter();
509                             Int_t ipL = iparticle->GetLastDaughter();   
510                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
511                         }
512 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
513                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
514                     }
515                 }
516                 if (pSelected[i] == -1) pSelected[i] = 0;
517                 if (!pSelected[i]) continue;
518                 // Count quarks only if you did not include fragmentation
519                 if (fFragmentation && kf <= 10) continue;
520
521                 nc++;
522 // Decision on tracking
523                 trackIt[i] = 0;
524 //
525 // Track final state particle
526                 if (ks == 1) trackIt[i] = 1;
527 // Track semi-stable particles
528                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
529 // Track particles selected by process if undecayed. 
530                 if (fForceDecay == kNoDecay) {
531                     if (ParentSelected(kf)) trackIt[i] = 1;
532                 } else {
533                     if (ParentSelected(kf)) trackIt[i] = 0;
534                 }
535                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
536 //
537 //
538
539             } // particle selection loop
540             if (nc > 0) {
541                 for (i = 0; i<np; i++) {
542                     if (!pSelected[i]) continue;
543                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
544                     kf = CheckPDGCode(iparticle->GetPdgCode());
545                     Int_t ks = iparticle->GetStatusCode();  
546                     p[0] = iparticle->Px();
547                     p[1] = iparticle->Py();
548                     p[2] = iparticle->Pz();
549                     p[3] = iparticle->Energy();
550                     
551                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
552                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
553                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
554                     
555                     Float_t tof   = kconv*iparticle->T();
556                     Int_t ipa     = iparticle->GetFirstMother()-1;
557                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
558  
559                     PushTrack(fTrackIt*trackIt[i], iparent, kf, 
560                               p[0], p[1], p[2], p[3], 
561                               origin[0], origin[1], origin[2], tof, 
562                               polar[0], polar[1], polar[2],
563                               kPPrimary, nt, 1., ks);
564                     pParent[i] = nt;
565                     KeepTrack(nt); 
566                 } //  PushTrack loop
567             }
568         } else {
569             nc = GenerateMB();
570         } // mb ?
571         
572         GetSubEventTime();
573
574         if (pParent)   delete[] pParent;
575         if (pSelected) delete[] pSelected;
576         if (trackIt)   delete[] trackIt;
577
578         if (nc > 0) {
579           switch (fCountMode) {
580           case kCountAll:
581             // printf(" Count all \n");
582             jev += nc;
583             break;
584           case kCountParents:
585             // printf(" Count parents \n");
586             jev += nParents;
587             break;
588           case kCountTrackables:
589             // printf(" Count trackable \n");
590             jev += nTkbles;
591             break;
592           }
593             if (jev >= fNpart || fNpart == -1) {
594                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
595                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
596
597                 fQ  += fPythia->GetVINT(51);
598                 fX1 += fPythia->GetVINT(41);
599                 fX2 += fPythia->GetVINT(42);
600                 fTrialsRun += fTrials;
601                 fNev++;
602                 MakeHeader();
603                 break;
604             }
605         }
606     } // event loop
607     SetHighWaterMark(nt);
608 //  adjust weight due to kinematic selection
609 //    AdjustWeights();
610 //  get cross-section
611     fXsection=fPythia->GetPARI(1);
612 }
613
614 Int_t  AliGenPythia::GenerateMB()
615 {
616 //
617 // Min Bias selection and other global selections
618 //
619     Int_t i, kf, nt, iparent;
620     Int_t nc = 0;
621     Float_t p[4];
622     Float_t polar[3]   =   {0,0,0};
623     Float_t origin[3]  =   {0,0,0};
624 //  converts from mm/c to s
625     const Float_t kconv=0.001/2.999792458e8;
626     
627
628     
629     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
630
631
632     Int_t* pParent = new Int_t[np];
633     for (i=0; i< np; i++) pParent[i] = -1;
634     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
635         TParticle* jet1 = (TParticle *) fParticles->At(6);
636         TParticle* jet2 = (TParticle *) fParticles->At(7);
637         if (!CheckTrigger(jet1, jet2)) return 0;
638     }
639     
640     for (i = 0; i < np; i++) {
641         Int_t trackIt = 0;
642         TParticle *  iparticle = (TParticle *) fParticles->At(i);
643         kf = CheckPDGCode(iparticle->GetPdgCode());
644         Int_t ks = iparticle->GetStatusCode();
645         Int_t km = iparticle->GetFirstMother();
646         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
647             (ks != 1) ||
648             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
649             nc++;
650             if (ks == 1) trackIt = 1;
651             Int_t ipa = iparticle->GetFirstMother()-1;
652             
653             iparent = (ipa > -1) ? pParent[ipa] : -1;
654             
655 //
656 // store track information
657             p[0] = iparticle->Px();
658             p[1] = iparticle->Py();
659             p[2] = iparticle->Pz();
660             p[3] = iparticle->Energy();
661
662             
663             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
664             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
665             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
666             
667             Float_t tof = fEventTime + kconv * iparticle->T();
668
669             PushTrack(fTrackIt*trackIt, iparent, kf, 
670                       p[0], p[1], p[2], p[3], 
671                       origin[0], origin[1], origin[2], tof, 
672                       polar[0], polar[1], polar[2],
673                       kPPrimary, nt, 1., ks);
674             //
675             // Special Treatment to store color-flow
676             //
677             if (ks == 3 || ks == 13 || ks == 14) {
678                 TParticle* particle = 0;
679                 if (fStack) {
680                     particle = fStack->Particle(nt);
681                 } else {
682                     particle = gAlice->Stack()->Particle(nt);
683                 }
684                 particle->SetFirstDaughter(fPythia->GetK(2, i));
685                 particle->SetLastDaughter(fPythia->GetK(3, i));         
686             }
687             
688             KeepTrack(nt);
689             pParent[i] = nt;
690             SetHighWaterMark(nt);
691             
692         } // select particle
693     } // particle loop 
694
695     if (pParent) delete[] pParent;
696     
697     printf("\n I've put %i particles on the stack \n",nc);
698     return 1;
699 }
700
701
702 void AliGenPythia::FinishRun()
703 {
704 // Print x-section summary
705     fPythia->Pystat(1);
706
707     if (fNev > 0.) {
708         fQ  /= fNev;
709         fX1 /= fNev;
710         fX2 /= fNev;    
711     }
712     
713     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
714     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
715 }
716
717 void AliGenPythia::AdjustWeights()
718 {
719 // Adjust the weights after generation of all events
720 //
721     if (gAlice) {
722         TParticle *part;
723         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
724         for (Int_t i=0; i<ntrack; i++) {
725             part= gAlice->GetMCApp()->Particle(i);
726             part->SetWeight(part->GetWeight()*fKineBias);
727         }
728     }
729 }
730     
731 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
732 {
733 // Treat protons as inside nuclei with mass numbers a1 and a2  
734
735     fAProjectile = a1;
736     fATarget     = a2;
737     fSetNuclei   = kTRUE;
738 }
739
740
741 void AliGenPythia::MakeHeader()
742 {
743   if (gAlice) {
744     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
745         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
746   }
747
748 // Builds the event header, to be called after each event
749     if (fHeader) delete fHeader;
750     fHeader = new AliGenPythiaEventHeader("Pythia");
751 //
752 // Event type  
753     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
754 //
755 // Number of trials
756     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
757 //
758 // Event Vertex 
759     fHeader->SetPrimaryVertex(fVertex);
760 //
761 // Jets that have triggered
762
763     if (fProcess == kPyJets)
764     {
765         Int_t ntrig, njet;
766         Float_t jets[4][10];
767         GetJets(njet, ntrig, jets);
768
769         
770         for (Int_t i = 0; i < ntrig; i++) {
771             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i], 
772                                                         jets[3][i]);
773         }
774     }
775 //
776 // Copy relevant information from external header, if present.
777 //
778     Float_t uqJet[4];
779     
780     if (fRL) {
781         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
782         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
783         {
784             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
785             
786             
787             exHeader->TriggerJet(i, uqJet);
788             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
789         }
790     }
791 //
792 // Store quenching parameters
793 //
794     if (fQuench){
795         Double_t z[4];
796         Double_t xp, yp;
797         if (fQuench == 1) {
798             // Pythia::Quench()
799             fPythia->GetQuenchingParameters(xp, yp, z);
800         } else {
801             // Pyquen
802             Double_t r1 = PARIMP.rb1;
803             Double_t r2 = PARIMP.rb2;
804             Double_t b  = PARIMP.b1;
805             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
806             Double_t phi = PARIMP.psib1;
807             xp = r * TMath::Cos(phi);
808             yp = r * TMath::Sin(phi);
809             
810         }
811             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
812             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
813         }
814     
815 //
816 //  Pass header to RunLoader
817 //
818     AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(fHeader);   
819 }
820         
821
822 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
823 {
824 // Check the kinematic trigger condition
825 //
826     Double_t eta[2];
827     eta[0] = jet1->Eta();
828     eta[1] = jet2->Eta();
829     Double_t phi[2];
830     phi[0] = jet1->Phi();
831     phi[1] = jet2->Phi();
832     Int_t    pdg[2]; 
833     pdg[0] = jet1->GetPdgCode();
834     pdg[1] = jet2->GetPdgCode();    
835     Bool_t   triggered = kFALSE;
836
837     if (fProcess == kPyJets) {
838         Int_t njets = 0;
839         Int_t ntrig = 0;
840         Float_t jets[4][10];
841 //
842 // Use Pythia clustering on parton level to determine jet axis
843 //
844         GetJets(njets, ntrig, jets);
845         
846         if (ntrig) triggered = kTRUE;
847 //
848     } else {
849         Int_t ij = 0;
850         Int_t ig = 1;
851         if (pdg[0] == kGamma) {
852             ij = 1;
853             ig = 0;
854         }
855         //Check eta range first...
856         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
857             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
858         {
859             //Eta is okay, now check phi range
860             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
861                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
862             {
863                 triggered = kTRUE;
864             }
865         }
866     }
867     return triggered;
868 }
869           
870 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
871 {
872 // Assignment operator
873     rhs.Copy(*this);
874     return *this;
875 }
876
877 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
878 {
879 //
880 // Load event into Pythia Common Block
881 //
882
883     Int_t npart = stack -> GetNprimary();
884     Int_t n0 = 0;
885     
886     if (!flag) {
887         (fPythia->GetPyjets())->N = npart;
888     } else {
889         n0 = (fPythia->GetPyjets())->N;
890         (fPythia->GetPyjets())->N = n0 + npart;
891     }
892     
893     
894     for (Int_t part = 0; part < npart; part++) {
895         TParticle *MPart = stack->Particle(part);
896         
897         Int_t kf     =  MPart->GetPdgCode();
898         Int_t ks     =  MPart->GetStatusCode();
899         Int_t idf    =  MPart->GetFirstDaughter();
900         Int_t idl    =  MPart->GetLastDaughter();
901         
902         if (reHadr) {
903             if (ks == 11 || ks == 12) {
904                 ks  -= 10;
905                 idf  = -1;
906                 idl  = -1;
907             }
908         }
909         
910         Float_t px = MPart->Px();
911         Float_t py = MPart->Py();
912         Float_t pz = MPart->Pz();
913         Float_t e  = MPart->Energy();
914         Float_t m  = MPart->GetCalcMass();
915         
916         
917         (fPythia->GetPyjets())->P[0][part+n0] = px;
918         (fPythia->GetPyjets())->P[1][part+n0] = py;
919         (fPythia->GetPyjets())->P[2][part+n0] = pz;
920         (fPythia->GetPyjets())->P[3][part+n0] = e;
921         (fPythia->GetPyjets())->P[4][part+n0] = m;
922         
923         (fPythia->GetPyjets())->K[1][part+n0] = kf;
924         (fPythia->GetPyjets())->K[0][part+n0] = ks;
925         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
926         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
927         (fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
928     }
929 }
930
931
932 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
933 {
934 //
935 //  Calls the Pythia jet finding algorithm to find jets in the current event
936 //
937 //
938 //
939 //  Save jets
940     Int_t n     = fPythia->GetN();
941
942 //
943 //  Run Jet Finder
944     fPythia->Pycell(njets);
945     Int_t i;
946     for (i = 0; i < njets; i++) {
947         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
948         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
949         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
950         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
951
952         jets[0][i] = px;
953         jets[1][i] = py;
954         jets[2][i] = pz;
955         jets[3][i] = e;
956     }
957 }
958
959
960
961 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
962 {
963 //
964 //  Calls the Pythia clustering algorithm to find jets in the current event
965 //
966     Int_t n     = fPythia->GetN();
967     nJets       = 0;
968     nJetsTrig   = 0;
969     if (fJetReconstruction == kCluster) {
970 //
971 //  Configure cluster algorithm
972 //    
973         fPythia->SetPARU(43, 2.);
974         fPythia->SetMSTU(41, 1);
975 //
976 //  Call cluster algorithm
977 //    
978         fPythia->Pyclus(nJets);
979 //
980 //  Loading jets from common block
981 //
982     } else {
983
984 //
985 //  Run Jet Finder
986         fPythia->Pycell(nJets);
987     }
988
989     Int_t i;
990     for (i = 0; i < nJets; i++) {
991         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
992         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
993         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
994         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
995         Float_t pt    = TMath::Sqrt(px * px + py * py);
996         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);  
997         Float_t theta = TMath::ATan2(pt,pz);
998         Float_t et    = e * TMath::Sin(theta);
999         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1000         if (
1001             eta > fEtaMinJet && eta < fEtaMaxJet && 
1002             phi > fPhiMinJet && phi < fPhiMaxJet &&
1003             et  > fEtMinJet  && et  < fEtMaxJet     
1004             ) 
1005         {
1006             jets[0][nJetsTrig] = px;
1007             jets[1][nJetsTrig] = py;
1008             jets[2][nJetsTrig] = pz;
1009             jets[3][nJetsTrig] = e;
1010             nJetsTrig++;
1011 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1012         } else {
1013 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1014         }
1015     }
1016 }
1017
1018 void AliGenPythia::GetSubEventTime()
1019 {
1020   // Calculates time of the next subevent
1021   fEventTime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
1022 }
1023
1024 #ifdef never
1025 void AliGenPythia::Streamer(TBuffer &R__b)
1026 {
1027    // Stream an object of class AliGenPythia.
1028
1029    if (R__b.IsReading()) {
1030       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1031       AliGenerator::Streamer(R__b);
1032       R__b >> (Int_t&)fProcess;
1033       R__b >> (Int_t&)fStrucFunc;
1034       R__b >> (Int_t&)fForceDecay;
1035       R__b >> fEnergyCMS;
1036       R__b >> fKineBias;
1037       R__b >> fTrials;
1038       fParentSelect.Streamer(R__b);
1039       fChildSelect.Streamer(R__b);
1040       R__b >> fXsection;
1041 //      (AliPythia::Instance())->Streamer(R__b);
1042       R__b >> fPtHardMin;
1043       R__b >> fPtHardMax;
1044 //      if (fDecayer) fDecayer->Streamer(R__b);
1045    } else {
1046       R__b.WriteVersion(AliGenPythia::IsA());
1047       AliGenerator::Streamer(R__b);
1048       R__b << (Int_t)fProcess;
1049       R__b << (Int_t)fStrucFunc;
1050       R__b << (Int_t)fForceDecay;
1051       R__b << fEnergyCMS;
1052       R__b << fKineBias;
1053       R__b << fTrials;
1054       fParentSelect.Streamer(R__b);
1055       fChildSelect.Streamer(R__b);
1056       R__b << fXsection;
1057 //      R__b << fPythia;
1058       R__b << fPtHardMin;
1059       R__b << fPtHardMax;
1060       //     fDecayer->Streamer(R__b);
1061    }
1062 }
1063 #endif
1064
1065
1066