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