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