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