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