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 :
17 : // Generator for muons according to kinematic parametrizations at ALICE
18 : // (not at the surface).
19 : // Origin: andrea.dainese@lnl.infn.it
20 :
21 : #include <TParticle.h>
22 : #include <TF1.h>
23 :
24 : #include "AliRun.h"
25 : #include "AliESDtrack.h"
26 : #include "AliESDVertex.h"
27 : #include "AliGenCosmicsParam.h"
28 :
29 6 : ClassImp(AliGenCosmicsParam)
30 :
31 : //-----------------------------------------------------------------------------
32 : AliGenCosmicsParam::AliGenCosmicsParam():
33 0 : AliGenerator(),
34 0 : fParamMI(kFALSE),
35 0 : fParamACORDE(kFALSE),
36 0 : fParamDataTPC(kTRUE),
37 0 : fYOrigin(600.),
38 0 : fMaxAngleWRTVertical(-99.),
39 0 : fBkG(0.),
40 0 : fTPC(kFALSE),
41 0 : fITS(kFALSE),
42 0 : fSPDinner(kFALSE),
43 0 : fSPDouter(kFALSE),
44 0 : fSDDinner(kFALSE),
45 0 : fSDDouter(kFALSE),
46 0 : fSSDinner(kFALSE),
47 0 : fSSDouter(kFALSE),
48 0 : fACORDE(kFALSE),
49 0 : fACORDE4ITS(kFALSE),
50 0 : fBottomScintillator(kFALSE)
51 0 : {
52 : //
53 : // Default constructor
54 : //
55 0 : SetNumberParticles(1);
56 0 : }
57 : //-----------------------------------------------------------------------------
58 : void AliGenCosmicsParam::Generate()
59 : {
60 : //
61 : // Generate one muon
62 : //
63 :
64 : //
65 0 : Float_t origin[3];
66 0 : Float_t p[3];
67 0 : Int_t nt;
68 : Double_t ptot=0,pt=0,angleWRTVertical=0;
69 : Bool_t okMom=kFALSE,okAngle=kFALSE;
70 : //
71 : Float_t rtrigger=1000.0,ztrigger=600.0;
72 0 : if(fTPC) { rtrigger=250.0; ztrigger=250.0; }
73 0 : if(fITS) { rtrigger=50.0; ztrigger=50.0; }
74 0 : if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }
75 0 : if(fSPDouter) { rtrigger=6.5; ztrigger=14.0; }
76 0 : if(fSDDinner) { rtrigger=14.0; ztrigger=21.0; }
77 0 : if(fSDDouter) { rtrigger=23.0; ztrigger=29.0; }
78 0 : if(fSSDinner) { rtrigger=37.0; ztrigger=42.0; }
79 0 : if(fSSDouter) { rtrigger=42.0; ztrigger=48.0; }
80 :
81 :
82 : // mu+ or mu-
83 : Float_t muMinusFraction = 4./9.; // mu+/mu- = 1.25
84 : Int_t ipart;
85 :
86 : Int_t trials=0;
87 : Int_t npart=0;
88 :
89 0 : while (npart<fNpart) {
90 :
91 0 : if(gRandom->Rndm()<muMinusFraction) {
92 : ipart = 13; // mu-
93 0 : } else {
94 : ipart = -13; // mu+
95 : }
96 :
97 0 : if(fParamACORDE) { // extracted from AliGenACORDE events
98 : // sample total momentum only once (to speed up)
99 0 : TF1 *dNdpACORDE = new TF1("dNdpACORDE","x/(1.+(x/12.8)*(x/12.8))^1.96",fPMin,fPMax);
100 0 : ptot = (Double_t)dNdpACORDE->GetRandom();
101 0 : delete dNdpACORDE;
102 : dNdpACORDE = 0;
103 0 : } else if(fParamDataTPC) { // extracted from cosmics in TPC (Summer 08)
104 : // sample total momentum only once (to speed up)
105 0 : TF1 *dNdpTPC = new TF1("dNdpTPC","x/(1.+(x/3.)*(x/3.))^1.",fPMin,fPMax);
106 0 : ptot = (Double_t)dNdpTPC->GetRandom();
107 0 : delete dNdpTPC;
108 : dNdpTPC = 0;
109 0 : }
110 :
111 : while(1) {
112 0 : trials++;
113 : // origin
114 0 : origin[0] = (fYOrigin*TMath::Tan(fMaxAngleWRTVertical)+rtrigger)*(-1.+2.*gRandom->Rndm());
115 0 : origin[1] = fYOrigin;
116 0 : origin[2] = (fYOrigin*TMath::Tan(fMaxAngleWRTVertical)+ztrigger)*(-1.+2.*gRandom->Rndm());
117 :
118 : // momentum
119 0 : while(1) {
120 : okMom=kFALSE; okAngle=kFALSE;
121 :
122 0 : if(fParamMI) { // parametrization by M.Ivanov of LEP cosmics data
123 0 : Float_t pref = 1. + gRandom->Exp(30.);
124 0 : p[1] = -pref;
125 0 : p[0] = gRandom->Gaus(0.0,0.2)*pref;
126 0 : p[2] = gRandom->Gaus(0.0,0.2)*pref;
127 0 : if(gRandom->Rndm()>0.9) {
128 0 : p[0] = gRandom->Gaus(0.0,0.4)*pref;
129 0 : p[2] = gRandom->Gaus(0.0,0.4)*pref;
130 0 : }
131 0 : ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
132 0 : pt=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
133 0 : } else if(fParamACORDE || fParamDataTPC) {
134 : Float_t theta,phi;
135 0 : while(1) {
136 0 : theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
137 0 : if(TMath::Abs(theta-0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
138 : }
139 : while(1) {
140 0 : phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
141 0 : if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
142 : }
143 0 : pt = ptot*TMath::Sin(theta);
144 0 : p[0] = pt*TMath::Cos(phi);
145 0 : p[1] = pt*TMath::Sin(phi);
146 0 : p[2] = ptot*TMath::Cos(theta);
147 0 : } else {
148 0 : AliFatal("Parametrization not set: use SetParamDataTPC, SetParamMI, or SetParamACORDE");
149 : }
150 :
151 :
152 : // check kinematic cuts
153 0 : if(TestBit(kMomentumRange)) {
154 0 : if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
155 : } else {
156 : okMom=kTRUE;
157 : }
158 :
159 0 : angleWRTVertical=TMath::ACos(TMath::Abs(p[1])/ptot); // acos(|py|/ptot)
160 0 : if(angleWRTVertical<fMaxAngleWRTVertical) okAngle=kTRUE;
161 :
162 0 : if(okAngle&&okMom) break;
163 : }
164 :
165 : // acceptance trigger
166 0 : if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
167 0 : if(fACORDE && !fBottomScintillator) {
168 0 : if(IntersectACORDE(ipart,origin,p)) break;
169 0 : } else if(!fACORDE && fBottomScintillator) {
170 0 : if(IntersectBottomScintillator(ipart,origin,p)) break;
171 0 : } else if(fACORDE && fBottomScintillator) {
172 0 : if(IntersectACORDE(ipart,origin,p) &&
173 0 : IntersectBottomScintillator(ipart,origin,p)) break;
174 : } else { // !fACORDE && !fBottomScintillator
175 : break;
176 : }
177 : }
178 : //
179 : }
180 :
181 0 : Float_t polarization[3]= {0,0,0};
182 0 : PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
183 0 : npart++;
184 : //printf("TRIALS %d\n",trials);
185 0 : }
186 :
187 : return;
188 0 : }
189 : //-----------------------------------------------------------------------------
190 : void AliGenCosmicsParam::Init()
191 : {
192 : //
193 : // Initialisation, check consistency of selected ranges
194 : //
195 0 : if(TestBit(kPtRange))
196 0 : AliFatal("You cannot set the pt range for this generator! Only momentum range");
197 : Double_t pmin=8.; // fParamACORDE
198 0 : if(fParamDataTPC) pmin=0.5;
199 0 : if(fPMin<pmin) {
200 0 : fPMin=pmin;
201 0 : if(TestBit(kMomentumRange))
202 0 : AliWarning(Form("Minimum momentum cannot be < %f GeV/c",pmin));
203 : }
204 0 : if(fMaxAngleWRTVertical<0.)
205 0 : AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");
206 :
207 0 : printf("************ AliGenCosmicsParam ****************\n");
208 0 : printf("***** Muons generated at Y = %f cm\n",fYOrigin);
209 0 : printf("************************************************\n");
210 :
211 : return;
212 0 : }
213 : //-----------------------------------------------------------------------------
214 : Bool_t AliGenCosmicsParam::IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
215 : Float_t o[3],Float_t p[3]) const
216 : {
217 : //
218 : // Intersection between muon and cylinder [-z,+z] with radius r
219 : //
220 :
221 0 : Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
222 0 : TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
223 0 : AliESDtrack track(&part);
224 0 : Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
225 0 : AliESDVertex origin(pos,sigma);
226 :
227 0 : track.RelateToVertex(&origin,fBkG,10000.);
228 :
229 0 : Float_t d0z0[2],covd0z0[3];
230 0 : track.GetImpactParameters(d0z0,covd0z0);
231 :
232 : // check rphi
233 0 : if(TMath::Abs(d0z0[0])>r) return kFALSE;
234 : // check z
235 0 : if(TMath::Abs(d0z0[1])>z) return kFALSE;
236 :
237 : /*
238 : if(TMath::Abs(fB)<0.01) { // NO FIELD
239 : Float_t drphi = TMath::Abs(o[1]-p[1]/p[0]*o[0])/
240 : TMath::Sqrt(p[1]*p[1]/p[0]/p[0]+1.);
241 : if(drphi>r) return kFALSE;
242 : Float_t dz = o[2]-p[2]/p[0]*o[0]+p[2]/p[0]*
243 : (p[1]*p[1]/p[0]/p[0]*o[0]-p[1]/p[0]*o[1])/(1.+p[1]*p[1]/p[0]/p[0]);
244 : if(TMath::Abs(dz)>z) return kFALSE;
245 : }
246 : */
247 :
248 0 : return kTRUE;
249 0 : }
250 : //-----------------------------------------------------------------------------
251 : Bool_t AliGenCosmicsParam::IntersectACORDE(Int_t pdg,
252 : Float_t o[3],Float_t p[3]) const
253 : {
254 : //
255 : // Intersection between muon and ACORDE (very rough)
256 : //
257 :
258 0 : Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
259 0 : TParticle part(pdg,0,0,0,0,0,-p[0],-p[1],-p[2],en,o[0],o[1],o[2],0);
260 0 : AliESDtrack track(&part);
261 :
262 : Float_t rACORDE=800.0,xACORDE=750.0/*250.0*/,zACORDE=500.0;
263 0 : if(fACORDE4ITS) { xACORDE=100.0; zACORDE=100.0; }
264 :
265 0 : Double_t planepoint[3]={0.,rACORDE,0.};
266 0 : Double_t planenorm[3]={0.,1.,0.};
267 :
268 0 : if(!track.Intersect(planepoint,planenorm,fBkG)) return kFALSE;
269 :
270 0 : Double_t xyz[3]={planepoint[0],planepoint[1],planepoint[2]};
271 : //printf("XYZ = %f %f %f\n",xyz[0],xyz[1],xyz[2]);
272 :
273 : // check global x
274 0 : if(TMath::Abs(xyz[0]) > xACORDE) return kFALSE;
275 : // check global z
276 0 : if(TMath::Abs(xyz[2]) > zACORDE) return kFALSE;
277 :
278 0 : return kTRUE;
279 0 : }
280 : //-----------------------------------------------------------------------------
281 : Bool_t AliGenCosmicsParam::IntersectBottomScintillator(Int_t pdg,
282 : Float_t o[3],Float_t p[3]) const
283 : {
284 : //
285 : // Intersection between muon and ACORDE (very rough)
286 : //
287 :
288 0 : Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
289 0 : TParticle part(pdg,0,0,0,0,0,-p[0],-p[1],-p[2],en,o[0],o[1],o[2],0);
290 0 : AliESDtrack track(&part);
291 :
292 : Double_t xSc=40.,ySc=-350.,zSc=40.;
293 :
294 0 : Double_t planepoint[3]={0.,ySc,0.};
295 0 : Double_t planenorm[3]={0.,1.,0.};
296 :
297 0 : if(!track.Intersect(planepoint,planenorm,fBkG)) return kFALSE;
298 :
299 0 : Double_t xyz[3]={planepoint[0],planepoint[1],planepoint[2]};
300 : //printf("XYZ = %f %f %f\n",xyz[0],xyz[1],xyz[2]);
301 :
302 : // check global x
303 0 : if(TMath::Abs(xyz[0]) > xSc) return kFALSE;
304 : // check global z
305 0 : if(TMath::Abs(xyz[2]) > zSc) return kFALSE;
306 :
307 0 : return kTRUE;
308 0 : }
309 : //-----------------------------------------------------------------------------
|
