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