New reader for the pedestal run and vdrift (Julian) and some bug fixing (Raphaelle)
[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"
61#include "AliESDtrack.h"
62#include "AliESDVertex.h"
63#include "AliGenFunction.h"
64
65ClassImp(AliGenFunction)
66
67//-----------------------------------------------------------------------------
68AliGenFunction::AliGenFunction():
69 AliGenerator(),
70 fBkG(0.),
71 fFMomentum(0), // momentum distribution function
72 fFPhi(0), // phi distribution function
73 fFTheta(0), // theta distribution function
74 fFPosition(0), // position distribution function
75 fFPdg(0), // pdg distribution function
76 fRefRadius(0), // reference radius to be crossed
77 fZmin(0), // minimal z at reference radius
78 fZmax(0), // z at reference radius
79 fMaxTrial(10000) // maximal number of attempts
80{
81 //
82 // Default constructor
83 //
84 SetNumberParticles(1);
85}
21eff5bc 86
87AliGenFunction::AliGenFunction(const AliGenFunction& func):
88 AliGenerator(),
89 fBkG(func.fBkG),
90 fFMomentum(func.fFMomentum), // momentum distribution function
91 fFPhi(func.fFPhi), // phi distribution function
92 fFTheta(func.fFTheta), // theta distribution function
93 fFPosition(func.fFPosition), // position distribution function
94 fFPdg(func.fFPdg), // pdg distribution function
95 fRefRadius(func.fRefRadius), // reference radius to be crossed
96 fZmin(func.fZmin), // minimal z at reference radius
97 fZmax(func.fZmax), // z at reference radius
98 fMaxTrial(10000) // maximal number of attempts
99{
100 // Copy constructor
101 SetNumberParticles(1);
102}
103
104AliGenFunction & AliGenFunction::operator=(const AliGenFunction& func)
105{
106 // Assigment operator
107 if(&func == this) return *this;
108 fBkG = func.fBkG;
109 fFMomentum = func.fFMomentum;
110 fFPhi = func.fFPhi;
111 fFTheta = func.fFTheta;
112 fFPosition = func.fFPosition;
113 fFPdg = func.fFPdg;
114 fRefRadius = func.fRefRadius;
115 fZmin = func.fZmin;
116 fZmax = func.fZmax;
117 fMaxTrial = func.fMaxTrial;
118 return *this;
119}
120
121
ccc5bf38 122//-----------------------------------------------------------------------------
123void AliGenFunction::Generate()
124{
125 //
126 // Generate one muon
127 //
128 Int_t naccepted =0;
129
130 for (Int_t ipart=0; ipart<fMaxTrial && naccepted<fNpart; ipart++){
131 //
132 //
133 //
134 Float_t mom[3];
135 Float_t posf[3];
136 Double_t pos[3];
137 Int_t pdg;
21eff5bc 138 Double_t ptot, pt, phi, theta;
ccc5bf38 139 //
140 ptot = fFMomentum->GetRandom();
141 phi = fFPhi->GetRandom();
142 theta = fFTheta->GetRandom();
143 pt = ptot*TMath::Sin(theta);
144 mom[0] = pt*TMath::Cos(phi);
145 mom[1] = pt*TMath::Sin(phi);
146 mom[2] = ptot*TMath::Cos(theta);
147
148 //
149 fFPosition->GetRandom3(pos[0],pos[1],pos[2]);
150 posf[0]=pos[0];
151 posf[1]=pos[1];
152 posf[2]=pos[2];
153 pdg = TMath::Nint(fFPdg->GetRandom());
154 if (!IntersectCylinder(fRefRadius,fZmin, fZmax, pdg, posf, mom)) continue;
155 //
156 //
157 Float_t polarization[3]= {0,0,0};
158 Int_t nt;
159 PushTrack(fTrackIt,-1,pdg,mom, posf, polarization,0,kPPrimary,nt);
160 naccepted++;
161 }
162
163
164 return;
165}
166//-----------------------------------------------------------------------------
167void AliGenFunction::Init()
168{
169 //
170 // Initialisation, check consistency of selected ranges
171 //
172 printf("************ AliGenFunction ****************\n");
173 printf("************************************************\n");
174 if (!fFMomentum){
175 AliFatal("Momentum distribution function not specified");
176 }
177 if (!fFPosition){
178 AliFatal("Position distribution function not specified");
179 }
180 if (!fFPdg){
181 AliFatal("PDG distribution function not specified");
182 }
183 if (fZmin==fZmax){
184 AliFatal("Z range not specified");
185 }
186
187 return;
188}
189
190void AliGenFunction::SetFunctions(TF1 * momentum, TF1 *fphi, TF1 *ftheta,TF3 * position, TF1* pdg){
191 //
192 //
193 //
194 fFMomentum = momentum;
195 fFPhi = fphi;
196 fFTheta = ftheta;
197 fFPosition = position;
198 fFPdg = pdg;
199}
200
201void AliGenFunction::SetCylinder(Double_t refR, Double_t zmin, Double_t zmax){
202 //
203 //
204 //
205 fRefRadius = refR; // reference radius to be crossed
206 fZmin = zmin; // minimal z at reference radius
207 fZmax = zmax; // maximal z at reference radius
208
209}
210
211
212//-----------------------------------------------------------------------------
213Bool_t AliGenFunction::IntersectCylinder(Float_t r,Float_t zmin, Float_t zmax,Int_t pdg,
214 Float_t o[3],Float_t p[3]) const
215{
216 //
217 // Intersection between muon and cylinder [-z,+z] with radius r
218 //
219 Double_t mass=0;
220 if (TDatabasePDG::Instance()->GetParticle(pdg)){
221 mass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
222 }
223
224 Float_t en = TMath::Sqrt(mass*mass+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
225 TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
226 AliESDtrack track(&part);
227 Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
228 AliESDVertex origin(pos,sigma);
229
230 track.RelateToVertex(&origin,fBkG,10000.);
231
232 Float_t d0z0[2],covd0z0[3];
233 track.GetImpactParameters(d0z0,covd0z0);
234
235 // check rphi
236 if(TMath::Abs(d0z0[0])>r) return kFALSE;
237 // check z
238 if(d0z0[1]>zmax) return kFALSE;
239 if(d0z0[1]<zmin) return kFALSE;
240 //
241 return kTRUE;
242}