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