TRD/AliTRDsim.cxx

   1
   2 /**************************************************************************
   3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *                                                                        *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  **************************************************************************/
  16
  17 /* $Id$ */
  18
  19 ///////////////////////////////////////////////////////////////////////////////
  20 //                                                                           //
  21 //  TRD simulation - multimodule (regular rad.)                              //
  22 //  after: M. CASTELLANO et al., COMP. PHYS. COMM. 51 (1988) 431             //
  23 //                             + COMP. PHYS. COMM. 61 (1990) 395             //
  24 //                                                                           //
  25 //   17.07.1998 - A.Andronic                                                 //
  26 //   08.12.1998 - simplified version                                         //
  27 //   11.07.2000 - Adapted code to aliroot environment (C.Blume)              //
  28 //   04.06.2004 - Momentum dependent parameters implemented (CBL)            //
  29 //                                                                           //
  30 ///////////////////////////////////////////////////////////////////////////////
  31
  32 #include <stdlib.h>
  33
  34 #include <TH1.h>
  35 #include <TRandom.h>
  36 #include <TMath.h>
  37 #include <TParticle.h>
  38
  39 #include "AliModule.h"
  40 #include "AliLog.h"
  41
  42 #include "AliTRDsim.h"
  43
  44 ClassImp(AliTRDsim)
  45
  46 //_____________________________________________________________________________
  47 AliTRDsim::AliTRDsim()
  48   :TObject()
  49   ,fNFoilsDim(0)
  50   ,fNFoils(0)
  51   ,fNFoilsUp(0)
  52   ,fFoilThick(0)
  53   ,fGapThick(0)
  54   ,fFoilDens(0)
  55   ,fGapDens(0)
  56   ,fFoilOmega(0)
  57   ,fGapOmega()
  58   ,fFoilZ(0)
  59   ,fGapZ(0)
  60   ,fFoilA(0)
  61   ,fGapA(0)
  62   ,fTemp(0)
  63   ,fSpNBins(0)
  64   ,fSpRange(0)
  65   ,fSpBinWidth(0)
  66   ,fSpLower(0)
  67   ,fSpUpper(0)
  68   ,fSigma(0)
  69   ,fSpectrum(0)
  70 {
  71   //
  72   // AliTRDsim default constructor
  73   //
  74
  75   Init();
  76
  77 }
  78
  79 //_____________________________________________________________________________
  80 AliTRDsim::AliTRDsim(AliModule *mod, Int_t foil, Int_t gap)
  81   :TObject()
  82   ,fNFoilsDim(0)
  83   ,fNFoils(0)
  84   ,fNFoilsUp(0)
  85   ,fFoilThick(0)
  86   ,fGapThick(0)
  87   ,fFoilDens(0)
  88   ,fGapDens(0)
  89   ,fFoilOmega(0)
  90   ,fGapOmega()
  91   ,fFoilZ(0)
  92   ,fGapZ(0)
  93   ,fFoilA(0)
  94   ,fGapA(0)
  95   ,fTemp(0)
  96   ,fSpNBins(0)
  97   ,fSpRange(0)
  98   ,fSpBinWidth(0)
  99   ,fSpLower(0)
 100   ,fSpUpper(0)
 101   ,fSigma(0)
 102   ,fSpectrum(0)
 103 {
 104   //
 105   // AliTRDsim constructor. Takes the material properties of the radiator
 106   // foils and the gas in the gaps from AliModule <mod>.
 107   // The default number of foils is 100 with a thickness of 20 mu. The
 108   // thickness of the gaps is 500 mu.
 109   //
 110
 111   Float_t aFoil;
 112   Float_t zFoil;
 113   Float_t rhoFoil;
 114
 115   Float_t aGap;
 116   Float_t zGap;
 117   Float_t rhoGap;
 118
 119   Float_t rad;
 120   Float_t abs;
 121
 122   Char_t  name[21];
 123
 124   Init();
 125
 126   mod->AliGetMaterial(foil,name,aFoil,zFoil,rhoFoil,rad,abs);
 127   mod->AliGetMaterial(gap ,name,aGap ,zGap ,rhoGap ,rad,abs);
 128
 129   fFoilDens  = rhoFoil;
 130   fFoilA     = aFoil;
 131   fFoilZ     = zFoil;
 132   fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA);
 133
 134   fGapDens   = rhoGap;
 135   fGapA      = aGap;
 136   fGapZ      = zGap;
 137   fGapOmega  = Omega(fGapDens ,fGapZ ,fGapA );
 138
 139 }
 140
 141 //_____________________________________________________________________________
 142 AliTRDsim::AliTRDsim(const AliTRDsim &s)
 143   :TObject(s)
 144   ,fNFoilsDim(s.fNFoilsDim)
 145   ,fNFoils(0)
 146   ,fNFoilsUp(0)
 147   ,fFoilThick(s.fFoilThick)
 148   ,fGapThick(s.fGapThick)
 149   ,fFoilDens(s.fFoilDens)
 150   ,fGapDens(s.fGapDens)
 151   ,fFoilOmega(s.fFoilOmega)
 152   ,fGapOmega(s.fGapOmega)
 153   ,fFoilZ(s.fFoilZ)
 154   ,fGapZ(s.fGapZ)
 155   ,fFoilA(s.fFoilA)
 156   ,fGapA(s.fGapA)
 157   ,fTemp(s.fTemp)
 158   ,fSpNBins(s.fSpNBins)
 159   ,fSpRange(s.fSpRange)
 160   ,fSpBinWidth(s.fSpBinWidth)
 161   ,fSpLower(s.fSpLower)
 162   ,fSpUpper(s.fSpUpper)
 163   ,fSigma(0)
 164   ,fSpectrum(0)
 165 {
 166   //
 167   // AliTRDsim copy constructor
 168   //
 169
 170   if (((AliTRDsim &) s).fNFoils) {
 171     delete [] ((AliTRDsim &) s).fNFoils;
 172   }
 173   ((AliTRDsim &) s).fNFoils   = new Int_t[fNFoilsDim];
 174   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
 175     ((AliTRDsim &) s).fNFoils[iFoil]   = fNFoils[iFoil];
 176   }
 177
 178   if (((AliTRDsim &) s).fNFoilsUp) {
 179     delete [] ((AliTRDsim &) s).fNFoilsUp;
 180   }
 181   ((AliTRDsim &) s).fNFoilsUp = new Double_t[fNFoilsDim];
 182   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
 183     ((AliTRDsim &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
 184   }
 185
 186   if (((AliTRDsim &) s).fSigma) {
 187     delete [] ((AliTRDsim &) s).fSigma;
 188   }
 189   ((AliTRDsim &) s).fSigma    = new Double_t[fSpNBins];
 190   for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
 191     ((AliTRDsim &) s).fSigma[iBin]     = fSigma[iBin];
 192   }
 193
 194   fSpectrum->Copy(*((AliTRDsim &) s).fSpectrum);
 195
 196 }
 197
 198 //_____________________________________________________________________________
 199 AliTRDsim::~AliTRDsim()
 200 {
 201   //
 202   // AliTRDsim destructor
 203   //
 204
 205   if (fSigma) {
 206     delete [] fSigma;
 207     fSigma    = 0;
 208   }
 209
 210   if (fNFoils) {
 211     delete [] fNFoils;
 212     fNFoils   = 0;
 213   }
 214
 215   if (fNFoilsUp) {
 216     delete [] fNFoilsUp;
 217     fNFoilsUp = 0;
 218   }
 219
 220 }
 221
 222 //_____________________________________________________________________________
 223 AliTRDsim &AliTRDsim::operator=(const AliTRDsim &s)
 224 {
 225   //
 226   // Assignment operator
 227   //
 228
 229   if (this != &s) ((AliTRDsim &) s).Copy(*this);
 230
 231   return *this;
 232
 233 }
 234
 235 //_____________________________________________________________________________
 236 void AliTRDsim::Copy(TObject &s) const
 237 {
 238   //
 239   // Copy function
 240   //
 241
 242   ((AliTRDsim &) s).fFoilThick  = fFoilThick;
 243   ((AliTRDsim &) s).fFoilDens   = fFoilDens;
 244   ((AliTRDsim &) s).fFoilOmega  = fFoilOmega;
 245   ((AliTRDsim &) s).fFoilZ      = fFoilZ;
 246   ((AliTRDsim &) s).fFoilA      = fFoilA;
 247   ((AliTRDsim &) s).fGapThick   = fGapThick;
 248   ((AliTRDsim &) s).fGapDens    = fGapDens;
 249   ((AliTRDsim &) s).fGapOmega   = fGapOmega;
 250   ((AliTRDsim &) s).fGapZ       = fGapZ;
 251   ((AliTRDsim &) s).fGapA       = fGapA;
 252   ((AliTRDsim &) s).fTemp       = fTemp;
 253   ((AliTRDsim &) s).fSpNBins    = fSpNBins;
 254   ((AliTRDsim &) s).fSpRange    = fSpRange;
 255   ((AliTRDsim &) s).fSpBinWidth = fSpBinWidth;
 256   ((AliTRDsim &) s).fSpLower    = fSpLower;
 257   ((AliTRDsim &) s).fSpUpper    = fSpUpper;
 258
 259   if (((AliTRDsim &) s).fNFoils) {
 260     delete [] ((AliTRDsim &) s).fNFoils;
 261   }
 262   ((AliTRDsim &) s).fNFoils   = new Int_t[fNFoilsDim];
 263   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
 264     ((AliTRDsim &) s).fNFoils[iFoil]   = fNFoils[iFoil];
 265   }
 266
 267   if (((AliTRDsim &) s).fNFoilsUp) {
 268     delete [] ((AliTRDsim &) s).fNFoilsUp;
 269   }
 270   ((AliTRDsim &) s).fNFoilsUp = new Double_t[fNFoilsDim];
 271   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
 272     ((AliTRDsim &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
 273   }
 274
 275   if (((AliTRDsim &) s).fSigma) {
 276     delete [] ((AliTRDsim &) s).fSigma;
 277   }
 278   ((AliTRDsim &) s).fSigma    = new Double_t[fSpNBins];
 279   for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
 280     ((AliTRDsim &) s).fSigma[iBin]     = fSigma[iBin];
 281   }
 282
 283   fSpectrum->Copy(*((AliTRDsim &) s).fSpectrum);
 284
 285 }
 286
 287 //_____________________________________________________________________________
 288 void AliTRDsim::Init()
 289 {
 290   //
 291   // Initialization
 292   // The default radiator are prolypropilene foils of 10 mu thickness
 293   // with gaps of 80 mu filled with N2.
 294   //
 295
 296   fNFoilsDim   = 7;
 297
 298   if (fNFoils) {
 299     delete [] fNFoils;
 300   }
 301   fNFoils      = new Int_t[fNFoilsDim];
 302   fNFoils[0]   = 170;
 303   fNFoils[1]   = 225;
 304   fNFoils[2]   = 275;
 305   fNFoils[3]   = 305;
 306   fNFoils[4]   = 325;
 307   fNFoils[5]   = 340;
 308   fNFoils[6]   = 350;
 309
 310   if (fNFoilsUp) {
 311     delete [] fNFoilsUp;
 312   }
 313   fNFoilsUp    = new Double_t[fNFoilsDim];
 314   fNFoilsUp[0] = 1.25;
 315   fNFoilsUp[1] = 1.75;
 316   fNFoilsUp[2] = 2.50;
 317   fNFoilsUp[3] = 3.50;
 318   fNFoilsUp[4] = 4.50;
 319   fNFoilsUp[5] = 5.50;
 320   fNFoilsUp[6] = 10000.0;
 321
 322   fFoilThick  = 0.0013;
 323   fFoilDens   = 0.92;
 324   fFoilZ      = 5.28571;
 325   fFoilA      = 10.4286;
 326   fFoilOmega  = Omega(fFoilDens,fFoilZ,fFoilA);
 327
 328   fGapThick   = 0.0060;
 329   fGapDens    = 0.00125;
 330   fGapZ       = 7.0;
 331   fGapA       = 14.00674;
 332   fGapOmega   = Omega(fGapDens ,fGapZ ,fGapA );
 333
 334   fTemp       = 293.16;
 335
 336   fSpNBins    = 200;
 337   fSpRange    = 100;
 338   fSpBinWidth = fSpRange / fSpNBins;
 339   fSpLower    = 1.0 - 0.5 * fSpBinWidth;
 340   fSpUpper    = fSpLower + fSpRange;
 341
 342   if (fSpectrum) delete fSpectrum;
 343   fSpectrum   = new TH1D("TRspectrum","TR spectrum",fSpNBins,fSpLower,fSpUpper);
 344   fSpectrum->SetDirectory(0);
 345
 346   // Set the sigma values
 347   SetSigma();
 348
 349 }
 350
 351 //_____________________________________________________________________________
 352 Int_t AliTRDsim::CreatePhotons(Int_t pdg, Float_t p
 353                              , Int_t &nPhoton, Float_t *ePhoton)
 354 {
 355   //
 356   // Create TRD photons for a charged particle of type <pdg> with the total
 357   // momentum <p>.
 358   // Number of produced TR photons:       <nPhoton>
 359   // Energies of the produced TR photons: <ePhoton>
 360   //
 361
 362   // PDG codes
 363   const Int_t kPdgEle  =  11;
 364   const Int_t kPdgMuon =  13;
 365   const Int_t kPdgPion = 211;
 366   const Int_t kPdgKaon = 321;
 367
 368   Float_t  mass        = 0;
 369   switch (TMath::Abs(pdg)) {
 370   case kPdgEle:
 371     mass      =  5.11e-4;
 372     break;
 373   case kPdgMuon:
 374     mass      =  0.10566;
 375     break;
 376   case kPdgPion:
 377     mass      =  0.13957;
 378     break;
 379   case kPdgKaon:
 380     mass      =  0.4937;
 381     break;
 382   default:
 383     return 0;
 384     break;
 385   };
 386
 387   // Calculate the TR photons
 388   return TrPhotons(p, mass, nPhoton, ePhoton);
 389
 390 }
 391
 392 //_____________________________________________________________________________
 393 Int_t AliTRDsim::TrPhotons(Float_t p, Float_t mass
 394                          , Int_t &nPhoton, Float_t *ePhoton)
 395 {
 396   //
 397   // Produces TR photons.
 398   //
 399
 400   const Double_t kAlpha  = 0.0072973;
 401   const Int_t    kSumMax = 30;
 402
 403   Double_t tau   = fGapThick / fFoilThick;
 404
 405   // Calculate gamma
 406   Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
 407
 408   // Select the number of foils corresponding to momentum
 409   Int_t    foils = SelectNFoils(p);
 410
 411   fSpectrum->Reset();
 412
 413   //
 414   // The TR spectrum
 415   //
 416
 417   Double_t csi1;
 418   Double_t csi2;
 419   Double_t rho1;
 420   Double_t rho2;
 421   Double_t sSigma;
 422   Double_t sum;
 423   Double_t nEqu;
 424   Double_t thetaN;
 425   Double_t aux;
 426   Double_t energyeV;
 427   Double_t energykeV;
 428
 429   for (Int_t iBin = 1; iBin <= fSpNBins; iBin++) {
 430
 431     energykeV = fSpectrum->GetBinCenter(iBin);
 432     energyeV  = energykeV * 1.e3;
 433     sSigma    = Sigma(energykeV);
 434
 435     csi1      = fFoilOmega / energyeV;
 436     csi2      = fGapOmega  / energyeV;
 437
 438     rho1      = 2.5 * energyeV * fFoilThick * 1.0e4
 439               * (1. / (gamma*gamma) + csi1*csi1);
 440     rho2      = 2.5 * energyeV * fFoilThick * 1.0e4
 441               * (1.0 / (gamma*gamma) + csi2 *csi2);
 442
 443     // Calculate the sum
 444     sum = 0;
 445     for (Int_t n = 1; n <= kSumMax; n++) {
 446       thetaN = (TMath::Pi() * 2.0 * n - (rho1 + tau * rho2)) / (1.0 + tau);
 447       if (thetaN < 0.0) {
 448         thetaN = 0.0;
 449       }
 450       aux    = 1.0 / (rho1 + thetaN) - 1.0 / (rho2 + thetaN);
 451       sum  += thetaN * (aux*aux) * (1.0 - TMath::Cos(rho1 + thetaN));
 452     }
 453
 454     // Equivalent number of foils
 455     nEqu = (1.0 - TMath::Exp(-foils * sSigma)) / (1.0 - TMath::Exp(-sSigma));
 456
 457     // dN / domega
 458     fSpectrum->SetBinContent(iBin,4.0 * kAlpha * nEqu * sum / (energykeV * (1.0 + tau)));
 459
 460   }
 461
 462   // <nTR> (binsize corr.)
 463   Float_t ntr = fSpBinWidth*fSpectrum->Integral();
 464   // Number of TR photons from Poisson distribution with mean <ntr>
 465   nPhoton     = gRandom->Poisson(ntr);
 466   // Energy of the TR photons
 467   for (Int_t iPhoton = 0; iPhoton < nPhoton; iPhoton++) {
 468     ePhoton[iPhoton] = fSpectrum->GetRandom();
 469   }
 470
 471   return 1;
 472
 473 }
 474
 475 //_____________________________________________________________________________
 476 void AliTRDsim::SetSigma()
 477 {
 478   //
 479   // Sets the absorbtion crosssection for the energies of the TR spectrum
 480   //
 481
 482   if (fSigma) {
 483     delete [] fSigma;
 484   }
 485   fSigma = new Double_t[fSpNBins];
 486
 487   for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
 488     Double_t energykeV = iBin * fSpBinWidth + 1.0;
 489     fSigma[iBin]       = Sigma(energykeV);
 490   }
 491
 492 }
 493
 494 //_____________________________________________________________________________
 495 Double_t AliTRDsim::Sigma(Double_t energykeV)
 496 {
 497   //
 498   // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
 499   //
 500
 501   // keV -> MeV
 502   Double_t energyMeV = energykeV * 0.001;
 503   if (energyMeV >= 0.001) {
 504     return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
 505            GetMuAi(energyMeV) * fGapDens  * fGapThick  * GetTemp());
 506   }
 507   else {
 508     return 1.0e6;
 509   }
 510
 511 }
 512
 513 //_____________________________________________________________________________
 514 Double_t AliTRDsim::GetMuPo(Double_t energyMeV)
 515 {
 516   //
 517   // Returns the photon absorbtion cross section for polypropylene
 518   //
 519
 520   const Int_t kN = 36;
 521
 522   Double_t mu[kN] = { 1.894E+03, 5.999E+02, 2.593E+02
 523                     , 7.743E+01, 3.242E+01, 1.643E+01
 524                     , 9.432E+00, 3.975E+00, 2.088E+00
 525                     , 7.452E-01, 4.315E-01, 2.706E-01
 526                     , 2.275E-01, 2.084E-01, 1.970E-01
 527                     , 1.823E-01, 1.719E-01, 1.534E-01
 528                     , 1.402E-01, 1.217E-01, 1.089E-01
 529                     , 9.947E-02, 9.198E-02, 8.078E-02
 530                     , 7.262E-02, 6.495E-02, 5.910E-02
 531                     , 5.064E-02, 4.045E-02, 3.444E-02
 532                     , 3.045E-02, 2.760E-02, 2.383E-02
 533                     , 2.145E-02, 1.819E-02, 1.658E-02 };
 534
 535   Double_t en[kN] = { 1.000E-03, 1.500E-03, 2.000E-03
 536                     , 3.000E-03, 4.000E-03, 5.000E-03
 537                     , 6.000E-03, 8.000E-03, 1.000E-02
 538                     , 1.500E-02, 2.000E-02, 3.000E-02
 539                     , 4.000E-02, 5.000E-02, 6.000E-02
 540                     , 8.000E-02, 1.000E-01, 1.500E-01
 541                     , 2.000E-01, 3.000E-01, 4.000E-01
 542                     , 5.000E-01, 6.000E-01, 8.000E-01
 543                     , 1.000E+00, 1.250E+00, 1.500E+00
 544                     , 2.000E+00, 3.000E+00, 4.000E+00
 545                     , 5.000E+00, 6.000E+00, 8.000E+00
 546                     , 1.000E+01, 1.500E+01, 2.000E+01 };
 547
 548   return Interpolate(energyMeV,en,mu,kN);
 549
 550 }
 551
 552 //_____________________________________________________________________________
 553 Double_t AliTRDsim::GetMuCO(Double_t energyMeV)
 554 {
 555   //
 556   // Returns the photon absorbtion cross section for CO2
 557   //
 558
 559   const Int_t kN = 36;
 560
 561   Double_t mu[kN] = { 0.39383E+04, 0.13166E+04, 0.58750E+03
 562                     , 0.18240E+03, 0.77996E+02, 0.40024E+02
 563                     , 0.23116E+02, 0.96997E+01, 0.49726E+01
 564                     , 0.15543E+01, 0.74915E+00, 0.34442E+00
 565                     , 0.24440E+00, 0.20589E+00, 0.18632E+00
 566                     , 0.16578E+00, 0.15394E+00, 0.13558E+00
 567                     , 0.12336E+00, 0.10678E+00, 0.95510E-01
 568                     , 0.87165E-01, 0.80587E-01, 0.70769E-01
 569                     , 0.63626E-01, 0.56894E-01, 0.51782E-01
 570                     , 0.44499E-01, 0.35839E-01, 0.30825E-01
 571                     , 0.27555E-01, 0.25269E-01, 0.22311E-01
 572                     , 0.20516E-01, 0.18184E-01, 0.17152E-01 };
 573
 574   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
 575                     , 0.30000E-02, 0.40000E-02, 0.50000E-02
 576                     , 0.60000E-02, 0.80000E-02, 0.10000E-01
 577                     , 0.15000E-01, 0.20000E-01, 0.30000E-01
 578                     , 0.40000E-01, 0.50000E-01, 0.60000E-01
 579                     , 0.80000E-01, 0.10000E+00, 0.15000E+00
 580                     , 0.20000E+00, 0.30000E+00, 0.40000E+00
 581                     , 0.50000E+00, 0.60000E+00, 0.80000E+00
 582                     , 0.10000E+01, 0.12500E+01, 0.15000E+01
 583                     , 0.20000E+01, 0.30000E+01, 0.40000E+01
 584                     , 0.50000E+01, 0.60000E+01, 0.80000E+01
 585                     , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
 586
 587   return Interpolate(energyMeV,en,mu,kN);
 588
 589 }
 590
 591 //_____________________________________________________________________________
 592 Double_t AliTRDsim::GetMuXe(Double_t energyMeV)
 593 {
 594   //
 595   // Returns the photon absorbtion cross section for xenon
 596   //
 597
 598   const Int_t kN = 48;
 599
 600   Double_t mu[kN] = { 9.413E+03, 8.151E+03, 7.035E+03
 601                     , 7.338E+03, 4.085E+03, 2.088E+03
 602                     , 7.780E+02, 3.787E+02, 2.408E+02
 603                     , 6.941E+02, 6.392E+02, 6.044E+02
 604                     , 8.181E+02, 7.579E+02, 6.991E+02
 605                     , 8.064E+02, 6.376E+02, 3.032E+02
 606                     , 1.690E+02, 5.743E+01, 2.652E+01
 607                     , 8.930E+00, 6.129E+00, 3.316E+01
 608                     , 2.270E+01, 1.272E+01, 7.825E+00
 609                     , 3.633E+00, 2.011E+00, 7.202E-01
 610                     , 3.760E-01, 1.797E-01, 1.223E-01
 611                     , 9.699E-02, 8.281E-02, 6.696E-02
 612                     , 5.785E-02, 5.054E-02, 4.594E-02
 613                     , 4.078E-02, 3.681E-02, 3.577E-02
 614                     , 3.583E-02, 3.634E-02, 3.797E-02
 615                     , 3.987E-02, 4.445E-02, 4.815E-02 };
 616
 617   Double_t en[kN] = { 1.00000E-03, 1.07191E-03, 1.14900E-03
 618                     , 1.14900E-03, 1.50000E-03, 2.00000E-03
 619                     , 3.00000E-03, 4.00000E-03, 4.78220E-03
 620                     , 4.78220E-03, 5.00000E-03, 5.10370E-03
 621                     , 5.10370E-03, 5.27536E-03, 5.45280E-03
 622                     , 5.45280E-03, 6.00000E-03, 8.00000E-03
 623                     , 1.00000E-02, 1.50000E-02, 2.00000E-02
 624                     , 3.00000E-02, 3.45614E-02, 3.45614E-02
 625                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 626                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 627                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 628                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 629                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 630                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 631                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 632                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 633
 634   return Interpolate(energyMeV,en,mu,kN);
 635
 636 }
 637
 638 //_____________________________________________________________________________
 639 Double_t AliTRDsim::GetMuBu(Double_t energyMeV)
 640 {
 641   //
 642   // Returns the photon absorbtion cross section for isobutane
 643   //
 644
 645   const Int_t kN = 36;
 646
 647   Double_t mu[kN] = { 0.38846E+03, 0.12291E+03, 0.53225E+02
 648                     , 0.16091E+02, 0.69114E+01, 0.36541E+01
 649                     , 0.22282E+01, 0.11149E+01, 0.72887E+00
 650                     , 0.45053E+00, 0.38167E+00, 0.33920E+00
 651                     , 0.32155E+00, 0.30949E+00, 0.29960E+00
 652                     , 0.28317E+00, 0.26937E+00, 0.24228E+00
 653                     , 0.22190E+00, 0.19289E+00, 0.17288E+00
 654                     , 0.15789E+00, 0.14602E+00, 0.12829E+00
 655                     , 0.11533E+00, 0.10310E+00, 0.93790E-01
 656                     , 0.80117E-01, 0.63330E-01, 0.53229E-01
 657                     , 0.46390E-01, 0.41425E-01, 0.34668E-01
 658                     , 0.30267E-01, 0.23910E-01, 0.20509E-01 };
 659
 660   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
 661                     , 0.30000E-02, 0.40000E-02, 0.50000E-02
 662                     , 0.60000E-02, 0.80000E-02, 0.10000E-01
 663                     , 0.15000E-01, 0.20000E-01, 0.30000E-01
 664                     , 0.40000E-01, 0.50000E-01, 0.60000E-01
 665                     , 0.80000E-01, 0.10000E+00, 0.15000E+00
 666                     , 0.20000E+00, 0.30000E+00, 0.40000E+00
 667                     , 0.50000E+00, 0.60000E+00, 0.80000E+00
 668                     , 0.10000E+01, 0.12500E+01, 0.15000E+01
 669                     , 0.20000E+01, 0.30000E+01, 0.40000E+01
 670                     , 0.50000E+01, 0.60000E+01, 0.80000E+01
 671                     , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
 672
 673   return Interpolate(energyMeV,en,mu,kN);
 674
 675 }
 676
 677 //_____________________________________________________________________________
 678 Double_t AliTRDsim::GetMuMy(Double_t energyMeV)
 679 {
 680   //
 681   // Returns the photon absorbtion cross section for mylar
 682   //
 683
 684   const Int_t kN = 36;
 685
 686   Double_t mu[kN] = { 2.911E+03, 9.536E+02, 4.206E+02
 687                     , 1.288E+02, 5.466E+01, 2.792E+01
 688                     , 1.608E+01, 6.750E+00, 3.481E+00
 689                     , 1.132E+00, 5.798E-01, 3.009E-01
 690                     , 2.304E-01, 2.020E-01, 1.868E-01
 691                     , 1.695E-01, 1.586E-01, 1.406E-01
 692                     , 1.282E-01, 1.111E-01, 9.947E-02
 693                     , 9.079E-02, 8.395E-02, 7.372E-02
 694                     , 6.628E-02, 5.927E-02, 5.395E-02
 695                     , 4.630E-02, 3.715E-02, 3.181E-02
 696                     , 2.829E-02, 2.582E-02, 2.257E-02
 697                     , 2.057E-02, 1.789E-02, 1.664E-02 };
 698
 699   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 700                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 701                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 702                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 703                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 704                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 705                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 706                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 707                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 708                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 709                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 710                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 711
 712   return Interpolate(energyMeV,en,mu,kN);
 713
 714 }
 715
 716 //_____________________________________________________________________________
 717 Double_t AliTRDsim::GetMuN2(Double_t energyMeV)
 718 {
 719   //
 720   // Returns the photon absorbtion cross section for nitrogen
 721   //
 722
 723   const Int_t kN = 36;
 724
 725   Double_t mu[kN] = { 3.311E+03, 1.083E+03, 4.769E+02
 726                     , 1.456E+02, 6.166E+01, 3.144E+01
 727                     , 1.809E+01, 7.562E+00, 3.879E+00
 728                     , 1.236E+00, 6.178E-01, 3.066E-01
 729                     , 2.288E-01, 1.980E-01, 1.817E-01
 730                     , 1.639E-01, 1.529E-01, 1.353E-01
 731                     , 1.233E-01, 1.068E-01, 9.557E-02
 732                     , 8.719E-02, 8.063E-02, 7.081E-02
 733                     , 6.364E-02, 5.693E-02, 5.180E-02
 734                     , 4.450E-02, 3.579E-02, 3.073E-02
 735                     , 2.742E-02, 2.511E-02, 2.209E-02
 736                     , 2.024E-02, 1.782E-02, 1.673E-02 };
 737
 738   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 739                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 740                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 741                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 742                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 743                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 744                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 745                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 746                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 747                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 748                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 749                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 750
 751   return Interpolate(energyMeV,en,mu,kN);
 752
 753 }
 754
 755 //_____________________________________________________________________________
 756 Double_t AliTRDsim::GetMuO2(Double_t energyMeV)
 757 {
 758   //
 759   // Returns the photon absorbtion cross section for oxygen
 760   //
 761
 762   const Int_t kN = 36;
 763
 764   Double_t mu[kN] = { 4.590E+03, 1.549E+03, 6.949E+02
 765                     , 2.171E+02, 9.315E+01, 4.790E+01
 766                     , 2.770E+01, 1.163E+01, 5.952E+00
 767                     , 1.836E+00, 8.651E-01, 3.779E-01
 768                     , 2.585E-01, 2.132E-01, 1.907E-01
 769                     , 1.678E-01, 1.551E-01, 1.361E-01
 770                     , 1.237E-01, 1.070E-01, 9.566E-02
 771                     , 8.729E-02, 8.070E-02, 7.087E-02
 772                     , 6.372E-02, 5.697E-02, 5.185E-02
 773                     , 4.459E-02, 3.597E-02, 3.100E-02
 774                     , 2.777E-02, 2.552E-02, 2.263E-02
 775                     , 2.089E-02, 1.866E-02, 1.770E-02 };
 776
 777   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 778                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 779                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 780                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 781                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 782                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 783                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 784                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 785                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 786                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 787                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 788                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 789
 790   return Interpolate(energyMeV,en,mu,kN);
 791
 792 }
 793
 794 //_____________________________________________________________________________
 795 Double_t AliTRDsim::GetMuHe(Double_t energyMeV)
 796 {
 797   //
 798   // Returns the photon absorbtion cross section for helium
 799   //
 800
 801   const Int_t kN = 36;
 802
 803   Double_t mu[kN] = { 6.084E+01, 1.676E+01, 6.863E+00
 804                     , 2.007E+00, 9.329E-01, 5.766E-01
 805                     , 4.195E-01, 2.933E-01, 2.476E-01
 806                     , 2.092E-01, 1.960E-01, 1.838E-01
 807                     , 1.763E-01, 1.703E-01, 1.651E-01
 808                     , 1.562E-01, 1.486E-01, 1.336E-01
 809                     , 1.224E-01, 1.064E-01, 9.535E-02
 810                     , 8.707E-02, 8.054E-02, 7.076E-02
 811                     , 6.362E-02, 5.688E-02, 5.173E-02
 812                     , 4.422E-02, 3.503E-02, 2.949E-02
 813                     , 2.577E-02, 2.307E-02, 1.940E-02
 814                     , 1.703E-02, 1.363E-02, 1.183E-02 };
 815
 816   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 817                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 818                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 819                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 820                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 821                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 822                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 823                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 824                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 825                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 826                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 827                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 828
 829   return Interpolate(energyMeV,en,mu,kN);
 830
 831 }
 832
 833 //_____________________________________________________________________________
 834 Double_t AliTRDsim::GetMuAi(Double_t energyMeV)
 835 {
 836   //
 837   // Returns the photon absorbtion cross section for air
 838   // Implemented by Oliver Busch
 839   //
 840
 841   const Int_t kN = 38;
 842
 843   Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
 844                       0.16143E+03, 0.14250E+03, 0.15722E+03,
 845                       0.77538E+02, 0.40099E+02, 0.23313E+02,
 846                       0.98816E+01, 0.51000E+01, 0.16079E+01,
 847                       0.77536E+00, 0.35282E+00, 0.24790E+00,
 848                       0.20750E+00, 0.18703E+00, 0.16589E+00,
 849                       0.15375E+00, 0.13530E+00, 0.12311E+00,
 850                       0.10654E+00, 0.95297E-01, 0.86939E-01,
 851                       0.80390E-01, 0.70596E-01, 0.63452E-01,
 852                       0.56754E-01, 0.51644E-01, 0.44382E-01,
 853                       0.35733E-01, 0.30721E-01, 0.27450E-01,
 854                       0.25171E-01, 0.22205E-01, 0.20399E-01,
 855                       0.18053E-01, 0.18057E-01 };
 856
 857
 858
 859   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
 860                       0.30000E-02, 0.32029E-02, 0.32029E-02,
 861                       0.40000E-02, 0.50000E-02, 0.60000E-02,
 862                       0.80000E-02, 0.10000E-01, 0.15000E-01,
 863                       0.20000E-01, 0.30000E-01, 0.40000E-01,
 864                       0.50000E-01, 0.60000E-01, 0.80000E-01,
 865                       0.10000E+00, 0.15000E+00, 0.20000E+00,
 866                       0.30000E+00, 0.40000E+00, 0.50000E+00,
 867                       0.60000E+00, 0.80000E+00, 0.10000E+01,
 868                       0.12500E+01, 0.15000E+01, 0.20000E+01,
 869                       0.30000E+01, 0.40000E+01, 0.50000E+01,
 870                       0.60000E+01, 0.80000E+01, 0.10000E+02,
 871                       0.15000E+02, 0.20000E+02 };
 872
 873   return Interpolate(energyMeV,en,mu,kN);
 874
 875 }
 876
 877 //_____________________________________________________________________________
 878 Double_t AliTRDsim::Interpolate(Double_t energyMeV
 879                               , Double_t *en, Double_t *mu, Int_t n)
 880 {
 881   //
 882   // Interpolates the photon absorbtion cross section
 883   // for a given energy <energyMeV>.
 884   //
 885
 886   Double_t de    = 0;
 887   Int_t    index = 0;
 888   Int_t    istat = Locate(en,n,energyMeV,index,de);
 889   if (istat == 0) {
 890     return (mu[index] - de * (mu[index]   - mu[index+1])
 891                            / (en[index+1] - en[index]  ));
 892   }
 893   else {
 894     return 0.0;
 895   }
 896
 897 }
 898
 899 //_____________________________________________________________________________
 900 Int_t AliTRDsim::Locate(Double_t *xv, Int_t n, Double_t xval
 901                       , Int_t &kl, Double_t &dx)
 902 {
 903   //
 904   // Locates a point (xval) in a 1-dim grid (xv(n))
 905   //
 906
 907   if (xval >= xv[n-1]) {
 908     return  1;
 909   }
 910   if (xval <  xv[0]) {
 911     return -1;
 912   }
 913
 914   Int_t km;
 915   Int_t kh = n - 1;
 916
 917   kl = 0;
 918   while (kh - kl > 1) {
 919     if (xval < xv[km = (kl+kh)/2]) {
 920       kh = km;
 921     }
 922     else {
 923       kl = km;
 924     }
 925   }
 926   if ((xval <  xv[kl])   ||
 927       (xval >  xv[kl+1]) ||
 928       (kl   >= n-1)) {
 929     AliError(Form("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
 930                  ,kl,xv[kl],xval,kl+1,xv[kl+1]));
 931     exit(1);
 932   }
 933
 934   dx = xval - xv[kl];
 935
 936   return 0;
 937
 938 }
 939
 940 //_____________________________________________________________________________
 941 Int_t AliTRDsim::SelectNFoils(Float_t p)
 942 {
 943   //
 944   // Selects the number of foils corresponding to the momentum
 945   //
 946
 947   Int_t foils = fNFoils[fNFoilsDim-1];
 948
 949   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
 950     if (p < fNFoilsUp[iFoil]) {
 951       foils = fNFoils[iFoil];
 952       break;
 953     }
 954   }
 955
 956   return foils;
 957
 958 }