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