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