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.),
45 // Default constructor
48 //-----------------------------------------------------------------------------
49 void AliGenCosmicsParam::Generate()
59 Double_t ptot=0,pt=0,angleWRTVertical=0;
60 Bool_t okMom=kFALSE,okAngle=kFALSE;
65 if(gRandom->Rndm()<0.5) ipart *= -1;
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();
80 origin[0] = fYOrigin*TMath::Tan(fMaxAngleWRTVertical)*(-1.+2.*gRandom->Rndm());
82 origin[2] = fYOrigin*TMath::Tan(fMaxAngleWRTVertical)*(-1.+2.*gRandom->Rndm());
86 okMom=kFALSE; okAngle=kFALSE;
88 if(fParamMI) { // parametrization by M.Ivanov of LEP cosmics data
89 Float_t pref = 1. + gRandom->Exp(30.);
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;
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
102 theta = gRandom->Gaus(0.5*TMath::Pi(),0.42);
103 if(TMath::Abs(theta-0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
106 phi = gRandom->Gaus(-0.5*TMath::Pi(),0.42);
107 if(TMath::Abs(phi+0.5*TMath::Pi())<fMaxAngleWRTVertical) break;
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);
114 AliFatal("Parametrization not set: use SetParamMI or SetParamACORDE");
118 // check kinematic cuts
119 if(TestBit(kMomentumRange)) {
120 if(ptot>fPMin && ptot<fPMax) okMom=kTRUE;
123 angleWRTVertical=TMath::ACos(TMath::Abs(p[1])/ptot); // acos(|py|/ptot)
124 if(angleWRTVertical<fMaxAngleWRTVertical) okAngle=kTRUE;
126 if(okAngle&&okMom) break;
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;
137 Float_t polarization[3]= {0,0,0};
138 PushTrack(fTrackIt,-1,ipart,p,origin,polarization,0,kPPrimary,nt);
140 //printf("TRIALS %d\n",trials);
145 //-----------------------------------------------------------------------------
146 void AliGenCosmicsParam::Init()
149 // Initialisation, check consistency of selected ranges
151 if(TestBit(kPtRange))
152 AliFatal("You cannot set the pt range for this generator! Only momentum range");
155 if(TestBit(kMomentumRange))
156 AliWarning("Minimum momentum cannot be < 8 GeV/c");
158 if(fMaxAngleWRTVertical<0.)
159 AliFatal("You must use SetMaxAngleWRTVertical() instead of SetThetaRange(), SetPhiRange()");
161 printf("************ AliGenCosmicsParam ****************\n");
162 printf("***** Muons generated at Y = %f cm\n",fYOrigin);
163 printf("************************************************\n");
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
172 // Intersection between muon and cylinder [-z,+z] with radius r
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);
181 track.RelateToVertex(&origin,fBkG,10000.);
183 Float_t d0z0[2],covd0z0[3];
184 track.GetImpactParameters(d0z0,covd0z0);
187 if(TMath::Abs(d0z0[0])>r) return kFALSE;
189 if(TMath::Abs(d0z0[1])>z) return kFALSE;
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;
204 //-----------------------------------------------------------------------------