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