This commit was generated by cvs2svn to compensate for changes in r2,
[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 /*
18 <img src="gif/AliGenParam.gif">
19 */
20 //End_Html
21
22//____________________________________________________________
23//____________________________________________________________
24AliGenParam::AliGenParam()
25 :AliGenerator()
26{
27 fPtPara = 0;
28 fYPara = 0;
29 fIpart = 0;
30 fAnalog = 1;
31}
32
33//____________________________________________________________
34AliGenParam::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//____________________________________________________________
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 /*
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//____________________________________________________________
123void 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
236Bool_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
245Bool_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