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