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