Use SetRandom() to initialize random number generator in constructor.
[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.27  2000/11/30 20:29:02  morsch
19 Initialise static variable sRandom in constructor: sRandom = fRandom;
20
21 Revision 1.26  2000/11/30 07:12:50  alibrary
22 Introducing new Rndm and QA classes
23
24 Revision 1.25  2000/10/18 19:11:27  hristov
25 Division by zero fixed
26
27 Revision 1.24  2000/09/18 10:41:35  morsch
28 Add possibility to use nuclear structure functions from PDF library V8.
29
30 Revision 1.23  2000/09/14 14:05:40  morsch
31 dito
32
33 Revision 1.22  2000/09/14 14:02:22  morsch
34 - Correct conversion from mm to cm when passing particle vertex to MC.
35 - Correct handling of fForceDecay == all.
36
37 Revision 1.21  2000/09/12 14:14:55  morsch
38 Call fDecayer->ForceDecay() at the beginning of Generate().
39
40 Revision 1.20  2000/09/06 14:29:33  morsch
41 Use AliPythia for event generation an AliDecayPythia for decays.
42 Correct handling of "nodecay" option
43
44 Revision 1.19  2000/07/11 18:24:56  fca
45 Coding convention corrections + few minor bug fixes
46
47 Revision 1.18  2000/06/30 12:40:34  morsch
48 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
49 Geant units (cm).
50
51 Revision 1.17  2000/06/09 20:34:07  morsch
52 All coding rule violations except RS3 corrected
53
54 Revision 1.16  2000/05/15 15:04:20  morsch
55 The full event is written for fNtrack = -1
56 Coding rule violations corrected.
57
58 Revision 1.15  2000/04/26 10:14:24  morsch
59 Particles array has one entry more than pythia particle list. Upper bound of
60 particle loop changed to np-1 (R. Guernane, AM)
61
62 Revision 1.14  2000/04/05 08:36:13  morsch
63 Check status code of particles in Pythia event
64 to avoid double counting as partonic state and final state particle.
65
66 Revision 1.13  1999/11/09 07:38:48  fca
67 Changes for compatibility with version 2.23 of ROOT
68
69 Revision 1.12  1999/11/03 17:43:20  fca
70 New version from G.Martinez & A.Morsch
71
72 Revision 1.11  1999/09/29 09:24:14  fca
73 Introduction of the Copyright and cvs Log
74 */
75
76 #include "AliGenPythia.h"
77 #include "AliDecayerPythia.h"
78 #include "AliRun.h"
79 #include "AliPythia.h"
80 #include "AliPDG.h"
81 #include <TParticle.h>
82 #include <TSystem.h>
83
84  ClassImp(AliGenPythia)
85
86 AliGenPythia::AliGenPythia()
87                  :AliGenerator()
88 {
89 // Default Constructor
90   fDecayer = new AliDecayerPythia();
91 }
92
93 AliGenPythia::AliGenPythia(Int_t npart)
94                  :AliGenerator(npart)
95 {
96 // default charm production at 5. 5 TeV
97 // semimuonic decay
98 // structure function GRVHO
99 //
100     fXsection  = 0.;
101     fNucA1=0;
102     fNucA2=0;
103     fParentSelect.Set(5);
104     fChildSelect.Set(5);
105     for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
106     SetProcess();
107     SetStrucFunc();
108     SetForceDecay();
109     SetPtHard();
110     SetEnergyCMS();
111     fDecayer = new AliDecayerPythia();
112     // Set random number generator 
113     SetRandom();
114 }
115
116 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
117 {
118 // copy constructor
119 }
120
121 AliGenPythia::~AliGenPythia()
122 {
123 // Destructor
124 }
125
126 void AliGenPythia::Init()
127 {
128 // Initialisation
129   SetMC(AliPythia::Instance());
130     fPythia=(AliPythia*) fgMCEvGen;
131 //
132     fParentWeight=1./Float_t(fNpart);
133 //
134 //  Forward Paramters to the AliPythia object
135     //    gSystem->Exec("ln -s $ALICE_ROOT/data/Decay.table fort.1");
136     //    fPythia->Pyupda(2,1);    
137     //    gSystem->Exec("rm fort.1");
138
139     fDecayer->SetForceDecay(fForceDecay);    
140     fDecayer->Init();
141
142
143     fPythia->SetCKIN(3,fPtHardMin);
144     fPythia->SetCKIN(4,fPtHardMax);    
145     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);  
146     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
147
148     //    fPythia->Pylist(0);
149     //    fPythia->Pystat(2);
150 //  Parent and Children Selection
151     switch (fProcess) 
152     {
153     case charm:
154
155         fParentSelect[0]=411;
156         fParentSelect[1]=421;
157         fParentSelect[2]=431;
158         fParentSelect[3]=4122;  
159         break;
160     case charm_unforced:
161
162         fParentSelect[0]=411;
163         fParentSelect[1]=421;
164         fParentSelect[2]=431;
165         fParentSelect[3]=4122;  
166         break;
167     case beauty:
168         fParentSelect[0]=511;
169         fParentSelect[1]=521;
170         fParentSelect[2]=531;
171         fParentSelect[3]=5122;  
172         break;
173     case beauty_unforced:
174         fParentSelect[0]=511;
175         fParentSelect[1]=521;
176         fParentSelect[2]=531;
177         fParentSelect[3]=5122;  
178         break;
179     case jpsi_chi:
180     case jpsi:
181         fParentSelect[0]=443;
182         break;
183     case mb:
184         break;
185     }
186
187     switch (fForceDecay) 
188     {
189     case semielectronic:
190     case dielectron:
191     case b_jpsi_dielectron:
192     case b_psip_dielectron:
193         fChildSelect[0]=kElectron;      
194         break;
195     case semimuonic:
196     case dimuon:
197     case b_jpsi_dimuon:
198     case b_psip_dimuon:
199     case pitomu:
200     case katomu:
201         fChildSelect[0]=kMuonMinus;
202         break;
203     case hadronicD:
204       fChildSelect[0]=kPiPlus;
205       fChildSelect[1]=kKPlus;
206       break;
207     case all:
208     case nodecay:
209       break;
210     }
211 }
212
213 void AliGenPythia::Generate()
214 {
215 // Generate one event
216     fDecayer->ForceDecay();
217
218     Float_t polar[3] =   {0,0,0};
219     Float_t origin[3]=   {0,0,0};
220     Float_t originP[3]= {0,0,0};
221     Float_t origin0[3]=  {0,0,0};
222     Float_t p[3], pP[4];
223 //    Float_t random[6];
224     static TClonesArray *particles;
225 //  converts from mm/c to s
226     const Float_t kconv=0.001/2.999792458e8;
227     
228     
229 //
230     Int_t nt=0;
231     Int_t ntP=0;
232     Int_t jev=0;
233     Int_t j, kf;
234
235     if(!particles) particles=new TClonesArray("TParticle",1000);
236     
237     fTrials=0;
238     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
239     if(fVertexSmear==kPerEvent) {
240         fPythia->SetMSTP(151,1);
241         for (j=0;j<3;j++) {
242             fPythia->SetPARP(151+j, fOsigma[j]/10.);
243         }
244     } else if (fVertexSmear==kPerTrack) {
245         fPythia->SetMSTP(151,0);
246     }
247     
248     while(1)
249     {
250         fPythia->Pyevnt();
251 //      fPythia->Pylist(1);
252         fTrials++;
253         fPythia->ImportParticles(particles,"All");
254         Int_t np = particles->GetEntriesFast();
255         printf("\n **************************************************%d\n",np);
256         Int_t nc=0;
257         if (np == 0 ) continue;
258         if (fProcess != mb) {
259             for (Int_t i = 0; i<np-1; i++) {
260                 TParticle *  iparticle = (TParticle *) particles->At(i);
261                 Int_t ks = iparticle->GetStatusCode();
262                 kf = CheckPDGCode(iparticle->GetPdgCode());
263                 if (ks==21) continue;
264
265                 fChildWeight=(fDecayer->GetPartialBranchingRatio(kf))*fParentWeight;      
266 //
267 // Parent
268                 if (ParentSelected(TMath::Abs(kf))) {
269                     if (KinematicSelection(iparticle)) {
270                         if (nc==0) {
271 //
272 // Store information concerning the hard scattering process
273 //
274                             Float_t massP  = fPythia->GetPARI(13);
275                             Float_t   ptP  = fPythia->GetPARI(17);
276                             Float_t    yP  = fPythia->GetPARI(37);
277                             Float_t  xmtP  = sqrt(ptP*ptP+massP*massP);
278                             Float_t    ty  = Float_t(TMath::TanH(yP));
279                             pP[0] = ptP;
280                             pP[1] = 0;
281                             pP[2] = xmtP*ty/sqrt(1.-ty*ty);
282                             pP[3] = massP;
283                             gAlice->SetTrack(0,-1,-1,
284                                              pP,originP,polar,
285                                              0,kPPrimary,ntP,fParentWeight);
286 //                                           0,"Hard Scat.",ntP,fParentWeight);
287                             gAlice->KeepTrack(ntP);
288                         }
289                         nc++;
290 //
291 // store parent track information
292                         p[0]=iparticle->Px();
293                         p[1]=iparticle->Py();
294                         p[2]=iparticle->Pz();
295                         origin[0]=origin0[0]+iparticle->Vx()/10.;
296                         origin[1]=origin0[1]+iparticle->Vy()/10.;
297                         origin[2]=origin0[2]+iparticle->Vz()/10.;
298
299                         Int_t ifch=iparticle->GetFirstDaughter();
300                         Int_t ilch=iparticle->GetLastDaughter();        
301
302                         if ((ifch !=0 && ilch !=0) || fForceDecay == nodecay) {
303                           Int_t trackit=0;
304                           if (fForceDecay == nodecay) trackit = 1;
305                             gAlice->SetTrack(trackit,ntP,kf,
306                                              p,origin,polar,
307                                              0,kPPrimary,nt,fParentWeight);
308                             gAlice->KeepTrack(nt);
309                             Int_t iparent = nt;
310 //
311 // Children         
312                             if (fForceDecay != nodecay) {
313                               for (j=ifch; j<=ilch; j++)
314                                 {
315                                   TParticle *  ichild = 
316                                     (TParticle *) particles->At(j-1);
317                                   kf = CheckPDGCode(ichild->GetPdgCode());
318 //
319 // 
320                                   if (ChildSelected(TMath::Abs(kf))) {
321                                     origin[0]=origin0[0]+ichild->Vx()/10.;
322                                     origin[1]=origin0[1]+ichild->Vy()/10.;
323                                     origin[2]=origin0[2]+ichild->Vz()/10.;              
324                                     p[0]=ichild->Px();
325                                     p[1]=ichild->Py();
326                                     p[2]=ichild->Pz();
327                                     Float_t tof=kconv*ichild->T();
328                                     gAlice->SetTrack(fTrackIt, iparent, kf,
329                                                      p,origin,polar,
330                                                      tof,kPDecay,nt,fChildWeight);
331                                     gAlice->KeepTrack(nt);
332                                   } // select child
333                                 } // child loop
334                             } 
335                         }
336                     } // kinematic selection
337                 } // select particle
338             } // particle loop
339         } else {
340             for (Int_t i = 0; i<np-1; i++) {
341                 TParticle *  iparticle = (TParticle *) particles->At(i);
342                 kf = CheckPDGCode(iparticle->GetPdgCode());
343                 Int_t ks = iparticle->GetStatusCode();
344
345                 if (ks==1 && kf!=0 && KinematicSelection(iparticle)) {
346                         nc++;
347 //
348 // store track information
349                         p[0]=iparticle->Px();
350                         p[1]=iparticle->Py();
351                         p[2]=iparticle->Pz();
352                         origin[0]=origin0[0]+iparticle->Vx()/10.;
353                         origin[1]=origin0[1]+iparticle->Vy()/10.;
354                         origin[2]=origin0[2]+iparticle->Vz()/10.;
355                         Float_t tof=kconv*iparticle->T();
356                         gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
357                                          tof,kPPrimary,nt);
358                         gAlice->KeepTrack(nt);
359                 } // select particle
360             } // particle loop 
361             printf("\n I've put %i particles on the stack \n",nc);
362         } // mb ?
363         if (nc > 0) {
364             jev+=nc;
365             if (jev >= fNpart || fNpart == -1) {
366                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
367                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
368 // Print x-section summary
369                 fPythia->Pystat(1);
370                 break;
371             }
372         }
373     } // event loop
374 //  adjust weight due to kinematic selection
375     AdjustWeights();
376 //  get cross-section
377     fXsection=fPythia->GetPARI(1);
378 }
379
380 Bool_t AliGenPythia::ParentSelected(Int_t ip)
381 {
382 // True if particle is in list of parent particles to be selected
383     for (Int_t i=0; i<5; i++)
384     {
385         if (fParentSelect[i]==ip) return kTRUE;
386     }
387     return kFALSE;
388 }
389
390 Bool_t AliGenPythia::ChildSelected(Int_t ip)
391 {
392 // True if particle is in list of decay products to be selected
393     if (fForceDecay == all) return kTRUE;
394     
395     for (Int_t i=0; i<5; i++)
396     {
397         if (fChildSelect[i]==ip) return kTRUE;
398     }
399     return kFALSE;
400 }
401
402 Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
403 {
404 // Perform kinematic selection
405     Float_t px=particle->Px();
406     Float_t py=particle->Py();
407     Float_t pz=particle->Pz();
408     Float_t  e=particle->Energy();
409
410 //
411 //  transverse momentum cut    
412     Float_t pt=TMath::Sqrt(px*px+py*py);
413     if (pt > fPtMax || pt < fPtMin) 
414     {
415 //      printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
416         return kFALSE;
417     }
418 //
419 // momentum cut
420     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
421     if (p > fPMax || p < fPMin) 
422     {
423 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
424         return kFALSE;
425     }
426     
427 //
428 // theta cut
429     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
430     if (theta > fThetaMax || theta < fThetaMin) 
431     {
432 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
433         return kFALSE;
434     }
435
436 //
437 // rapidity cut
438     if (e==pz) {
439       return kFALSE;
440     }
441     else {
442       Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
443       if (y > fYMax || y < fYMin)
444         {
445 //      printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
446           return kFALSE;
447         }
448     }
449
450 //
451 // phi cut
452     Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
453     if (phi > fPhiMax || phi < fPhiMin)
454     {
455 //      printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
456         return kFALSE;
457     }
458
459     return kTRUE;
460 }
461 void AliGenPythia::AdjustWeights()
462 {
463 // Adjust the weights after generation of all events
464 //
465     TClonesArray *partArray = gAlice->Particles();
466     TParticle *part;
467     Int_t ntrack=gAlice->GetNtrack();
468     for (Int_t i=0; i<ntrack; i++) {
469         part= (TParticle*) partArray->UncheckedAt(i);
470         part->SetWeight(part->GetWeight()*fKineBias);
471     }
472 }
473
474 Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
475 {
476 //
477 //  If the particle is in a diffractive state, then take action accordingly
478   switch (pdgcode) {
479   case 91:
480     return 92;
481   case 110:
482     //rho_diff0 -- difficult to translate, return rho0
483     return 113;
484   case 210:
485     //pi_diffr+ -- change to pi+
486     return 211;
487   case 220:
488     //omega_di0 -- change to omega0
489     return 223;
490   case 330:
491     //phi_diff0 -- return phi0
492     return 333;
493   case 440:
494     //J/psi_di0 -- return J/psi
495     return 443;
496   case 2110:
497     //n_diffr -- return neutron
498     return 2112;
499   case 2210:
500     //p_diffr+ -- return proton
501     return 2212;
502   }
503   //non diffractive state -- return code unchanged
504   return pdgcode;
505 }
506
507     
508 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
509 {
510 // Treat protons as inside nuclei with mass numbers a1 and a2  
511     fNucA1 = a1;
512     fNucA2 = a2;
513 }
514         
515           
516 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
517 {
518 // Assignment operator
519     return *this;
520 }
521
522
523
524 void AliGenPythia::Streamer(TBuffer &R__b)
525 {
526    // Stream an object of class AliGenPythia.
527
528    if (R__b.IsReading()) {
529       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
530       AliGenerator::Streamer(R__b);
531       R__b >> (Int_t&)fProcess;
532       R__b >> (Int_t&)fStrucFunc;
533       R__b >> (Int_t&)fForceDecay;
534       R__b >> fEnergyCMS;
535       R__b >> fKineBias;
536       R__b >> fTrials;
537       fParentSelect.Streamer(R__b);
538       fChildSelect.Streamer(R__b);
539       R__b >> fXsection;
540 //      (AliPythia::Instance())->Streamer(R__b);
541       R__b >> fPtHardMin;
542       R__b >> fPtHardMax;
543 //      if (fDecayer) fDecayer->Streamer(R__b);
544    } else {
545       R__b.WriteVersion(AliGenPythia::IsA());
546       AliGenerator::Streamer(R__b);
547       R__b << (Int_t)fProcess;
548       R__b << (Int_t)fStrucFunc;
549       R__b << (Int_t)fForceDecay;
550       R__b << fEnergyCMS;
551       R__b << fKineBias;
552       R__b << fTrials;
553       fParentSelect.Streamer(R__b);
554       fChildSelect.Streamer(R__b);
555       R__b << fXsection;
556 //      R__b << fPythia;
557       R__b << fPtHardMin;
558       R__b << fPtHardMax;
559       //     fDecayer->Streamer(R__b);
560    }
561 }
562
563
564
565 #ifndef WIN32
566 #define pyr    pyr_
567 #define pyrset pyrset_
568 #define pyrget pyrget_
569 #else
570 #define pyr    PYR
571 #define pyrset PYRSET
572 #define pyrget PYRGET
573 #endif
574
575 extern "C" {
576   Double_t pyr(Int_t*) {return sRandom->Rndm();}
577   void pyrset(Int_t*,Int_t*) {}
578   void pyrget(Int_t*,Int_t*) {}
579 }