newly implemented model for slow nucleon production
[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 fParamDataTPC(kTRUE),
37 fYOrigin(600.),
38 fMaxAngleWRTVertical(-99.),
39 fBkG(0.),
40 fTPC(kFALSE),
41 fITS(kFALSE),
42 fSPDinner(kFALSE),
43 fSPDouter(kFALSE),
44 fSDDinner(kFALSE),
45 fSDDouter(kFALSE),
46 fSSDinner(kFALSE),
47 fSSDouter(kFALSE),
48 fACORDE(kFALSE),
49 fACORDE4ITS(kFALSE),
50 fBottomScintillator(kFALSE)
51 {
52   //
53   // Default constructor
54   //
55   SetNumberParticles(1);
56 }
57 //-----------------------------------------------------------------------------
58 void 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   //
71   Float_t rtrigger=1000.0,ztrigger=600.0;
72   if(fTPC)      { rtrigger=250.0; ztrigger=250.0; }
73   if(fITS)      { rtrigger=50.0; ztrigger=50.0; }
74   if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }
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; }
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   while (npart<fNpart) {
90
91     if(gRandom->Rndm()<muMinusFraction) {
92       ipart = 13; // mu-
93     } else {
94       ipart = -13; // mu+
95     }
96
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;
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;
109     }
110
111     while(1) {
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]);
133         } else if(fParamACORDE || fParamDataTPC) {
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 {
148           AliFatal("Parametrization not set: use SetParamDataTPC, SetParamMI, or SetParamACORDE");
149         }
150         
151         
152         // check kinematic cuts
153         if(TestBit(kMomentumRange)) {
154           if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
155         } else {
156           okMom=kTRUE;
157         }
158
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       }
164
165       // acceptance trigger
166       if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
167         if(fACORDE && !fBottomScintillator) {
168           if(IntersectACORDE(ipart,origin,p)) break;
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;
176         }
177       }
178       //
179     }
180
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   }
186
187   return;
188 }
189 //-----------------------------------------------------------------------------
190 void 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");
197   Double_t pmin=8.; // fParamACORDE
198   if(fParamDataTPC) pmin=0.5;
199   if(fPMin<pmin) { 
200     fPMin=pmin; 
201     if(TestBit(kMomentumRange)) 
202       AliWarning(Form("Minimum momentum cannot be < %f GeV/c",pmin)); 
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 //-----------------------------------------------------------------------------
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   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 //-----------------------------------------------------------------------------
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   Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
259   TParticle part(pdg,0,0,0,0,0,-p[0],-p[1],-p[2],en,o[0],o[1],o[2],0);
260   AliESDtrack track(&part);
261
262   Float_t rACORDE=800.0,xACORDE=750.0/*250.0*/,zACORDE=500.0;
263   if(fACORDE4ITS) { xACORDE=100.0; zACORDE=100.0; }
264
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]);
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 //-----------------------------------------------------------------------------
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   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 //-----------------------------------------------------------------------------