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