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