]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDsim.cxx
Avoiding relative path to include files
[u/mrichter/AliRoot.git] / TRD / AliTRDsim.cxx
index 82b68e25fd7e450f97de2aac7ee66d9a5d721dec..8dbc34f2b150c838d817700486e214aa76683095 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.7  2000/12/20 13:00:45  cblume
-Modifications for the HP-compiler
-
-Revision 1.6  2000/12/12 10:20:10  cblume
-Initialize fSepctrum = 0 in ctors
-
-Revision 1.5  2000/10/15 23:40:01  cblume
-Remove AliTRDconst
-
-Revision 1.4  2000/10/06 16:49:46  cblume
-Made Getters const
-
-Revision 1.3.2.1  2000/09/18 13:45:30  cblume
-New class AliTRDsim that simulates TR photons
-
-Revision 1.2  1999/09/29 09:24:35  fca
-Introduction of the Copyright and cvs Log
-
-*/
+/* $Id$ */
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
@@ -44,6 +24,7 @@ Introduction of the Copyright and cvs Log
 //   17.07.1998 - A.Andronic                                                 //
 //   08.12.1998 - simplified version                                         //
 //   11.07.2000 - Adapted code to aliroot environment (C.Blume)              //
+//   04.06.2004 - Momentum dependent parameters implemented (CBL)            //
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
@@ -69,6 +50,8 @@ AliTRDsim::AliTRDsim():TObject()
 
   fSpectrum = 0;
   fSigma    = 0;
+  fNFoils   = 0;
+  fNFoilsUp = 0;
 
   Init();
 
@@ -91,6 +74,8 @@ AliTRDsim::AliTRDsim(AliModule *mod, Int_t foil, Int_t gap)
 
   fSpectrum = 0;
   fSigma    = 0;
+  fNFoils   = 0;
+  fNFoilsUp = 0;
 
   Init();
 
@@ -110,7 +95,7 @@ AliTRDsim::AliTRDsim(AliModule *mod, Int_t foil, Int_t gap)
 }
 
 //_____________________________________________________________________________
-AliTRDsim::AliTRDsim(const AliTRDsim &s)
+AliTRDsim::AliTRDsim(const AliTRDsim &s):TObject(s)
 {
   //
   // AliTRDsim copy constructor
@@ -127,8 +112,10 @@ AliTRDsim::~AliTRDsim()
   // AliTRDsim destructor
   //
 
-  if (fSpectrum) delete fSpectrum;
-  if (fSigma)    delete fSigma;
+  //  if (fSpectrum) delete fSpectrum;
+  if (fSigma)    delete [] fSigma;
+  if (fNFoils)   delete [] fNFoils;
+  if (fNFoilsUp) delete [] fNFoilsUp;
 
 }
 
@@ -151,7 +138,6 @@ void AliTRDsim::Copy(TObject &s)
   // Copy function
   //
 
-  ((AliTRDsim &) s).fNFoils     = fNFoils;
   ((AliTRDsim &) s).fFoilThick  = fFoilThick;
   ((AliTRDsim &) s).fFoilDens   = fFoilDens;
   ((AliTRDsim &) s).fFoilOmega  = fFoilOmega;
@@ -169,7 +155,19 @@ void AliTRDsim::Copy(TObject &s)
   ((AliTRDsim &) s).fSpLower    = fSpLower;
   ((AliTRDsim &) s).fSpUpper    = fSpUpper;
 
