http://savannah.cern.ch/bugs/?103960
[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 {
74   // Default constructor
75 }
76 //____________________________________________________________
77 AliGenParam::AliGenParam(Int_t npart, const AliGenLib * Library,  Int_t param, const char* tname)
78     :AliGenMC(npart),
79      fPtParaFunc(Library->GetPt(param, tname)),
80      fYParaFunc (Library->GetY (param, tname)),
81      fIpParaFunc(Library->GetIp(param, tname)),
82    fV2ParaFunc(Library->GetV2(param, tname)),
83      fPtPara(0),
84      fYPara(0),
85    fV2Para(0),
86    fdNdPhi(0),
87      fParam(param),
88      fdNdy0(0.),
89      fYWgt(0.),
90      fPtWgt(0.),
91      fBias(0.),
92      fTrials(0),
93      fDeltaPt(0.01),
94      fSelectAll(kFALSE),
95    fDecayer(0),
96    fForceConv(kFALSE)
97 {
98   // Constructor using number of particles parameterisation id and library
99     fName = "Param";
100     fTitle= "Particle Generator using pT and y parameterisation";
101     fAnalog = kAnalog;
102     SetForceDecay();
103 }
104 //____________________________________________________________
105 AliGenParam::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),
110   fV2ParaFunc(0),
111     fPtPara(0),
112     fYPara(0),
113   fV2Para(0),
114   fdNdPhi(0),
115     fParam(param),
116     fdNdy0(0.),
117     fYWgt(0.),
118     fPtWgt(0.),
119     fBias(0.),
120     fTrials(0),
121     fDeltaPt(0.01),
122     fSelectAll(kFALSE),
123   fDecayer(0),
124   fForceConv(kFALSE)
125 {
126   // Constructor using parameterisation id and number of particles
127   //
128     fName = name;
129     fTitle= "Particle Generator using pT and y parameterisation";
130       
131     AliGenLib* pLibrary = new AliGenMUONlib();
132     fPtParaFunc = pLibrary->GetPt(param, tname);
133     fYParaFunc  = pLibrary->GetY (param, tname);
134     fIpParaFunc = pLibrary->GetIp(param, tname);
135   fV2ParaFunc = pLibrary->GetV2(param, tname);
136     
137     fAnalog = kAnalog;
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(); 
146 }
147 //____________________________________________________________
148
149 AliGenParam::AliGenParam(Int_t npart, Int_t param,
150                          Double_t (*PtPara) (const Double_t*, const Double_t*),
151                          Double_t (*YPara ) (const Double_t* ,const Double_t*),
152                          Double_t (*V2Para) (const Double_t* ,const Double_t*),
153                          Int_t    (*IpPara) (TRandom *))                 
154     :AliGenMC(npart),
155      
156      fPtParaFunc(PtPara),
157      fYParaFunc(YPara),
158      fIpParaFunc(IpPara),
159    fV2ParaFunc(V2Para),
160      fPtPara(0),
161      fYPara(0),
162    fV2Para(0),
163    fdNdPhi(0),
164      fParam(param),
165      fdNdy0(0.),
166      fYWgt(0.),
167      fPtWgt(0.),
168      fBias(0.),
169      fTrials(0),
170      fDeltaPt(0.01),
171      fSelectAll(kFALSE),
172   fDecayer(0),
173   fForceConv(kFALSE)
174 {
175   // Constructor
176   // Gines Martinez 1/10/99 
177     fName = "Param";
178     fTitle= "Particle Generator using pT and y parameterisation";
179
180     fAnalog = kAnalog;
181     fChildSelect.Set(5);
182     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
183     SetForceDecay();
184     SetCutOnChild();
185     SetChildMomentumRange();
186     SetChildPtRange();
187     SetChildPhiRange();
188     SetChildThetaRange();  
189 }
190
191 //____________________________________________________________
192 AliGenParam::~AliGenParam()
193 {
194   // Destructor
195     delete  fPtPara;
196     delete  fYPara;
197   delete  fV2Para;
198   delete  fdNdPhi;
199 }
200
201 //-------------------------------------------------------------------
202 TVector3 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
218 double 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
225 double 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
232 double 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
260 double 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
270 Int_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
315 //____________________________________________________________
316 void AliGenParam::Init()
317 {
318   // Initialisation
319
320     if (gMC) fDecayer = gMC->GetDecayer();
321   //Begin_Html
322   /*
323     <img src="picts/AliGenParam.gif">
324   */
325   //End_Html
326     char name[256];
327     snprintf(name, 256, "pt-parameterisation for %s", GetName());
328     
329     if (fPtPara) fPtPara->Delete();
330     fPtPara = new TF1(name, fPtParaFunc, fPtMin, fPtMax,0);
331     gROOT->GetListOfFunctions()->Remove(fPtPara);
332   //  Set representation precision to 10 MeV
333     Int_t npx= Int_t((fPtMax - fPtMin) / fDeltaPt);
334     
335     fPtPara->SetNpx(npx);
336
337     snprintf(name, 256, "y-parameterisation  for %s", GetName());
338     if (fYPara) fYPara->Delete();
339     fYPara  = new TF1(name, fYParaFunc, fYMin, fYMax, 0);
340     gROOT->GetListOfFunctions()->Remove(fYPara);
341
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?
355     
356     snprintf(name, 256, "pt-for-%s", GetName());
357     TF1 ptPara(name ,fPtParaFunc, 0, 15, 0);
358     snprintf(name, 256, "y-for-%s", GetName());
359     TF1 yPara(name, fYParaFunc, -6, 6, 0);
360
361   //
362   // dN/dy| y=0
363     Double_t y1=0;
364     Double_t y2=0;
365     
366     fdNdy0=fYParaFunc(&y1,&y2);
367   //
368   // Integral over generation region
369 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
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);
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
378   Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();    //TR: should probably be done differently in case of anisotropic phi...
379     if (fAnalog == kAnalog) {
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     }
388   //
389   // particle decay related initialization
390     fDecayer->SetForceDecay(fForceDecay);
391     fDecayer->Init();
392
393   //
394     AliGenMC::Init();
395 }
396
397 //____________________________________________________________
398 void AliGenParam::Generate()
399 {
400   // 
401   // Generate 1 event (see Generate(Int_t ntimes) for details
402   //
403   GenerateN(1);
404 }
405 //____________________________________________________________
406 void AliGenParam::GenerateN(Int_t ntimes)
407 {
408   //
409   // Generate ntimes*'npart' light and heavy mesons (J/Psi, upsilon or phi, Pion,
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
420   fDecayer->SetForceDecay(fForceDecay);
421   fDecayer->Init();
422
423   //
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)
426   Float_t time0;              // Time0 of the generated parent particle
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], 
430           och[3];             // Momentum, polarisation and origin of the children particles from lujet
431   Double_t ty, xmt;
432   Int_t nt, i, j;
433   Float_t  wgtp, wgtch;
434   Double_t dummy;
435   static TClonesArray *particles;
436   //
437   if(!particles) particles = new TClonesArray("TParticle",1000);
438   
439   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
440   //
441   Float_t random[6];
442  
443   // Calculating vertex position per event
444   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
445   time0 = fTimeOrigin;
446   if(fVertexSmear==kPerEvent) {
447       Vertex();
448       for (j=0;j<3;j++) origin0[j]=fVertex[j];
449       time0 = fTime;
450   }
451   
452   Int_t ipa=0;
453   
454   // Generating fNpart particles
455   fNprimaries = 0;
456   
457   Int_t nGen = fNpart*ntimes;
458   while (ipa<nGen) {
459       while(1) {
460       //
461       // particle type 
462           Int_t iPart = fIpParaFunc(fRandom);
463           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
464           TParticlePDG *particle = pDataBase->GetParticle(iPart);
465           Float_t am = particle->Mass();
466           
467           Rndm(random,2);
468       //
469       // y
470           ty = TMath::TanH(fYPara->GetRandom());
471       //
472       // pT
473           if (fAnalog == kAnalog) {
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);
484           if (TMath::Abs(ty)==1.) {
485               ty=0.;
486               Fatal("AliGenParam", 
487                     "Division by 0: Please check you rapidity range !");
488           }
489       //
490       // phi
491           //      if(!ipa)
492           //phi=fEvPlane; //align first particle of each event with event plane
493           //else{
494         double v2 = fV2Para->Eval(pt);
495         fdNdPhi->SetParameter(0,v2);
496         fdNdPhi->SetParameter(1,fEvPlane);
497         phi=fdNdPhi->GetRandom();
498         //     }
499           
500           pl=xmt*ty/sqrt((1.-ty)*(1.+ty));
501           theta=TMath::ATan2(pt,pl);
502       // Cut on theta
503           if(theta<fThetaMin || theta>fThetaMax) continue;
504           ptot=TMath::Sqrt(pt*pt+pl*pl);
505       // Cut on momentum
506           if(ptot<fPMin || ptot>fPMax) continue;
507       //
508           p[0]=pt*TMath::Cos(phi);
509           p[1]=pt*TMath::Sin(phi);
510           p[2]=pl;
511           if(fVertexSmear==kPerTrack) {
512               Rndm(random,6);
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               }
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]));
522           }
523           
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       //
531           Bool_t decayed = kFALSE;
532           
533
534           if (fForceDecay != kNoDecay) {
535         // Using lujet to decay particle
536               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
537               TLorentzVector pmom(p[0], p[1], p[2], energy);
538               fDecayer->Decay(iPart,&pmom);
539         //
540         // select decay particles
541               Int_t np=fDecayer->ImportParticles(particles);
542
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
551               //  Selecting  GeometryAcceptance for particles fPdgCodeParticleforAcceptanceCut;
552               if (fGeometryAcceptance) 
553                 if (!CheckAcceptanceGeometry(np,particles)) continue;
554               Int_t ncsel=0;
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               
566               if (np >1) {
567                   decayed = kTRUE;
568                   TParticle* iparticle =  0;
569                   Int_t ipF, ipL;
570                   for (i = 1; i<np ; i++) {
571                       trackIt[i] = 1;
572                       iparticle = (TParticle *) particles->At(i);
573                       Int_t kf = iparticle->GetPdgCode();
574                       Int_t ks = iparticle->GetStatusCode();
575             // flagged particle
576
577                       if (pFlag[i] == 1) {
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
584             // flag decay products of particles with long life-time (c tau > .3 mum)                  
585                       
586                       if (ks != 1) { 
587               //                          TParticlePDG *particle = pDataBase->GetParticle(kf);
588                           
589                           Double_t lifeTime = fDecayer->GetLifetime(kf);
590               //                          Double_t mass     = particle->Mass();
591               //                          Double_t width    = particle->Width();
592                           if (lifeTime > (Double_t) fMaxLifeTime) {
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 ?
601             //
602             // children
603                       
604                       if ((ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll || fSelectAll) && trackIt[i])
605                       {
606                           if (fCutOnChild) {
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;
613                                   ncsel++;
614                               } else {
615                                   ncsel=-1;
616                                   break;
617                               } // child kine cuts
618                           } else {
619                               pSelected[i]  = 1;
620                               ncsel++;
621                           } // if child selection
622                       } // select muon
623                   } // decay particle loop
624               } // if decay products
625               
626               Int_t iparent;
627               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
628                   ipa++;
629           //
630           // Parent
631                   
632                   
633                   PushTrack(0, -1, iPart, p, origin0, polar, time0, kPPrimary, nt, wgtp, ((decayed)? 11 : 1));
634                   pParent[0] = nt;
635                   KeepTrack(nt); 
636                   fNprimaries++;
637                   
638           //
639           // Decay Products
640           //              
641                   for (i = 1; i < np; i++) {
642                       if (pSelected[i]) {
643                           TParticle* iparticle = (TParticle *) particles->At(i);
644                           Int_t kf   = iparticle->GetPdgCode();
645                           Int_t ksc  = iparticle->GetStatusCode();
646                           Int_t jpa  = iparticle->GetFirstMother()-1;
647                           
648                           och[0] = origin0[0]+iparticle->Vx();
649                           och[1] = origin0[1]+iparticle->Vy();
650                           och[2] = origin0[2]+iparticle->Vz();
651                           pc[0]  = iparticle->Px();
652                           pc[1]  = iparticle->Py();
653                           pc[2]  = iparticle->Pz();
654                           
655                           if (jpa > -1) {
656                               iparent = pParent[jpa];
657                           } else {
658                               iparent = -1;
659                           }
660                          
661                           PushTrack(fTrackIt * trackIt[i], iparent, kf,
662                                            pc, och, polar,
663                                            time0 + iparticle->T(), kPDecay, nt, wgtch, ksc);
664                           pParent[i] = nt;
665                           KeepTrack(nt); 
666                           fNprimaries++;
667                       } // Selected
668                   } // Particle loop 
669               }  // Decays by Lujet
670               particles->Clear();
671               if (pFlag)      delete[] pFlag;
672               if (pParent)    delete[] pParent;
673               if (pSelected)  delete[] pSelected;          
674               if (trackIt)    delete[] trackIt;
675           } // kinematic selection
676           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
677           {
678             gAlice->GetMCApp()->
679                 PushTrack(fTrackIt,-1,iPart,p,origin0,polar,time0,kPPrimary,nt,wgtp, 1);
680             ipa++; 
681             fNprimaries++;
682           }
683           break;
684     } // while
685   } // event loop
686   
687   SetHighWaterMark(nt);
688
689   AliGenEventHeader* header = new AliGenEventHeader("PARAM");
690   header->SetPrimaryVertex(fVertex);
691   header->SetInteractionTime(fTime);
692   header->SetNProduced(fNprimaries);
693   AddHeader(header);
694 }
695 //____________________________________________________________________________________
696 Float_t AliGenParam::GetRelativeArea(Float_t ptMin, Float_t ptMax, Float_t yMin, Float_t yMax, Float_t phiMin, Float_t phiMax)
697 {
698   //
699   // Normalisation for selected kinematic region
700   //
701 #if ROOT_VERSION_CODE < ROOT_VERSION(5,99,0)
702   Float_t ratio =  
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)   *
705     (phiMax-phiMin)/360.;
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
712   return TMath::Abs(ratio);
713 }
714
715 //____________________________________________________________________________________
716
717 void AliGenParam::Draw( const char * /*opt*/)
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
732
733
734