1 /**************************************************************************
2 * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Generator for muons according to kinematic parametrizations at ALICE
18 // (not at the surface).
19 // Origin: andrea.dainese@lnl.infn.it
21 #include <TParticle.h>
25 #include "AliESDtrack.h"
26 #include "AliESDVertex.h"
27 #include "AliGenCosmicsParam.h"
29 ClassImp(AliGenCosmicsParam)
31 //-----------------------------------------------------------------------------
32 AliGenCosmicsParam::AliGenCosmicsParam():
37 fMaxAngleWRTVertical(-99.),
51 // Default constructor
53 SetNumberParticles(1);
55 //-----------------------------------------------------------------------------
56 void AliGenCosmicsParam::Generate()
66 Double_t ptot=0,pt=0,angleWRTVertical=0;
67 Bool_t okMom=kFALSE,okAngle=kFALSE;
69 Float_t rtrigger=1000.0,ztrigger=600.0;
70 if(fTPC) { rtrigger=250.0; ztrigger=250.0; }
71 if(fITS) { rtrigger=50.0; ztrigger=50.0; }
72 if(fSPDinner) { rtrigger=3.5; ztrigger=14.0; }
73 if(fSPDouter) { rtrigger=6.5; ztrigger=14.0; }
74 if(fSDDinner) { rtrigger=14.0; ztrigger=21.0; }
75 if(fSDDouter) { rtrigger=23.0; ztrigger=29.0; }
76 if(fSSDinner) { rtrigger=37.0; ztrigger=42.0; }
77 if(fSSDouter) { rtrigger=42.0; ztrigger=48.0; }
81 Float_t muMinusFraction = 4./9.; // mu+/mu- = 1.25
87 while (npart<fNpart) {
89 if(gRandom->Rndm()<muMinusFraction) {
95 if(fParamACORDE) { // extracted from AliGenACORDE events
96 // sample total momentum only once (to speed up)
97 TF1 *dNdpACORDE = new TF1("dNdpACORDE","x/(1.+(x/12.8)*(x/12.8))^1.96",fPMin,fPMax);
98 ptot = (Double_t)dNdpACORDE->GetRandom();
106 origin[0] = (fYOrigin*TMath::Tan(fMaxAngleWRTVertical)+rtrigger)*(-1.+2.*gRandom->Rndm());
107 origin[1] = fYOrigin;
108 origin[2] = (fYOrigin*TMath::Tan(fMaxAngleWRTVertical)+ztrigger)*(-1.+2.*gRandom->Rndm());
112 okMom=kFALSE; okAngle=kFALSE;
114 if(fParamMI) { // parametrization by M.Ivanov of LEP cosmics data
115 Float_t pref = 1. + gRandom->Exp(30.);
117 p[0] = gRandom->Gaus(0.0,0.2)*pref;
118 p[2] = gRandom->Gaus(0.0,0.2)*pref;
119 if(gRandom->Rndm()>0.9) {
120 p[0] = gRandom->Gaus(0.0,0.4)*pref;
121 p[2] = gRandom->Gaus(0.0,0.4)*pref;
123 ptot=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
124 pt=TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
125 } else if(fParamACORDE) { // extracted from AliGenACORDE events
128 theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
129 if(TMath::Abs(theta-0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
132 phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
133 if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
135 pt = ptot*TMath::Sin(theta);
136 p[0] = pt*TMath::Cos(phi);
137 p[1] = pt*TMath::Sin(phi);
138 p[2] = ptot*TMath::Cos(theta);
140 AliFatal("Parametrization not set: use SetParamMI or SetParamACORDE");
144 // check kinematic cuts
145 if(TestBit(kMomentumRange)) {
146 if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
151 angleWRTVertical=TMath::ACos(TMath::Abs(p[1])/ptot); // acos(|py|/ptot)
152 if(angleWRTVertical<fMaxAngleWRTVertical) okAngle=kTRUE;
154 if(okAngle&&okMom) break;
157 // acceptance trigger
158 if(IntersectCylinder(rtrigger,ztrigger,ipart,origin,p)) {
162 if(IntersectACORDE(ipart,origin,p)) break;
168 Float_t polarization[3]= {0,0,0};
169 PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
171 //printf("TRIALS %d\n",trials);
176 //-----------------------------------------------------------------------------
177 void AliGenCosmicsParam::Init()
180 // Initialisation, check consistency of selected ranges
182 if(TestBit(kPtRange))
183 AliFatal("You cannot set the pt range for this generator! Only momentum range");
186 if(TestBit(kMomentumRange))
187 AliWarning("Minimum momentum cannot be < 8 GeV/c");
189 if(fMaxAngleWRTVertical<0.)
190 AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");
192 printf("************ AliGenCosmicsParam ****************\n");
193 printf("***** Muons generated at Y = %f cm\n",fYOrigin);
194 printf("************************************************\n");
198 //-----------------------------------------------------------------------------
199 Bool_t AliGenCosmicsParam::IntersectCylinder(Float_t r,Float_t z,Int_t pdg,
200 Float_t o[3],Float_t p[3]) const
203 // Intersection between muon and cylinder [-z,+z] with radius r
206 Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
207 TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
208 AliESDtrack track(&part);
209 Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
210 AliESDVertex origin(pos,sigma);
212 track.RelateToVertex(&origin,fBkG,10000.);
214 Float_t d0z0[2],covd0z0[3];
215 track.GetImpactParameters(d0z0,covd0z0);
218 if(TMath::Abs(d0z0[0])>r) return kFALSE;
220 if(TMath::Abs(d0z0[1])>z) return kFALSE;
223 if(TMath::Abs(fB)<0.01) { // NO FIELD
224 Float_t drphi = TMath::Abs(o[1]-p[1]/p[0]*o[0])/
225 TMath::Sqrt(p[1]*p[1]/p[0]/p[0]+1.);
226 if(drphi>r) return kFALSE;
227 Float_t dz = o[2]-p[2]/p[0]*o[0]+p[2]/p[0]*
228 (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]);
229 if(TMath::Abs(dz)>z) return kFALSE;
235 //-----------------------------------------------------------------------------
236 Bool_t AliGenCosmicsParam::IntersectACORDE(Int_t pdg,
237 Float_t o[3],Float_t p[3]) const
240 // Intersection between muon and ACORDE (very rough)
243 Float_t en = TMath::Sqrt(0.105*0.105+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
244 TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
245 AliESDtrack track(&part);
247 Float_t rACORDE=800.0,xACORDE=750.0,zACORDE=500.0;
248 if(fACORDE4ITS) { xACORDE=100.0; zACORDE=100.0; }
251 track.GetXYZAt(rACORDE,fBkG,xyz);
254 if(TMath::Abs(xyz[0]) > xACORDE) return kFALSE;
256 if(TMath::Abs(xyz[2]) > zACORDE) return kFALSE;
260 //-----------------------------------------------------------------------------