-  if (((AliTRDsim &) s).fSigma) delete ((AliTRDsim &) s).fSigma;
+  if (((AliTRDsim &) s).fNFoils) delete [] ((AliTRDsim &) s).fNFoils;
+  ((AliTRDsim &) s).fNFoils   = new Int_t[fNFoilsDim];
+  for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+    ((AliTRDsim &) s).fNFoils[iFoil] = fNFoils[iFoil];
+  }  
+
+  if (((AliTRDsim &) s).fNFoilsUp) delete [] ((AliTRDsim &) s).fNFoilsUp;
+  ((AliTRDsim &) s).fNFoilsUp = new Double_t[fNFoilsDim];
+  for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+    ((AliTRDsim &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
+  }  
+
+  if (((AliTRDsim &) s).fSigma)  delete [] ((AliTRDsim &) s).fSigma;
   ((AliTRDsim &) s).fSigma = new Double_t[fSpNBins];
   for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
     ((AliTRDsim &) s).fSigma[iBin] = fSigma[iBin];
@@ -184,23 +182,42 @@ void AliTRDsim::Init()
 {
   //
   // Initialization 
-  // The default radiator are 100 prolypropilene foils of 20 mu thickness
-  // with gaps of 500 mu filled with CO2.
-  //      
+  // The default radiator are prolypropilene foils of 10 mu thickness
+  // with gaps of 80 mu filled with N2.
   // 
 
-  fNFoils     = 100;
-
-  fFoilThick  = 0.0020;
+  fNFoilsDim   = 7;
+
+  if (fNFoils)   delete [] fNFoils;
+  fNFoils      = new Int_t[fNFoilsDim];
+  fNFoils[0]   = 170;
+  fNFoils[1]   = 250;
+  fNFoils[2]   = 310;
+  fNFoils[3]   = 380;
+  fNFoils[4]   = 430;
+  fNFoils[5]   = 490;
+  fNFoils[6]   = 550;
+
+  if (fNFoilsUp) delete [] fNFoilsUp;
+  fNFoilsUp    = new Double_t[fNFoilsDim];
+  fNFoilsUp[0] = 1.25;
+  fNFoilsUp[1] = 1.75;
+  fNFoilsUp[2] = 2.50;
+  fNFoilsUp[3] = 3.50;
+  fNFoilsUp[4] = 4.50;
+  fNFoilsUp[5] = 5.50;
+  fNFoilsUp[6] = 10000.0;
+
+  fFoilThick  = 0.0013;
   fFoilDens   = 0.92;   
   fFoilZ      = 5.28571;
   fFoilA      = 10.4286;
   fFoilOmega  = Omega(fFoilDens,fFoilZ,fFoilA);
 
-  fGapThick   = 0.0500;
-  fGapDens    = 0.001977;  
-  fGapZ       = 7.45455;
-  fGapA       = 14.9091;
+  fGapThick   = 0.0060;
+  fGapDens    = 0.00125;  
+  fGapZ       = 7.0;
+  fGapA       = 14.00674;
   fGapOmega   = Omega(fGapDens ,fGapZ ,fGapA );
 
   fTemp       = 293.16;
@@ -213,6 +230,7 @@ void AliTRDsim::Init()
 
   if (fSpectrum) delete fSpectrum;
   fSpectrum   = new TH1D("TRspectrum","TR spectrum",fSpNBins,fSpLower,fSpUpper);
+  fSpectrum->SetDirectory(0);
 
   // Set the sigma values 
   SetSigma();
@@ -255,16 +273,14 @@ Int_t AliTRDsim::CreatePhotons(Int_t pdg, Float_t p
     break;
   };
 
-  // Calculate gamma
-  Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
-
   // Calculate the TR photons
-  return TrPhotons(gamma, nPhoton, ePhoton);
+  return TrPhotons(p, mass, nPhoton, ePhoton);
 
 }
 
 //_____________________________________________________________________________
-Int_t AliTRDsim::TrPhotons(Double_t gamma, Int_t &nPhoton, Float_t *ePhoton)
+Int_t AliTRDsim::TrPhotons(Float_t p, Float_t mass
+                         , Int_t &nPhoton, Float_t *ePhoton)
 {
   //
   // Produces TR photons.
@@ -275,6 +291,12 @@ Int_t AliTRDsim::TrPhotons(Double_t gamma, Int_t &nPhoton, Float_t *ePhoton)
 
   Double_t kappa = fGapThick / fFoilThick;
 
+  // Calculate gamma
+  Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
+
+  // Select the number of foils corresponding to momentum
+  Int_t    foils = SelectNFoils(p);
+
   fSpectrum->Reset();
 
   // The TR spectrum
@@ -303,7 +325,7 @@ Int_t AliTRDsim::TrPhotons(Double_t gamma, Int_t &nPhoton, Float_t *ePhoton)
     }
 
     // Absorbtion
-    Double_t conv      = 1.0 - TMath::Exp(-fNFoils * fSigma[iBin]);
+    Double_t conv      = 1.0 - TMath::Exp(-foils * fSigma[iBin]);
 
     // eV -> keV
     Float_t  energykeV = energyeV * 0.001;
@@ -337,7 +359,7 @@ void AliTRDsim::SetSigma()
   // Sets the absorbtion crosssection for the energies of the TR spectrum
   //
 
-  if (fSigma) delete fSigma;
+  if (fSigma) delete [] fSigma;
   fSigma = new Double_t[fSpNBins];
   for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
     Double_t energykeV = iBin * fSpBinWidth + 1.0;
@@ -354,14 +376,11 @@ Double_t AliTRDsim::Sigma(Double_t energykeV)
   // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
   //
 
-  // Gas at 0 C
-  const Double_t kTemp0 = 273.16;
-
   // keV -> MeV
   Double_t energyMeV = energykeV * 0.001;
   if (energyMeV >= 0.001) {
-    return(GetMuPo(energyMeV) * fFoilDens * fFoilThick + 
-           GetMuCO(energyMeV) * fGapDens  * fGapThick  * fTemp/kTemp0);
+    return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
+           GetMuAi(energyMeV) * fGapDens  * fGapThick  * GetTemp());
   }
   else {
     return 1e6;
@@ -689,6 +708,50 @@ Double_t AliTRDsim::GetMuHe(Double_t energyMeV)
 
 }
 
