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