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