+//_____________________________________________________________________________
+Double_t AliTRDsim::GetMuAi(Double_t energyMeV)
+{
+  //
+  // Returns the photon absorbtion cross section for air
+  // Implemented by Oliver Busch
+  //
+
+  const Int_t kN = 38;
+
+  Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
+                      0.16143E+03, 0.14250E+03, 0.15722E+03,
+                      0.77538E+02, 0.40099E+02, 0.23313E+02,
+                      0.98816E+01, 0.51000E+01, 0.16079E+01,
+                      0.77536E+00, 0.35282E+00, 0.24790E+00,
+                      0.20750E+00, 0.18703E+00, 0.16589E+00,
+                      0.15375E+00, 0.13530E+00, 0.12311E+00,
+                      0.10654E+00, 0.95297E-01, 0.86939E-01,
+                      0.80390E-01, 0.70596E-01, 0.63452E-01,
+                      0.56754E-01, 0.51644E-01, 0.44382E-01,
+                      0.35733E-01, 0.30721E-01, 0.27450E-01,
+                      0.25171E-01, 0.22205E-01, 0.20399E-01,
+                      0.18053E-01, 0.18057E-01 };
+
+
+
+  Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
+                      0.30000E-02, 0.32029E-02, 0.32029E-02,
+                      0.40000E-02, 0.50000E-02, 0.60000E-02,
+                      0.80000E-02, 0.10000E-01, 0.15000E-01,
+                      0.20000E-01, 0.30000E-01, 0.40000E-01,
+                      0.50000E-01, 0.60000E-01, 0.80000E-01,
+                      0.10000E+00, 0.15000E+00, 0.20000E+00,
+                      0.30000E+00, 0.40000E+00, 0.50000E+00,
+                      0.60000E+00, 0.80000E+00, 0.10000E+01,
+                      0.12500E+01, 0.15000E+01, 0.20000E+01,
+                      0.30000E+01, 0.40000E+01, 0.50000E+01,
+                      0.60000E+01, 0.80000E+01, 0.10000E+02,
+                      0.15000E+02, 0.20000E+02 };
+
+  return Interpolate(energyMeV,en,mu,kN);
+
+}
+
 //_____________________________________________________________________________
 Double_t AliTRDsim::Interpolate(Double_t energyMeV
                               , Double_t *en, Double_t *mu, Int_t n)
@@ -741,3 +804,23 @@ Int_t AliTRDsim::Locate(Double_t *xv, Int_t n, Double_t xval
   return 0;
 
 }
+
+//_____________________________________________________________________________
+Int_t AliTRDsim::SelectNFoils(Float_t p)
+{
+  //
+  // Selects the number of foils corresponding to the momentum
+  //
+
+  Int_t foils = fNFoils[fNFoilsDim-1];
+
+  for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
+    if (p < fNFoilsUp[iFoil]) {
+      foils = fNFoils[iFoil];
+      break;
+    }
+  }
+
+  return foils;
+
+}