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