New version from G.Martinez & A.Morsch
[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.11  1999/09/29 09:24:14  fca
19 Introduction of the Copyright and cvs Log
20 */
21
22 #include "AliGenerator.h"
23 #include "AliGenPythia.h"
24 #include "AliRun.h"
25 #include "AliPythia.h"
26 #include <TDirectory.h>
27 #include <TFile.h>
28 #include <TTree.h>
29 #include <stdlib.h>
30 #include <AliPythia.h>
31 #include <TParticle.h>
32 //#include <GParticle.h>
33  ClassImp(AliGenPythia)
34
35 AliGenPythia::AliGenPythia()
36                  :AliGenerator()
37 {
38 }
39
40 AliGenPythia::AliGenPythia(Int_t npart)
41                  :AliGenerator(npart)
42 {
43 // default charm production at 5. 5 TeV
44 // semimuonic decay
45 // structure function GRVHO
46 //
47     fXsection  = 0.;
48     fParentSelect.Set(5);
49     fChildSelect.Set(5);
50     for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0;
51     SetProcess();
52     SetStrucFunc();
53     SetForceDecay();
54     SetPtHard();
55     SetEnergyCMS();
56 }
57
58 AliGenPythia::~AliGenPythia()
59 {
60 }
61
62 void AliGenPythia::Init()
63 {
64     SetMC(new AliPythia());
65     fPythia=(AliPythia*) fgMCEvGen;
66 //
67     fParentWeight=1./Float_t(fNpart);
68 //
69 //  Forward Paramters to the AliPythia object
70     fPythia->DefineParticles();
71     fPythia->SetCKIN(3,fPtHardMin);
72     fPythia->SetCKIN(4,fPtHardMax);    
73     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
74     fPythia->ForceDecay(fForceDecay);
75     fPythia->LuList(0);
76     fPythia->PyStat(2);
77 //  Parent and Children Selection
78     switch (fProcess) 
79     {
80     case charm:
81
82         fParentSelect[0]=411;
83         fParentSelect[1]=421;
84         fParentSelect[2]=431;
85         fParentSelect[3]=4122;  
86         break;
87     case charm_unforced:
88
89         fParentSelect[0]=411;
90         fParentSelect[1]=421;
91         fParentSelect[2]=431;
92         fParentSelect[3]=4122;  
93         break;
94     case beauty:
95         fParentSelect[0]=511;
96         fParentSelect[1]=521;
97         fParentSelect[2]=531;
98         fParentSelect[3]=5122;  
99         break;
100     case beauty_unforced:
101         fParentSelect[0]=511;
102         fParentSelect[1]=521;
103         fParentSelect[2]=531;
104         fParentSelect[3]=5122;  
105         break;
106     case jpsi_chi:
107     case jpsi:
108         fParentSelect[0]=443;
109         break;
110     case mb:
111         break;
112     }
113
114     switch (fForceDecay) 
115     {
116     case semielectronic:
117     case dielectron:
118     case b_jpsi_dielectron:
119     case b_psip_dielectron:
120         fChildSelect[0]=11;     
121         break;
122     case semimuonic:
123     case dimuon:
124     case b_jpsi_dimuon:
125     case b_psip_dimuon:
126     case pitomu:
127     case katomu:
128         fChildSelect[0]=13;
129         break;
130     case all:
131     case nodecay:
132         break;
133     }
134 }
135
136 void AliGenPythia::Generate()
137 {
138
139     Float_t polar[3] =   {0,0,0};
140     Float_t origin[3]=   {0,0,0};
141     Float_t origin_p[3]= {0,0,0};
142     Float_t origin0[3]=  {0,0,0};
143     Float_t p[3], p_p[4], random[6];
144     static TClonesArray *particles;
145 //  converts from mm/c to s
146     const Float_t kconv=0.001/2.999792458e8;
147     
148     
149 //
150     Int_t nt=0;
151     Int_t nt_p=0;
152     Int_t jev=0;
153     Int_t j, kf;
154
155     if(!particles) particles=new TClonesArray("TParticle",1000);
156     
157     fTrials=0;
158     for (j=0;j<3;j++) origin0[j]=fOrigin[j];
159     if(fVertexSmear==perEvent) {
160         gMC->Rndm(random,6);
161         for (j=0;j<3;j++) {
162             origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
163                 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
164             fPythia->SetMSTP(151,0);
165         }
166     } else if (fVertexSmear==perTrack) {
167         fPythia->SetMSTP(151,0);
168         for (j=0;j<3;j++) {
169             fPythia->SetPARP(151+j, fOsigma[j]*10.);
170         }
171     }
172     while(1)
173     {
174         fPythia->PyEvnt();
175         fTrials++;
176         fPythia->ImportParticles(particles,"All");
177         Int_t np = particles->GetEntriesFast();
178         printf("\n **************************************************%d\n",np);
179         Int_t nc=0;
180         if (np == 0 ) continue;
181         if (fProcess != mb) {
182             for (Int_t i = 0; i<np; i++) {
183                 TParticle *  iparticle = (TParticle *) particles->At(i);
184                 kf = CheckPDGCode(iparticle->GetPdgCode());
185                 fChildWeight=(fPythia->GetBraPart(kf))*fParentWeight;     
186 //
187 // Parent
188                 if (ParentSelected(TMath::Abs(kf))) {
189                     if (KinematicSelection(iparticle)) {
190                         if (nc==0) {
191 //
192 // Store information concerning the hard scattering process
193 //
194                             Float_t mass_p = fPythia->GetPARI(13);
195                             Float_t   pt_p = fPythia->GetPARI(17);
196                             Float_t    y_p = fPythia->GetPARI(37);
197                             Float_t  xmt_p = sqrt(pt_p*pt_p+mass_p*mass_p);
198                             Float_t     ty = Float_t(TMath::TanH(y_p));
199                             p_p[0] = pt_p;
200                             p_p[1] = 0;
201                             p_p[2] = xmt_p*ty/sqrt(1.-ty*ty);
202                             p_p[3] = mass_p;
203                             gAlice->SetTrack(0,-1,-1,
204                                              p_p,origin_p,polar,
205                                              0,"Hard Scat.",nt_p,fParentWeight);
206                             gAlice->KeepTrack(nt_p);
207                         }
208                         nc++;
209 //
210 // store parent track information
211                         p[0]=iparticle->Px();
212                         p[1]=iparticle->Py();
213                         p[2]=iparticle->Pz();
214                         origin[0]=origin0[0]+iparticle->Vx()/10;
215                         origin[1]=origin0[1]+iparticle->Vy()/10;
216                         origin[2]=origin0[2]+iparticle->Vz()/10;
217
218                         Int_t ifch=iparticle->GetFirstDaughter();
219                         Int_t ilch=iparticle->GetLastDaughter();        
220                         if (ifch !=0 && ilch !=0) {
221                             gAlice->SetTrack(0,nt_p,kf,
222                                              p,origin,polar,
223                                              0,"Primary",nt,fParentWeight);
224                             gAlice->KeepTrack(nt);
225                             Int_t iparent = nt;
226 //
227 // Children         
228
229                             for (j=ifch; j<=ilch; j++)
230                             {
231                                 TParticle *  ichild = 
232                                     (TParticle *) particles->At(j-1);
233                                 kf = CheckPDGCode(ichild->GetPdgCode());
234 //
235 // 
236                                 if (ChildSelected(TMath::Abs(kf))) {
237                                     origin[0]=ichild->Vx();
238                                     origin[1]=ichild->Vy();
239                                     origin[2]=ichild->Vz();             
240                                     p[0]=ichild->Px();
241                                     p[1]=ichild->Py();
242                                     p[2]=ichild->Pz();
243                                     Float_t tof=kconv*ichild->T();
244                                     gAlice->SetTrack(fTrackIt, iparent, kf,
245                                                      p,origin,polar,
246                                                      tof,"Decay",nt,fChildWeight);
247                                     gAlice->KeepTrack(nt);
248                                 } // select child
249                             } // child loop
250                         }
251                     } // kinematic selection
252                 } // select particle
253             } // particle loop
254         } else {
255             for (Int_t i = 0; i<np; i++) {
256                 TParticle *  iparticle = (TParticle *) particles->At(i);
257                 kf = CheckPDGCode(iparticle->GetPdgCode());
258                 Int_t ks = iparticle->GetStatusCode();
259                 if (ks==1 && kf!=0 && KinematicSelection(iparticle)) {
260                         nc++;
261 //
262 // store track information
263                         p[0]=iparticle->Px();
264                         p[1]=iparticle->Py();
265                         p[2]=iparticle->Pz();
266                         origin[0]=origin0[0]+iparticle->Vx()/10;
267                         origin[1]=origin0[1]+iparticle->Vy()/10;
268                         origin[2]=origin0[2]+iparticle->Vz()/10;
269                         Float_t tof=kconv*iparticle->T();
270                         gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar,
271                                          tof,"Primary",nt);
272                         gAlice->KeepTrack(nt);
273                 } // select particle
274             } // particle loop 
275             printf("\n I've put %i particles on the stack \n",nc);
276         } // mb ?
277         if (nc > 0) {
278             jev+=nc;
279             if (jev >= fNpart) {
280                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
281                 printf("\n Trials: %i\n",fTrials);
282                 break;
283             }
284         }
285     } // event loop
286 //  adjust weight due to kinematic selection
287     AdjustWeights();
288 //  get cross-section
289     fXsection=fPythia->GetPARI(1);
290 }
291
292 Bool_t AliGenPythia::ParentSelected(Int_t ip)
293 {
294     for (Int_t i=0; i<5; i++)
295     {
296         if (fParentSelect[i]==ip) return kTRUE;
297     }
298     return kFALSE;
299 }
300
301 Bool_t AliGenPythia::ChildSelected(Int_t ip)
302 {
303     for (Int_t i=0; i<5; i++)
304     {
305         if (fChildSelect[i]==ip) return kTRUE;
306     }
307     return kFALSE;
308 }
309
310 Bool_t AliGenPythia::KinematicSelection(TParticle *particle)
311 {
312     Float_t px=particle->Px();
313     Float_t py=particle->Py();
314     Float_t pz=particle->Pz();
315     Float_t  e=particle->Energy();
316
317 //
318 //  transverse momentum cut    
319     Float_t pt=TMath::Sqrt(px*px+py*py);
320     if (pt > fPtMax || pt < fPtMin) 
321     {
322 //      printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
323         return kFALSE;
324     }
325 //
326 // momentum cut
327     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
328     if (p > fPMax || p < fPMin) 
329     {
330 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
331         return kFALSE;
332     }
333     
334 //
335 // theta cut
336     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
337     if (theta > fThetaMax || theta < fThetaMin) 
338     {
339 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
340         return kFALSE;
341     }
342
343 //
344 // rapidity cut
345     Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
346     if (y > fYMax || y < fYMin)
347     {
348 //      printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
349         return kFALSE;
350     }
351
352 //
353 // phi cut
354     Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)))+TMath::Pi();
355     if (phi > fPhiMax || phi < fPhiMin)
356     {
357 //      printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
358         return kFALSE;
359     }
360
361     return kTRUE;
362 }
363 void AliGenPythia::AdjustWeights()
364 {
365     TClonesArray *PartArray = gAlice->Particles();
366     TParticle *Part;
367     Int_t ntrack=gAlice->GetNtrack();
368     for (Int_t i=0; i<ntrack; i++) {
369         Part= (TParticle*) PartArray->UncheckedAt(i);
370         Part->SetWeight(Part->GetWeight()*fKineBias);
371     }
372 }
373
374 Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode)
375 {
376 //
377 //  If the particle is in a diffractive state, then take action accordigly
378   switch (pdgcode) {
379   case 110:
380     //rho_diff0 -- difficult to translate, return rho0
381     return 113;
382   case 210:
383     //pi_diffr+ -- change to pi+
384     return 211;
385   case 220:
386     //omega_di0 -- change to omega0
387     return 223;
388   case 330:
389     //phi_diff0 -- return phi0
390     return 333;
391   case 440:
392     //J/psi_di0 -- return J/psi
393     return 443;
394   case 2110:
395     //n_diffr -- return neutron
396     return 2112;
397   case 2210:
398     //p_diffr+ -- return proton
399     return 2212;
400   }
401   //non diffractive state -- return code unchanged
402   return pdgcode;
403 }
404                   
405
406
407
408
409
410
411
412