doxy: install THtml converter
[u/mrichter/AliRoot.git] / EVGEN / AliGenFunction.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2007, 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 // Generator for particles according generic functions
17 //  
18 //  TF1 *   fFMomentum;           // momentum distribution function inGeV
19 //  TF1 *   fFPhi;                // phi distribution function in rad
20 //  TF1 *   fFTheta;              // theta distribution function in rad
21 //  TF3 *   fFPosition;           // position distribution function in cm
22 //  TF1 *   fFPdg;                // pdg distribution function  
23 //  We assume that the moment, postion and PDG code of particles are independent  
24 //  Only tracks/particle crossing the reference radius at given z range
25 //
26 // Origin: marian.ivanov@cern.ch
27
28
29 /*
30   Example generic function for cosmic generation:
31  
32 //
33         AliGenFunction *generCosmic = new AliGenFunction;
34         generCosmic->SetNumberParticles(100);
35         TF1 *fmom   = new TF1("fmom","TMath::Exp(-x/8)",0,30);
36         //
37         TF1 *fphi   = new TF1("fphi","TMath::Gaus(x,-TMath::Pi()/2,0.3)",-3.14,3.14);
38         TF1 *ftheta = new TF1("ftheta","TMath::Gaus(x,TMath::Pi()/2,0.3)",-3.14,3.14);
39         TF3 *fpos = new TF3("fpos","1+(x+y+z)*0",-2000,2000,700,701,-2000,2000);
40         TF1 * fpdg= new TF1("fpdg","1*(abs(x-13)<0.1)+1*(abs(x+13)<0.1)",-300,300);
41         fpdg->SetNpx(10000);  // neccessary - number of bins for generation
42         generCosmic->SetVertexSmear(kPerEvent);
43         generCosmic->SetFunctions(fmom,fphi,ftheta,fpos,fpdg);
44         generCosmic->SetCylinder(375,-375,375);
45         generCosmic->SetBkG(0.2);
46         gener->AddGenerator(generCosmic,"generCosmic",1);
47
48
49   
50 */
51
52
53
54
55 #include <TParticle.h>
56 #include <TF1.h>
57 #include <TF3.h>
58 #include <TDatabasePDG.h>
59
60 #include "AliRun.h"
61 #include "AliLog.h"
62 #include "AliESDtrack.h"
63 #include "AliESDVertex.h"
64 #include "AliGenFunction.h"
65 #include "AliGenEventHeader.h"
66
67 ClassImp(AliGenFunction)
68
69 //-----------------------------------------------------------------------------
70 AliGenFunction::AliGenFunction():
71   AliGenerator(),
72   fBkG(0.),  
73   fFMomentum(0),           // momentum distribution function 
74   fFPhi(0),                // phi distribution function
75   fFTheta(0),              // theta distribution function
76   fFPosition(0),           // position distribution function 
77   fFPdg(0),                // pdg distribution function  
78   fRefRadius(0),           // reference radius to be crossed
79   fZmin(0),                // minimal z at reference radius
80   fZmax(0),                // z at reference radius
81   fMaxTrial(10000)         // maximal number of attempts
82 {
83   //
84   // Default constructor
85   //
86   SetNumberParticles(1);
87 }
88
89 AliGenFunction::AliGenFunction(const AliGenFunction& func):
90   AliGenerator(),
91   fBkG(func.fBkG),  
92   fFMomentum(func.fFMomentum),           // momentum distribution function 
93   fFPhi(func.fFPhi),                     // phi distribution function
94   fFTheta(func.fFTheta),                 // theta distribution function
95   fFPosition(func.fFPosition),           // position distribution function 
96   fFPdg(func.fFPdg),                     // pdg distribution function  
97   fRefRadius(func.fRefRadius),           // reference radius to be crossed
98   fZmin(func.fZmin),                     // minimal z at reference radius
99   fZmax(func.fZmax),                     // z at reference radius
100   fMaxTrial(10000)                       // maximal number of attempts
101 {
102     // Copy constructor
103     SetNumberParticles(1);
104 }
105
106 AliGenFunction & AliGenFunction::operator=(const AliGenFunction& func)
107 {
108     // Assigment operator
109       if(&func == this) return *this;
110       fBkG       = func.fBkG;
111       fFMomentum = func.fFMomentum; 
112       fFPhi      = func.fFPhi;      
113       fFTheta    = func.fFTheta;    
114       fFPosition = func.fFPosition; 
115       fFPdg      = func.fFPdg;      
116       fRefRadius = func.fRefRadius;      
117       fZmin      = func.fZmin;                
118       fZmax      = func.fZmax;                
119       fMaxTrial  = func.fMaxTrial;
120       return *this;
121 }
122
123
124 //-----------------------------------------------------------------------------
125 void AliGenFunction::Generate()
126 {
127   //
128   // Generate one muon
129   //
130   Int_t naccepted =0;
131   
132   for (Int_t ipart=0; ipart<fMaxTrial && naccepted<fNpart; ipart++){
133     //
134     //
135     //
136     Float_t mom[3];
137     Float_t posf[3];
138     Double_t pos[3];
139     Int_t pdg;
140     Double_t ptot, pt,  phi, theta; 
141     //
142     ptot     = fFMomentum->GetRandom();
143     phi      = fFPhi->GetRandom();
144     theta    = fFTheta->GetRandom();
145     pt     = ptot*TMath::Sin(theta);
146     mom[0] = pt*TMath::Cos(phi); 
147     mom[1] = pt*TMath::Sin(phi); 
148     mom[2] = ptot*TMath::Cos(theta);
149
150     //
151     fFPosition->GetRandom3(pos[0],pos[1],pos[2]);
152     posf[0]=pos[0];
153     posf[1]=pos[1];
154     posf[2]=pos[2];
155     pdg = TMath::Nint(fFPdg->GetRandom());
156     if (!IntersectCylinder(fRefRadius,fZmin, fZmax, pdg, posf, mom)) continue;
157     //
158     //
159     Float_t polarization[3]= {0,0,0};
160     Int_t nt;
161     PushTrack(fTrackIt,-1,pdg,mom, posf, polarization,0,kPPrimary,nt);
162     naccepted++;
163   }
164
165   AliGenEventHeader* header = new AliGenEventHeader("THn");
166   gAlice->SetGenEventHeader(header);
167   return;
168 }
169 //-----------------------------------------------------------------------------
170 void AliGenFunction::Init()
171 {
172   // 
173   // Initialisation, check consistency of selected ranges
174   //
175   printf("************ AliGenFunction ****************\n");
176   printf("************************************************\n");
177   if (!fFMomentum){
178     AliFatal("Momentum distribution function not specified");
179   }
180   if (!fFPosition){
181     AliFatal("Position distribution function not specified");
182   }
183   if (!fFPdg){
184     AliFatal("PDG distribution function not specified");
185   }
186   if (fZmin==fZmax){
187     AliFatal("Z range not specified");
188   }
189
190   return;
191 }
192
193 void AliGenFunction::SetFunctions(TF1 * momentum, TF1 *fphi, TF1 *ftheta,TF3 * position, TF1* pdg){
194   //
195   // Set the function
196   //
197   fFMomentum = momentum;
198   fFPhi = fphi;
199   fFTheta = ftheta;
200   fFPosition = position;
201   fFPdg = pdg;
202 }
203
204 void AliGenFunction::SetCylinder(Double_t refR, Double_t zmin, Double_t zmax){
205   //
206   // Set the cylinder geometry
207   //
208   fRefRadius = refR;          // reference radius to be crossed
209   fZmin = zmin;               // minimal z at reference radius
210   fZmax = zmax;               // maximal z at reference radius
211
212 }
213
214
215 //-----------------------------------------------------------------------------
216 Bool_t AliGenFunction::IntersectCylinder(Float_t r,Float_t zmin, Float_t zmax,Int_t pdg,
217                                          Float_t o[3],Float_t p[3]) const
218 {
219   //
220   // Intersection between muon and cylinder [-z,+z] with radius r
221   //
222   Double_t mass=0;
223   if (TDatabasePDG::Instance()->GetParticle(pdg)){
224     mass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
225   }
226
227   Float_t en = TMath::Sqrt(mass*mass+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
228   TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
229   AliESDtrack track(&part);
230   Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
231   AliESDVertex origin(pos,sigma);
232
233   track.RelateToVertex(&origin,fBkG,10000.);
234
235   Float_t d0z0[2],covd0z0[3];
236   track.GetImpactParameters(d0z0,covd0z0);
237
238   // check rphi 
239   if(TMath::Abs(d0z0[0])>r) return kFALSE;
240   // check z
241   if(d0z0[1]>zmax) return kFALSE;
242   if(d0z0[1]<zmin) return kFALSE;
243   //
244   return kTRUE;
245 }