http://savannah.cern.ch/bugs/?103960
[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
218double AliGenParam::ScreenFunc1(double d){
219 if(d>1)
220 return 21.12-4.184*log(d+0.952);
221 else
222 return 20.867-3.242*d+0.652*d*d;
223}
224
225double AliGenParam::ScreenFunc2(double d){
226 if(d>1)
227 return 21.12-4.184*log(d+0.952);
228 else
229 return 20.209-1.93*d-0.086*d*d;
230}
231
232double AliGenParam::EnergyFraction(double Z, double E){
233 double e0=0.000511/E;
234 double aZ=Z/137.036;
235 double dmin=ScreenVar(Z,e0,0.5);
236 double fcZ=(aZ*aZ)*(1/(1+aZ*aZ)+0.20206-0.0368*aZ*aZ+0.0083*aZ*aZ*aZ);
237 double Fz=8*log(Z)/3;
238 if(E>0.05)
239 Fz+=8*fcZ;
240 double dmax=exp((42.24-Fz)/8.368)-0.952;
241 if(!(dmax>dmin))
242 return fRandom->Uniform(e0,0.5);
243
244 double e1=0.5-0.5*sqrt(1-dmin/dmax);
245 double emin=TMath::Max(e0,e1);
246
247 double normval=1/(0.5*(ScreenFunc1(dmin)-0.5*Fz)+0.1666667*(ScreenFunc2(dmin)-0.5*Fz));
248 while(true){
249 double y=fRandom->Uniform(emin,1);
250 double eps=y*(0.38453+y*(0.10234+y*(0.026072+y*(0.014367-0.027313*y)))); //inverse to the enveloping cumulative probability distribution
251 double val=fRandom->Uniform(0,2.12*(eps*eps-eps)+1.53); //enveloping probability density
252 double d=ScreenVar(Z,e0,eps);
253 double bh=((eps*eps)+(1-eps)*(1-eps))*(ScreenFunc1(d)-0.5*Fz)+0.6666667*eps*(1-eps)*(ScreenFunc2(d)-0.5*Fz);
254 bh*=normval;
255 if(val<bh)
256 return eps;
257 }
258}
259
260double AliGenParam::PolarAngle(double E){
261 float rand[3];
262 AliRndm rndm;
263 rndm.Rndm(rand,3);
264 double u=-8*log(rand[1]*rand[2])/5;
265 if(!(rand[0]<9.0/36))
266 u/=3;
267 return u*0.000511/E;
268}
269
270Int_t AliGenParam::ForceGammaConversion(TClonesArray *particles, Int_t nPart)
271{
272 //based on: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node27.html
273 // and: http://geant4.cern.ch/G4UsersDocuments/UsersGuides/PhysicsReferenceManual/html/node58.html
274 Int_t nPartNew=nPart;
275 for(int iPart=0; iPart<nPart; iPart++){
276 TParticle *gamma = (TParticle *) particles->At(iPart);
277 if(gamma->GetPdgCode()!=22) continue;
278
279 TVector3 gammaV3(gamma->Px(),gamma->Py(),gamma->Pz());
280 Float_t az=fRandom->Uniform(TMath::Pi()*2);
281 double frac=EnergyFraction(1,gamma->Energy());
282 double Ee1=frac*gamma->Energy();
283 double Ee2=(1-frac)*gamma->Energy();
284 double Pe1=Ee1;//sqrt(Ee1*Ee1-0.000511*0.000511);
285 double Pe2=Ee2;//sqrt(Ee2*Ee2-0.000511*0.000511);
286
287 TVector3 rotAxis(OrthogonalVector(gammaV3));
288 rotAxis.Rotate(az,gammaV3);
289 TVector3 e1V3(gammaV3);
290 e1V3.Rotate(PolarAngle(Ee1),rotAxis);
291 e1V3=e1V3.Unit();
292 e1V3*=Pe1;
293 TVector3 e2V3(gammaV3);
294 e2V3.Rotate(-PolarAngle(Ee2),rotAxis);
295 e2V3=e2V3.Unit();
296 e2V3*=Pe2;
297 // gamma = new TParticle(*gamma);
298 // particles->RemoveAt(iPart);
299 gamma->SetFirstDaughter(nPartNew+1);
300 gamma->SetLastDaughter(nPartNew+2);
301 // new((*particles)[iPart]) TParticle(*gamma);
302 // delete gamma;
303
304 TLorentzVector vtx;
305 gamma->ProductionVertex(vtx);
306 new((*particles)[nPartNew]) TParticle(11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e1V3,Ee1), vtx);
307 nPartNew++;
308 new((*particles)[nPartNew]) TParticle(-11, gamma->GetStatusCode(), iPart+1, -1, 0, 0, TLorentzVector(e2V3,Ee2), vtx);
309 nPartNew++;
310 }
311 // particles->Compress();
312 return particles->GetEntriesFast();
313}
314
fe4da5cc 315//____________________________________________________________
316void AliGenParam::Init()
317{
6078e216 318 // Initialisation
fa480368 319
2904363f 320 if (gMC) fDecayer = gMC->GetDecayer();
fe4da5cc 321 //Begin_Html
322 /*
1439f98e 323 <img src="picts/AliGenParam.gif">
fe4da5cc 324 */
325 //End_Html
c5c3e4b0 326 char name[256];
8f8399ab 327 snprintf(name, 256, "pt-parameterisation for %s", GetName());
c5c3e4b0 328
2ad38d56 329 if (fPtPara) fPtPara->Delete();
c5c3e4b0 330 fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
828fff0d 331 gROOT->GetListOfFunctions()->Remove(fPtPara);
6078e216 332 // Set representation precision to 10 MeV
c5c3e4b0 333 Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
749070b6 334
335 fPtPara->SetNpx(npx);
c5c3e4b0 336
8f8399ab 337 snprintf(name, 256, "y-parameterisation for %s", GetName());
2ad38d56 338 if (fYPara) fYPara->Delete();
c5c3e4b0 339 fYPara = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
828fff0d 340 gROOT->GetListOfFunctions()->Remove(fYPara);
341
6078e216 342 snprintf(name, 256, "v2-parameterisation for %s", GetName());
343 if (fV2Para) fV2Para->Delete();
344 fV2Para = new TF1(name, fV2ParaFunc, fPtMin, fPtMax, 0);
345 // fV2Para = new TF1(name, "2*[0]/(1+TMath::Exp([1]*([2]-x)))-[0]", fPtMin, fPtMax);
346 // fV2Para->SetParameter(0, 0.236910);
347 // fV2Para->SetParameter(1, 1.71122);
348 // fV2Para->SetParameter(2, 0.0827617);
349 //gROOT->GetListOfFunctions()->Remove(fV2Para); //TR: necessary?
350
351 snprintf(name, 256, "dNdPhi for %s", GetName());
352 if (fdNdPhi) fdNdPhi->Delete();
353 fdNdPhi = new TF1(name, "1+2*[0]*TMath::Cos(2*(x-[1]))", fPhiMin, fPhiMax);
354 //gROOT->GetListOfFunctions()->Remove(fdNdPhi); //TR: necessary?
749070b6 355
8f8399ab 356 snprintf(name, 256, "pt-for-%s", GetName());
c5c3e4b0 357 TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
8f8399ab 358 snprintf(name, 256, "y-for-%s", GetName());
c5c3e4b0 359 TF1 yPara(name, fYParaFunc, -6, 6, 0);
b7601ac4 360
6078e216 361 //
362 // dN/dy| y=0
749070b6 363 Double_t y1=0;
364 Double_t y2=0;
365
366 fdNdy0=fYParaFunc(&y1,&y2);
6078e216 367 //
368 // Integral over generation region
faf00b0d 369#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
389a9124 370 Float_t intYS = yPara.Integral(fYMin, fYMax,(Double_t*) 0x0,1.e-6);
371 Float_t intPt0 = ptPara.Integral(0,15,(Double_t *) 0x0,1.e-6);
372 Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,(Double_t*) 0x0,1.e-6);
faf00b0d 373#else
374 Float_t intYS = yPara.Integral(fYMin, fYMax,1.e-6);
375 Float_t intPt0 = ptPara.Integral(0,15,1.e-6);
376 Float_t intPtS = ptPara.Integral(fPtMin,fPtMax,1.e-6);
377#endif
6078e216 378 Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi(); //TR: should probably be done differently in case of anisotropic phi...
34f60c01 379 if (fAnalog == kAnalog) {
749070b6 380 fYWgt = intYS/fdNdy0;
381 fPtWgt = intPtS/intPt0;
382 fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
383 } else {
384 fYWgt = intYS/fdNdy0;
385 fPtWgt = (fPtMax-fPtMin)/intPt0;
386 fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
387 }
6078e216 388 //
389 // particle decay related initialization
00d6ce7d 390 fDecayer->SetForceDecay(fForceDecay);
fa480368 391 fDecayer->Init();
00d6ce7d 392
6078e216 393 //
fff02fee 394 AliGenMC::Init();
fe4da5cc 395}
396
397//____________________________________________________________
398void AliGenParam::Generate()
399{
d5615b85 400 //
401 // Generate 1 event (see Generate(Int_t ntimes) for details
6078e216 402 //
d5615b85 403 GenerateN(1);
404}
405//____________________________________________________________
406void AliGenParam::GenerateN(Int_t ntimes)
407{
408 //
409 // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion,
6078e216 410 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
411 // antineutrons in the the desired theta, phi and momentum windows;
412 // Gaussian smearing on the vertex is done if selected.
413 // The decay of heavy mesons is done using lujet,
414 // and the childern particle are tracked by GEANT
415 // However, light mesons are directly tracked by GEANT
416 // setting fForceDecay = nodecay (SetForceDecay(nodecay))
417 //
418 //
419 // Reinitialize decayer
2ad38d56 420 fDecayer->SetForceDecay(fForceDecay);
421 fDecayer->Init();
422
6078e216 423 //
28337bc1 424 Float_t polar[3]= {0,0,0}; // Polarisation of the parent particle (for GEANT tracking)
425 Float_t origin0[3]; // Origin of the generated parent particle (for GEANT tracking)
21391258 426 Float_t time0; // Time0 of the generated parent particle
28337bc1 427 Float_t pt, pl, ptot; // Transverse, logitudinal and total momenta of the parent particle
428 Float_t phi, theta; // Phi and theta spherical angles of the parent particle momentum
429 Float_t p[3], pc[3],
fff02fee 430 och[3]; // Momentum, polarisation and origin of the children particles from lujet
2d7a47be 431 Double_t ty, xmt;
fff02fee 432 Int_t nt, i, j;
fe4da5cc 433 Float_t wgtp, wgtch;
434 Double_t dummy;
09fd3ea2 435 static TClonesArray *particles;
fe4da5cc 436 //
fff02fee 437 if(!particles) particles = new TClonesArray("TParticle",1000);
8b31bfac 438
349be858 439 TDatabasePDG *pDataBase = TDatabasePDG::Instance();
fe4da5cc 440 //
441 Float_t random[6];
28337bc1 442
6078e216 443 // Calculating vertex position per event
fe4da5cc 444 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
21391258 445 time0 = fTimeOrigin;
aee8290b 446 if(fVertexSmear==kPerEvent) {
d9ea0e3b 447 Vertex();
448 for (j=0;j<3;j++) origin0[j]=fVertex[j];
21391258 449 time0 = fTime;
fe4da5cc 450 }
d9ea0e3b 451
21aaa175 452 Int_t ipa=0;
fff02fee 453
6078e216 454 // Generating fNpart particles
7cf36cae 455 fNprimaries = 0;
456
d5615b85 457 Int_t nGen = fNpart*ntimes;
458 while (ipa<nGen) {
cc5d764c 459 while(1) {
6078e216 460 //
461 // particle type
65fb704d 462 Int_t iPart = fIpParaFunc(fRandom);
fa480368 463 fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;
675e9664 464 TParticlePDG *particle = pDataBase->GetParticle(iPart);
fa480368 465 Float_t am = particle->Mass();
8b31bfac 466
65fb704d 467 Rndm(random,2);
6078e216 468 //
469 // y
2d7a47be 470 ty = TMath::TanH(fYPara->GetRandom());
6078e216 471 //
472 // pT
34f60c01 473 if (fAnalog == kAnalog) {
fe4da5cc 474 pt=fPtPara->GetRandom();
475 wgtp=fParentWeight;
476 wgtch=fChildWeight;
477 } else {
478 pt=fPtMin+random[1]*(fPtMax-fPtMin);
479 Double_t ptd=pt;
480 wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
481 wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
482 }
483 xmt=sqrt(pt*pt+am*am);
2d7a47be 484 if (TMath::Abs(ty)==1.) {
485 ty=0.;
486 Fatal("AliGenParam",
487 "Division by 0: Please check you rapidity range !");
488 }
6078e216 489 //
490 // phi
4ae1c9f0 491 // if(!ipa)
492 //phi=fEvPlane; //align first particle of each event with event plane
493 //else{
6078e216 494 double v2 = fV2Para->Eval(pt);
495 fdNdPhi->SetParameter(0,v2);
496 fdNdPhi->SetParameter(1,fEvPlane);
497 phi=fdNdPhi->GetRandom();
4ae1c9f0 498 // }
2d7a47be 499
60e55aee 500 pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
fe4da5cc 501 theta=TMath::ATan2(pt,pl);
6078e216 502 // Cut on theta
fe4da5cc 503 if(theta<fThetaMin || theta>fThetaMax) continue;
504 ptot=TMath::Sqrt(pt*pt+pl*pl);
6078e216 505 // Cut on momentum
fe4da5cc 506 if(ptot<fPMin || ptot>fPMax) continue;
6078e216 507 //
fe4da5cc 508 p[0]=pt*TMath::Cos(phi);
509 p[1]=pt*TMath::Sin(phi);
510 p[2]=pl;
aee8290b 511 if(fVertexSmear==kPerTrack) {
65fb704d 512 Rndm(random,6);
fe4da5cc 513 for (j=0;j<3;j++) {
514 origin0[j]=
515 fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
516 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
517 }
21391258 518 Rndm(random,2);
519 time0 = fTimeOrigin + fOsigma[2]/TMath::Ccgs()*
520 TMath::Cos(2*random[0]*TMath::Pi())*
521 TMath::Sqrt(-2*TMath::Log(random[1]));
fe4da5cc 522 }
28337bc1 523
6078e216 524 // Looking at fForceDecay :
525 // if fForceDecay != none Primary particle decays using
526 // AliPythia and children are tracked by GEANT
527 //
528 // if fForceDecay == none Primary particle is tracked by GEANT
529 // (In the latest, make sure that GEANT actually does all the decays you want)
530 //
cd716030 531 Bool_t decayed = kFALSE;
532
fff02fee 533
34f60c01 534 if (fForceDecay != kNoDecay) {
6078e216 535 // Using lujet to decay particle
886b6f73 536 Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
fa480368 537 TLorentzVector pmom(p[0], p[1], p[2], energy);
538 fDecayer->Decay(iPart,&pmom);
6078e216 539 //
540 // select decay particles
fa480368 541 Int_t np=fDecayer->ImportParticles(particles);
1242532d 542
4ae1c9f0 543 if(fForceConv) np=ForceGammaConversion(particles,np);
544
545 // for(int iPart=0; iPart<np; iPart++){
546 // TParticle *gamma = (TParticle *) particles->At(iPart);
547 // printf("%i %i:", iPart, gamma->GetPdgCode());
548 // printf("%i %i %i|",gamma->GetFirstMother(),gamma->GetFirstDaughter(),gamma->GetLastDaughter());
549 // }
550
1242532d 551 // Selecting GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
552 if (fGeometryAcceptance)
553 if (!CheckAcceptanceGeometry(np,particles)) continue;
886b6f73 554 Int_t ncsel=0;
fff02fee 555 Int_t* pFlag = new Int_t[np];
556 Int_t* pParent = new Int_t[np];
557 Int_t* pSelected = new Int_t[np];
558 Int_t* trackIt = new Int_t[np];
559
560 for (i=0; i<np; i++) {
561 pFlag[i] = 0;
562 pSelected[i] = 0;
563 pParent[i] = -1;
564 }
565
6ba00c52 566 if (np >1) {
cd716030 567 decayed = kTRUE;
08bffa4d 568 TParticle* iparticle = 0;
fff02fee 569 Int_t ipF, ipL;
570 for (i = 1; i<np ; i++) {
571 trackIt[i] = 1;
6ba00c52 572 iparticle = (TParticle *) particles->At(i);
573 Int_t kf = iparticle->GetPdgCode();
fff02fee 574 Int_t ks = iparticle->GetStatusCode();
6078e216 575 // flagged particle
fff02fee 576
577 if (pFlag[i] == 1) {
fff02fee 578 ipF = iparticle->GetFirstDaughter();
579 ipL = iparticle->GetLastDaughter();
580 if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
581 continue;
582 }
583
6078e216 584 // flag decay products of particles with long life-time (c tau > .3 mum)
fff02fee 585
586 if (ks != 1) {
6078e216 587 // TParticlePDG *particle = pDataBase->GetParticle(kf);
fff02fee 588
589 Double_t lifeTime = fDecayer->GetLifetime(kf);
6078e216 590 // Double_t mass = particle->Mass();
591 // Double_t width = particle->Width();
47fc6bd5 592 if (lifeTime > (Double_t) fMaxLifeTime) {
fff02fee 593 ipF = iparticle->GetFirstDaughter();
594 ipL = iparticle->GetLastDaughter();
595 if (ipF > 0) for (j=ipF-1; j<ipL; j++) pFlag[j]=1;
596 } else{
597 trackIt[i] = 0;
598 pSelected[i] = 1;
599 }
600 } // ks==1 ?
6078e216 601 //
602 // children
47fc6bd5 603
18e09c20 604 if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
6ba00c52 605 {
6ba00c52 606 if (fCutOnChild) {
fff02fee 607 pc[0]=iparticle->Px();
608 pc[1]=iparticle->Py();
609 pc[2]=iparticle->Pz();
610 Bool_t childok = KinematicSelection(iparticle, 1);
611 if(childok) {
612 pSelected[i] = 1;
6ba00c52 613 ncsel++;
614 } else {
615 ncsel=-1;
616 break;
617 } // child kine cuts
618 } else {
fff02fee 619 pSelected[i] = 1;
886b6f73 620 ncsel++;
6ba00c52 621 } // if child selection
622 } // select muon
623 } // decay particle loop
624 } // if decay products
625
886b6f73 626 Int_t iparent;
627 if ((fCutOnChild && ncsel >0) || !fCutOnChild){
628 ipa++;
6078e216 629 //
630 // Parent
cd716030 631
632
21391258 633 PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
fff02fee 634 pParent[0] = nt;
a99cf51f 635 KeepTrack(nt);
7cf36cae 636 fNprimaries++;
637
6078e216 638 //
639 // Decay Products
640 //
fff02fee 641 for (i = 1; i < np; i++) {
642 if (pSelected[i]) {
643 TParticle* iparticle = (TParticle *) particles->At(i);
ea79897e 644 Int_t kf = iparticle->GetPdgCode();
cd716030 645 Int_t ksc = iparticle->GetStatusCode();
ea79897e 646 Int_t jpa = iparticle->GetFirstMother()-1;
fff02fee 647
774ceaaf 648 och[0] = origin0[0]+iparticle->Vx();
649 och[1] = origin0[1]+iparticle->Vy();
650 och[2] = origin0[2]+iparticle->Vz();
fff02fee 651 pc[0] = iparticle->Px();
652 pc[1] = iparticle->Py();
653 pc[2] = iparticle->Pz();
654
ea79897e 655 if (jpa > -1) {
656 iparent = pParent[jpa];
fff02fee 657 } else {
658 iparent = -1;
659 }
1242532d 660
ea79897e 661 PushTrack(fTrackIt * trackIt[i], iparent, kf,
fff02fee 662 pc, och, polar,
21391258 663 time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
fff02fee 664 pParent[i] = nt;
a99cf51f 665 KeepTrack(nt);
7cf36cae 666 fNprimaries++;
8b31bfac 667 } // Selected
668 } // Particle loop
28337bc1 669 } // Decays by Lujet
8b31bfac 670 particles->Clear();
fff02fee 671 if (pFlag) delete[] pFlag;
672 if (pParent) delete[] pParent;
673 if (pSelected) delete[] pSelected;
8b31bfac 674 if (trackIt) delete[] trackIt;
21aaa175 675 } // kinematic selection
28337bc1 676 else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
677 {
5d12ce38 678 gAlice->GetMCApp()->
21391258 679 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
28337bc1 680 ipa++;
7cf36cae 681 fNprimaries++;
28337bc1 682 }
fe4da5cc 683 break;
21aaa175 684 } // while
fe4da5cc 685 } // event loop
7cf36cae 686
a99cf51f 687 SetHighWaterMark(nt);
7cf36cae 688
689 AliGenEventHeader* header = new AliGenEventHeader("PARAM");
690 header->SetPrimaryVertex(fVertex);
21391258 691 header->SetInteractionTime(fTime);
7cf36cae 692 header->SetNProduced(fNprimaries);
693 AddHeader(header);
fe4da5cc 694}
2ad38d56 695//____________________________________________________________________________________
696Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
697{
6078e216 698 //
699 // Normalisation for selected kinematic region
700 //
faf00b0d 701#if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
2ad38d56 702 Float_t ratio =
389a9124 703 fPtPara->Integral(ptMin,ptMax,(Double_t *)0,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),(Double_t *)0,1.e-6) *
704 fYPara->Integral(yMin,yMax,(Double_t *)0,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),(Double_t *)0,1.e-6) *
2ad38d56 705 (phiMax-phiMin)/360.;
faf00b0d 706#else
707 Float_t ratio =
708 fPtPara->Integral(ptMin,ptMax,1.e-6) / fPtPara->Integral( fPtPara->GetXmin(), fPtPara->GetXmax(),1.e-6) *
709 fYPara->Integral(yMin,yMax,1.e-6)/fYPara->Integral(fYPara->GetXmin(),fYPara->GetXmax(),1.e-6) *
710 (phiMax-phiMin)/360.;
711#endif
2ad38d56 712 return TMath::Abs(ratio);
713}
714
715//____________________________________________________________________________________
fe4da5cc 716
dc1d768c 717void AliGenParam::Draw( const char * /*opt*/)
5bd39445 718{
719 //
720 // Draw the pT and y Distributions
721 //
722 TCanvas *c0 = new TCanvas("c0","Canvas 0",400,10,600,700);
723 c0->Divide(2,1);
724 c0->cd(1);
725 fPtPara->Draw();
726 fPtPara->GetHistogram()->SetXTitle("p_{T} (GeV)");
727 c0->cd(2);
728 fYPara->Draw();
729 fYPara->GetHistogram()->SetXTitle("y");
730}
731
f87cfe57 732
fe4da5cc 733
734