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