Pass status code in SetTrack.
[u/mrichter/AliRoot.git] / EVGEN / 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 /*
17 $Log$
18 Revision 1.53  2002/03/25 14:51:13  morsch
19 New stack-fill and count options introduced (N. Carrer).
20
21 Revision 1.51  2002/03/06 08:46:57  morsch
22 - Loop until np-1
23 - delete dyn. alloc. arrays (N. Carrer)
24
25 Revision 1.50  2002/03/03 13:48:50  morsch
26 Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
27 NLO calculations (Nicola Carrer).
28
29 Revision 1.49  2002/02/08 16:50:50  morsch
30 Add name and title in constructor.
31
32 Revision 1.48  2001/12/20 11:44:28  morsch
33 Add kinematic bias for direct gamma production.
34
35 Revision 1.47  2001/12/19 14:45:00  morsch
36 Store number of trials in header.
37
38 Revision 1.46  2001/12/19 10:36:19  morsch
39 Add possibility if jet kinematic biasing.
40
41 Revision 1.45  2001/11/28 08:06:52  morsch
42 Use fMaxLifeTime parameter.
43
44 Revision 1.44  2001/11/27 13:13:07  morsch
45 Maximum lifetime for long-lived particles to be put on the stack is parameter.
46 It can be set via SetMaximumLifetime(..).
47
48 Revision 1.43  2001/10/21 18:35:56  hristov
49 Several pointers were set to zero in the default constructors to avoid memory management problems
50
51 Revision 1.42  2001/10/15 08:21:55  morsch
52 Vertex truncation settings moved to AliGenMC.
53
54 Revision 1.41  2001/10/08 08:45:42  morsch
55 Possibility of vertex cut added.
56
57 Revision 1.40  2001/09/25 11:30:23  morsch
58 Pass event vertex to header.
59
60 Revision 1.39  2001/07/27 17:09:36  morsch
61 Use local SetTrack, KeepTrack and SetHighWaterMark methods
62 to delegate either to local stack or to stack owned by AliRun.
63 (Piotr Skowronski, A.M.)
64
65 Revision 1.38  2001/07/13 10:58:54  morsch
66 - Some coded moved to AliGenMC
67 - Improved handling of secondary vertices.
68
69 Revision 1.37  2001/06/28 11:17:28  morsch
70 SetEventListRange setter added. Events in specified range are listed for
71 debugging. (Yuri Kharlov)
72
73 Revision 1.36  2001/03/30 07:05:49  morsch
74 Final print-out in finish run.
75 Write parton system for jet-production (preliminary solution).
76
77 Revision 1.35  2001/03/09 13:03:40  morsch
78 Process_t and Struc_Func_t moved to AliPythia.h
79
80 Revision 1.34  2001/02/14 15:50:40  hristov
81 The last particle in event marked using SetHighWaterMark
82
83 Revision 1.33  2001/01/30 09:23:12  hristov
84 Streamers removed (R.Brun)
85
86 Revision 1.32  2001/01/26 19:55:51  hristov
87 Major upgrade of AliRoot code
88
89 Revision 1.31  2001/01/17 10:54:31  hristov
90 Better protection against FPE
91
92 Revision 1.30  2000/12/18 08:55:35  morsch
93 Make AliPythia dependent generartors work with new scheme of random number generation
94
95 Revision 1.29  2000/12/04 11:22:03  morsch
96 Init of sRandom as in 1.15
97
98 Revision 1.28  2000/12/02 11:41:39  morsch
99 Use SetRandom() to initialize random number generator in constructor.
100
101 Revision 1.27  2000/11/30 20:29:02  morsch
102 Initialise static variable sRandom in constructor: sRandom = fRandom;
103
104 Revision 1.26  2000/11/30 07:12:50  alibrary
105 Introducing new Rndm and QA classes
106
107 Revision 1.25  2000/10/18 19:11:27  hristov
108 Division by zero fixed
109
110 Revision 1.24  2000/09/18 10:41:35  morsch
111 Add possibility to use nuclear structure functions from PDF library V8.
112
113 Revision 1.23  2000/09/14 14:05:40  morsch
114 dito
115
116 Revision 1.22  2000/09/14 14:02:22  morsch
117 - Correct conversion from mm to cm when passing particle vertex to MC.
118 - Correct handling of fForceDecay == all.
119
120 Revision 1.21  2000/09/12 14:14:55  morsch
121 Call fDecayer->ForceDecay() at the beginning of Generate().
122
123 Revision 1.20  2000/09/06 14:29:33  morsch
124 Use AliPythia for event generation an AliDecayPythia for decays.
125 Correct handling of "nodecay" option
126
127 Revision 1.19  2000/07/11 18:24:56  fca
128 Coding convention corrections + few minor bug fixes
129
130 Revision 1.18  2000/06/30 12:40:34  morsch
131 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
132 Geant units (cm).
133
134 Revision 1.17  2000/06/09 20:34:07  morsch
135 All coding rule violations except RS3 corrected
136
137 Revision 1.16  2000/05/15 15:04:20  morsch
138 The full event is written for fNtrack = -1
139 Coding rule violations corrected.
140
141 Revision 1.15  2000/04/26 10:14:24  morsch
142 Particles array has one entry more than pythia particle list. Upper bound of
143 particle loop changed to np-1 (R. Guernane, AM)
144
145 Revision 1.14  2000/04/05 08:36:13  morsch
146 Check status code of particles in Pythia event
147 to avoid double counting as partonic state and final state particle.
148
149 Revision 1.13  1999/11/09 07:38:48  fca
150 Changes for compatibility with version 2.23 of ROOT
151
152 Revision 1.12  1999/11/03 17:43:20  fca
153 New version from G.Martinez & A.Morsch
154
155 Revision 1.11  1999/09/29 09:24:14  fca
156 Introduction of the Copyright and cvs Log
157 */
158
159 #include "AliGenPythia.h"
160 #include "AliGenPythiaEventHeader.h"
161 #include "AliDecayerPythia.h"
162 #include "AliRun.h"
163 #include "AliPythia.h"
164 #include "AliPDG.h"
165 #include <TParticle.h>
166 #include <TSystem.h>
167
168  ClassImp(AliGenPythia)
169
170 AliGenPythia::AliGenPythia()
171                  :AliGenMC()
172 {
173 // Default Constructor
174   fParticles = 0;
175   fPythia    = 0;
176   fDecayer   = new AliDecayerPythia();
177   SetEventListRange();
178   SetJetPhiRange();
179   SetJetEtaRange();
180 }
181
182 AliGenPythia::AliGenPythia(Int_t npart)
183                  :AliGenMC(npart)
184 {
185 // default charm production at 5. 5 TeV
186 // semimuonic decay
187 // structure function GRVHO
188 //
189     fName = "Pythia";
190     fTitle= "Particle Generator using PYTHIA";
191     fXsection  = 0.;
192     fNucA1=0;
193     fNucA2=0;
194     SetProcess();
195     SetStrucFunc();
196     SetForceDecay();
197     SetPtHard();
198     SetEnergyCMS();
199     fDecayer = new AliDecayerPythia();
200     // Set random number generator 
201     sRandom=fRandom;
202     fFlavorSelect   = 0;
203     // Produced particles  
204     fParticles = new TClonesArray("TParticle",1000);
205     fEventVertex.Set(3);
206     SetEventListRange();
207     SetJetPhiRange();
208     SetJetEtaRange();
209     // Options determining what to keep in the stack (Heavy flavour generation)
210     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
211     fFeedDownOpt = kTRUE;             // allow feed down from higher family
212     // Fragmentation on/off
213     fFragmentation = kTRUE;
214     // Default counting mode
215     fCountMode = kCountAll;
216 }
217
218 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
219 {
220 // copy constructor
221 }
222
223 AliGenPythia::~AliGenPythia()
224 {
225 // Destructor
226 }
227
228 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
229 {
230   // Set a range of event numbers, for which a table
231   // of generated particle will be printed
232   fDebugEventFirst = eventFirst;
233   fDebugEventLast  = eventLast;
234   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
235 }
236
237 void AliGenPythia::Init()
238 {
239 // Initialisation
240   SetMC(AliPythia::Instance());
241     fPythia=(AliPythia*) fgMCEvGen;
242 //
243     fParentWeight=1./Float_t(fNpart);
244 //
245 //  Forward Paramters to the AliPythia object
246     fDecayer->SetForceDecay(fForceDecay);    
247     fDecayer->Init();
248
249
250     fPythia->SetCKIN(3,fPtHardMin);
251     fPythia->SetCKIN(4,fPtHardMax);    
252     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
253     // Fragmentation?
254     if (fFragmentation) {
255       fPythia->SetMSTP(111,1);
256     } else {
257       fPythia->SetMSTP(111,0);
258     }
259
260     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
261
262     //    fPythia->Pylist(0);
263     //    fPythia->Pystat(2);
264 //  Parent and Children Selection
265     switch (fProcess) 
266     {
267     case kPyCharm:
268     case kPyCharmUnforced:
269     case kPyCharmPbMNR:
270         fParentSelect[0] =   411;
271         fParentSelect[1] =   421;
272         fParentSelect[2] =   431;
273         fParentSelect[3] =  4122;
274         fFlavorSelect    =  4;  
275         break;
276     case kPyD0PbMNR:
277         fParentSelect[0] =   421;
278         fFlavorSelect    =  4;  
279         break;
280     case kPyBeauty:
281         fParentSelect[0]=  511;
282         fParentSelect[1]=  521;
283         fParentSelect[2]=  531;
284         fParentSelect[3]= 5122;
285         fParentSelect[4]= 5132;
286         fParentSelect[5]= 5232;
287         fParentSelect[6]= 5332;
288         fFlavorSelect   = 5;    
289         break;
290     case kPyBeautyUnforced:
291         fParentSelect[0] =  511;
292         fParentSelect[1] =  521;
293         fParentSelect[2] =  531;
294         fParentSelect[3] = 5122;
295         fParentSelect[4] = 5132;
296         fParentSelect[5] = 5232;
297         fParentSelect[6] = 5332;
298         fFlavorSelect    = 5;   
299         break;
300     case kPyJpsiChi:
301     case kPyJpsi:
302         fParentSelect[0] = 443;
303         break;
304     case kPyMb:
305     case kPyJets:
306     case kPyDirectGamma:
307         break;
308     }
309     AliGenMC::Init();
310 }
311
312 void AliGenPythia::Generate()
313 {
314 // Generate one event
315     fDecayer->ForceDecay();
316
317     Float_t polar[3]   =   {0,0,0};
318     Float_t origin[3]  =   {0,0,0};
319     Float_t p[3];
320 //  converts from mm/c to s
321     const Float_t kconv=0.001/2.999792458e8;
322 //
323     Int_t nt=0;
324     Int_t jev=0;
325     Int_t j, kf;
326     fTrials=0;
327
328     //  Set collision vertex position 
329     if(fVertexSmear==kPerEvent) {
330         fPythia->SetMSTP(151,1);
331         for (j=0;j<3;j++) {
332             fPythia->SetPARP(151+j, fOsigma[j]*10.);
333         }
334     } else if (fVertexSmear==kPerTrack) {
335         fPythia->SetMSTP(151,0);
336     }
337 //  event loop    
338     while(1)
339     {
340         fPythia->Pyevnt();
341         if (gAlice->GetEvNumber()>=fDebugEventFirst &&
342             gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
343         fTrials++;
344         
345         fPythia->ImportParticles(fParticles,"All");
346
347 //
348 //
349 //
350         Int_t i;
351         
352         Int_t np = fParticles->GetEntriesFast();
353         if (np == 0 ) continue;
354 // Get event vertex and discard the event if the Z coord. is too big    
355         TParticle *iparticle = (TParticle *) fParticles->At(0);
356         Float_t distz = iparticle->Vz()/10.;
357         if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
358 //
359         fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
360         fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
361         fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
362 //
363         Int_t* pParent   = new Int_t[np];
364         Int_t* pSelected = new Int_t[np];
365         Int_t* trackIt   = new Int_t[np];
366         for (i=0; i< np; i++) {
367             pParent[i]   = -1;
368             pSelected[i] =  0;
369             trackIt[i]   =  0;
370         }
371         // printf("\n **************************************************%d\n",np);
372         Int_t nc = 0;        // Total n. of selected particles
373         Int_t nParents = 0;  // Selected parents
374         Int_t nTkbles = 0;   // Trackable particles
375         if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
376             
377             for (i = 0; i<np; i++) {
378                 iparticle = (TParticle *) fParticles->At(i);
379                 Int_t ks = iparticle->GetStatusCode();
380                 kf = CheckPDGCode(iparticle->GetPdgCode());
381 // No initial state partons
382                 if (ks==21) continue;
383 //
384 // Heavy Flavor Selection
385 //
386                 // quark ?
387                 kf = TMath::Abs(kf);
388                 Int_t kfl = kf;
389                 // meson ?
390                 if  (kfl > 10) kfl/=100;
391                 // baryon
392                 if (kfl > 10) kfl/=10;
393                 if (kfl > 10) kfl/=10;
394
395                 Int_t ipa = iparticle->GetFirstMother()-1;
396                 Int_t kfMo = 0;
397                 
398                 if (ipa > -1) {
399                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
400                     kfMo = TMath::Abs(mother->GetPdgCode());
401                 }
402 //              printf("\n particle (all)  %d %d %d", i, pSelected[i], kf);
403                 // What to keep in Stack?
404                 Bool_t flavorOK = kFALSE;
405                 Bool_t selectOK = kFALSE;
406                 if (fFeedDownOpt) {
407                   if (kfl >= fFlavorSelect) flavorOK = kTRUE;
408                 } else {
409                   if (kfl > fFlavorSelect) {
410                     nc = -1;
411                     break;
412                   }
413                   if (kfl == fFlavorSelect) flavorOK = kTRUE;
414                 }
415                 switch (fStackFillOpt) {
416                 case kFlavorSelection:
417                   selectOK = kTRUE;
418                   break;
419                 case kParentSelection:
420                   if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
421                   break;
422                 }
423                 if (flavorOK && selectOK) { 
424 //
425 // Heavy flavor hadron or quark
426 //
427 // Kinematic seletion on final state heavy flavor mesons
428                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
429                     {
430                       continue;
431                     }
432                     pSelected[i] = 1;
433                     if (ParentSelected(kf)) ++nParents; // Update parent count
434 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
435                 } else {
436 // Kinematic seletion on decay products
437                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
438                         && !KinematicSelection(iparticle, 1))
439                     {
440                       continue;
441                     }
442 //
443 // Decay products 
444 // Select if mother was selected and is not tracked
445
446                     if (pSelected[ipa] && 
447                         !trackIt[ipa]  &&     // mother will be  tracked ?
448                         kfMo !=  5 &&         // mother is b-quark, don't store fragments          
449                         kfMo !=  4 &&         // mother is c-quark, don't store fragments 
450                         kf   != 92)           // don't store string
451                     {
452 //
453 // Semi-stable or de-selected: diselect decay products:
454 // 
455 //
456                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
457                         {
458                             Int_t ipF = iparticle->GetFirstDaughter();
459                             Int_t ipL = iparticle->GetLastDaughter();   
460                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
461                         }
462 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
463                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
464                     }
465                 }
466                 if (pSelected[i] == -1) pSelected[i] = 0;
467                 if (!pSelected[i]) continue;
468                 // Count quarks only if you did not include fragmentation
469                 if (fFragmentation && kf <= 10) continue;
470                 nc++;
471 // Decision on tracking
472                 trackIt[i] = 0;
473 //
474 // Track final state particle
475                 if (ks == 1) trackIt[i] = 1;
476 // Track semi-stable particles
477                 if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
478 // Track particles selected by process if undecayed. 
479                 if (fForceDecay == kNoDecay) {
480                     if (ParentSelected(kf)) trackIt[i] = 1;
481                 } else {
482                     if (ParentSelected(kf)) trackIt[i] = 0;
483                 }
484                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
485 //
486 //
487
488             } // particle selection loop
489             if (nc > 0) {
490                 for (i = 0; i<np; i++) {
491                     if (!pSelected[i]) continue;
492                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
493                     kf = CheckPDGCode(iparticle->GetPdgCode());
494                     Int_t ks = iparticle->GetStatusCode();  
495                     p[0] = iparticle->Px();
496                     p[1] = iparticle->Py();
497                     p[2] = iparticle->Pz();
498                     origin[0] = fOrigin[0]+iparticle->Vx()/10.;
499                     origin[1] = fOrigin[1]+iparticle->Vy()/10.;
500                     origin[2] = fOrigin[2]+iparticle->Vz()/10.;
501                     Float_t tof   = kconv*iparticle->T();
502                     Int_t ipa     = iparticle->GetFirstMother()-1;
503                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
504                     SetTrack(fTrackIt*trackIt[i] ,
505                                      iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
506                     pParent[i] = nt;
507                     KeepTrack(nt); 
508                 } //  SetTrack loop
509             }
510         } else {
511             nc = GenerateMB();
512         } // mb ?
513
514         if (pParent)   delete[] pParent;
515         if (pSelected) delete[] pSelected;
516         if (trackIt)   delete[] trackIt;
517
518         if (nc > 0) {
519           switch (fCountMode) {
520           case kCountAll:
521             // printf(" Count all \n");
522             jev += nc;
523             break;
524           case kCountParents:
525             // printf(" Count parents \n");
526             jev += nParents;
527             break;
528           case kCountTrackables:
529             // printf(" Count trackable \n");
530             jev += nTkbles;
531             break;
532           }
533             if (jev >= fNpart || fNpart == -1) {
534                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
535                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
536                 MakeHeader();
537                 break;
538             }
539         }
540     } // event loop
541     SetHighWaterMark(nt);
542 //  adjust weight due to kinematic selection
543     AdjustWeights();
544 //  get cross-section
545     fXsection=fPythia->GetPARI(1);
546 }
547
548 Int_t  AliGenPythia::GenerateMB()
549 {
550 //
551 // Min Bias selection and other global selections
552 //
553     Int_t i, kf, nt, iparent;
554     Int_t nc = 0;
555     Float_t p[3];
556     Float_t polar[3]   =   {0,0,0};
557     Float_t origin[3]  =   {0,0,0};
558 //  converts from mm/c to s
559     const Float_t kconv=0.001/2.999792458e8;
560     
561     Int_t np = fParticles->GetEntriesFast();
562     Int_t* pParent = new Int_t[np];
563     for (i=0; i< np; i++) pParent[i] = -1;
564     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
565         TParticle* jet1 = (TParticle *) fParticles->At(6);
566         TParticle* jet2 = (TParticle *) fParticles->At(7);
567         if (!CheckTrigger(jet1, jet2)) return 0;
568     }
569     
570     for (i = 0; i<np; i++) {
571         Int_t trackIt = 0;
572         TParticle *  iparticle = (TParticle *) fParticles->At(i);
573         kf = CheckPDGCode(iparticle->GetPdgCode());
574         Int_t ks = iparticle->GetStatusCode();
575         Int_t km = iparticle->GetFirstMother();
576 //      printf("\n Particle: %d %d %d", i, kf, ks);
577         
578         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
579             (ks != 1) ||
580             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
581             nc++;
582             if (ks == 1) trackIt = 1;
583             Int_t ipa = iparticle->GetFirstMother()-1;
584
585             iparent = (ipa > -1) ? pParent[ipa] : -1;
586
587 //
588 // store track information
589             p[0] = iparticle->Px();
590             p[1] = iparticle->Py();
591             p[2] = iparticle->Pz();
592             origin[0] = fOrigin[0]+iparticle->Vx()/10.;
593             origin[1] = fOrigin[1]+iparticle->Vy()/10.;
594             origin[2] = fOrigin[2]+iparticle->Vz()/10.;
595             Float_t tof=kconv*iparticle->T();
596             SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
597                          tof, kPPrimary, nt, 1., ks);
598             KeepTrack(nt);
599             pParent[i] = nt;
600         } // select particle
601     } // particle loop 
602
603     if (pParent) delete[] pParent;
604     
605     printf("\n I've put %i particles on the stack \n",nc);
606     return nc;
607 }
608
609
610 void AliGenPythia::FinishRun()
611 {
612 // Print x-section summary
613     fPythia->Pystat(1);
614 }
615
616 void AliGenPythia::AdjustWeights()
617 {
618 // Adjust the weights after generation of all events
619 //
620     TParticle *part;
621     Int_t ntrack=gAlice->GetNtrack();
622     for (Int_t i=0; i<ntrack; i++) {
623         part= gAlice->Particle(i);
624         part->SetWeight(part->GetWeight()*fKineBias);
625     }
626 }
627     
628 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
629 {
630 // Treat protons as inside nuclei with mass numbers a1 and a2  
631     fNucA1 = a1;
632     fNucA2 = a2;
633 }
634
635
636 void AliGenPythia::MakeHeader()
637 {
638 // Builds the event header, to be called after each event
639     AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
640     ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
641     ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
642     header->SetPrimaryVertex(fEventVertex);
643     gAlice->SetGenEventHeader(header);
644 }
645         
646
647 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) const
648 {
649 // Check the kinematic trigger condition
650 //
651     Double_t eta[2];
652     eta[0] = jet1->Eta();
653     eta[1] = jet2->Eta();
654     Double_t phi[2];
655     phi[0] = jet1->Phi();
656     phi[1] = jet2->Phi();
657     Int_t    pdg[2]; 
658     pdg[0] = jet1->GetPdgCode();
659     pdg[1] = jet2->GetPdgCode();    
660     Bool_t   triggered = kFALSE;
661
662     if (fProcess == kPyJets) {
663         //Check eta range first...
664         if ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet) ||
665             (eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet))
666         {
667             //Eta is okay, now check phi range
668             if ((phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet) ||
669                 (phi[1] < fPhiMaxJet && phi[1] > fPhiMinJet))
670             {
671                 triggered = kTRUE;
672             }
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     
694     
695     return triggered;
696 }
697           
698 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
699 {
700 // Assignment operator
701     return *this;
702 }
703
704
705
706 #ifdef never
707 void AliGenPythia::Streamer(TBuffer &R__b)
708 {
709    // Stream an object of class AliGenPythia.
710
711    if (R__b.IsReading()) {
712       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
713       AliGenerator::Streamer(R__b);
714       R__b >> (Int_t&)fProcess;
715       R__b >> (Int_t&)fStrucFunc;
716       R__b >> (Int_t&)fForceDecay;
717       R__b >> fEnergyCMS;
718       R__b >> fKineBias;
719       R__b >> fTrials;
720       fParentSelect.Streamer(R__b);
721       fChildSelect.Streamer(R__b);
722       R__b >> fXsection;
723 //      (AliPythia::Instance())->Streamer(R__b);
724       R__b >> fPtHardMin;
725       R__b >> fPtHardMax;
726 //      if (fDecayer) fDecayer->Streamer(R__b);
727    } else {
728       R__b.WriteVersion(AliGenPythia::IsA());
729       AliGenerator::Streamer(R__b);
730       R__b << (Int_t)fProcess;
731       R__b << (Int_t)fStrucFunc;
732       R__b << (Int_t)fForceDecay;
733       R__b << fEnergyCMS;
734       R__b << fKineBias;
735       R__b << fTrials;
736       fParentSelect.Streamer(R__b);
737       fChildSelect.Streamer(R__b);
738       R__b << fXsection;
739 //      R__b << fPythia;
740       R__b << fPtHardMin;
741       R__b << fPtHardMax;
742       //     fDecayer->Streamer(R__b);
743    }
744 }
745 #endif
746