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