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