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