]>
Commit | Line | Data |
---|---|---|
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 | ||
88cb7938 | 16 | /* $Id$ */ |
6d4dd317 | 17 | |
18 | // Class to generate particles from using paramtrized pT and y distributions. | |
19 | // Distributions are obtained from pointer to object of type AliGenLib. | |
20 | // (For example AliGenMUONlib) | |
21 | // Decays are performed using Pythia. | |
22 | // andreas.morsch@cern.ch | |
23 | ||
ac3faee4 | 24 | #include <TCanvas.h> |
7ca4655f | 25 | #include <TClonesArray.h> |
fa480368 | 26 | #include <TDatabasePDG.h> |
ac3faee4 | 27 | #include <TF1.h> |
28 | #include <TH1F.h> | |
fa480368 | 29 | #include <TLorentzVector.h> |
ac3faee4 | 30 | #include <TMath.h> |
31 | #include <TParticle.h> | |
32 | #include <TParticlePDG.h> | |
828fff0d | 33 | #include <TROOT.h> |
ac3faee4 | 34 | #include <TVirtualMC.h> |
fa480368 | 35 | |
ac3faee4 | 36 | #include "AliDecayer.h" |
37 | #include "AliGenMUONlib.h" | |
38 | #include "AliGenParam.h" | |
5d12ce38 | 39 | #include "AliMC.h" |
ac3faee4 | 40 | #include "AliRun.h" |
7cf36cae | 41 | #include "AliGenEventHeader.h" |
fe4da5cc | 42 | |
43 | ClassImp(AliGenParam) | |
44 | ||
45 | //------------------------------------------------------------ | |
46 | ||
6078e216 | 47 | //Begin_Html |
48 | /* | |
1439f98e | 49 | <img src="picts/AliGenParam.gif"> |
6078e216 | 50 | */ |
51 | //End_Html | |
fe4da5cc | 52 | |
fe4da5cc | 53 | //____________________________________________________________ |
6078e216 | 54 | AliGenParam::AliGenParam() |
55 | : fPtParaFunc(0), | |
1c56e311 | 56 | fYParaFunc(0), |
57 | fIpParaFunc(0), | |
6078e216 | 58 | fV2ParaFunc(0), |
1c56e311 | 59 | fPtPara(0), |
60 | fYPara(0), | |
6078e216 | 61 | fV2Para(0), |
62 | fdNdPhi(0), | |
1c56e311 | 63 | fParam(0), |
64 | fdNdy0(0.), | |
65 | fYWgt(0.), | |
66 | fPtWgt(0.), | |
67 | fBias(0.), | |
68 | fTrials(0), | |
69 | fDeltaPt(0.01), | |
18e09c20 | 70 | fSelectAll(kFALSE), |
4ae1c9f0 | 71 | fDecayer(0), |
72 | fForceConv(kFALSE) | |
fe4da5cc | 73 | { |
6078e216 | 74 | // Default constructor |
b22ee262 | 75 | } |
2ad38d56 | 76 | //____________________________________________________________ |
f2c13e03 | 77 | AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library, Int_t param, const char* tname) |
1c56e311 | 78 | :AliGenMC(npart), |
79 | fPtParaFunc(Library->GetPt(param, tname)), | |
80 | fYParaFunc (Library->GetY (param, tname)), | |
81 | fIpParaFunc(Library->GetIp(param, tname)), | |
6078e216 | 82 | fV2ParaFunc(Library->GetV2(param, tname)), |
1c56e311 | 83 | fPtPara(0), |
84 | fYPara(0), | |
6078e216 | 85 | fV2Para(0), |
86 | fdNdPhi(0), | |
1c56e311 | 87 | fParam(param), |
88 | fdNdy0(0.), | |
89 | fYWgt(0.), | |
90 | fPtWgt(0.), | |
91 | fBias(0.), | |
92 | fTrials(0), | |
93 | fDeltaPt(0.01), | |
18e09c20 | 94 | fSelectAll(kFALSE), |
4ae1c9f0 | 95 | fDecayer(0), |
96 | fForceConv(kFALSE) | |
b22ee262 | 97 | { |
6078e216 | 98 | // Constructor using number of particles parameterisation id and library |
8b31bfac | 99 | fName = "Param"; |
100 | fTitle= "Particle Generator using pT and y parameterisation"; | |
34f60c01 | 101 | fAnalog = kAnalog; |
b22ee262 | 102 | SetForceDecay(); |
fe4da5cc | 103 | } |
fe4da5cc | 104 | //____________________________________________________________ |
1c56e311 | 105 | AliGenParam::AliGenParam(Int_t npart, Int_t param, const char* tname, const char* name): |
106 | AliGenMC(npart), | |
107 | fPtParaFunc(0), | |
108 | fYParaFunc (0), | |
109 | fIpParaFunc(0), | |
6078e216 | 110 | fV2ParaFunc(0), |
1c56e311 | 111 | fPtPara(0), |
112 | fYPara(0), | |
6078e216 | 113 | fV2Para(0), |
114 | fdNdPhi(0), | |
1c56e311 | 115 | fParam(param), |
116 | fdNdy0(0.), | |
117 | fYWgt(0.), | |
118 | fPtWgt(0.), | |
119 | fBias(0.), | |
120 | fTrials(0), | |
121 | fDeltaPt(0.01), | |
18e09c20 | 122 | fSelectAll(kFALSE), |
4ae1c9f0 | 123 | fDecayer(0), |
124 | fForceConv(kFALSE) | |
fe4da5cc | 125 | { |
6078e216 | 126 | // Constructor using parameterisation id and number of particles |
127 | // | |
1c56e311 | 128 | fName = name; |
129 | fTitle= "Particle Generator using pT and y parameterisation"; | |
8b31bfac | 130 | |
675e9664 | 131 | AliGenLib* pLibrary = new AliGenMUONlib(); |
675e9664 | 132 | fPtParaFunc = pLibrary->GetPt(param, tname); |
133 | fYParaFunc = pLibrary->GetY (param, tname); | |
134 | fIpParaFunc = pLibrary->GetIp(param, tname); | |
6078e216 | 135 | fV2ParaFunc = pLibrary->GetV2(param, tname); |
b22ee262 | 136 | |
34f60c01 | 137 | fAnalog = kAnalog; |
b22ee262 | 138 | fChildSelect.Set(5); |
139 | for (Int_t i=0; i<5; i++) fChildSelect[i]=0; | |
140 | SetForceDecay(); | |
141 | SetCutOnChild(); | |
142 | SetChildMomentumRange(); | |
143 | SetChildPtRange(); | |
144 | SetChildPhiRange(); | |
145 | SetChildThetaRange(); | |
fe4da5cc | 146 | } |
2ad38d56 | 147 | //____________________________________________________________ |
fe4da5cc | 148 | |
34f60c01 | 149 | AliGenParam::AliGenParam(Int_t npart, Int_t param, |
75e0cc59 | 150 | Double_t (*PtPara) (const Double_t*, const Double_t*), |
151 | Double_t (*YPara ) (const Double_t* ,const Double_t*), | |
6078e216 | 152 | Double_t (*V2Para) (const Double_t* ,const Double_t*), |
65fb704d | 153 | Int_t (*IpPara) (TRandom *)) |
1c56e311 | 154 | :AliGenMC(npart), |
155 | ||
156 | fPtParaFunc(PtPara), | |
157 | fYParaFunc(YPara), | |
158 | fIpParaFunc(IpPara), | |
6078e216 | 159 | fV2ParaFunc(V2Para), |
1c56e311 | 160 | fPtPara(0), |
161 | fYPara(0), | |
6078e216 | 162 | fV2Para(0), |
163 | fdNdPhi(0), | |
1c56e311 | 164 | fParam(param), |
165 | fdNdy0(0.), | |
166 | fYWgt(0.), | |
167 | fPtWgt(0.), | |
168 | fBias(0.), | |
169 | fTrials(0), | |
170 | fDeltaPt(0.01), | |
18e09c20 | 171 | fSelectAll(kFALSE), |
4ae1c9f0 | 172 | fDecayer(0), |
173 | fForceConv(kFALSE) | |
886b6f73 | 174 | { |
6078e216 | 175 | // Constructor |
176 | // Gines Martinez 1/10/99 | |
8b31bfac | 177 | fName = "Param"; |
178 | fTitle= "Particle Generator using pT and y parameterisation"; | |
179 | ||
34f60c01 | 180 | fAnalog = kAnalog; |
886b6f73 | 181 | fChildSelect.Set(5); |
182 | for (Int_t i=0; i<5; i++) fChildSelect[i]=0; | |
183 | SetForceDecay(); | |
184 | SetCutOnChild(); | |
749070b6 | 185 | SetChildMomentumRange(); |
186 | SetChildPtRange(); | |
187 | SetChildPhiRange(); | |
188 | SetChildThetaRange(); | |
886b6f73 | 189 | } |
190 | ||
fe4da5cc | 191 | //____________________________________________________________ |
192 | AliGenParam::~AliGenParam() | |
193 | { | |
6078e216 | 194 | // Destructor |
fe4da5cc | 195 | delete fPtPara; |
196 | delete fYPara; | |
6078e216 | 197 | delete fV2Para; |
198 | delete fdNdPhi; | |
fe4da5cc | 199 | } |
200 | ||
4ae1c9f0 | 201 | //------------------------------------------------------------------- |
202 | TVector3 AliGenParam::OrthogonalVector(TVector3 &inVec){ | |
203 | double abc[]={inVec.x(), inVec.y(), inVec.z()}; | |
204 | double xyz[]={1,1,1}; | |
205 | int solvDim=0; | |
206 | double tmp=abc[0]; | |
207 | for(int i=0; i<3; i++) | |
208 | if(abs(abc[i])>tmp){ | |
209 | solvDim=i; | |
210 | tmp=abs(abc[i]); | |
211 | } | |
212 | xyz[solvDim]=(-abc[(1+solvDim)%3]-abc[(2+solvDim)%3])/abc[(0+solvDim)%3]; | |
213 | ||
214 | TVector3 res(xyz[0],xyz[1],xyz[2]); | |
215 | return res; | |
216 | } | |
217 | ||
71443190 | 218 | void AliGenParam::RotateVector(Double_t *pin, Double_t *pout, Double_t costheta, Double_t sintheta, |
219 | Double_t cosphi, Double_t sinphi) | |
220 | { | |
221 | // Perform rotation | |
222 | pout[0] = pin[0]*costheta*cosphi-pin[1]*sinphi+pin[2]*sintheta*cosphi; | |
223 | pout[1] = pin[0]*costheta*sinphi+pin[1]*cosphi+pin[2]*sintheta*sinphi; | |
224 | pout[2] = -1.0 * pin[0] * sintheta + pin[2] * costheta; | |
225 | return; | |
226 | } | |
227 | ||
71443190 | 228 | double AliGenParam::ScreenFunction1(double screenVariable){ |
229 | if(screenVariable>1) | |
230 | return 42.24 - 8.368 * log(screenVariable + 0.952); | |
4ae1c9f0 | 231 | else |
71443190 | 232 | return 42.392 - screenVariable * (7.796 - 1.961 * screenVariable); |
4ae1c9f0 | 233 | } |
234 | ||
71443190 | 235 | double AliGenParam::ScreenFunction2(double screenVariable){ |
236 | if(screenVariable>1) | |
237 | return 42.24 - 8.368 * log(screenVariable + 0.952); | |
4ae1c9f0 | 238 | else |
71443190 | 239 | return 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable); |
4ae1c9f0 | 240 | } |
241 | ||
71443190 | 242 | double AliGenParam::RandomEnergyFraction(double Z, double photonEnergy){ |
4ae1c9f0 | 243 | double aZ=Z/137.036; |
71443190 | 244 | double epsilon ; |
245 | double epsilon0Local = 0.000511 / photonEnergy ; | |
246 | ||
247 | // Do it fast if photon energy < 2. MeV | |
248 | if (photonEnergy < 0.002 ) | |
249 | { | |
250 | epsilon = epsilon0Local + (0.5 - epsilon0Local) * fRandom->Rndm(); | |
251 | } | |
252 | else | |
253 | { | |
254 | double fZ = 8*log(Z)/3; | |
4ae1c9f0 | 255 | double fcZ=(aZ*aZ)*(1/(1+aZ*aZ)+0.20206-0.0368*aZ*aZ+0.0083*aZ*aZ*aZ); |
71443190 | 256 | if (photonEnergy > 0.050) fZ += 8*fcZ; |
257 | ||
258 | // Limits of the screening variable | |
259 | double screenFactor = 136. * epsilon0Local / pow (Z,1/3); | |
260 | double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ; | |
261 | double screenMin = std::min(4.*screenFactor,screenMax) ; | |
262 | ||
263 | // Limits of the energy sampling | |
264 | double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ; | |
265 | double epsilonMin = std::max(epsilon0Local,epsilon1); | |
266 | double epsilonRange = 0.5 - epsilonMin ; | |
267 | ||
268 | // Sample the energy rate of the created electron (or positron) | |
269 | double screen; | |
270 | double gReject ; | |
271 | ||
272 | double f10 = ScreenFunction1(screenMin) - fZ; | |
273 | double f20 = ScreenFunction2(screenMin) - fZ; | |
274 | double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.); | |
275 | double normF2 = std::max(1.5 * f20,0.); | |
276 | ||
277 | do | |
278 | { | |
279 | if (normF1 / (normF1 + normF2) > fRandom->Rndm() ) | |
280 | { | |
281 | epsilon = 0.5 - epsilonRange * pow(fRandom->Rndm(), 0.333333) ; | |
282 | screen = screenFactor / (epsilon * (1. - epsilon)); | |
283 | gReject = (ScreenFunction1(screen) - fZ) / f10 ; | |
284 | } | |
285 | else | |
286 | { | |
287 | epsilon = epsilonMin + epsilonRange * fRandom->Rndm(); | |
288 | screen = screenFactor / (epsilon * (1 - epsilon)); | |
289 | gReject = (ScreenFunction2(screen) - fZ) / f20 ; | |
290 | } | |
291 | } while ( gReject < fRandom->Rndm() ); | |
292 | ||
293 | } // End of epsilon sampling | |
294 | return epsilon; | |
295 | } | |
296 | ||
297 | double AliGenParam::RandomPolarAngle(){ | |
298 | double u; | |
299 | const double a1 = 0.625; | |
300 | double a2 = 3. * a1; | |
301 | // double d = 27. ; | |
302 | ||
303 | // if (9. / (9. + d) > fRandom->Rndm()) | |
304 | if (0.25 > fRandom->Rndm()) | |
305 | { | |
306 | u = - log(fRandom->Rndm() * fRandom->Rndm()) / a1 ; | |
307 | } | |
308 | else | |
309 | { | |
310 | u = - log(fRandom->Rndm() * fRandom->Rndm()) / a2 ; | |
311 | } | |
312 | return u*0.000511; | |
313 | } | |
4ae1c9f0 | 314 | |
71443190 | 315 | Double_t AliGenParam::RandomMass(Double_t mh){ |
4ae1c9f0 | 316 | while(true){ |
71443190 | 317 | double y=fRandom->Rndm(); |
318 | double mee=2*0.000511*TMath::Power(2*0.000511/mh,-y); //inverse of the enveloping cumulative distribution | |
6a8b015a | 319 | double apxkw=2.0/3.0/137.036/TMath::Pi()/mee; //enveloping probability density |
320 | double val=fRandom->Uniform(0,apxkw); | |
71443190 | 321 | double kw=apxkw*sqrt(1-4*0.000511*0.000511/mee/mee)*(1+2*0.000511*0.000511/mee/mee)*1*1*TMath::Power(1-mee*mee/mh/mh,3); |
322 | if(val<kw) | |
323 | return mee; | |
4ae1c9f0 | 324 | } |
325 | } | |
326 | ||
71443190 | 327 | Int_t AliGenParam::VirtualGammaPairProduction(TClonesArray *particles, Int_t nPart) |
328 | { | |
329 | Int_t nPartNew=nPart; | |
330 | for(int iPart=0; iPart<nPart; iPart++){ | |
331 | TParticle *gamma = (TParticle *) particles->At(iPart); | |
6a8b015a | 332 | if(gamma->GetPdgCode()!=220000) continue; |
333 | if(gamma->Pt()<0.002941) continue; //approximation of kw in AliGenEMlib is 0 below 0.002941 | |
71443190 | 334 | double mass=RandomMass(gamma->Pt()); |
335 | ||
336 | // lepton pair kinematics in virtual photon rest frame | |
337 | double Ee=mass/2; | |
338 | double Pe=TMath::Sqrt((Ee+0.000511)*(Ee-0.000511)); | |
339 | ||
340 | double costheta = (2.0 * gRandom->Rndm()) - 1.; | |
341 | double sintheta = TMath::Sqrt((1. + costheta) * (1. - costheta)); | |
342 | double phi = 2.0 * TMath::ACos(-1.) * gRandom->Rndm(); | |
343 | double sinphi = TMath::Sin(phi); | |
344 | double cosphi = TMath::Cos(phi); | |
345 | ||
346 | // momentum vectors of leptons in virtual photon rest frame | |
347 | Double_t pProd1[3] = {Pe * sintheta * cosphi, | |
348 | Pe * sintheta * sinphi, | |
349 | Pe * costheta}; | |
350 | ||
351 | Double_t pProd2[3] = {-1.0 * Pe * sintheta * cosphi, | |
352 | -1.0 * Pe * sintheta * sinphi, | |
353 | -1.0 * Pe * costheta}; | |
354 | ||
355 | // lepton 4-vectors in properly rotated virtual photon rest frame | |
356 | Double_t pRot1[3] = {0.}; | |
357 | RotateVector(pProd1, pRot1, costheta, -sintheta, -cosphi, -sinphi); | |
358 | Double_t pRot2[3] = {0.}; | |
359 | RotateVector(pProd2, pRot2, costheta, -sintheta, -cosphi, -sinphi); | |
360 | ||
361 | TLorentzVector e1V4(pRot1[0],pRot1[1],pRot1[2],Ee); | |
362 | TLorentzVector e2V4(pRot2[0],pRot2[1],pRot2[2],Ee); | |
363 | ||
364 | TVector3 boost(gamma->Px(),gamma->Py(),gamma->Pz()); | |
365 | boost*=1/sqrt(gamma->P()*gamma->P()+mass*mass); | |
366 | e1V4.Boost(boost); | |
367 | e2V4.Boost(boost); | |
368 | ||
369 | TLorentzVector vtx; | |
370 | gamma->ProductionVertex(vtx); | |
371 | new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e1V4, vtx); | |
372 | nPartNew++; | |
373 | new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, e2V4, vtx); | |
374 | nPartNew++; | |
375 | } | |
6a8b015a | 376 | return nPartNew; |
4ae1c9f0 | 377 | } |
378 | ||
379 | Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart) | |
380 | { | |
381 | //based on: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node27.html | |
382 | // and: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node58.html | |
71443190 | 383 | // and: G4LivermoreGammaConversionModel.cc |
4ae1c9f0 | 384 | Int_t nPartNew=nPart; |
385 | for(int iPart=0; iPart<nPart; iPart++){ | |
386 | TParticle *gamma = (TParticle *) particles->At(iPart); | |
387 | if(gamma->GetPdgCode()!=22) continue; | |
71443190 | 388 | if(gamma->Energy()<0.001022) continue; |
4ae1c9f0 | 389 | TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz()); |
71443190 | 390 | double frac=RandomEnergyFraction(1,gamma->Energy()); |
4ae1c9f0 | 391 | double Ee1=frac*gamma->Energy(); |
392 | double Ee2=(1-frac)*gamma->Energy(); | |
71443190 | 393 | double Pe1=sqrt((Ee1+0.000511)*(Ee1-0.000511)); |
394 | double Pe2=sqrt((Ee2+0.000511)*(Ee2-0.000511)); | |
4ae1c9f0 | 395 | |
396 | TVector3 rotAxis(OrthogonalVector(gammaV3)); | |
71443190 | 397 | Float_t az=fRandom->Uniform(TMath::Pi()*2); |
4ae1c9f0 | 398 | rotAxis.Rotate(az,gammaV3); |
399 | TVector3 e1V3(gammaV3); | |
71443190 | 400 | double u=RandomPolarAngle(); |
401 | e1V3.Rotate(u/Ee1,rotAxis); | |
4ae1c9f0 | 402 | e1V3=e1V3.Unit(); |
403 | e1V3*=Pe1; | |
404 | TVector3 e2V3(gammaV3); | |
71443190 | 405 | e2V3.Rotate(-u/Ee2,rotAxis); |
4ae1c9f0 | 406 | e2V3=e2V3.Unit(); |
407 | e2V3*=Pe2; | |
408 | // gamma = new TParticle(*gamma); | |
409 | // particles->RemoveAt(iPart); | |
410 | gamma->SetFirstDaughter(nPartNew+1); | |
411 | gamma->SetLastDaughter(nPartNew+2); | |
412 | // new((*particles)[iPart]) TParticle(*gamma); | |
413 | // delete gamma; | |
414 | ||
415 | TLorentzVector vtx; | |
416 | gamma->ProductionVertex(vtx); | |
417 | new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e1V3,Ee1), vtx); | |
418 | nPartNew++; | |
419 | new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx); | |
420 | nPartNew++; | |
421 | } | |
6a8b015a | 422 | return nPartNew; |
4ae1c9f0 | 423 | } |
424 | ||
fe4da5cc | 425 | //____________________________________________________________ |
426 | void AliGenParam::Init() | |
427 | { | |
6078e216 | 428 | // Initialisation |
fa480368 | 429 | |
2942f542 | 430 | if (TVirtualMC::GetMC()) fDecayer = TVirtualMC::GetMC()->GetDecayer(); |
fe4da5cc | 431 | //Begin_Html |
432 | /* | |
1439f98e | 433 | <img src="picts/AliGenParam.gif"> |
fe4da5cc | 434 | */ |
435 | //End_Html | |
c5c3e4b0 | 436 | char name[256]; |
8f8399ab | 437 | snprintf(name, 256, "pt-parameterisation for %s", GetName()); |
c5c3e4b0 | 438 | |
2ad38d56 | 439 | if (fPtPara) fPtPara->Delete(); |
c5c3e4b0 | 440 | fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0); |
828fff0d | 441 | gROOT->GetListOfFunctions()->Remove(fPtPara); |
6078e216 | 442 | // Set representation precision to 10 MeV |
c5c3e4b0 | 443 | Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt); |
749070b6 | 444 | |
445 | fPtPara->SetNpx(npx); | |
c5c3e4b0 | 446 | |
8f8399ab | 447 | snprintf(name, 256, "y-parameterisation for %s", GetName()); |
2ad38d56 | 448 | if (fYPara) fYPara->Delete(); |
c5c3e4b0 | 449 | fYPara = new TF1(name, fYParaFunc, fYMin, fYMax, 0); |
828fff0d | 450 | gROOT->GetListOfFunctions()->Remove(fYPara); |
451 | ||
6078e216 | 452 | snprintf(name, 256, "v2-parameterisation for %s", GetName()); |
453 | if (fV2Para) fV2Para->Delete(); | |
454 | fV2Para = new TF1(name, fV2ParaFunc, fPtMin, fPtMax, 0); | |
455 | // fV2Para = new TF1(name, "2*[0]/(1+TMath::Exp([1]*([2]-x)))-[0]", fPtMin, fPtMax); | |
456 | // fV2Para->SetParameter(0, 0.236910); | |
457 | // fV2Para->SetParameter(1, 1.71122); | |
458 | // fV2Para->SetParameter(2, 0.0827617); | |
459 | //gROOT->GetListOfFunctions()->Remove(fV2Para); //TR: necessary? | |
460 | ||
461 | snprintf(name, 256, "dNdPhi for %s", GetName()); | |
462 | if (fdNdPhi) fdNdPhi->Delete(); | |
463 | fdNdPhi = new TF1(name, "1+2*[0]*TMath::Cos(2*(x-[1]))", fPhiMin, fPhiMax); | |
464 | //gROOT->GetListOfFunctions()->Remove(fdNdPhi); //TR: necessary? | |
749070b6 | 465 | |
8f8399ab | 466 | snprintf(name, 256, "pt-for-%s", GetName()); |
c5c3e4b0 | 467 | TF1 ptPara(name ,fPtParaFunc, 0, 15, 0); |
8f8399ab | 468 | snprintf(name, 256, "y-for-%s", GetName()); |
c5c3e4b0 | 469 | TF1 yPara(name, fYParaFunc, -6, 6, 0); |
b7601ac4 | 470 | |
6078e216 | 471 | // |
472 | // dN/dy| y=0 | |
749070b6 | 473 | Double_t y1=0; |
474 | Double_t y2=0; | |
475 | ||
476 | fdNdy0=fYParaFunc(&y1,&y2); | |
6078e216 | 477 | // |
478 | // Integral over generation region | |
faf00b0d | 479 | #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0) |
389a9124 | 480 | Float_t intYS = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6); |
481 | Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6); | |
482 | Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6); | |
faf00b0d | 483 | #else |
484 | Float_t intYS = yPara.Integral(fYMin, fYMax,1.e-6); | |
485 | Float_t intPt0 = ptPara.Integral(0,15,1.e-6); | |
486 | Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,1.e-6); | |
487 | #endif | |
6078e216 | 488 | Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi(); //TR: should probably be done differently in case of anisotropic phi... |
34f60c01 | 489 | if (fAnalog == kAnalog) { |
749070b6 | 490 | fYWgt = intYS/fdNdy0; |
491 | fPtWgt = intPtS/intPt0; | |
492 | fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart; | |
493 | } else { | |
494 | fYWgt = intYS/fdNdy0; | |
495 | fPtWgt = (fPtMax-fPtMin)/intPt0; | |
496 | fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart; | |
497 | } | |
6078e216 | 498 | // |
499 | // particle decay related initialization | |
00d6ce7d | 500 | fDecayer->SetForceDecay(fForceDecay); |
fa480368 | 501 | fDecayer->Init(); |
00d6ce7d | 502 | |
6078e216 | 503 | // |
fff02fee | 504 | AliGenMC::Init(); |
fe4da5cc | 505 | } |
506 | ||
507 | //____________________________________________________________ | |
508 | void AliGenParam::Generate() | |
509 | { | |
d5615b85 | 510 | // |
511 | // Generate 1 event (see Generate(Int_t ntimes) for details | |
6078e216 | 512 | // |
d5615b85 | 513 | GenerateN(1); |
514 | } | |
515 | //____________________________________________________________ | |
516 | void AliGenParam::GenerateN(Int_t ntimes) | |
517 | { | |
518 | // | |
519 | // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion, | |
6078e216 | 520 | // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and |
521 | // antineutrons in the the desired theta, phi and momentum windows; | |
522 | // Gaussian smearing on the vertex is done if selected. | |
523 | // The decay of heavy mesons is done using lujet, | |
524 | // and the childern particle are tracked by GEANT | |
525 | // However, light mesons are directly tracked by GEANT | |
526 | // setting fForceDecay = nodecay (SetForceDecay(nodecay)) | |
527 | // | |
528 | // | |
529 | // Reinitialize decayer | |
2ad38d56 | 530 | fDecayer->SetForceDecay(fForceDecay); |
531 | fDecayer->Init(); | |
532 | ||
6078e216 | 533 | // |
28337bc1 | 534 | Float_t polar[3]= {0,0,0}; // Polarisation of the parent particle (for GEANT tracking) |
535 | Float_t origin0[3]; // Origin of the generated parent particle (for GEANT tracking) | |
21391258 | 536 | Float_t time0; // Time0 of the generated parent particle |
28337bc1 | 537 | Float_t pt, pl, ptot; // Transverse, logitudinal and total momenta of the parent particle |
538 | Float_t phi, theta; // Phi and theta spherical angles of the parent particle momentum | |
539 | Float_t p[3], pc[3], | |
fff02fee | 540 | och[3]; // Momentum, polarisation and origin of the children particles from lujet |
2d7a47be | 541 | Double_t ty, xmt; |
fff02fee | 542 | Int_t nt, i, j; |
fe4da5cc | 543 | Float_t wgtp, wgtch; |
544 | Double_t dummy; | |
09fd3ea2 | 545 | static TClonesArray *particles; |
fe4da5cc | 546 | // |
fff02fee | 547 | if(!particles) particles = new TClonesArray("TParticle",1000); |
8b31bfac | 548 | |
349be858 | 549 | TDatabasePDG *pDataBase = TDatabasePDG::Instance(); |
fe4da5cc | 550 | // |
551 | Float_t random[6]; | |
28337bc1 | 552 | |
6078e216 | 553 | // Calculating vertex position per event |
fe4da5cc | 554 | for (j=0;j<3;j++) origin0[j]=fOrigin[j]; |
21391258 | 555 | time0 = fTimeOrigin; |
aee8290b | 556 | if(fVertexSmear==kPerEvent) { |
d9ea0e3b | 557 | Vertex(); |
558 | for (j=0;j<3;j++) origin0[j]=fVertex[j]; | |
21391258 | 559 | time0 = fTime; |
fe4da5cc | 560 | } |
d9ea0e3b | 561 | |
21aaa175 | 562 | Int_t ipa=0; |
fff02fee | 563 | |
6078e216 | 564 | // Generating fNpart particles |
7cf36cae | 565 | fNprimaries = 0; |
566 | ||
d5615b85 | 567 | Int_t nGen = fNpart*ntimes; |
568 | while (ipa<nGen) { | |
cc5d764c | 569 | while(1) { |
6078e216 | 570 | // |
571 | // particle type | |
65fb704d | 572 | Int_t iPart = fIpParaFunc(fRandom); |
71443190 | 573 | Int_t iTemp = iPart; |
574 | ||
575 | // custom pdg codes for to destinguish direct photons | |
6a8b015a | 576 | if(iPart==220000) iPart=22; |
71443190 | 577 | |
fa480368 | 578 | fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight; |
675e9664 | 579 | TParticlePDG *particle = pDataBase->GetParticle(iPart); |
fa480368 | 580 | Float_t am = particle->Mass(); |
8b31bfac | 581 | |
65fb704d | 582 | Rndm(random,2); |
6078e216 | 583 | // |
584 | // y | |
2d7a47be | 585 | ty = TMath::TanH(fYPara->GetRandom()); |
71443190 | 586 | |
6078e216 | 587 | // |
588 | // pT | |
34f60c01 | 589 | if (fAnalog == kAnalog) { |
fe4da5cc | 590 | pt=fPtPara->GetRandom(); |
591 | wgtp=fParentWeight; | |
592 | wgtch=fChildWeight; | |
593 | } else { | |
594 | pt=fPtMin+random[1]*(fPtMax-fPtMin); | |
595 | Double_t ptd=pt; | |
596 | wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy); | |
597 | wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy); | |
598 | } | |
599 | xmt=sqrt(pt*pt+am*am); | |
2d7a47be | 600 | if (TMath::Abs(ty)==1.) { |
601 | ty=0.; | |
602 | Fatal("AliGenParam", | |
603 | "Division by 0: Please check you rapidity range !"); | |
604 | } | |
6078e216 | 605 | // |
606 | // phi | |
4ae1c9f0 | 607 | // if(!ipa) |
608 | //phi=fEvPlane; //align first particle of each event with event plane | |
609 | //else{ | |
6078e216 | 610 | double v2 = fV2Para->Eval(pt); |
611 | fdNdPhi->SetParameter(0,v2); | |
612 | fdNdPhi->SetParameter(1,fEvPlane); | |
613 | phi=fdNdPhi->GetRandom(); | |
4ae1c9f0 | 614 | // } |
2d7a47be | 615 | |
60e55aee | 616 | pl=xmt*ty/sqrt((1.-ty)*(1.+ty)); |
fe4da5cc | 617 | theta=TMath::ATan2(pt,pl); |
6078e216 | 618 | // Cut on theta |
fe4da5cc | 619 | if(theta<fThetaMin || theta>fThetaMax) continue; |
620 | ptot=TMath::Sqrt(pt*pt+pl*pl); | |
6078e216 | 621 | // Cut on momentum |
fe4da5cc | 622 | if(ptot<fPMin || ptot>fPMax) continue; |
6078e216 | 623 | // |
fe4da5cc | 624 | p[0]=pt*TMath::Cos(phi); |
625 | p[1]=pt*TMath::Sin(phi); | |
626 | p[2]=pl; | |
71443190 | 627 | |
aee8290b | 628 | if(fVertexSmear==kPerTrack) { |
65fb704d | 629 | Rndm(random,6); |
fe4da5cc | 630 | for (j=0;j<3;j++) { |
631 | origin0[j]= | |
632 | fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* | |
633 | TMath::Sqrt(-2*TMath::Log(random[2*j+1])); | |
634 | } | |
21391258 | 635 | Rndm(random,2); |
636 | time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()* | |
637 | TMath::Cos(2*random[0]*TMath::Pi())* | |
638 | TMath::Sqrt(-2*TMath::Log(random[1])); | |
fe4da5cc | 639 | } |
28337bc1 | 640 | |
6078e216 | 641 | // Looking at fForceDecay : |
642 | // if fForceDecay != none Primary particle decays using | |
643 | // AliPythia and children are tracked by GEANT | |
644 | // | |
645 | // if fForceDecay == none Primary particle is tracked by GEANT | |
646 | // (In the latest, make sure that GEANT actually does all the decays you want) | |
647 | // | |
cd716030 | 648 | Bool_t decayed = kFALSE; |
649 | ||
fff02fee | 650 | |
34f60c01 | 651 | if (fForceDecay != kNoDecay) { |
6078e216 | 652 | // Using lujet to decay particle |
886b6f73 | 653 | Float_t energy=TMath::Sqrt(ptot*ptot+am*am); |
fa480368 | 654 | TLorentzVector pmom(p[0], p[1], p[2], energy); |
655 | fDecayer->Decay(iPart,&pmom); | |
6078e216 | 656 | // |
657 | // select decay particles | |
fa480368 | 658 | Int_t np=fDecayer->ImportParticles(particles); |
1242532d | 659 | |
71443190 | 660 | iPart=iTemp; |
6a8b015a | 661 | if(iPart==220000){ |
662 | TParticle *gamma = (TParticle *)particles->At(0); | |
663 | gamma->SetPdgCode(iPart); | |
664 | np=VirtualGammaPairProduction(particles,np); | |
71443190 | 665 | } |
6a8b015a | 666 | if(fForceConv) np=ForceGammaConversion(particles,np); |
4ae1c9f0 | 667 | |
1242532d | 668 | // Selecting GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut; |
669 | if (fGeometryAcceptance) | |
670 | if (!CheckAcceptanceGeometry(np,particles)) continue; | |
886b6f73 | 671 | Int_t ncsel=0; |
fff02fee | 672 | Int_t* pFlag = new Int_t[np]; |
673 | Int_t* pParent = new Int_t[np]; | |
674 | Int_t* pSelected = new Int_t[np]; | |
675 | Int_t* trackIt = new Int_t[np]; | |
676 | ||
677 | for (i=0; i<np; i++) { | |
678 | pFlag[i] = 0; | |
679 | pSelected[i] = 0; | |
680 | pParent[i] = -1; | |
681 | } | |
682 | ||
6ba00c52 | 683 | if (np >1) { |
cd716030 | 684 | decayed = kTRUE; |
08bffa4d | 685 | TParticle* iparticle = 0; |
fff02fee | 686 | Int_t ipF, ipL; |
687 | for (i = 1; i<np ; i++) { | |
688 | trackIt[i] = 1; | |
6ba00c52 | 689 | iparticle = (TParticle *) particles->At(i); |
690 | Int_t kf = iparticle->GetPdgCode(); | |
fff02fee | 691 | Int_t ks = iparticle->GetStatusCode(); |
6078e216 | 692 | // flagged particle |
fff02fee | 693 | |
694 | if (pFlag[i] == 1) { | |
fff02fee | 695 | ipF = iparticle->GetFirstDaughter(); |
696 | ipL = iparticle->GetLastDaughter(); | |
697 | if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1; | |
698 | continue; | |
699 | } | |
700 | ||
6078e216 | 701 | // flag decay products of particles with long life-time (c tau > .3 mum) |
fff02fee | 702 | |
703 | if (ks != 1) { | |
6078e216 | 704 | // TParticlePDG *particle = pDataBase->GetParticle(kf); |
fff02fee | 705 | |
706 | Double_t lifeTime = fDecayer->GetLifetime(kf); | |
6078e216 | 707 | // Double_t mass = particle->Mass(); |
708 | // Double_t width = particle->Width(); | |
47fc6bd5 | 709 | if (lifeTime > (Double_t) fMaxLifeTime) { |
fff02fee | 710 | ipF = iparticle->GetFirstDaughter(); |
711 | ipL = iparticle->GetLastDaughter(); | |
712 | if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1; | |
713 | } else{ | |
714 | trackIt[i] = 0; | |
715 | pSelected[i] = 1; | |
716 | } | |
717 | } // ks==1 ? | |
6078e216 | 718 | // |
719 | // children | |
47fc6bd5 | 720 | |
18e09c20 | 721 | if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i]) |
6ba00c52 | 722 | { |
6ba00c52 | 723 | if (fCutOnChild) { |
fff02fee | 724 | pc[0]=iparticle->Px(); |
725 | pc[1]=iparticle->Py(); | |
726 | pc[2]=iparticle->Pz(); | |
727 | Bool_t childok = KinematicSelection(iparticle, 1); | |
728 | if(childok) { | |
729 | pSelected[i] = 1; | |
6ba00c52 | 730 | ncsel++; |
731 | } else { | |
732 | ncsel=-1; | |
733 | break; | |
734 | } // child kine cuts | |
735 | } else { | |
fff02fee | 736 | pSelected[i] = 1; |
886b6f73 | 737 | ncsel++; |
6ba00c52 | 738 | } // if child selection |
739 | } // select muon | |
740 | } // decay particle loop | |
741 | } // if decay products | |
742 | ||
886b6f73 | 743 | Int_t iparent; |
744 | if ((fCutOnChild && ncsel >0) || !fCutOnChild){ | |
745 | ipa++; | |
6078e216 | 746 | // |
747 | // Parent | |
cd716030 | 748 | |
749 | ||
21391258 | 750 | PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1)); |
fff02fee | 751 | pParent[0] = nt; |
a99cf51f | 752 | KeepTrack(nt); |
7cf36cae | 753 | fNprimaries++; |
754 | ||
6078e216 | 755 | // |
756 | // Decay Products | |
757 | // | |
fff02fee | 758 | for (i = 1; i < np; i++) { |
759 | if (pSelected[i]) { | |
760 | TParticle* iparticle = (TParticle *) particles->At(i); | |
ea79897e | 761 | Int_t kf = iparticle->GetPdgCode(); |
cd716030 | 762 | Int_t ksc = iparticle->GetStatusCode(); |
ea79897e | 763 | Int_t jpa = iparticle->GetFirstMother()-1; |
fff02fee | 764 | |
774ceaaf | 765 | och[0] = origin0[0]+iparticle->Vx(); |
766 | och[1] = origin0[1]+iparticle->Vy(); | |
767 | och[2] = origin0[2]+iparticle->Vz(); | |
fff02fee | 768 | pc[0] = iparticle->Px(); |
769 | pc[1] = iparticle->Py(); | |
770 | pc[2] = iparticle->Pz(); | |
771 | ||
ea79897e | 772 | if (jpa > -1) { |
773 | iparent = pParent[jpa]; | |
fff02fee | 774 | } else { |
775 | iparent = -1; | |
776 | } | |
1242532d | 777 | |
ea79897e | 778 | PushTrack(fTrackIt * trackIt[i], iparent, kf, |
fff02fee | 779 | pc, och, polar, |
21391258 | 780 | time0 + iparticle->T(), kPDecay, nt, wgtch, ksc); |
fff02fee | 781 | pParent[i] = nt; |
a99cf51f | 782 | KeepTrack(nt); |
7cf36cae | 783 | fNprimaries++; |
8b31bfac | 784 | } // Selected |
785 | } // Particle loop | |
28337bc1 | 786 | } // Decays by Lujet |
8b31bfac | 787 | particles->Clear(); |
fff02fee | 788 | if (pFlag) delete[] pFlag; |
789 | if (pParent) delete[] pParent; | |
790 | if (pSelected) delete[] pSelected; | |
8b31bfac | 791 | if (trackIt) delete[] trackIt; |
21aaa175 | 792 | } // kinematic selection |
28337bc1 | 793 | else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons) |
794 | { | |
5d12ce38 | 795 | gAlice->GetMCApp()-> |
21391258 | 796 | PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1); |
28337bc1 | 797 | ipa++; |
7cf36cae | 798 | fNprimaries++; |
28337bc1 | 799 | } |
fe4da5cc | 800 | break; |
21aaa175 | 801 | } // while |
fe4da5cc | 802 | } // event loop |
7cf36cae | 803 | |
a99cf51f | 804 | SetHighWaterMark(nt); |
7cf36cae | 805 | |
806 | AliGenEventHeader* header = new AliGenEventHeader("PARAM"); | |
807 | header->SetPrimaryVertex(fVertex); | |
21391258 | 808 | header->SetInteractionTime(fTime); |
7cf36cae | 809 | header->SetNProduced(fNprimaries); |
810 | AddHeader(header); | |
fe4da5cc | 811 | } |
2ad38d56 | 812 | //____________________________________________________________________________________ |
813 | Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax) | |
814 | { | |
6078e216 | 815 | // |
816 | // Normalisation for selected kinematic region | |
817 | // | |
faf00b0d | 818 | #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0) |
2ad38d56 | 819 | Float_t ratio = |
389a9124 | 820 | fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) * |
821 | fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6) * | |
2ad38d56 | 822 | (phiMax-phiMin)/360.; |
faf00b0d | 823 | #else |
824 | Float_t ratio = | |
825 | fPtPara->Integral(ptMin,ptMax,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),1.e-6) * | |
826 | fYPara->Integral(yMin,yMax,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),1.e-6) * | |
827 | (phiMax-phiMin)/360.; | |
828 | #endif | |
2ad38d56 | 829 | return TMath::Abs(ratio); |
830 | } | |
831 | ||
832 | //____________________________________________________________________________________ | |
fe4da5cc | 833 | |
dc1d768c | 834 | void AliGenParam::Draw( const char * /*opt*/) |
5bd39445 | 835 | { |
836 | // | |
837 | // Draw the pT and y Distributions | |
838 | // | |
839 | TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700); | |
840 | c0->Divide(2,1); | |
841 | c0->cd(1); | |
842 | fPtPara->Draw(); | |
843 | fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)"); | |
844 | c0->cd(2); | |
845 | fYPara->Draw(); | |
846 | fYPara->GetHistogram()->SetXTitle("y"); | |
847 | } | |
848 | ||
f87cfe57 | 849 | |
fe4da5cc | 850 | |
851 |