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