effc++ warnings
[u/mrichter/AliRoot.git] / EVGEN / AliGenCosmicsParam.cxx
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 ClassImp(AliGenCosmicsParam)
30
31 //-----------------------------------------------------------------------------
32 AliGenCosmicsParam::AliGenCosmicsParam():
33 AliGenerator(),
34 fParamMI(kFALSE),
35 fParamACORDE(kFALSE),
36 fYOrigin(600.),
37 fMaxAngleWRTVertical(-99.),
38 fBkG(0.),
39 fTPC(kFALSE),
40 fITS(kFALSE),
41 fSPDinner(kFALSE),
42 fSPDouter(kFALSE),
43 fSDDinner(kFALSE),
44 fSDDouter(kFALSE),
45 fSSDinner(kFALSE),
46 fSSDouter(kFALSE),
47 fACORDE(kFALSE),
48 fACORDE4ITS(kFALSE),
49 fBottomScintillator(kFALSE)
50 {
51   //
52   // Default constructor
53   //
54   SetNumberParticles(1);
55 }
56 //-----------------------------------------------------------------------------
57 void 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   //
70   Float_t rtrigger=1000.0,ztrigger=600.0;
71   if(fTPC)      { rtrigger=250.0; ztrigger=250.0; }
72   if(fITS)      { rtrigger=50.0; ztrigger=50.0; }
73   if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }
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; }
79
80
81   // mu+ or mu-
82   Float_t muMinusFraction = 4./9.; // mu+/mu- = 1.25
83   Int_t ipart;
84
85   Int_t trials=0;
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     }
95
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     }
103
104     while(1) {
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");
142         }
143         
144         
145         // check kinematic cuts
146         if(TestBit(kMomentumRange)) {
147           if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
148         } else {
149           okMom=kTRUE;
150         }
151
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       }
157
158       // acceptance trigger
159       if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
160         if(fACORDE && !fBottomScintillator) {
161           if(IntersectACORDE(ipart,origin,p)) break;
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;
169         }
170       }
171       //
172     }
173
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   }
179
180   return;
181 }
182 //-----------------------------------------------------------------------------
183 void 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 //-----------------------------------------------------------------------------
205 Bool_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 //-----------------------------------------------------------------------------
242 Bool_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]);
250   TParticle part(pdg,0,0,0,0,0,-p[0],-p[1],-p[2],en,o[0],o[1],o[2],0);
251   AliESDtrack track(&part);
252
253   Float_t rACORDE=800.0,xACORDE=750.0/*250.0*/,zACORDE=500.0;
254   if(fACORDE4ITS) { xACORDE=100.0; zACORDE=100.0; }
255
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]);
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 //-----------------------------------------------------------------------------
272 Bool_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 //-----------------------------------------------------------------------------