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