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