First version, generation of cosmic muons on the surface
authorgamez <gamez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 7 Oct 2002 11:25:47 +0000 (11:25 +0000)
committergamez <gamez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 7 Oct 2002 11:25:47 +0000 (11:25 +0000)
CRT/AliGenCRT.cxx [new file with mode: 0644]
CRT/AliGenCRT.h [new file with mode: 0644]

diff --git a/CRT/AliGenCRT.cxx b/CRT/AliGenCRT.cxx
new file mode 100644 (file)
index 0000000..915c925
--- /dev/null
@@ -0,0 +1,548 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Log$
+
+*/
+
+/////////////////////////////////////////////////////////////////////////////
+//
+//  Contain parametrizations to generate atmospheric muons, and also
+//  to generate single muons and muon bundles at surface level.
+//
+//Begin_Html
+/*
+<img src="picts/AliGenCRTClass.gif">
+</pre>
+<br clear=left>
+<font size=+2 color=red>
+<p>The responsible person for this module is
+<a href="mailto:Enrique.Gamez.Flores@cern.ch">Enrique Gamez</a>.
+</font>
+<pre>
+*/
+//End_Html
+//
+/////////////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+
+#include "AliRun.h"
+#include "AliConst.h"
+#include "AliPDG.h"
+#include "AliMCProcess.h"
+
+#include "AliGenCRT.h"
+
+ClassImp(AliGenCRT)
+
+//_____________________________________________________________________________
+AliGenCRT::AliGenCRT() : AliGenerator(-1)
+{
+  //
+  // Default ctor.
+  //
+  fIpart = 0;
+
+  fXwidth=0;
+  fNx=1;
+  fZwidth=0;
+  fNz=1;
+  fMuonGrid = kFALSE;
+
+
+  // Set the origin above the vertex, on the surface.
+  fOrigin[0] = 0.;
+  fOrigin[1] = 0.;
+  fOrigin[2] = 0.;
+
+  fZenithMin = 0.; 
+  fZenithMax = 0.;
+
+  fAzimuthMin = 0.;
+  fAzimuthMax = 0.;
+
+  fPResolution = 1.;
+  fPRange = 0.;
+
+  fMomentumDist = 0;
+  fZenithDist = 0;
+
+  fPDist = 0;
+
+  fAp = 0;
+  fUnfoldedMomentumDist = 0;
+
+  fCRModeName = 0;
+}
+
+//_____________________________________________________________________________
+AliGenCRT::AliGenCRT(Int_t npart) 
+  : AliGenerator(npart)
+{
+  //
+  // Standard ctor.
+  //
+  fName = "CRT";
+  fTitle = "Cosmic Muons generator";
+
+  // Generate Muon- by default
+  fIpart = kMuonMinus;
+
+  fXwidth=0;
+  fNx=1;
+  fZwidth=0;
+  fNz=1;
+  fMuonGrid = kFALSE;
+
+  // Set the origin above the vertex, on the surface.
+  fOrigin[0] = 0.;
+  fOrigin[1] = AliCRTConstants::fgDepth; // At the surface by default.
+  fOrigin[2] = 0.;
+
+  fZenithMin = 0.; // Generate veritcals by default.
+  fZenithMax = 0.;
+
+  fAzimuthMin = 0.;
+  fAzimuthMax = 0.;
+
+  fPResolution = 1.; // 1 GeV by default.
+  fPRange = 0.;      // 0 GeV by default.
+
+  fMomentumDist = 0;
+  fZenithDist = 0;
+
+  fPDist = 0;
+
+  fAp = 0;
+  fUnfoldedMomentumDist = 0;
+
+  fCRModeName = 0;
+}
+
+//_____________________________________________________________________________
+AliGenCRT::AliGenCRT(const AliGenCRT& gen) : AliGenerator()
+{
+  //
+  // Copy ctor.
+  //
+  gen.Copy(*this);
+}
+
+//_____________________________________________________________________________
+AliGenCRT& AliGenCRT::operator= (const AliGenCRT& gen)
+{
+  //
+  // Asingment ctor.
+  //
+  gen.Copy(*this);
+  return *this;
+}
+
+//_____________________________________________________________________________
+AliGenCRT::~AliGenCRT()
+{
+  //
+  // Default dtor.
+  //
+  if ( fAp )           delete fAp;
+  if ( fMomentumDist ) delete fMomentumDist;
+  if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
+  if ( fZenithDist )   delete fZenithDist;
+  if ( fPDist )        delete fPDist;
+}
+
+//_____________________________________________________________________________
+void AliGenCRT::Generate()
+{
+  //
+  // Generate on one trigger
+  //
+
+  Float_t polar[3]= {0,0,0}; // Polarization parameters
+  //
+  Float_t origin[3];
+  Float_t p[3];
+  Int_t i, j, nt;
+  Float_t pmom, pt;
+  Float_t zenith, azimuth;
+  //
+  Float_t random[6];
+
+
+  // Check if you need to use a distribution function for the momentum
+  if ( fMomentumDist ) {
+    pmom = - this->GetMomentum(); // Always downwards.
+
+  } else {
+    pmom = -fPMin;
+
+  }
+  
+  zenith = this->GetZenithAngle(pmom);  // In degrees
+  
+  // Take the azimuth random.
+  Rndm(random, 2);
+  azimuth = (fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0]); // In degrees
+  
+  origin[0] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
+    TMath::Sin(azimuth*kDegrad);
+  origin[1] = AliCRTConstants::fgDepth;
+  origin[2] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
+    TMath::Cos(azimuth*kDegrad);
+  
+  if ( fVertexSmear == kPerEvent ) {
+    Rndm(random,6);
+    for (j=0;j<3;j++) {
+      if ( j == 1 ) {
+       // Don't smear the vertical position.
+       origin[j] = AliCRTConstants::fgDepth;
+      } else {
+       origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
+         TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
+      }
+    }
+  }
+
+
+  if ( fCRMode == kSingleMuons ) {
+
+    //
+    // Generate kinematics a la AliGenBox
+    //
+
+    for(i=0;i<fNpart;i++) {
+      Rndm(random,3);
+
+      if(TestBit(kMomentumRange)) {
+       pmom = -( fPMin + random[1]*(fPMax - fPMin) ); // always downwards.
+       pt=pmom*TMath::Sin(zenith*kDegrad);
+      } else {
+
+       pt= -(fPtMin+random[1]*(fPtMax-fPtMin)); // always downwards.
+       pmom=pt/TMath::Sin(zenith*kDegrad);
+      }
+      
+      p[0] = pt*TMath::Sin(azimuth*kDegrad);
+      p[1] = pmom*TMath::Cos(zenith*kDegrad);
+      p[2] = pt*TMath::Cos(azimuth*kDegrad);
+      
+      if(fVertexSmear==kPerTrack) {
+       Rndm(random,6);
+       for (j=0;j<3;j++) {
+         if ( j == 1 ) {
+           origin[j] = AliCRTConstants::fgDepth;
+         } else {
+           origin[j]=fOrigin[j]+fOsigma[j]*
+             TMath::Cos(2*random[2*j]*TMath::Pi())*
+             TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
+         }
+       }
+      }
+
+    }
+    
+    // End of generation a la AliGenBox
+    
+  } else if ( fCRMode == kMuonBundle ) {
+
+    //
+    // Generate kinematics a la AliGenScan
+    //
+    Float_t dx,dz;
+    //
+    if (fNx > 0) {
+      dx=fXwidth/fNx;
+    } else {
+      dx=1e10;
+    }
+    
+    if (fNz > 0) {
+      dz=fZwidth/fNz;
+    } else {
+      dz=1e10;
+    }
+
+    for (Int_t ix=0; ix<fNx; ix++) {
+      for (Int_t iz=0; iz<fNz; iz++){
+       Rndm(random,6);
+       origin[0]+=ix*dx+2*(random[0]-0.5)*fOsigma[0];
+       origin[2]+=iz*dz+2*(random[1]-0.5)*fOsigma[2];       
+
+       pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downward.
+       p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
+       p[1] = TMath::Cos(zenith*kDegrad)*pmom;
+       p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
+       
+      }
+    }
+
+    // End of generation a la AliGenScan
+
+  } else if ( fCRMode == kMuonFlux ) {
+    
+    // Always generate the things downwards.
+    p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
+    p[1] = TMath::Cos(zenith*kDegrad)*pmom;
+    p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
+
+  } else {
+    Error("Generate", "Generation Mode unknown!\n");
+  }
+
+  // Put the track on the stack.
+  SetTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
+
+
+}
+
+//_____________________________________________________________________________
+void AliGenCRT::Init()
+{
+  //
+  // Initialize some internal methods.
+  //
+
+  // Determine some specific data members.
+  fPRange = TMath::Abs(fPMax-fPMin);
+
+  if ( fCRMode == kSingleMuons ) {
+    fCRModeName = new TString("Single Muons");
+    // Initialisation, check consistency of selected ranges
+    if(TestBit(kPtRange)&&TestBit(kMomentumRange)) 
+      Fatal("Init","You should not set the momentum range and the pt range!");
+    
+    if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange))) 
+      Fatal("Init","You should set either the momentum or the pt range!");
+    
+  } else if ( fCRMode == kMuonBundle ) {
+    fCRModeName = new TString("Muon Bundles");
+
+  } else if ( fCRMode == kMuonFlux ) {
+    fCRModeName = new TString("Muon Fluxes");
+    // Initialize the ditribution functions.
+    this->InitMomentumGeneration();
+    this->InitZenithalAngleGeneration();
+    
+  } else {
+    Fatal("Generate", "Generation Mode unknown!\n");
+
+  }
+
+}
+
+//____________________________________________________________________________
+void AliGenCRT::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
+{
+  //
+  // Define the grid
+  // This data shuold be used for Muon bundles generation.
+  //
+  fXwidth=xwidth;
+  fNx=nx;
+  fZwidth=zwidth;
+  fNz=nz;
+
+  // Print a message  about the use, if the Mode has not been set, or
+  // it has to a different Mode.
+  if ( fCRMode != kMuonBundle ) {
+    Warning("SetRange","You have been specified a grid to generate muon bundles, but seems that you haven't choose this generation mode, or you have already select a different one");
+    fMuonGrid = kTRUE;
+  }
+}
+
+//____________________________________________________________________________
+void AliGenCRT::InitApWeightFactors()
+{
+  //
+  // This factors will be  to correct the zenithal angle distribution
+  // acording the momentum
+
+  //
+  // Fill the array for the flux zenith angle dependence.
+  // at the index 0 of fAp[] will be the "factor" if we have a muon
+  // of 0 GeV.
+  Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
+  Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
+
+  // Get the information from the momentum
+  Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
+  // Initialize the Array of floats to hold the a(p) factors.
+  fAp = new TArrayF(pEnd);
+  
+  Int_t index = 0;
+
+  for (Int_t i = 0; i < pEnd; i++ ) {
+    Float_t currentP = ((Float_t)i)*fPResolution;
+    if ( currentP < p[1] )                          index = 0;
+    else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
+    else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
+    else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
+    else if ( currentP >= p[4] )                    index = 4;
+
+    Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
+                 (p[index+1] - p[index]) + a[index];
+    fAp->AddAt(ap, i);
+  }
+
+}
+
+//___________________________________________________________________________
+void AliGenCRT::InitMomentumGeneration()
+{
+  //
+  // Initialize a funtion to generate the momentum randomly
+  // acording this function.
+  //
+
+  // Check if we nned to initialize the function
+  if ( fPMin != fPMax ) {
+
+    // Check also if the function have been defined yet.
+    if ( !fMomentumDist ) {
+
+      // If not, use this function
+      const char* y      = "log10(x)";
+      
+      const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
+      const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
+      const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
+      const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
+      
+      const char* h = "%s + %s + %s + %s";
+      const char* flux = "pow(10., %s)";
+      const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
+      const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
+      
+      char buffer1[1024];
+      char buffer2[1024];
+      char buffer3[1024];
+      char buffer4[1024];
+      char buffer5[1024];
+      char buffer6[1024];
+      char buffer7[1024];
+
+      sprintf(buffer1, h1Coef, y, y, y, y, y, y);
+      sprintf(buffer2, h2Coef, y, y, y, y, y, y);
+      sprintf(buffer3, h3Coef, y, y, y, y, y, y);
+      sprintf(buffer4, s2Coef, y, y, y, y, y, y);
+      
+      sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
+      
+      sprintf(buffer6, flux, buffer5);
+      
+      fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
+      sprintf(buffer7, normalizedFlux, buffer5);
+      fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
+      for (Int_t i = 0; i < 4; i++ ) {
+       fMomentumDist->SetParName(i, paramNames[i]);
+       fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
+      }
+      
+      fMomentumDist->SetParameter(0, 0.133);
+      fMomentumDist->SetParameter(1, -2.521);
+      fMomentumDist->SetParameter(2, -5.78);
+      fMomentumDist->SetParameter(3, -2.11);
+
+      fUnfoldedMomentumDist->SetParameter(0, 0.133);
+      fUnfoldedMomentumDist->SetParameter(1, -2.521);
+      fUnfoldedMomentumDist->SetParameter(2, -5.78);
+      fUnfoldedMomentumDist->SetParameter(3, -2.11);
+      
+    }
+
+  }
+
+}
+
+//____________________________________________________________________________
+void AliGenCRT::InitZenithalAngleGeneration()
+{
+  //
+  // Initalize a distribution function for the zenith angle.
+  // This angle will be obtained randomly acording this function.
+  // The generated angles  will been in degrees.
+
+  // Check if we need to create the function.
+  if ( fZenithMin != fZenithMax ) {
+
+    // Check also if another function have been defined.
+    if ( !fZenithDist ) {
+      
+      // initialize the momentum dependent coefficients, a(p) 
+      InitApWeightFactors();
+      
+      // Define the standard function.
+      char* zenithalDisributionFunction = "1 + [0]*(1 - cos(x*3.14159265358979312/180))";
+
+      Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
+      fPDist = new TClonesArray("TF1", pEnd);
+      TClonesArray &angle = *fPDist;
+      for ( Int_t i = 0; i < pEnd; i++ ) {
+       // Fill the distribution
+       TF1* zenith = new(angle[i]) TF1("zenith",zenithalDisributionFunction, fZenithMin, fZenithMax);
+
+       // Fix the momentum dependent coefficients
+       zenith->SetParName(0, "a(p)");
+       zenith->SetParameter(0, fAp->At(i));
+
+      }
+
+    } 
+
+  }
+
+}
+
+//____________________________________________________________________________
+const Float_t AliGenCRT::GetZenithAngle(Float_t mom)
+{
+
+  Float_t zenith = 0.;
+  // Check if you need to generate a constant zenith angle.
+  if ( !fZenithDist ) {
+    // Check if you have defined an array of momentum functions
+    if ( fPDist ) {
+      Int_t pIndex = TMath::Abs(TMath::Nint(mom));
+      TF1* zenithAngle = (TF1*)fPDist->UncheckedAt(pIndex);
+      zenith = zenithAngle->GetRandom();
+      return zenith;
+    } else {
+
+      if ( fCRMode != kMuonFlux ) {
+       // If you aren't generating muons obeying any ditribution
+       // only generate a flat zenith angle, acording the input settings
+       Float_t random[2];
+       Rndm(random, 2);
+       zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
+
+      } else {
+       // Even if you are generating muons acording some distribution,
+       // but you don't want to ...
+       zenith = fZenithMin;
+
+      }
+
+    }
+  } else {
+    zenith = fZenithDist->GetRandom();
+  }
+     
+  return zenith;
+  
+}
+
+//____________________________________________________________________________
diff --git a/CRT/AliGenCRT.h b/CRT/AliGenCRT.h
new file mode 100644 (file)
index 0000000..fd51220
--- /dev/null
@@ -0,0 +1,84 @@
+#ifndef ALIGENCRT_H
+#define ALIGENCRT_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+#include <TF2.h>
+#include <TClonesArray.h>
+
+#include "AliGenerator.h"
+#include "AliCRTConstants.h"
+
+class AliGenCRT : public AliGenerator {
+
+public:
+  AliGenCRT();
+  AliGenCRT(Int_t npart);
+  AliGenCRT(const AliGenCRT& gen);
+  AliGenCRT& operator= (const AliGenCRT& gen);
+  virtual ~AliGenCRT();
+
+  // Overloaded methods.
+  virtual void Init();
+  virtual void Generate();
+  virtual void SetPart(Int_t part) {fIpart = part;}
+
+  // Local methods.
+  void SetMode(ECRMode mode) {fCRMode = mode;}
+  const TString*  GetMode() const {return fCRModeName;}
+  
+  void SetZenithalAngleRange(Float_t min,Float_t max=0) {fZenithMin=min;fZenithMax=max;}
+  void SetAzimuthalAngleRange(Float_t min, Float_t max=0) {fAzimuthMin=min;fAzimuthMax=max;}
+  
+  void SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth);
+  const Float_t GetMomentumResolution() const {return fPResolution;}
+  
+  void SetMomentumDistrubutionFunction(TF1 *func) {fMomentumDist=func;}
+  void SetZenithalDistributionFunction(TF1 *func) {fZenithDist = func;}
+  void SetMomentumResolution(Float_t res=1.) {fPResolution=res;}
+  
+  const Float_t GetMomentum() {return fMomentumDist->GetRandom();}
+  const Float_t GetZenithAngle(Float_t mom=0.);
+
+  // The following methods are for testing pourpuses
+  TF1* GetMomentumDistibution() {return fMomentumDist;}
+  TF1* GetUnfoldedDistribution() {return fUnfoldedMomentumDist;}
+
+protected:
+  void InitApWeightFactors();
+  void InitMomentumGeneration();
+  void InitZenithalAngleGeneration();
+
+  Int_t    fIpart;              // Particle type.
+  ECRMode  fCRMode;             // Cosmic muons generation method flag
+  TString* fCRModeName;         // Cosmic muons generation mode name
+
+  Float_t  fXwidth;             // X width of the grid
+  Int_t    fNx;                 // Number of divisions in x
+  Float_t  fZwidth;             // Z widht of the  grid
+  Int_t    fNz;                 // Number of divisions in z
+  Bool_t   fMuonGrid;           // Flag for method (Muon-bundles) checkout
+
+  Float_t  fZenithMin;          // Minimum zenithal angle.
+  Float_t  fZenithMax;          // Maximum zenithal angle.
+
+  Float_t  fAzimuthMin;         // Minimum azimuthal angle.
+  Float_t  fAzimuthMax;         // Maximum azimuthal angle.
+
+  Float_t  fPRange;             // Cosmic muon momentum range width in GeVs.
+  Float_t  fPResolution;        // Momentum resolution in GeVs.
+
+  TArrayF* fAp;                 // a(p) weught factors for the ang. dist.
+
+  TF1*     fMomentumDist;       // Function to generate the momentum dist.
+  TF1*     fUnfoldedMomentumDist;
+  TF1*     fZenithDist;         // Function to generate the zenith angle dist.
+
+  TClonesArray* fPDist;         // Array of fZenithDist, to be used by a(p).
+
+private:
+  ClassDef(AliGenCRT, 0) // Generator for AliCRT class
+};
+#endif // ALIGENCRT_H