doxy: TPC/stressTest/testSparse converted
[u/mrichter/AliRoot.git] / EVGEN / AliGenFunction.cxx
CommitLineData
ccc5bf38 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>
48c27d61 58#include <TDatabasePDG.h>
ccc5bf38 59
60#include "AliRun.h"
bc094e42 61#include "AliLog.h"
ccc5bf38 62#include "AliESDtrack.h"
63#include "AliESDVertex.h"
64#include "AliGenFunction.h"
fee5f683 65#include "AliGenEventHeader.h"
ccc5bf38 66
67ClassImp(AliGenFunction)
68
69//-----------------------------------------------------------------------------
70AliGenFunction::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}
21eff5bc 88
89AliGenFunction::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
106AliGenFunction & 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
ccc5bf38 124//-----------------------------------------------------------------------------
125void 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;
21eff5bc 140 Double_t ptot, pt, phi, theta;
ccc5bf38 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
fee5f683 165 AliGenEventHeader* header = new AliGenEventHeader("THn");
166 gAlice->SetGenEventHeader(header);
ccc5bf38 167 return;
168}
169//-----------------------------------------------------------------------------
170void 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
193void AliGenFunction::SetFunctions(TF1 * momentum, TF1 *fphi, TF1 *ftheta,TF3 * position, TF1* pdg){
194 //
bc094e42 195 // Set the function
ccc5bf38 196 //
197 fFMomentum = momentum;
198 fFPhi = fphi;
199 fFTheta = ftheta;
200 fFPosition = position;
201 fFPdg = pdg;
202}
203
204void AliGenFunction::SetCylinder(Double_t refR, Double_t zmin, Double_t zmax){
205 //
bc094e42 206 // Set the cylinder geometry
ccc5bf38 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//-----------------------------------------------------------------------------
216Bool_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}