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