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