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