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