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