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