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