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