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