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