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