Some rationalisation of the documentation. In particular pictures are all now in...
[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="picts/AliGenParam.gif">
19   */
20   //End_Html
21
22 //____________________________________________________________
23 //____________________________________________________________
24 AliGenParam::AliGenParam()
25                  :AliGenerator()
26 {
27   fPtPara = 0;
28   fYPara  = 0;
29   fParam  = jpsi_p;
30   fAnalog = analog;
31 }
32
33 //____________________________________________________________
34 AliGenParam::AliGenParam(Int_t npart, Param_t param) 
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(param);
43   fYParaFunc  = AliGenMUONlib::GetY(param);
44   fIpParaFunc = AliGenMUONlib::GetIp(param);
45   
46   fPtPara = 0;
47   fYPara  = 0;
48   fParam  = param;
49   fAnalog = analog;
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="picts/AliGenParam.gif">
72   */
73   //End_Html
74  
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
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
88   Float_t IntYS  = YPara ->Integral(fYMin, fYMax);
89   Float_t IntPt0 = PtPara->Integral(0,15);
90   Float_t IntPtS = PtPara->Integral(fPtMin,fPtMax);
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;
121     case pitomu:
122         fChildSelect[0]=13;
123         break;
124     case katomu:
125         fChildSelect[0]=13;
126         break;
127     }
128
129 }
130
131 //____________________________________________________________
132 void 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);
217 //        fPythia->LuList(1);
218           
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;
236                   gAlice->SetTrack(fTrackIt,iparent,fPythia->GetGeantCode(kf),
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
247 Bool_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
256 Bool_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