cafe5933b82b73c1912161fb5f82e605c3b7ea59
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
1 #include "AliGenParam.h"
2 #include "AliGenMUONlib.h"
3 #include "AliRun.h"
4 #include "AliPythia.h"
5 #include <TDirectory.h>
6 #include <TFile.h>
7 #include <TTree.h>
8 #include <stdlib.h>
9 #include <TParticle.h>
10
11 ClassImp(AliGenParam)
12
13 //------------------------------------------------------------
14
15   //Begin_Html
16   /*
17     <img src="picts/AliGenParam.gif">
18   */
19   //End_Html
20
21 //____________________________________________________________
22 //____________________________________________________________
23 AliGenParam::AliGenParam()
24                  :AliGenerator()
25 {
26   fPtPara = 0;
27   fYPara  = 0;
28   fParam  = jpsi_p;
29   fAnalog = analog;
30 }
31
32 //____________________________________________________________
33 AliGenParam::AliGenParam(Int_t npart, Param_t param) 
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";
41   fPtParaFunc = AliGenMUONlib::GetPt(param);
42   fYParaFunc  = AliGenMUONlib::GetY(param);
43   fIpParaFunc = AliGenMUONlib::GetIp(param);
44   
45   fPtPara = 0;
46   fYPara  = 0;
47   fParam  = param;
48   fAnalog = analog;
49   fChildSelect.Set(5);
50   for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
51   ForceDecay();
52 }
53
54 //____________________________________________________________
55 AliGenParam::~AliGenParam()
56 {
57     delete  fPtPara;
58     delete  fYPara;
59 }
60
61 //____________________________________________________________
62 void AliGenParam::Init()
63 {
64     SetMC(new AliPythia());
65     fPythia= (AliPythia*) fgMCEvGen;
66
67 //  End of the test !!!
68   //Begin_Html
69   /*
70     <img src="picts/AliGenParam.gif">
71   */
72   //End_Html
73  
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
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
87   Float_t IntYS  = YPara ->Integral(fYMin, fYMax);
88   Float_t IntPt0 = PtPara->Integral(0,15);
89   Float_t IntPtS = PtPara->Integral(fPtMin,fPtMax);
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;
120     case pitomu:
121         fChildSelect[0]=13;
122         break;
123     case katomu:
124         fChildSelect[0]=13;
125         break;
126     }
127
128 }
129
130 //____________________________________________________________
131 void 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
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;
148   static TClonesArray *particles;
149   //
150   if(!particles) particles=new TClonesArray("TParticle",1000);
151   //
152   Float_t random[6];
153   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
154   if(fVertexSmear==perEvent) {
155       gMC->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           gMC->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               gMC->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 //
206 // parent
207
208           gAlice->
209               SetTrack(0,-1,Ipart,p,origin,polar,0,"Primary",nt,wgtp);
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);
216 //        fPythia->LuList(1);
217           
218 //
219 // select muons
220           fPythia->ImportParticles(particles) ;
221           Int_t np = particles->GetEntriesFast();
222           for (Int_t i = 0; i<np; i++) {
223               TParticle *  iparticle = (TParticle *) particles->At(i);
224               Int_t kf = iparticle->GetPdgCode();
225 //
226 // children
227               if (ChildSelected(TMath::Abs(kf)))
228               {
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,
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
246 Bool_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
255 Bool_t AliGenParam::KinematicSelection(TParticle *particle)
256 {
257     Float_t px=particle->Px();
258     Float_t py=particle->Py();
259     Float_t pz=particle->Pz();
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