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