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