Change the name of the default its files from _60 to _5
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
CommitLineData
fe4da5cc 1#include "AliGenParam.h"
2#include "AliGenMUONlib.h"
3#include "AliRun.h"
fe4da5cc 4#include "AliPythia.h"
5#include <TDirectory.h>
6#include <TFile.h>
7#include <TTree.h>
8#include <stdlib.h>
1578254f 9#include <TParticle.h>
fe4da5cc 10
11ClassImp(AliGenParam)
12
13//------------------------------------------------------------
14
15 //Begin_Html
16 /*
1439f98e 17 <img src="picts/AliGenParam.gif">
fe4da5cc 18 */
19 //End_Html
20
21//____________________________________________________________
22//____________________________________________________________
23AliGenParam::AliGenParam()
24 :AliGenerator()
25{
26 fPtPara = 0;
27 fYPara = 0;
b7601ac4 28 fParam = jpsi_p;
29 fAnalog = analog;
fe4da5cc 30}
31
32//____________________________________________________________
b7601ac4 33AliGenParam::AliGenParam(Int_t npart, Param_t param)
fe4da5cc 34// Double_t (*PtPara)(Double_t*, Double_t*),
35// Double_t (*YPara) (Double_t* ,Double_t*))
36 :AliGenerator(npart)
37{
38 //
39 // fName="HMESONpara";
40 // fTitle="Heavy Mesons Parametrisation";
b7601ac4 41 fPtParaFunc = AliGenMUONlib::GetPt(param);
42 fYParaFunc = AliGenMUONlib::GetY(param);
43 fIpParaFunc = AliGenMUONlib::GetIp(param);
fe4da5cc 44
45 fPtPara = 0;
46 fYPara = 0;
b7601ac4 47 fParam = param;
48 fAnalog = analog;
fe4da5cc 49 fChildSelect.Set(5);
50 for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
51 ForceDecay();
52}
53
54//____________________________________________________________
55AliGenParam::~AliGenParam()
56{
57 delete fPtPara;
58 delete fYPara;
59}
60
61//____________________________________________________________
62void AliGenParam::Init()
63{
64 SetMC(new AliPythia());
65 fPythia= (AliPythia*) fgMCEvGen;
66
67// End of the test !!!
68 //Begin_Html
69 /*
1439f98e 70 <img src="picts/AliGenParam.gif">
fe4da5cc 71 */
72 //End_Html
73
b7601ac4 74 fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
75 fYPara = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
76 TF1* PtPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
77 TF1* YPara = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
78
fe4da5cc 79//
80// dN/dy| y=0
81 Double_t y1=0;
82 Double_t y2=0;
83
84 fdNdy0=fYParaFunc(&y1,&y2);
85//
86// Integral over generation region
b7601ac4 87 Float_t IntYS = YPara ->Integral(fYMin, fYMax);
88 Float_t IntPt0 = PtPara->Integral(0,15);
89 Float_t IntPtS = PtPara->Integral(fPtMin,fPtMax);
fe4da5cc 90 Float_t PhiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
91 if (fAnalog) {
92 fYWgt = IntYS/fdNdy0;
93 fPtWgt = IntPtS/IntPt0;
94 fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
95 } else {
96 fYWgt = IntYS/fdNdy0;
97 fPtWgt = (fPtMax-fPtMin)/IntPt0;
98 fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
99 }
100//
101// particle decay related initialization
102 fPythia->DefineParticles();
103// semimuonic decays of charm and beauty
104 fPythia->ForceDecay(fForceDecay);
105//
106 switch (fForceDecay)
107 {
108 case semielectronic:
109 case dielectron:
110 case b_jpsi_dielectron:
111 case b_psip_dielectron:
112 fChildSelect[0]=11;
113 break;
114 case semimuonic:
115 case dimuon:
116 case b_jpsi_dimuon:
117 case b_psip_dimuon:
118 fChildSelect[0]=13;
119 break;
b7601ac4 120 case pitomu:
121 fChildSelect[0]=13;
122 break;
123 case katomu:
124 fChildSelect[0]=13;
125 break;
fe4da5cc 126 }
127
128}
129
130//____________________________________________________________
131void AliGenParam::Generate()
132{
133// Generate 'npart' of heavy mesons (J/Psi, upsilon or phi) in the
134// the desired theta, phi and momentum windows; Gaussian smearing
135// on the vertex is done if selected
136
fe4da5cc 137
138 Float_t polar[3]= {0,0,0};
139 //
140 Float_t origin[3], origin0[3];
141 Float_t pt, pl, ptot;
142 Float_t phi, theta;
143 Float_t p[3];
144 Float_t ty, xmt;
145 Int_t i, nt, j;
146 Float_t wgtp, wgtch;
147 Double_t dummy;
09fd3ea2 148 static TClonesArray *particles;
fe4da5cc 149 //
09fd3ea2 150 if(!particles) particles=new TClonesArray("TParticle",1000);
fe4da5cc 151 //
152 Float_t random[6];
153 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
154 if(fVertexSmear==perEvent) {
cfce8870 155 gMC->Rndm(random,6);
fe4da5cc 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);
cfce8870 168 gMC->Rndm(random,2);
fe4da5cc 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) {
cfce8870 197 gMC->Rndm(random,6);
fe4da5cc 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
fe4da5cc 205//
206// parent
207
208 gAlice->
1578254f 209 SetTrack(0,-1,Ipart,p,origin,polar,0,"Primary",nt,wgtp);
fe4da5cc 210 Int_t iparent=nt;
211//
212// use lujet to decay particle
213
214 Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
215 fPythia->DecayParticle(Ipart,energy,theta,phi);
b7601ac4 216// fPythia->LuList(1);
217
fe4da5cc 218//
219// select muons
09fd3ea2 220 fPythia->ImportParticles(particles) ;
fe4da5cc 221 Int_t np = particles->GetEntriesFast();
222 for (Int_t i = 0; i<np; i++) {
1578254f 223 TParticle * iparticle = (TParticle *) particles->At(i);
224 Int_t kf = iparticle->GetPdgCode();
fe4da5cc 225//
226// children
227 if (ChildSelected(TMath::Abs(kf)))
228 {
1578254f 229 p[0]=iparticle->Px();
230 p[1]=iparticle->Py();
231 p[2]=iparticle->Pz();
232 origin[0]=origin0[0]+iparticle->Vx()/10;
233 origin[1]=origin0[1]+iparticle->Vy()/10;
234 origin[2]=origin0[2]+iparticle->Vz()/10;
235 gAlice->SetTrack(fTrackIt,iparent,kf,
fe4da5cc 236 p,origin,polar,
237 0,"Decay",nt,wgtch);
238 gAlice->KeepTrack(nt);
239 } // select muon
240 } // decay particle loop
241 break;
242 } // kinematic selection
243 } // event loop
244}
245
246Bool_t AliGenParam::ChildSelected(Int_t ip)
247{
248 for (Int_t i=0; i<5; i++)
249 {
250 if (fChildSelect[i]==ip) return kTRUE;
251 }
252 return kFALSE;
253}
254
1578254f 255Bool_t AliGenParam::KinematicSelection(TParticle *particle)
fe4da5cc 256{
1578254f 257 Float_t px=particle->Px();
258 Float_t py=particle->Py();
259 Float_t pz=particle->Pz();
fe4da5cc 260//
261// momentum cut
262 Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
263 if (p > fPMax || p < fPMin)
264 {
265// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
266 return kFALSE;
267 }
268 Float_t pt=TMath::Sqrt(px*px+py*py);
269
270//
271// theta cut
272 Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
273 if (theta > fThetaMax || theta < fThetaMin)
274 {
275// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
276 return kFALSE;
277 }
278
279 return kTRUE;
280}
281
282
283
284