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