doxy: MUON macros, refs, src links, no coll graph
[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),
309923cb 72 fForceConv(kFALSE),
73 fKeepParent(kFALSE),
74 fKeepIfOneChildSelected(kFALSE)
fe4da5cc 75{
6078e216 76 // Default constructor
b22ee262 77}
2ad38d56 78//____________________________________________________________
f2c13e03 79AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library, Int_t param, const char* tname)
1c56e311 80 :AliGenMC(npart),
81 fPtParaFunc(Library->GetPt(param, tname)),
82 fYParaFunc (Library->GetY (param, tname)),
83 fIpParaFunc(Library->GetIp(param, tname)),
6078e216 84 fV2ParaFunc(Library->GetV2(param, tname)),
1c56e311 85 fPtPara(0),
86 fYPara(0),
6078e216 87 fV2Para(0),
88 fdNdPhi(0),
1c56e311 89 fParam(param),
90 fdNdy0(0.),
91 fYWgt(0.),
92 fPtWgt(0.),
93 fBias(0.),
94 fTrials(0),
95 fDeltaPt(0.01),
18e09c20 96 fSelectAll(kFALSE),
4ae1c9f0 97 fDecayer(0),
309923cb 98 fForceConv(kFALSE),
99 fKeepParent(kFALSE),
100 fKeepIfOneChildSelected(kFALSE)
b22ee262 101{
6078e216 102 // Constructor using number of particles parameterisation id and library
8b31bfac 103 fName = "Param";
104 fTitle= "Particle Generator using pT and y parameterisation";
34f60c01 105 fAnalog = kAnalog;
b22ee262 106 SetForceDecay();
fe4da5cc 107}
fe4da5cc 108//____________________________________________________________
1c56e311 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),
6078e216 114 fV2ParaFunc(0),
1c56e311 115 fPtPara(0),
116 fYPara(0),
6078e216 117 fV2Para(0),
118 fdNdPhi(0),
1c56e311 119 fParam(param),
120 fdNdy0(0.),
121 fYWgt(0.),
122 fPtWgt(0.),
123 fBias(0.),
124 fTrials(0),
125 fDeltaPt(0.01),
18e09c20 126 fSelectAll(kFALSE),
4ae1c9f0 127 fDecayer(0),
309923cb 128 fForceConv(kFALSE),
129 fKeepParent(kFALSE),
130 fKeepIfOneChildSelected(kFALSE)
fe4da5cc 131{
6078e216 132 // Constructor using parameterisation id and number of particles
133 //
1c56e311 134 fName = name;
135 fTitle= "Particle Generator using pT and y parameterisation";
8b31bfac 136
675e9664 137 AliGenLib* pLibrary = new AliGenMUONlib();
675e9664 138 fPtParaFunc = pLibrary->GetPt(param, tname);
139 fYParaFunc = pLibrary->GetY (param, tname);
140 fIpParaFunc = pLibrary->GetIp(param, tname);
6078e216 141 fV2ParaFunc = pLibrary->GetV2(param, tname);
b22ee262 142
34f60c01 143 fAnalog = kAnalog;
b22ee262 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();
fe4da5cc 152}
2ad38d56 153//____________________________________________________________
fe4da5cc 154
34f60c01 155AliGenParam::AliGenParam(Int_t npart, Int_t param,
75e0cc59 156 Double_t (*PtPara) (const Double_t*, const Double_t*),
157 Double_t (*YPara ) (const Double_t* ,const Double_t*),
6078e216 158 Double_t (*V2Para) (const Double_t* ,const Double_t*),
65fb704d 159 Int_t (*IpPara) (TRandom *))
1c56e311 160 :AliGenMC(npart),
161
162 fPtParaFunc(PtPara),
163 fYParaFunc(YPara),
164 fIpParaFunc(IpPara),
6078e216 165 fV2ParaFunc(V2Para),
1c56e311 166 fPtPara(0),
167 fYPara(0),
6078e216 168 fV2Para(0),
169 fdNdPhi(0),
1c56e311 170 fParam(param),
171 fdNdy0(0.),
172 fYWgt(0.),
173 fPtWgt(0.),
174 fBias(0.),
175 fTrials(0),
176 fDeltaPt(0.01),
18e09c20 177 fSelectAll(kFALSE),
4ae1c9f0 178 fDecayer(0),
309923cb 179 fForceConv(kFALSE),
180 fKeepParent(kFALSE),
181 fKeepIfOneChildSelected(kFALSE)
886b6f73 182{
6078e216 183 // Constructor
184 // Gines Martinez 1/10/99
8b31bfac 185 fName = "Param";
186 fTitle= "Particle Generator using pT and y parameterisation";
187
34f60c01 188 fAnalog = kAnalog;
886b6f73 189 fChildSelect.Set(5);
190 for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
191 SetForceDecay();
192 SetCutOnChild();
749070b6 193 SetChildMomentumRange();
194 SetChildPtRange();
195 SetChildPhiRange();
196 SetChildThetaRange();
886b6f73 197}
198
fe4da5cc 199//____________________________________________________________
200AliGenParam::~AliGenParam()
201{
6078e216 202 // Destructor
fe4da5cc 203 delete fPtPara;
204 delete fYPara;
6078e216 205 delete fV2Para;
206 delete fdNdPhi;
fe4da5cc 207}
208
4ae1c9f0 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++)
06e32752 216 if(fabs(abc[i])>tmp){
4ae1c9f0 217 solvDim=i;
06e32752 218 tmp=fabs(abc[i]);
4ae1c9f0 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
71443190 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
71443190 236double AliGenParam::ScreenFunction1(double screenVariable){
237 if(screenVariable>1)
238 return 42.24 - 8.368 * log(screenVariable + 0.952);
4ae1c9f0 239 else
71443190 240 return 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
4ae1c9f0 241}
242
71443190 243double AliGenParam::ScreenFunction2(double screenVariable){
244 if(screenVariable>1)
245 return 42.24 - 8.368 * log(screenVariable + 0.952);
4ae1c9f0 246 else
71443190 247 return 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
4ae1c9f0 248}
249
71443190 250double AliGenParam::RandomEnergyFraction(double Z, double photonEnergy){
4ae1c9f0 251 double aZ=Z/137.036;
71443190 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;
4ae1c9f0 263 double fcZ=(aZ*aZ)*(1/(1+aZ*aZ)+0.20206-0.0368*aZ*aZ+0.0083*aZ*aZ*aZ);
71443190 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}
4ae1c9f0 322
71443190 323Double_t AliGenParam::RandomMass(Double_t mh){
4ae1c9f0 324 while(true){
71443190 325 double y=fRandom->Rndm();
326 double mee=2*0.000511*TMath::Power(2*0.000511/mh,-y); //inverse of the enveloping cumulative distribution
6a8b015a 327 double apxkw=2.0/3.0/137.036/TMath::Pi()/mee; //enveloping probability density
328 double val=fRandom->Uniform(0,apxkw);
71443190 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;
4ae1c9f0 332 }
333}
334
71443190 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);
6a8b015a 340 if(gamma->GetPdgCode()!=220000) continue;
341 if(gamma->Pt()<0.002941) continue; //approximation of kw in AliGenEMlib is 0 below 0.002941
71443190 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 }
6a8b015a 384 return nPartNew;
4ae1c9f0 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
71443190 391 // and: G4LivermoreGammaConversionModel.cc
4ae1c9f0 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;
71443190 396 if(gamma->Energy()<0.001022) continue;
4ae1c9f0 397 TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
71443190 398 double frac=RandomEnergyFraction(1,gamma->Energy());
4ae1c9f0 399 double Ee1=frac*gamma->Energy();
400 double Ee2=(1-frac)*gamma->Energy();
71443190 401 double Pe1=sqrt((Ee1+0.000511)*(Ee1-0.000511));
402 double Pe2=sqrt((Ee2+0.000511)*(Ee2-0.000511));
4ae1c9f0 403
404 TVector3 rotAxis(OrthogonalVector(gammaV3));
71443190 405 Float_t az=fRandom->Uniform(TMath::Pi()*2);
4ae1c9f0 406 rotAxis.Rotate(az,gammaV3);
407 TVector3 e1V3(gammaV3);
71443190 408 double u=RandomPolarAngle();
409 e1V3.Rotate(u/Ee1,rotAxis);
4ae1c9f0 410 e1V3=e1V3.Unit();
411 e1V3*=Pe1;
412 TVector3 e2V3(gammaV3);
71443190 413 e2V3.Rotate(-u/Ee2,rotAxis);
4ae1c9f0 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 }
6a8b015a 430 return nPartNew;
4ae1c9f0 431}
432
fe4da5cc 433//____________________________________________________________
434void AliGenParam::Init()
435{
6078e216 436 // Initialisation
fa480368 437
2942f542 438 if (TVirtualMC::GetMC()) fDecayer = TVirtualMC::GetMC()->GetDecayer();
fe4da5cc 439 //Begin_Html
440 /*
1439f98e 441 <img src="picts/AliGenParam.gif">
fe4da5cc 442 */
443 //End_Html
c5c3e4b0 444 char name[256];
8f8399ab 445 snprintf(name, 256, "pt-parameterisation for %s", GetName());
c5c3e4b0 446
2ad38d56 447 if (fPtPara) fPtPara->Delete();
c5c3e4b0 448 fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
828fff0d 449 gROOT->GetListOfFunctions()->Remove(fPtPara);
6078e216 450 // Set representation precision to 10 MeV
c5c3e4b0 451 Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
749070b6 452
453 fPtPara->SetNpx(npx);
c5c3e4b0 454
8f8399ab 455 snprintf(name, 256, "y-parameterisation for %s", GetName());
2ad38d56 456 if (fYPara) fYPara->Delete();
c5c3e4b0 457 fYPara = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
828fff0d 458 gROOT->GetListOfFunctions()->Remove(fYPara);
459
6078e216 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?
749070b6 473
8f8399ab 474 snprintf(name, 256, "pt-for-%s", GetName());
c5c3e4b0 475 TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
8f8399ab 476 snprintf(name, 256, "y-for-%s", GetName());
c5c3e4b0 477 TF1 yPara(name, fYParaFunc, -6, 6, 0);
b7601ac4 478
6078e216 479 //
480 // dN/dy| y=0
749070b6 481 Double_t y1=0;
482 Double_t y2=0;
483
484 fdNdy0=fYParaFunc(&y1,&y2);
6078e216 485 //
486 // Integral over generation region
faf00b0d 487#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
389a9124 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);
faf00b0d 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
6078e216 496 Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi(); //TR: should probably be done differently in case of anisotropic phi...
34f60c01 497 if (fAnalog == kAnalog) {
749070b6 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 }
6078e216 506 //
507 // particle decay related initialization
00d6ce7d 508 fDecayer->SetForceDecay(fForceDecay);
fa480368 509 fDecayer->Init();
00d6ce7d 510
6078e216 511 //
fff02fee 512 AliGenMC::Init();
fe4da5cc 513}
514
515//____________________________________________________________
516void AliGenParam::Generate()
517{
d5615b85 518 //
519 // Generate 1 event (see Generate(Int_t ntimes) for details
6078e216 520 //
d5615b85 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,
6078e216 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
2ad38d56 538 fDecayer->SetForceDecay(fForceDecay);
539 fDecayer->Init();
540
6078e216 541 //
28337bc1 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)
21391258 544 Float_t time0; // Time0 of the generated parent particle
28337bc1 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],
fff02fee 548 och[3]; // Momentum, polarisation and origin of the children particles from lujet
2d7a47be 549 Double_t ty, xmt;
fff02fee 550 Int_t nt, i, j;
fe4da5cc 551 Float_t wgtp, wgtch;
552 Double_t dummy;
09fd3ea2 553 static TClonesArray *particles;
fe4da5cc 554 //
fff02fee 555 if(!particles) particles = new TClonesArray("TParticle",1000);
8b31bfac 556
349be858 557 TDatabasePDG *pDataBase = TDatabasePDG::Instance();
fe4da5cc 558 //
559 Float_t random[6];
28337bc1 560
6078e216 561 // Calculating vertex position per event
fe4da5cc 562 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
21391258 563 time0 = fTimeOrigin;
aee8290b 564 if(fVertexSmear==kPerEvent) {
d9ea0e3b 565 Vertex();
566 for (j=0;j<3;j++) origin0[j]=fVertex[j];
21391258 567 time0 = fTime;
fe4da5cc 568 }
d9ea0e3b 569
21aaa175 570 Int_t ipa=0;
fff02fee 571
6078e216 572 // Generating fNpart particles
7cf36cae 573 fNprimaries = 0;
574
d5615b85 575 Int_t nGen = fNpart*ntimes;
576 while (ipa<nGen) {
cc5d764c 577 while(1) {
6078e216 578 //
579 // particle type
65fb704d 580 Int_t iPart = fIpParaFunc(fRandom);
71443190 581 Int_t iTemp = iPart;
582
583 // custom pdg codes for to destinguish direct photons
6a8b015a 584 if(iPart==220000) iPart=22;
71443190 585
fa480368 586 fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
675e9664 587 TParticlePDG *particle = pDataBase->GetParticle(iPart);
fa480368 588 Float_t am = particle->Mass();
8b31bfac 589
65fb704d 590 Rndm(random,2);
6078e216 591 //
592 // y
2d7a47be 593 ty = TMath::TanH(fYPara->GetRandom());
71443190 594
6078e216 595 //
596 // pT
34f60c01 597 if (fAnalog == kAnalog) {
fe4da5cc 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);
2d7a47be 608 if (TMath::Abs(ty)==1.) {
609 ty=0.;
610 Fatal("AliGenParam",
611 "Division by 0: Please check you rapidity range !");
612 }
6078e216 613 //
614 // phi
4ae1c9f0 615 // if(!ipa)
616 //phi=fEvPlane; //align first particle of each event with event plane
617 //else{
6078e216 618 double v2 = fV2Para->Eval(pt);
619 fdNdPhi->SetParameter(0,v2);
620 fdNdPhi->SetParameter(1,fEvPlane);
621 phi=fdNdPhi->GetRandom();
4ae1c9f0 622 // }
2d7a47be 623
60e55aee 624 pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
fe4da5cc 625 theta=TMath::ATan2(pt,pl);
6078e216 626 // Cut on theta
fe4da5cc 627 if(theta<fThetaMin || theta>fThetaMax) continue;
628 ptot=TMath::Sqrt(pt*pt+pl*pl);
6078e216 629 // Cut on momentum
fe4da5cc 630 if(ptot<fPMin || ptot>fPMax) continue;
6078e216 631 //
fe4da5cc 632 p[0]=pt*TMath::Cos(phi);
633 p[1]=pt*TMath::Sin(phi);
634 p[2]=pl;
71443190 635
aee8290b 636 if(fVertexSmear==kPerTrack) {
65fb704d 637 Rndm(random,6);
fe4da5cc 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 }
21391258 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]));
fe4da5cc 647 }
28337bc1 648
6078e216 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 //
cd716030 656 Bool_t decayed = kFALSE;
657
fff02fee 658
34f60c01 659 if (fForceDecay != kNoDecay) {
6078e216 660 // Using lujet to decay particle
886b6f73 661 Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
fa480368 662 TLorentzVector pmom(p[0], p[1], p[2], energy);
663 fDecayer->Decay(iPart,&pmom);
6078e216 664 //
665 // select decay particles
fa480368 666 Int_t np=fDecayer->ImportParticles(particles);
1242532d 667
71443190 668 iPart=iTemp;
6a8b015a 669 if(iPart==220000){
670 TParticle *gamma = (TParticle *)particles->At(0);
671 gamma->SetPdgCode(iPart);
672 np=VirtualGammaPairProduction(particles,np);
71443190 673 }
6a8b015a 674 if(fForceConv) np=ForceGammaConversion(particles,np);
4ae1c9f0 675
1242532d 676 // Selecting GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
677 if (fGeometryAcceptance)
678 if (!CheckAcceptanceGeometry(np,particles)) continue;
886b6f73 679 Int_t ncsel=0;
fff02fee 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
6ba00c52 691 if (np >1) {
cd716030 692 decayed = kTRUE;
08bffa4d 693 TParticle* iparticle = 0;
fff02fee 694 Int_t ipF, ipL;
695 for (i = 1; i<np ; i++) {
696 trackIt[i] = 1;
6ba00c52 697 iparticle = (TParticle *) particles->At(i);
698 Int_t kf = iparticle->GetPdgCode();
fff02fee 699 Int_t ks = iparticle->GetStatusCode();
6078e216 700 // flagged particle
fff02fee 701
702 if (pFlag[i] == 1) {
fff02fee 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
6078e216 709 // flag decay products of particles with long life-time (c tau > .3 mum)
fff02fee 710
711 if (ks != 1) {
6078e216 712 // TParticlePDG *particle = pDataBase->GetParticle(kf);
fff02fee 713
714 Double_t lifeTime = fDecayer->GetLifetime(kf);
6078e216 715 // Double_t mass = particle->Mass();
716 // Double_t width = particle->Width();
47fc6bd5 717 if (lifeTime > (Double_t) fMaxLifeTime) {
fff02fee 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 ?
6078e216 726 //
727 // children
47fc6bd5 728
18e09c20 729 if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
6ba00c52 730 {
6ba00c52 731 if (fCutOnChild) {
fff02fee 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;
6ba00c52 738 ncsel++;
739 } else {
309923cb 740 if(!fKeepIfOneChildSelected){
741 ncsel=-1;
742 break;
743 }
6ba00c52 744 } // child kine cuts
745 } else {
fff02fee 746 pSelected[i] = 1;
886b6f73 747 ncsel++;
6ba00c52 748 } // if child selection
749 } // select muon
750 } // decay particle loop
751 } // if decay products
752
886b6f73 753 Int_t iparent;
309923cb 754
755 if (fKeepParent || (fCutOnChild && ncsel >0) || !fCutOnChild){
6078e216 756 //
757 // Parent
cd716030 758
759
21391258 760 PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
fff02fee 761 pParent[0] = nt;
a99cf51f 762 KeepTrack(nt);
7cf36cae 763 fNprimaries++;
764
309923cb 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
6078e216 770 //
771 // Decay Products
772 //
fff02fee 773 for (i = 1; i < np; i++) {
774 if (pSelected[i]) {
775 TParticle* iparticle = (TParticle *) particles->At(i);
ea79897e 776 Int_t kf = iparticle->GetPdgCode();
cd716030 777 Int_t ksc = iparticle->GetStatusCode();
ea79897e 778 Int_t jpa = iparticle->GetFirstMother()-1;
fff02fee 779
774ceaaf 780 och[0] = origin0[0]+iparticle->Vx();
781 och[1] = origin0[1]+iparticle->Vy();
782 och[2] = origin0[2]+iparticle->Vz();
fff02fee 783 pc[0] = iparticle->Px();
784 pc[1] = iparticle->Py();
785 pc[2] = iparticle->Pz();
786
ea79897e 787 if (jpa > -1) {
788 iparent = pParent[jpa];
fff02fee 789 } else {
790 iparent = -1;
791 }
1242532d 792
ea79897e 793 PushTrack(fTrackIt * trackIt[i], iparent, kf,
fff02fee 794 pc, och, polar,
21391258 795 time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
fff02fee 796 pParent[i] = nt;
a99cf51f 797 KeepTrack(nt);
7cf36cae 798 fNprimaries++;
8b31bfac 799 } // Selected
800 } // Particle loop
28337bc1 801 } // Decays by Lujet
8b31bfac 802 particles->Clear();
fff02fee 803 if (pFlag) delete[] pFlag;
804 if (pParent) delete[] pParent;
805 if (pSelected) delete[] pSelected;
8b31bfac 806 if (trackIt) delete[] trackIt;
21aaa175 807 } // kinematic selection
28337bc1 808 else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
809 {
5d12ce38 810 gAlice->GetMCApp()->
21391258 811 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
28337bc1 812 ipa++;
7cf36cae 813 fNprimaries++;
28337bc1 814 }
fe4da5cc 815 break;
21aaa175 816 } // while
fe4da5cc 817 } // event loop
7cf36cae 818
a99cf51f 819 SetHighWaterMark(nt);
7cf36cae 820
821 AliGenEventHeader* header = new AliGenEventHeader("PARAM");
822 header->SetPrimaryVertex(fVertex);
21391258 823 header->SetInteractionTime(fTime);
7cf36cae 824 header->SetNProduced(fNprimaries);
825 AddHeader(header);
fe4da5cc 826}
2ad38d56 827//____________________________________________________________________________________
828Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
829{
6078e216 830 //
831 // Normalisation for selected kinematic region
832 //
faf00b0d 833#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
2ad38d56 834 Float_t ratio =
389a9124 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) *
2ad38d56 837 (phiMax-phiMin)/360.;
faf00b0d 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
2ad38d56 844 return TMath::Abs(ratio);
845}
846
847//____________________________________________________________________________________
fe4da5cc 848
dc1d768c 849void AliGenParam::Draw( const char * /*opt*/)
5bd39445 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
f87cfe57 864
fe4da5cc 865
866