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