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
228Double_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
233double 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 240double 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 247double 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
302double 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 320Double_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 332Int_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
383Int_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//____________________________________________________________
432void 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//____________________________________________________________
514void AliGenParam::Generate()
515{
d5615b85 516 //
517 // Generate 1 event (see Generate(Int_t ntimes) for details
6078e216 518 //
d5615b85 519 GenerateN(1);
520}
521//____________________________________________________________
522void 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//____________________________________________________________________________________
820Float_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 841void 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