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