EVGEN/AliGenFunction.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 // Generator for particles according generic functions
  17 //
  18 //  TF1 *   fFMomentum;           // momentum distribution function inGeV
  19 //  TF1 *   fFPhi;                // phi distribution function in rad
  20 //  TF1 *   fFTheta;              // theta distribution function in rad
  21 //  TF3 *   fFPosition;           // position distribution function in cm
  22 //  TF1 *   fFPdg;                // pdg distribution function
  23 //  We assume that the moment, postion and PDG code of particles are independent
  24 //  Only tracks/particle crossing the reference radius at given z range
  25 //
  26 // Origin: marian.ivanov@cern.ch
  27
  28
  29 /*
  30   Example generic function for cosmic generation:
  31
  32 //
  33         AliGenFunction *generCosmic = new AliGenFunction;
  34         generCosmic->SetNumberParticles(100);
  35         TF1 *fmom   = new TF1("fmom","TMath::Exp(-x/8)",0,30);
  36         //
  37         TF1 *fphi   = new TF1("fphi","TMath::Gaus(x,-TMath::Pi()/2,0.3)",-3.14,3.14);
  38         TF1 *ftheta = new TF1("ftheta","TMath::Gaus(x,TMath::Pi()/2,0.3)",-3.14,3.14);
  39         TF3 *fpos = new TF3("fpos","1+(x+y+z)*0",-2000,2000,700,701,-2000,2000);
  40         TF1 * fpdg= new TF1("fpdg","1*(abs(x-13)<0.1)+1*(abs(x+13)<0.1)",-300,300);
  41         fpdg->SetNpx(10000);  // neccessary - number of bins for generation
  42         generCosmic->SetVertexSmear(kPerEvent);
  43         generCosmic->SetFunctions(fmom,fphi,ftheta,fpos,fpdg);
  44         generCosmic->SetCylinder(375,-375,375);
  45         generCosmic->SetBkG(0.2);
  46         gener->AddGenerator(generCosmic,"generCosmic",1);
  47
  48
  49
  50 */
  51
  52
  53
  54
  55 #include <TParticle.h>
  56 #include <TF1.h>
  57 #include <TF3.h>
  58 #include <TDatabasePDG.h>
  59
  60 #include "AliRun.h"
  61 #include "AliLog.h"
  62 #include "AliESDtrack.h"
  63 #include "AliESDVertex.h"
  64 #include "AliGenFunction.h"
  65 #include "AliGenEventHeader.h"
  66
  67 ClassImp(AliGenFunction)
  68
  69 //-----------------------------------------------------------------------------
  70 AliGenFunction::AliGenFunction():
  71   AliGenerator(),
  72   fBkG(0.),
  73   fFMomentum(0),           // momentum distribution function
  74   fFPhi(0),                // phi distribution function
  75   fFTheta(0),              // theta distribution function
  76   fFPosition(0),           // position distribution function
  77   fFPdg(0),                // pdg distribution function
  78   fRefRadius(0),           // reference radius to be crossed
  79   fZmin(0),                // minimal z at reference radius
  80   fZmax(0),                // z at reference radius
  81   fMaxTrial(10000)         // maximal number of attempts
  82 {
  83   //
  84   // Default constructor
  85   //
  86   SetNumberParticles(1);
  87 }
  88
  89 AliGenFunction::AliGenFunction(const AliGenFunction& func):
  90   AliGenerator(),
  91   fBkG(func.fBkG),
  92   fFMomentum(func.fFMomentum),           // momentum distribution function
  93   fFPhi(func.fFPhi),                     // phi distribution function
  94   fFTheta(func.fFTheta),                 // theta distribution function
  95   fFPosition(func.fFPosition),           // position distribution function
  96   fFPdg(func.fFPdg),                     // pdg distribution function
  97   fRefRadius(func.fRefRadius),           // reference radius to be crossed
  98   fZmin(func.fZmin),                     // minimal z at reference radius
  99   fZmax(func.fZmax),                     // z at reference radius
 100   fMaxTrial(10000)                       // maximal number of attempts
 101 {
 102     // Copy constructor
 103     SetNumberParticles(1);
 104 }
 105
 106 AliGenFunction & AliGenFunction::operator=(const AliGenFunction& func)
 107 {
 108     // Assigment operator
 109       if(&func == this) return *this;
 110       fBkG       = func.fBkG;
 111       fFMomentum = func.fFMomentum;
 112       fFPhi      = func.fFPhi;
 113       fFTheta    = func.fFTheta;
 114       fFPosition = func.fFPosition;
 115       fFPdg      = func.fFPdg;
 116       fRefRadius = func.fRefRadius;
 117       fZmin      = func.fZmin;
 118       fZmax      = func.fZmax;
 119       fMaxTrial  = func.fMaxTrial;
 120       return *this;
 121 }
 122
 123
 124 //-----------------------------------------------------------------------------
 125 void AliGenFunction::Generate()
 126 {
 127   //
 128   // Generate one muon
 129   //
 130   Int_t naccepted =0;
 131
 132   for (Int_t ipart=0; ipart<fMaxTrial && naccepted<fNpart; ipart++){
 133     //
 134     //
 135     //
 136     Float_t mom[3];
 137     Float_t posf[3];
 138     Double_t pos[3];
 139     Int_t pdg;
 140     Double_t ptot, pt,  phi, theta;
 141     //
 142     ptot     = fFMomentum->GetRandom();
 143     phi      = fFPhi->GetRandom();
 144     theta    = fFTheta->GetRandom();
 145     pt     = ptot*TMath::Sin(theta);
 146     mom[0] = pt*TMath::Cos(phi);
 147     mom[1] = pt*TMath::Sin(phi);
 148     mom[2] = ptot*TMath::Cos(theta);
 149
 150     //
 151     fFPosition->GetRandom3(pos[0],pos[1],pos[2]);
 152     posf[0]=pos[0];
 153     posf[1]=pos[1];
 154     posf[2]=pos[2];
 155     pdg = TMath::Nint(fFPdg->GetRandom());
 156     if (!IntersectCylinder(fRefRadius,fZmin, fZmax, pdg, posf, mom)) continue;
 157     //
 158     //
 159     Float_t polarization[3]= {0,0,0};
 160     Int_t nt;
 161     PushTrack(fTrackIt,-1,pdg,mom, posf, polarization,0,kPPrimary,nt);
 162     naccepted++;
 163   }
 164
 165   AliGenEventHeader* header = new AliGenEventHeader("THn");
 166   gAlice->SetGenEventHeader(header);
 167   return;
 168 }
 169 //-----------------------------------------------------------------------------
 170 void AliGenFunction::Init()
 171 {
 172   //
 173   // Initialisation, check consistency of selected ranges
 174   //
 175   printf("************ AliGenFunction ****************\n");
 176   printf("************************************************\n");
 177   if (!fFMomentum){
 178     AliFatal("Momentum distribution function not specified");
 179   }
 180   if (!fFPosition){
 181     AliFatal("Position distribution function not specified");
 182   }
 183   if (!fFPdg){
 184     AliFatal("PDG distribution function not specified");
 185   }
 186   if (fZmin==fZmax){
 187     AliFatal("Z range not specified");
 188   }
 189
 190   return;
 191 }
 192
 193 void AliGenFunction::SetFunctions(TF1 * momentum, TF1 *fphi, TF1 *ftheta,TF3 * position, TF1* pdg){
 194   //
 195   // Set the function
 196   //
 197   fFMomentum = momentum;
 198   fFPhi = fphi;
 199   fFTheta = ftheta;
 200   fFPosition = position;
 201   fFPdg = pdg;
 202 }
 203
 204 void AliGenFunction::SetCylinder(Double_t refR, Double_t zmin, Double_t zmax){
 205   //
 206   // Set the cylinder geometry
 207   //
 208   fRefRadius = refR;          // reference radius to be crossed
 209   fZmin = zmin;               // minimal z at reference radius
 210   fZmax = zmax;               // maximal z at reference radius
 211
 212 }
 213
 214
 215 //-----------------------------------------------------------------------------
 216 Bool_t AliGenFunction::IntersectCylinder(Float_t r,Float_t zmin, Float_t zmax,Int_t pdg,
 217                                          Float_t o[3],Float_t p[3]) const
 218 {
 219   //
 220   // Intersection between muon and cylinder [-z,+z] with radius r
 221   //
 222   Double_t mass=0;
 223   if (TDatabasePDG::Instance()->GetParticle(pdg)){
 224     mass = TDatabasePDG::Instance()->GetParticle(13)->Mass();
 225   }
 226
 227   Float_t en = TMath::Sqrt(mass*mass+p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
 228   TParticle part(pdg,0,0,0,0,0,p[0],p[1],p[2],en,o[0],o[1],o[2],0);
 229   AliESDtrack track(&part);
 230   Double_t pos[3]={0.,0.,0.},sigma[3]={0.,0.,0.};
 231   AliESDVertex origin(pos,sigma);
 232
 233   track.RelateToVertex(&origin,fBkG,10000.);
 234
 235   Float_t d0z0[2],covd0z0[3];
 236   track.GetImpactParameters(d0z0,covd0z0);
 237
 238   // check rphi
 239   if(TMath::Abs(d0z0[0])>r) return kFALSE;
 240   // check z
 241   if(d0z0[1]>zmax) return kFALSE;
 242   if(d0z0[1]<zmin) return kFALSE;
 243   //
 244   return kTRUE;
 245 }