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 //-----------------------------------------------------------------------------