Put only first generation decay products on the stack.
[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 /*
17 $Log$
18 Revision 1.29  2001/03/27 10:58:41  morsch
19 Initialize decayer before generation. Important if run inside cocktail.
20
21 Revision 1.28  2001/03/09 13:01:41  morsch
22 - enum constants for paramterisation type (particle family) moved to AliGen*lib.h
23 - use AliGenGSIlib::kUpsilon, AliGenPHOSlib::kEtaPrime to access the constants
24
25 Revision 1.27  2001/02/02 15:21:10  morsch
26 Set high water mark after last particle.
27 Use Vertex() method for Vertex.
28
29 Revision 1.26  2000/12/21 16:24:06  morsch
30 Coding convention clean-up
31
32 Revision 1.25  2000/11/30 07:12:50  alibrary
33 Introducing new Rndm and QA classes
34
35 Revision 1.24  2000/10/18 19:11:27  hristov
36 Division by zero fixed
37
38 Revision 1.23  2000/10/02 21:28:06  fca
39 Removal of useless dependecies via forward declarations
40
41 Revision 1.22  2000/09/12 14:14:55  morsch
42 Call fDecayer->ForceDecay() at the beginning of Generate().
43
44 Revision 1.21  2000/09/08 15:39:01  morsch
45 Handle the case fForceDecay=all during the generation, i.e. select all secondaries.
46
47 Revision 1.20  2000/09/06 14:35:44  morsch
48 Use AliDecayerPythia for particle decays.
49
50 Revision 1.19  2000/07/11 18:24:56  fca
51 Coding convention corrections + few minor bug fixes
52
53 Revision 1.18  2000/06/29 21:08:27  morsch
54 All paramatrisation libraries derive from the pure virtual base class AliGenLib.
55 This allows to pass a pointer to a library directly to AliGenParam and avoids the
56 use of function pointers in Config.C.
57
58 Revision 1.17  2000/06/09 20:33:30  morsch
59 All coding rule violations except RS3 corrected
60
61 Revision 1.16  2000/05/02 07:51:31  morsch
62 - Control precision of pT sampling TF1::SetNpx(..)
63 - Correct initialisation of child-cuts in all constructors.
64 - Most coding rule violations corrected.
65
66 Revision 1.15  2000/04/03 15:42:12  morsch
67 Cuts on primary particles are separated from those on the decay products. Methods
68 SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.
69
70 Revision 1.14  1999/11/09 07:38:48  fca
71 Changes for compatibility with version 2.23 of ROOT
72
73 Revision 1.13  1999/11/04 11:30:31  fca
74 Correct the logics for SetForceDecay
75
76 Revision 1.12  1999/11/03 17:43:20  fca
77 New version from G.Martinez & A.Morsch
78
79 Revision 1.11  1999/09/29 09:24:14  fca
80 Introduction of the Copyright and cvs Log
81
82 */
83
84 #include "AliGenParam.h"
85 #include "AliDecayerPythia.h"
86 #include "AliGenMUONlib.h"
87 #include "AliRun.h"
88 #include <TParticle.h>
89 #include <TParticlePDG.h>
90 #include <TDatabasePDG.h>
91 #include <TLorentzVector.h>
92
93 #include <TF1.h>
94
95 ClassImp(AliGenParam)
96
97 //------------------------------------------------------------
98
99   //Begin_Html
100   /*
101     <img src="picts/AliGenParam.gif">
102   */
103   //End_Html
104
105 //____________________________________________________________
106 //____________________________________________________________
107 AliGenParam::AliGenParam()
108 {
109 // Deafault constructor
110     fPtPara = 0;
111     fYPara  = 0;
112     fParam  = 0;
113     fAnalog = kAnalog;
114     SetCutOnChild();
115     SetChildMomentumRange();
116     SetChildPtRange();
117     SetChildPhiRange();
118     SetChildThetaRange();  
119     SetDeltaPt();
120 //
121 // Set random number generator   
122     sRandom = fRandom;
123 }
124
125 AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Int_t param, char* tname):AliGenerator(npart)
126 {
127 // Constructor using number of particles parameterisation id and library
128     
129     fPtParaFunc = Library->GetPt(param, tname);
130     fYParaFunc  = Library->GetY (param, tname);
131     fIpParaFunc = Library->GetIp(param, tname);
132     
133     fPtPara = 0;
134     fYPara  = 0;
135     fParam  = param;
136     fAnalog = kAnalog;
137     fChildSelect.Set(5);
138     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
139     SetForceDecay();
140     SetCutOnChild();
141     SetChildMomentumRange();
142     SetChildPtRange();
143     SetChildPhiRange();
144     SetChildThetaRange(); 
145     SetDeltaPt(); 
146 //
147 // Set random number generator   
148     sRandom = fRandom;
149 }
150
151 //____________________________________________________________
152
153 AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenerator(npart)
154 {
155 // Constructor using parameterisation id and number of particles
156 //      
157     AliGenLib* pLibrary = new AliGenMUONlib();
158  
159     fPtParaFunc = pLibrary->GetPt(param, tname);
160     fYParaFunc  = pLibrary->GetY (param, tname);
161     fIpParaFunc = pLibrary->GetIp(param, tname);
162     
163     fPtPara = 0;
164     fYPara  = 0;
165     fParam  = param;
166     fAnalog = kAnalog;
167     fChildSelect.Set(5);
168     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
169     SetForceDecay();
170     SetCutOnChild();
171     SetChildMomentumRange();
172     SetChildPtRange();
173     SetChildPhiRange();
174     SetChildThetaRange(); 
175     SetDeltaPt(); 
176 }
177
178 AliGenParam::AliGenParam(Int_t npart, Int_t param,
179                          Double_t (*PtPara) (Double_t*, Double_t*),
180                          Double_t (*YPara ) (Double_t* ,Double_t*),
181                          Int_t    (*IpPara) (TRandom *))                 
182     :AliGenerator(npart)
183 {
184 // Constructor
185 // Gines Martinez 1/10/99 
186     fPtParaFunc = PtPara; 
187     fYParaFunc  = YPara;  
188     fIpParaFunc = IpPara;
189 //  
190     fPtPara = 0;
191     fYPara  = 0;
192     fParam  = param;
193     fAnalog = kAnalog;
194     fChildSelect.Set(5);
195     for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
196     SetForceDecay();
197     SetCutOnChild();
198     SetChildMomentumRange();
199     SetChildPtRange();
200     SetChildPhiRange();
201     SetChildThetaRange();  
202     SetDeltaPt();
203 }
204
205
206 AliGenParam::AliGenParam(const AliGenParam & Paramd)
207 {
208 // copy constructor
209 }
210
211 //____________________________________________________________
212 AliGenParam::~AliGenParam()
213 {
214 // Destructor
215     delete  fPtPara;
216     delete  fYPara;
217 }
218
219 //____________________________________________________________
220 void AliGenParam::Init()
221 {
222 // Initialisation
223
224     fDecayer = new AliDecayerPythia();
225   //Begin_Html
226   /*
227     <img src="picts/AliGenParam.gif">
228   */
229   //End_Html
230  
231     fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
232 //  Set representation precision to 10 MeV
233     Int_t npx= Int_t((fPtMax-fPtMin)/fDeltaPt);
234     
235     fPtPara->SetNpx(npx);
236     
237     fYPara  = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
238     TF1* ptPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
239     TF1* yPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
240
241 //
242 // dN/dy| y=0
243     Double_t y1=0;
244     Double_t y2=0;
245     
246     fdNdy0=fYParaFunc(&y1,&y2);
247 //
248 // Integral over generation region
249     Float_t intYS  = yPara ->Integral(fYMin, fYMax);
250     Float_t intPt0 = ptPara->Integral(0,15);
251     Float_t intPtS = ptPara->Integral(fPtMin,fPtMax);
252     Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
253     if (fAnalog == kAnalog) {
254         fYWgt  = intYS/fdNdy0;
255         fPtWgt = intPtS/intPt0;
256         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
257     } else {
258         fYWgt = intYS/fdNdy0;
259         fPtWgt = (fPtMax-fPtMin)/intPt0;
260         fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
261     }
262 //
263 // particle decay related initialization
264     fDecayer->SetForceDecay(fForceDecay);
265     fDecayer->Init();
266
267 //
268     switch (fForceDecay) 
269     {
270     case kSemiElectronic:
271     case kDiElectron:
272     case kBJpsiDiElectron:
273     case kBPsiPrimeDiElectron:
274         fChildSelect[0]=11;     
275         break;
276     case kSemiMuonic:
277     case kDiMuon:
278     case kBJpsiDiMuon:
279     case kBPsiPrimeDiMuon:
280         fChildSelect[0]=13;
281         break;
282     case kPiToMu:
283         fChildSelect[0]=13;
284         break;
285     case kKaToMu:
286         fChildSelect[0]=13;
287         break;
288     case kHadronicD:
289 // Implement me !!
290         break;
291     case kNoDecay:
292         break;
293     case kAll:
294         break;
295     }
296 }
297
298 //____________________________________________________________
299 void AliGenParam::Generate()
300 {
301 //
302 // Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
303 // Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
304 // antineutrons in the the desired theta, phi and momentum windows; 
305 // Gaussian smearing on the vertex is done if selected. 
306 // The decay of heavy mesons is done using lujet, 
307 //    and the childern particle are tracked by GEANT
308 // However, light mesons are directly tracked by GEANT 
309 // setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
310 //
311 //
312 //  Reinitialize decayer
313     fDecayer->Init();
314 //
315   Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
316   Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
317   Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
318   Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
319   Float_t p[3], pc[3], 
320           och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
321   Float_t ty, xmt;
322   Int_t nt, i, j, kfch[10];
323   Float_t  wgtp, wgtch;
324   Double_t dummy;
325   static TClonesArray *particles;
326   //
327   if(!particles) particles=new TClonesArray("TParticle",1000);
328
329   static TDatabasePDG *pDataBase = new TDatabasePDG();
330   if(!pDataBase) pDataBase = new TDatabasePDG();
331   //
332   Float_t random[6];
333  
334 // Calculating vertex position per event
335   for (j=0;j<3;j++) origin0[j]=fOrigin[j];
336   if(fVertexSmear==kPerEvent) {
337 //      Rndm(random,6);
338 //      for (j=0;j<3;j++) {
339 //        origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
340 //            TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
341 //            TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
342 //      }
343 //    }
344       Vertex();
345       for (j=0;j<3;j++) origin0[j]=fVertex[j];
346   }
347   
348   Int_t ipa=0;
349 // Generating fNpart particles
350   while (ipa<fNpart) {
351       while(1) {
352 //
353 // particle type
354           Int_t iPart = fIpParaFunc(fRandom);
355           fChildWeight=(fDecayer->GetPartialBranchingRatio(iPart))*fParentWeight;          
356           TParticlePDG *particle = pDataBase->GetParticle(iPart);
357           Float_t am = particle->Mass();
358
359           Rndm(random,2);
360 //
361 // phi
362           phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
363 //
364 // y
365           ty=Float_t(TMath::TanH(fYPara->GetRandom()));
366 //
367 // pT
368           if (fAnalog == kAnalog) {
369               pt=fPtPara->GetRandom();
370               wgtp=fParentWeight;
371               wgtch=fChildWeight;
372           } else {
373               pt=fPtMin+random[1]*(fPtMax-fPtMin);
374               Double_t ptd=pt;
375               wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
376               wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
377           }
378           xmt=sqrt(pt*pt+am*am);
379       if (TMath::Abs(ty)==1) ty=0;
380           pl=xmt*ty/sqrt(1.-ty*ty);
381           theta=TMath::ATan2(pt,pl);
382 // Cut on theta
383           if(theta<fThetaMin || theta>fThetaMax) continue;
384           ptot=TMath::Sqrt(pt*pt+pl*pl);
385 // Cut on momentum
386           if(ptot<fPMin || ptot>fPMax) continue;
387           p[0]=pt*TMath::Cos(phi);
388           p[1]=pt*TMath::Sin(phi);
389           p[2]=pl;
390           if(fVertexSmear==kPerTrack) {
391               Rndm(random,6);
392               for (j=0;j<3;j++) {
393                   origin0[j]=
394                       fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
395                       TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
396               }
397           }
398           
399 // Looking at fForceDecay : 
400 // if fForceDecay != none Primary particle decays using 
401 // AliPythia and children are tracked by GEANT
402 //
403 // if fForceDecay == none Primary particle is tracked by GEANT 
404 // (In the latest, make sure that GEANT actually does all the decays you want)    
405 //
406           if (fForceDecay != kNoDecay) {
407 // Using lujet to decay particle
408               Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
409               TLorentzVector pmom(p[0], p[1], p[2], energy);
410               fDecayer->Decay(iPart,&pmom);
411 //
412 // select decay particles
413               Int_t np=fDecayer->ImportParticles(particles);
414               Int_t ncsel=0;
415
416               if (np >1) {
417                   TParticle* iparticle =  (TParticle *) particles->At(0);
418                   Int_t ipF = iparticle->GetFirstDaughter();
419                   Int_t ipL = iparticle->GetLastDaughter();           
420                   for (i = ipF-1; i<ipL; i++) {
421                       iparticle = (TParticle *) particles->At(i);
422                       Int_t kf = iparticle->GetPdgCode();
423 //
424 // children
425                       if (ChildSelected(TMath::Abs(kf)) || fForceDecay == kAll)
426                       {
427                           pc[0]=iparticle->Px();
428                           pc[1]=iparticle->Py();
429                           pc[2]=iparticle->Pz();
430                           och[0]=origin0[0]+iparticle->Vx()/10;
431                           och[1]=origin0[1]+iparticle->Vy()/10;
432                           och[2]=origin0[2]+iparticle->Vz()/10;
433                           if (fCutOnChild) {
434                               Float_t ptChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
435                               Float_t pChild=TMath::Sqrt(ptChild*ptChild+pc[2]*pc[2]);
436                               Float_t thetaChild=TMath::ATan2(ptChild,pc[2]);
437                               Float_t phiChild=TMath::ATan2(pc[1],pc[0]);
438                               Bool_t childok = 
439                                   ((ptChild    > fChildPtMin    && ptChild    <fChildPtMax)      &&
440                                    (pChild     > fChildPMin     && pChild     <fChildPMax)       &&
441                                    (thetaChild > fChildThetaMin && thetaChild <fChildThetaMax)   &&
442                                    (phiChild   > fChildPhiMin   && phiChild   <fChildPhiMax));
443                               if(childok)
444                               {
445                                   pch[ncsel][0]=pc[0];
446                                   pch[ncsel][1]=pc[1];
447                                   pch[ncsel][2]=pc[2];
448                                   kfch[ncsel]=kf;
449                                   ncsel++;
450                               } else {
451                                   ncsel=-1;
452                                   break;
453                               } // child kine cuts
454                           } else {
455                               pch[ncsel][0]=pc[0];
456                               pch[ncsel][1]=pc[1];
457                               pch[ncsel][2]=pc[2];
458                               kfch[ncsel]=kf;
459                               ncsel++;
460                           } // if child selection
461                       } // select muon
462                   } // decay particle loop
463               } // if decay products
464               
465               Int_t iparent;
466               if ((fCutOnChild && ncsel >0) || !fCutOnChild){
467                   ipa++;
468 //
469 // parent
470                   gAlice->
471                       SetTrack(0,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
472                   iparent=nt;
473                   gAlice->KeepTrack(nt); 
474                   for (i=0; i< ncsel; i++) {
475                       gAlice->SetTrack(fTrackIt,iparent,kfch[i],
476                                        &pch[i][0],och,polar,
477                                        0,kPDecay,nt,wgtch);
478                       gAlice->KeepTrack(nt); 
479                   }
480               }  // Decays by Lujet
481           } // kinematic selection
482           else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
483           {
484             gAlice->
485                 SetTrack(fTrackIt,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp);
486             ipa++; 
487           }
488           break;
489     } // while
490   } // event loop
491   gAlice->SetHighWaterMark(nt);
492   
493 }
494
495 Bool_t AliGenParam::ChildSelected(Int_t ip)
496 {
497 // True if particle is in list of selected children
498     for (Int_t i=0; i<5; i++)
499     {
500         if (fChildSelect[i]==ip) return kTRUE;
501     }
502     return kFALSE;
503 }
504
505 Bool_t AliGenParam::KinematicSelection(TParticle *particle)
506 {
507 // Perform kinematic cuts
508     Float_t px=particle->Px();
509     Float_t py=particle->Py();
510     Float_t pz=particle->Pz();
511 //
512 // momentum cut
513     Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
514     if (p > fPMax || p < fPMin) 
515     {
516 //      printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
517         return kFALSE;
518     }
519     Float_t pt=TMath::Sqrt(px*px+py*py);
520     
521 //
522 // theta cut
523     Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
524     if (theta > fThetaMax || theta < fThetaMin) 
525     {
526 //      printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
527         return kFALSE;
528     }
529
530     return kTRUE;
531 }
532
533
534 AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
535 {
536 // Assignment operator
537     return *this;
538 }
539
540
541