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