Use local SetTrack, KeepTrack and SetHighWaterMark methods
[u/mrichter/AliRoot.git] / EVGEN / 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.38  2001/07/13 10:58:54  morsch
19 - Some coded moved to AliGenMC
20 - Improved handling of secondary vertices.
21
22 Revision 1.37  2001/06/28 11:17:28  morsch
23 SetEventListRange setter added. Events in specified range are listed for
24 debugging. (Yuri Kharlov)
25
26 Revision 1.36  2001/03/30 07:05:49  morsch
27 Final print-out in finish run.
28 Write parton system for jet-production (preliminary solution).
29
30 Revision 1.35  2001/03/09 13:03:40  morsch
31 Process_t and Struc_Func_t moved to AliPythia.h
32
33 Revision 1.34  2001/02/14 15:50:40  hristov
34 The last particle in event marked using SetHighWaterMark
35
36 Revision 1.33  2001/01/30 09:23:12  hristov
37 Streamers removed (R.Brun)
38
39 Revision 1.32  2001/01/26 19:55:51  hristov
40 Major upgrade of AliRoot code
41
42 Revision 1.31  2001/01/17 10:54:31  hristov
43 Better protection against FPE
44
45 Revision 1.30  2000/12/18 08:55:35  morsch
46 Make AliPythia dependent generartors work with new scheme of random number generation
47
48 Revision 1.29  2000/12/04 11:22:03  morsch
49 Init of sRandom as in 1.15
50
51 Revision 1.28  2000/12/02 11:41:39  morsch
52 Use SetRandom() to initialize random number generator in constructor.
53
54 Revision 1.27  2000/11/30 20:29:02  morsch
55 Initialise static variable sRandom in constructor: sRandom = fRandom;
56
57 Revision 1.26  2000/11/30 07:12:50  alibrary
58 Introducing new Rndm and QA classes
59
60 Revision 1.25  2000/10/18 19:11:27  hristov
61 Division by zero fixed
62
63 Revision 1.24  2000/09/18 10:41:35  morsch
64 Add possibility to use nuclear structure functions from PDF library V8.
65
66 Revision 1.23  2000/09/14 14:05:40  morsch
67 dito
68
69 Revision 1.22  2000/09/14 14:02:22  morsch
70 - Correct conversion from mm to cm when passing particle vertex to MC.
71 - Correct handling of fForceDecay == all.
72
73 Revision 1.21  2000/09/12 14:14:55  morsch
74 Call fDecayer->ForceDecay() at the beginning of Generate().
75
76 Revision 1.20  2000/09/06 14:29:33  morsch
77 Use AliPythia for event generation an AliDecayPythia for decays.
78 Correct handling of "nodecay" option
79
80 Revision 1.19  2000/07/11 18:24:56  fca
81 Coding convention corrections + few minor bug fixes
82
83 Revision 1.18  2000/06/30 12:40:34  morsch
84 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
85 Geant units (cm).
86
87 Revision 1.17  2000/06/09 20:34:07  morsch
88 All coding rule violations except RS3 corrected
89
90 Revision 1.16  2000/05/15 15:04:20  morsch
91 The full event is written for fNtrack = -1
92 Coding rule violations corrected.
93
94 Revision 1.15  2000/04/26 10:14:24  morsch
95 Particles array has one entry more than pythia particle list. Upper bound of
96 particle loop changed to np-1 (R. Guernane, AM)
97
98 Revision 1.14  2000/04/05 08:36:13  morsch
99 Check status code of particles in Pythia event
100 to avoid double counting as partonic state and final state particle.
101
102 Revision 1.13  1999/11/09 07:38:48  fca
103 Changes for compatibility with version 2.23 of ROOT
104
105 Revision 1.12  1999/11/03 17:43:20  fca
106 New version from G.Martinez & A.Morsch
107
108 Revision 1.11  1999/09/29 09:24:14  fca
109 Introduction of the Copyright and cvs Log
110 */
111
112 #include "AliGenPythia.h"
113 #include "AliGenPythiaEventHeader.h"
114 #include "AliDecayerPythia.h"
115 #include "AliRun.h"
116 #include "AliPythia.h"
117 #include "AliPDG.h"
118 #include <TParticle.h>
119 #include <TSystem.h>
120
121  ClassImp(AliGenPythia)
122
123 AliGenPythia::AliGenPythia()
124                  :AliGenMC()
125 {
126 // Default Constructor
127   fDecayer = new AliDecayerPythia();
128   SetEventListRange();
129 }
130
131 AliGenPythia::AliGenPythia(Int_t npart)
132                  :AliGenMC(npart)
133 {
134 // default charm production at 5. 5 TeV
135 // semimuonic decay
136 // structure function GRVHO
137 //
138     fXsection  = 0.;
139     fNucA1=0;
140     fNucA2=0;
141     SetProcess();
142     SetStrucFunc();
143     SetForceDecay();
144     SetPtHard();
145     SetEnergyCMS();
146     fDecayer = new AliDecayerPythia();
147     // Set random number generator 
148     sRandom=fRandom;
149     SetEventListRange();
150     fFlavorSelect   = 0;
151     // Produced particles  
152     fParticles = new TClonesArray("TParticle",1000);
153 }
154
155 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
156 {
157 // copy constructor
158 }
159
160 AliGenPythia::~AliGenPythia()
161 {
162 // Destructor
163 }
164
165 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
166 {
167   // Set a range of event numbers, for which a table
168   // of generated particle will be printed
169   fDebugEventFirst = eventFirst;
170   fDebugEventLast  = eventLast;
171   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
172 }
173
174 void AliGenPythia::Init()
175 {
176 // Initialisation
177   SetMC(AliPythia::Instance());
178     fPythia=(AliPythia*) fgMCEvGen;
179 //
180     fParentWeight=1./Float_t(fNpart);
181 //
182 //  Forward Paramters to the AliPythia object
183     fDecayer->SetForceDecay(fForceDecay);    
184     fDecayer->Init();
185
186
187     fPythia->SetCKIN(3,fPtHardMin);
188     fPythia->SetCKIN(4,fPtHardMax);    
189     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
190     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
191
192     //    fPythia->Pylist(0);
193     //    fPythia->Pystat(2);
194 //  Parent and Children Selection
195     switch (fProcess) 
196     {
197     case kPyCharm:
198         fParentSelect[0] =  411;
199         fParentSelect[1] =  421;
200         fParentSelect[2] =  431;
201         fParentSelect[3] = 4122;        
202         fFlavorSelect    = 4;
203         break;
204     case kPyCharmUnforced:
205         fParentSelect[0] =  411;
206         fParentSelect[1] =  421;
207         fParentSelect[2] =  431;
208         fParentSelect[3]=  4122;
209         fFlavorSelect    = 4;   
210         break;
211     case kPyBeauty:
212         fParentSelect[0]=  511;
213         fParentSelect[1]=  521;
214         fParentSelect[2]=  531;
215         fParentSelect[3]= 5122;
216         fParentSelect[4]= 5132;
217         fParentSelect[5]= 5232;
218         fParentSelect[6]= 5332;
219         fFlavorSelect   = 5;    
220         break;
221     case kPyBeautyUnforced:
222         fParentSelect[0] =  511;
223         fParentSelect[1] =  521;
224         fParentSelect[2] =  531;
225         fParentSelect[3] = 5122;
226         fParentSelect[4] = 5132;
227         fParentSelect[5] = 5232;
228         fParentSelect[6] = 5332;
229         fFlavorSelect    = 5;   
230         break;
231     case kPyJpsiChi:
232     case kPyJpsi:
233         fParentSelect[0] = 443;
234         break;
235     case kPyMb:
236     case kPyJets:
237     case kPyDirectGamma:
238         break;
239     }
240     AliGenMC::Init();
241 }
242
243 void AliGenPythia::Generate()
244 {
245 // Generate one event
246     fDecayer->ForceDecay();
247
248     Float_t polar[3]   =   {0,0,0};
249     Float_t origin[3]  =   {0,0,0};
250     Float_t origin0[3] =   {0,0,0};
251     Float_t p[3];
252 //  converts from mm/c to s
253     const Float_t kconv=0.001/2.999792458e8;
254 //
255     Int_t nt=0;
256     Int_t jev=0;
257     Int_t j, kf;
258     fTrials=0;
259
260 //  Set collision vertex position 
261     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
262     if(fVertexSmear==kPerEvent) {
263         fPythia->SetMSTP(151,1);
264         for (j=0;j<3;j++) {
265             fPythia->SetPARP(151+j, fOsigma[j]*10.);
266         }
267     } else if (fVertexSmear==kPerTrack) {
268         fPythia->SetMSTP(151,0);
269     }
270 //  event loop    
271     while(1)
272     {
273         fPythia->Pyevnt();
274         if (gAlice->GetEvNumber()>=fDebugEventFirst &&
275             gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
276         fTrials++;
277         fPythia->ImportParticles(fParticles,"All");
278
279 //
280 //
281 //
282         Int_t i;
283         
284         Int_t np = fParticles->GetEntriesFast();
285         Int_t* pParent   = new Int_t[np];
286         Int_t* pSelected = new Int_t[np];
287         Int_t* trackIt   = new Int_t[np];
288         for (i=0; i< np-1; i++) {
289             pParent[i]   = -1;
290             pSelected[i] =  0;
291         }
292         printf("\n **************************************************%d\n",np);
293         Int_t nc = 0;
294         if (np == 0 ) continue;
295         if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
296             for (i = 0; i<np-1; i++) {
297                 TParticle *  iparticle = (TParticle *) fParticles->At(i);
298                 Int_t ks = iparticle->GetStatusCode();
299                 kf = CheckPDGCode(iparticle->GetPdgCode());
300 // No initial state partons
301                 if (ks==21) continue;
302 //
303 // Heavy Flavor Selection
304 //
305                 // quark ?
306                 kf = TMath::Abs(kf);
307                 Int_t kfl = kf;
308                 // meson ?
309                 if  (kfl > 10) kfl/=100;
310                 // baryon
311                 if (kfl > 10) kfl/=10;
312                 if (kfl > 10) kfl/=10;
313
314                 Int_t ipa = iparticle->GetFirstMother()-1;
315                 Int_t kfMo = 0;
316                 
317                 if (ipa > -1) {
318                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
319                     kfMo = TMath::Abs(mother->GetPdgCode());
320                 }
321 //              printf("\n particle (all)  %d %d %d", i, pSelected[i], kf);
322                 if (kfl >= fFlavorSelect) { 
323 //
324 // Heavy flavor hadron or quark
325 //
326 // Kinematic seletion on final state heavy flavor mesons
327                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) 
328                     {
329                         nc = -1;
330                         break;
331                     }
332                     pSelected[i] = 1;
333 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
334                 } else {
335 // Kinematic seletion on decay products
336                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) 
337                         && !KinematicSelection(iparticle, 1))
338                     {
339                         nc = -1;
340                         break;
341                     }
342 //
343 // Decay products 
344 // Select if mother was selected and is not tracked
345
346                     if (pSelected[ipa] && 
347                         !trackIt[ipa]  &&     // mother will be  tracked ?
348                         kfMo !=  5 &&         // mother is b-quark, don't store fragments          
349                         kfMo !=  4 &&         // mother is c-quark, don't store fragments 
350                         kf   != 92)           // don't store string
351                     {
352 //
353 // Semi-stable or de-selected: diselect decay products:
354 // 
355 //
356                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > 10e-15)
357                         {
358                             Int_t ipF = iparticle->GetFirstDaughter();
359                             Int_t ipL = iparticle->GetLastDaughter();   
360                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
361                         }
362 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
363                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
364                     }
365                 }
366                 if (pSelected[i] == -1) pSelected[i] = 0;
367                 if (!pSelected[i]) continue;
368                 nc++;
369 // Decision on tracking
370                 trackIt[i] = 0;
371 //
372 // Track final state particle
373                 if (ks == 1) trackIt[i] = 1;
374 // Track semi-stable particles
375                 if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15))  trackIt[i] = 1;
376 // Track particles selected by process if undecayed. 
377                 if (fForceDecay == kNoDecay) {
378                     if (ParentSelected(kf)) trackIt[i] = 1;
379                 } else {
380                     if (ParentSelected(kf)) trackIt[i] = 0;
381                 }
382 //
383 //
384
385             } // particle selection loop
386             if (nc > -1) {
387                 for (i = 0; i<np-1; i++) {
388                     if (!pSelected[i]) continue;
389                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
390                     kf = CheckPDGCode(iparticle->GetPdgCode());
391                     p[0]=iparticle->Px();
392                     p[1]=iparticle->Py();
393                     p[2]=iparticle->Pz();
394                     origin[0]=origin0[0]+iparticle->Vx()/10.;
395                     origin[1]=origin0[1]+iparticle->Vy()/10.;
396                     origin[2]=origin0[2]+iparticle->Vz()/10.;
397                     Float_t tof=kconv*iparticle->T();
398                     Int_t ipa = iparticle->GetFirstMother()-1;
399                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
400 //                  printf("\n SetTrack %d %d %d %d", i, iparent, kf, trackIt[i]);
401                     SetTrack(fTrackIt*trackIt[i] ,
402                                      iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
403                     pParent[i] = nt;
404                     KeepTrack(nt); 
405                 } //  SetTrack loop
406             }
407         } else {
408             nc = GenerateMB();
409         } // mb ?
410
411         if (nc > 0) {
412             jev+=nc;
413             if (jev >= fNpart || fNpart == -1) {
414                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
415                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
416                 break;
417             }
418         }
419     } // event loop
420     SetHighWaterMark(nt);
421 //  adjust weight due to kinematic selection
422     AdjustWeights();
423 //  get cross-section
424     fXsection=fPythia->GetPARI(1);
425 }
426
427 Int_t  AliGenPythia::GenerateMB()
428 {
429     Int_t i, kf, nt, iparent;
430     Int_t nc = 0;
431     Float_t p[3];
432     Float_t polar[3]   =   {0,0,0};
433     Float_t origin[3]  =   {0,0,0};
434     Float_t origin0[3] =   {0,0,0};
435 //  converts from mm/c to s
436     const Float_t kconv=0.001/2.999792458e8;
437     
438     Int_t np = fParticles->GetEntriesFast();
439     Int_t* pParent = new Int_t[np];
440     for (i=0; i< np-1; i++) pParent[i] = -1;
441
442     for (i = 0; i<np-1; i++) {
443         Int_t trackIt = 0;
444         TParticle *  iparticle = (TParticle *) fParticles->At(i);
445         kf = CheckPDGCode(iparticle->GetPdgCode());
446         Int_t ks = iparticle->GetStatusCode();
447         Int_t km = iparticle->GetFirstMother();
448 //      printf("\n Particle: %d %d %d", i, kf, ks);
449         
450         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
451             (ks != 1) ||
452             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
453             nc++;
454             if (ks == 1) trackIt = 1;
455             Int_t ipa = iparticle->GetFirstMother()-1;
456
457             iparent = (ipa > -1) ? pParent[ipa] : -1;
458
459 //
460 // store track information
461             p[0]=iparticle->Px();
462             p[1]=iparticle->Py();
463             p[2]=iparticle->Pz();
464             origin[0]=origin0[0]+iparticle->Vx()/10.;
465             origin[1]=origin0[1]+iparticle->Vy()/10.;
466             origin[2]=origin0[2]+iparticle->Vz()/10.;
467             Float_t tof=kconv*iparticle->T();
468             SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
469                              tof, kPPrimary, nt);
470             KeepTrack(nt);
471             pParent[i] = nt;
472         } // select particle
473     } // particle loop 
474
475     if (pParent) delete[] pParent;
476     
477     printf("\n I've put %i particles on the stack \n",nc);
478     MakeHeader();
479     return nc;
480 }
481
482
483 void AliGenPythia::FinishRun()
484 {
485 // Print x-section summary
486     fPythia->Pystat(1);
487 }
488
489 void AliGenPythia::AdjustWeights()
490 {
491 // Adjust the weights after generation of all events
492 //
493     TParticle *part;
494     Int_t ntrack=gAlice->GetNtrack();
495     for (Int_t i=0; i<ntrack; i++) {
496         part= gAlice->Particle(i);
497         part->SetWeight(part->GetWeight()*fKineBias);
498     }
499 }
500     
501 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
502 {
503 // Treat protons as inside nuclei with mass numbers a1 and a2  
504     fNucA1 = a1;
505     fNucA2 = a2;
506 }
507
508
509 void AliGenPythia::MakeHeader()
510 {
511 // Builds the event header, to be called after each event
512     AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
513    ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
514    gAlice->SetGenEventHeader(header);
515 }
516         
517           
518 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
519 {
520 // Assignment operator
521     return *this;
522 }
523
524
525
526 #ifdef never
527 void AliGenPythia::Streamer(TBuffer &R__b)
528 {
529    // Stream an object of class AliGenPythia.
530
531    if (R__b.IsReading()) {
532       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
533       AliGenerator::Streamer(R__b);
534       R__b >> (Int_t&)fProcess;
535       R__b >> (Int_t&)fStrucFunc;
536       R__b >> (Int_t&)fForceDecay;
537       R__b >> fEnergyCMS;
538       R__b >> fKineBias;
539       R__b >> fTrials;
540       fParentSelect.Streamer(R__b);
541       fChildSelect.Streamer(R__b);
542       R__b >> fXsection;
543 //      (AliPythia::Instance())->Streamer(R__b);
544       R__b >> fPtHardMin;
545       R__b >> fPtHardMax;
546 //      if (fDecayer) fDecayer->Streamer(R__b);
547    } else {
548       R__b.WriteVersion(AliGenPythia::IsA());
549       AliGenerator::Streamer(R__b);
550       R__b << (Int_t)fProcess;
551       R__b << (Int_t)fStrucFunc;
552       R__b << (Int_t)fForceDecay;
553       R__b << fEnergyCMS;
554       R__b << fKineBias;
555       R__b << fTrials;
556       fParentSelect.Streamer(R__b);
557       fChildSelect.Streamer(R__b);
558       R__b << fXsection;
559 //      R__b << fPythia;
560       R__b << fPtHardMin;
561       R__b << fPtHardMax;
562       //     fDecayer->Streamer(R__b);
563    }
564 }
565 #endif
566