This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
1 #include "AliGenParam.h"
2 #include "AliGenMUONlib.h"
3 #include "AliRun.h"
4 #include "TGeant3.h"
5 #include "AliPythia.h"
6 #include <TDirectory.h>
7 #include <TFile.h>
8 #include <TTree.h>
9 #include <stdlib.h>
10 #include <TMCParticle.h>
11
12 ClassImp(AliGenParam)
13
14 //------------------------------------------------------------
15
16   //Begin_Html
17   /*
18     <img src="gif/AliGenParam.gif">
19   */
20   //End_Html
21
22 //____________________________________________________________
23 //____________________________________________________________
24 AliGenParam::AliGenParam()
25                  :AliGenerator()
26 {
27   fPtPara = 0;
28   fYPara  = 0;
29   fIpart  = 0;
30   fAnalog = 1;
31 }
32
33 //____________________________________________________________
34 AliGenParam::AliGenParam(Int_t npart, Int_t Ipart) 
35 //                                 Double_t (*PtPara)(Double_t*, Double_t*), 
36 //                                 Double_t (*YPara) (Double_t* ,Double_t*))
37     :AliGenerator(npart)
38 {
39   //
40   //  fName="HMESONpara";
41   //  fTitle="Heavy Mesons Parametrisation";
42   fPtParaFunc = AliGenMUONlib::GetPt(Ipart);
43   fYParaFunc  = AliGenMUONlib::GetY(Ipart);
44   fIpParaFunc = AliGenMUONlib::GetIp(Ipart);
45   
46   fPtPara = 0;
47   fYPara  = 0;
48   fIpart  = Ipart;
49   fAnalog = 1;
50   fChildSelect.Set(5);
51   for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
52   ForceDecay();
53 }
54
55 //____________________________________________________________
56 AliGenParam::~AliGenParam()
57 {
58     delete  fPtPara;
59     delete  fYPara;
60 }
61
62 //____________________________________________________________
63 void AliGenParam::Init()
64 {
65     SetMC(new AliPythia());
66     fPythia= (AliPythia*) fgMCEvGen;
67
68 //  End of the test !!!
69   //Begin_Html
70   /*
71     <img src="gif/AliGenParam.gif">
72   */
73   //End_Html
74  
75   fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
76   fYPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
77 //
78 // dN/dy| y=0
79   Double_t y1=0;
80   Double_t y2=0;
81   
82   fdNdy0=fYParaFunc(&y1,&y2);
83 //
84 // Integral over generation region
85   Float_t IntYS  = fYPara ->Integral(fYMin, fYMax);
86   Float_t IntPt0 = fPtPara->Integral(0,15);
87   Float_t IntPtS = fPtPara->Integral(fPtMin,fPtMax);
88   Float_t PhiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
89   if (fAnalog) {
90      fYWgt  = IntYS/fdNdy0;
91      fPtWgt = IntPtS/IntPt0;
92      fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
93   } else {
94       fYWgt = IntYS/fdNdy0;
95       fPtWgt = (fPtMax-fPtMin)/IntPt0;
96       fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
97   }
98 //
99 // particle decay related initialization
100   fPythia->DefineParticles();
101 // semimuonic decays of charm and beauty
102   fPythia->ForceDecay(fForceDecay);
103 //
104     switch (fForceDecay) 
105     {
106     case semielectronic:
107     case dielectron:
108     case b_jpsi_dielectron:
109     case b_psip_dielectron:
110         fChildSelect[0]=11;     
111         break;
112     case semimuonic:
113     case dimuon:
114     case b_jpsi_dimuon:
115     case b_psip_dimuon:
116         fChildSelect[0]=13;
117         break;
118     }
119
120 }
121
122 //____________________________________________________________
123 void AliGenParam::Generate()
124 {
125 // Generate 'npart' of heavy mesons (J/Psi, upsilon or phi) in the
126 // the desired theta, phi and momentum windows; Gaussian smearing 
127 // on the vertex is done if selected
128
129   AliMC* pMC = AliMC::GetMC();
130
131   Float_t polar[3]= {0,0,0};
132   //
133   Float_t origin[3], origin0[3];
134   Float_t pt, pl, ptot;
135   Float_t phi, theta;
136   Float_t p[3];
137   Float_t ty, xmt;
138   Int_t i, nt, j;
139   Float_t  wgtp, wgtch;
140   Double_t dummy;
141   //
142   //
143   Float_t random[6];
144   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
145   if(fVertexSmear==perEvent) {
146       pMC->Rndm(random,6);
147       for (j=0;j<3;j++) {
148           origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
149               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
150       }
151   }
152   for(i=0;i<fNpart;i++) {
153       while(1) {
154 //
155 // particle type
156           Int_t Ipart = fIpParaFunc();
157           fChildWeight=(fPythia->GetBraPart(Ipart))*fParentWeight;        
158           Float_t am=fPythia->GetPMAS(fPythia->LuComp(Ipart),1);
159           pMC->Rndm(random,2);
160 //
161 // phi
162           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
163 //
164 // y
165           ty=Float_t(TMath::TanH(fYPara->GetRandom()));
166 //
167 // pT
168           if (fAnalog) {
169               pt=fPtPara->GetRandom();
170               wgtp=fParentWeight;
171               wgtch=fChildWeight;
172           } else {
173               pt=fPtMin+random[1]*(fPtMax-fPtMin);
174               Double_t ptd=pt;
175               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
176               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
177           }
178           xmt=sqrt(pt*pt+am*am);
179           pl=xmt*ty/sqrt(1.-ty*ty);
180           theta=TMath::ATan2(pt,pl);
181           if(theta<fThetaMin || theta>fThetaMax) continue;
182           ptot=TMath::Sqrt(pt*pt+pl*pl);
183           if(ptot<fPMin || ptot>fPMax) continue;
184           p[0]=pt*TMath::Cos(phi);
185           p[1]=pt*TMath::Sin(phi);
186           p[2]=pl;
187           if(fVertexSmear==perTrack) {
188               pMC->Rndm(random,6);
189               for (j=0;j<3;j++) {
190                   origin0[j]=
191                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
192                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
193               }
194           }
195           
196           Int_t kg=fPythia->GetGeantCode(Ipart);
197 //
198 // parent
199
200           gAlice->
201               SetTrack(0,-1,kg,p,origin,polar,0,"Primary",nt,wgtp);
202           Int_t iparent=nt;
203 //
204 // use lujet to decay particle
205
206           Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
207           fPythia->DecayParticle(Ipart,energy,theta,phi);
208 //
209 // select muons
210           TObjArray* particles = fPythia->GetPrimaries() ;
211           Int_t np = particles->GetEntriesFast();
212           for (Int_t i = 0; i<np; i++) {
213               TMCParticle *  iparticle = (TMCParticle *) particles->At(i);
214               Int_t kf = iparticle->GetKF();
215 //
216 // children
217               if (ChildSelected(TMath::Abs(kf)))
218               {
219                   p[0]=iparticle->GetPx();
220                   p[1]=iparticle->GetPy();
221                   p[2]=iparticle->GetPz();
222                   origin[0]=origin0[0]+iparticle->GetVx()/10;
223                   origin[1]=origin0[1]+iparticle->GetVy()/10;
224                   origin[2]=origin0[2]+iparticle->GetVz()/10;
225                   gAlice->SetTrack(1,iparent,fPythia->GetGeantCode(kf),
226                                    p,origin,polar,
227                                    0,"Decay",nt,wgtch);
228                   gAlice->KeepTrack(nt);
229               } // select muon
230           } // decay particle loop
231           break;
232       } // kinematic selection
233   } // event loop
234 }
235
236 Bool_t AliGenParam::ChildSelected(Int_t ip)
237 {
238     for (Int_t i=0; i<5; i++)
239     {
240         if (fChildSelect[i]==ip) return kTRUE;
241     }
242     return kFALSE;
243 }
244
245 Bool_t AliGenParam::KinematicSelection(TMCParticle *particle)
246 {
247     Float_t px=particle->GetPx();
248     Float_t py=particle->GetPy();
249     Float_t pz=particle->GetPz();
250 //
251 // momentum cut
252     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
253     if (p > fPMax || p < fPMin) 
254     {
255 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
256         return kFALSE;
257     }
258     Float_t pt=TMath::Sqrt(px*px+py*py);
259     
260 //
261 // theta cut
262     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
263     if (theta > fThetaMax || theta < fThetaMin) 
264     {
265 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
266         return kFALSE;
267     }
268
269     return kTRUE;
270 }
271
272
273
274