Changes for compatibility with version 2.23 of ROOT
[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.13  1999/11/04 11:30:31  fca
19 Correct the logics for SetForceDecay
20
21 Revision 1.12  1999/11/03 17:43:20  fca
22 New version from G.Martinez & A.Morsch
23
24 Revision 1.11  1999/09/29 09:24:14  fca
25 Introduction of the Copyright and cvs Log
26
27 */
28
29 #include "AliGenParam.h"
30 #include "AliGenMUONlib.h"
31 #include "AliGenPHOSlib.h"
32 #include "AliRun.h"
33 #include "AliPythia.h"
34 #include <TDirectory.h>
35 #include <TFile.h>
36 #include <TTree.h>
37 #include <stdlib.h>
38 #include <TParticle.h>
39
40 ClassImp(AliGenParam)
41
42 //------------------------------------------------------------
43
44   //Begin_Html
45   /*
46     <img src="picts/AliGenParam.gif">
47   */
48   //End_Html
49
50 //____________________________________________________________
51 //____________________________________________________________
52 AliGenParam::AliGenParam()
53                  :AliGenerator()
54 {
55   fPtPara = 0;
56   fYPara  = 0;
57   fParam  = jpsi_p;
58   fAnalog = analog;
59   SetCutOnChild();
60 }
61
62 //____________________________________________________________
63
64 AliGenParam::AliGenParam(Int_t npart, Param_t param) :AliGenerator(npart)
65 {
66   //
67   //  fName="HMESONpara";
68   //  fTitle="Heavy Mesons Parametrisation";
69   fPtParaFunc = AliGenMUONlib::GetPt(param);
70   fYParaFunc  = AliGenMUONlib::GetY(param);
71   fIpParaFunc = AliGenMUONlib::GetIp(param);
72   
73   fPtPara = 0;
74   fYPara  = 0;
75   fParam  = param;
76   fAnalog = analog;
77   fChildSelect.Set(5);
78   for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
79   SetForceDecay();
80   SetCutOnChild();
81 }
82
83 AliGenParam::AliGenParam(Int_t npart, Param_t param,
84                          Double_t (*PtPara) (Double_t*, Double_t*),
85                          Double_t (*YPara ) (Double_t* ,Double_t*),
86                          Int_t    (*IpPara) ())                  
87     :AliGenerator(npart)
88 {
89 // Gines Martinez 1/10/99 
90     fPtParaFunc = PtPara; 
91     fYParaFunc  = YPara;  
92     fIpParaFunc = IpPara;
93 //  
94     fPtPara = 0;
95     fYPara  = 0;
96     fParam  = param;
97     fAnalog = analog;
98     fChildSelect.Set(5);
99     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
100     SetForceDecay();
101     SetCutOnChild();
102 }
103
104 //____________________________________________________________
105 AliGenParam::~AliGenParam()
106 {
107     delete  fPtPara;
108     delete  fYPara;
109 }
110
111 //____________________________________________________________
112 void AliGenParam::Init()
113 {
114     SetMC(new AliPythia());
115     fPythia= (AliPythia*) fgMCEvGen;
116
117 //  End of the test !!!
118   //Begin_Html
119   /*
120     <img src="picts/AliGenParam.gif">
121   */
122   //End_Html
123  
124   fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
125   fYPara  = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
126   TF1* PtPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
127   TF1* YPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
128
129 //
130 // dN/dy| y=0
131   Double_t y1=0;
132   Double_t y2=0;
133   
134   fdNdy0=fYParaFunc(&y1,&y2);
135 //
136 // Integral over generation region
137   Float_t IntYS  = YPara ->Integral(fYMin, fYMax);
138   Float_t IntPt0 = PtPara->Integral(0,15);
139   Float_t IntPtS = PtPara->Integral(fPtMin,fPtMax);
140   Float_t PhiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
141   if (fAnalog == analog) {
142      fYWgt  = IntYS/fdNdy0;
143      fPtWgt = IntPtS/IntPt0;
144      fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
145   } else {
146       fYWgt = IntYS/fdNdy0;
147       fPtWgt = (fPtMax-fPtMin)/IntPt0;
148       fParentWeight = fYWgt*fPtWgt*PhiWgt/fNpart;
149   }
150 //
151 // particle decay related initialization
152   fPythia->DefineParticles();
153 // semimuonic decays of charm and beauty
154   fPythia->ForceDecay(fForceDecay);
155 //
156     switch (fForceDecay) 
157     {
158     case semielectronic:
159     case dielectron:
160     case b_jpsi_dielectron:
161     case b_psip_dielectron:
162         fChildSelect[0]=11;     
163         break;
164     case semimuonic:
165     case dimuon:
166     case b_jpsi_dimuon:
167     case b_psip_dimuon:
168         fChildSelect[0]=13;
169         break;
170     case pitomu:
171         fChildSelect[0]=13;
172         break;
173     case katomu:
174         fChildSelect[0]=13;
175         break;
176     case nodecay:
177         break;
178     case all:
179         break;
180     }
181 }
182
183 //____________________________________________________________
184 void AliGenParam::Generate()
185 {
186 //
187 // Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
188 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
189 // antineutrons in the the desired theta, phi and momentum windows; 
190 // Gaussian smearing on the vertex is done if selected. 
191 // The decay of heavy mesons is done using lujet, and the childern particle are tracked by GEANT
192 // However, light mesons are directly tracked by GEANT setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
193 //
194
195 // printf("Generate !!!!!!!!!!!!!\n");
196
197   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
198   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
199   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
200   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
201   Float_t p[3], pc[3], 
202           och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
203   Float_t ty, xmt;
204   Int_t nt, i, j, kfch[10];
205   Float_t  wgtp, wgtch;
206   Double_t dummy;
207   static TClonesArray *particles;
208   //
209   if(!particles) particles=new TClonesArray("TParticle",1000);
210   //
211   Float_t random[6];
212  
213 // Calculating vertex position per event
214   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
215   if(fVertexSmear==perEvent) {
216       gMC->Rndm(random,6);
217       for (j=0;j<3;j++) {
218           origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
219               TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
220       }
221   }
222   Int_t ipa=0;
223 // Generating fNpart particles
224   while (ipa<fNpart) {
225     while(1) {
226 //
227 // particle type
228           Int_t Ipart = fIpParaFunc();
229           fChildWeight=(fPythia->GetBraPart(Ipart))*fParentWeight;        
230           Float_t am=fPythia->GetPMAS(fPythia->Lucomp(Ipart),1);
231           gMC->Rndm(random,2);
232 //
233 // phi
234           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
235 //
236 // y
237           ty=Float_t(TMath::TanH(fYPara->GetRandom()));
238 //
239 // pT
240           if (fAnalog == analog) {
241               pt=fPtPara->GetRandom();
242               wgtp=fParentWeight;
243               wgtch=fChildWeight;
244           } else {
245               pt=fPtMin+random[1]*(fPtMax-fPtMin);
246               Double_t ptd=pt;
247               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
248               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
249           }
250           xmt=sqrt(pt*pt+am*am);
251           pl=xmt*ty/sqrt(1.-ty*ty);
252           theta=TMath::ATan2(pt,pl);
253           if(theta<fThetaMin || theta>fThetaMax) continue;
254           ptot=TMath::Sqrt(pt*pt+pl*pl);
255           if(ptot<fPMin || ptot>fPMax) continue;
256           p[0]=pt*TMath::Cos(phi);
257           p[1]=pt*TMath::Sin(phi);
258           p[2]=pl;
259           if(fVertexSmear==perTrack) {
260               gMC->Rndm(random,6);
261               for (j=0;j<3;j++) {
262                   origin0[j]=
263                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
264                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
265               }
266           }
267           
268 // Looking at fForceDecay : 
269 //                          if fForceDecay != none Primary particle decays using 
270 //                             AliPythia  and children are tracked by GEANT
271 //
272 //                          if fForceDecay == none Primary particle is tracked by GEANT 
273 //                          (In the latest, make sure that GEANT actually does all the decays you want)   
274 //
275           if (fForceDecay != nodecay) {
276 // Using lujet to decay particle
277               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
278               fPythia->DecayParticle(Ipart,energy,theta,phi);
279               //          fPythia->LuList(1);
280               //printf("origin0 %f %f %f\n",origin0[0],origin0[1],origin0[2]);
281               //printf("fCutOnChild %d \n",fCutOnChild);
282 //
283 // select muons
284               Int_t np=fPythia->ImportParticles(particles,"All");
285               //printf("np     %d \n",np);
286               Int_t ncsel=0;
287               for (i = 1; i<np; i++) {
288                   TParticle *  iparticle = (TParticle *) particles->At(i);
289                   Int_t kf = iparticle->GetPdgCode();
290 //
291 // children
292                   if (ChildSelected(TMath::Abs(kf)))
293                   {
294                       pc[0]=iparticle->Px();
295                       pc[1]=iparticle->Py();
296                       pc[2]=iparticle->Pz();
297                       och[0]=origin0[0]+iparticle->Vx()/10;
298                       och[1]=origin0[1]+iparticle->Vy()/10;
299                       och[2]=origin0[2]+iparticle->Vz()/10;
300                       if (fCutOnChild) {
301                           Float_t PtChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
302                           Float_t PChild=TMath::Sqrt(PtChild*PtChild+pc[2]*pc[2]);
303                           Float_t ThetaChild=TMath::ATan2(PtChild,pc[2]);
304                           Float_t PhiChild=TMath::ATan2(pc[1],pc[0])+TMath::Pi();
305                           Bool_t childok = 
306                               ((PtChild   > fPtMin   && PtChild   <fPtMax)      &&
307                                (PChild    > fPMin    && PChild    <fPMax)       &&
308                                (ThetaChild>fThetaMin && ThetaChild<fThetaMax)   &&
309                                (PhiChild  >  fPhiMin && PhiChild  <fPhiMax));
310                           if(childok)
311                           {
312                               pch[ncsel][0]=pc[0];
313                               pch[ncsel][1]=pc[1];
314                               pch[ncsel][2]=pc[2];
315                               kfch[ncsel]=kf;
316                               ncsel++;
317                           } else {
318                               ncsel=-1;
319                               break;
320                           } // child kine cuts
321                       } else {
322                           pch[ncsel][0]=pc[0];
323                           pch[ncsel][1]=pc[1];
324                           pch[ncsel][2]=pc[2];
325                           kfch[ncsel]=kf;
326                           ncsel++;
327                       } // if child selection
328                   } // select muon
329               } // decay particle loop
330               Int_t iparent;
331               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
332                   ipa++;
333 //
334 // parent
335                   gAlice->
336                       SetTrack(0,-1,Ipart,p,origin0,polar,0,"Primary",nt,wgtp);
337                   iparent=nt;
338
339                   for (i=0; i< ncsel; i++) {
340                       gAlice->SetTrack(fTrackIt,iparent,kfch[i],
341                                        &pch[i][0],och,polar,
342                                        0,"Decay",nt,wgtch);
343                       gAlice->KeepTrack(nt); 
344                   }
345               }  // Decays by Lujet
346           } // kinematic selection
347           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
348           {
349             gAlice->
350                 SetTrack(fTrackIt,-1,Ipart,p,origin0,polar,0,"Primary",nt,wgtp);
351             ipa++; 
352           }
353           break;
354     } // while
355   } // event loop
356 }
357
358 Bool_t AliGenParam::ChildSelected(Int_t ip)
359 {
360     for (Int_t i=0; i<5; i++)
361     {
362         if (fChildSelect[i]==ip) return kTRUE;
363     }
364     return kFALSE;
365 }
366
367 Bool_t AliGenParam::KinematicSelection(TParticle *particle)
368 {
369     Float_t px=particle->Px();
370     Float_t py=particle->Py();
371     Float_t pz=particle->Pz();
372 //
373 // momentum cut
374     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
375     if (p > fPMax || p < fPMin) 
376     {
377 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
378         return kFALSE;
379     }
380     Float_t pt=TMath::Sqrt(px*px+py*py);
381     
382 //
383 // theta cut
384     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
385     if (theta > fThetaMax || theta < fThetaMin) 
386     {
387 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
388         return kFALSE;
389     }
390
391     return kTRUE;
392 }
393
394
395
396