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   Double_t t = gMC->TrackTime();
 498
 499   // Counter for TR analysed in custom code (e < 15keV)
 500   nPhoton = 0;
 501
 502   for (Int_t iPhoton = 0; iPhoton < nPhCand; iPhoton++) {
 503
 504     // Energy of the TR photon
 505     Double_t e = fSpectrum->GetRandom();
 506
 507     // Put TR photon on particle stack
 508     if (e > 15.0) {
 509
 510       e *= 1.0e-6; // Convert it to GeV
 511
 512       Int_t phtrack;
 513       stack->PushTrack(1                 // Must be 1
 514                       ,track             // Identifier of the parent track, -1 for a primary
 515                       ,22                // Particle code.
 516                       ,px*e              // 4 momentum (The photon is generated on the same
 517                       ,py*e              // direction as the parent. For irregular radiator one
 518                       ,pz*e              // can calculate also the angle but this is a secondary
 519                       ,e                 // order effect)
 520                       ,x,y,z,t           // 4 vertex
 521                       ,0.0,0.0,0.0       // Polarisation
 522                       ,kPFeedBackPhoton  // Production mechanism (there is no TR in G3 so one
 523                                          // has to make some convention)
 524                       ,phtrack           // On output the number of the track stored
 525                       ,1.0
 526                       ,1);
 527
 528     }
 529     // Custom treatment of TR photons
 530     else {
 531
 532       ePhoton[nPhoton++] = e;
 533
 534     }
 535
 536   }
 537
 538   return 1;
 539
 540 }
 541
 542 //_____________________________________________________________________________
 543 void AliTRDsimTR::SetSigma()
 544 {
 545   //
 546   // Sets the absorbtion crosssection for the energies of the TR spectrum
 547   //
 548
 549   if (fSigma) {
 550     delete [] fSigma;
 551   }
 552   fSigma = new Double_t[fSpNBins];
 553
 554   for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
 555     Double_t energykeV = iBin * fSpBinWidth + 1.0;
 556     fSigma[iBin]       = Sigma(energykeV);
 557   }
 558
 559 }
 560
 561 //_____________________________________________________________________________
 562 Double_t AliTRDsimTR::Sigma(Double_t energykeV)
 563 {
 564   //
 565   // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
 566   //
 567
 568   // keV -> MeV
 569   Double_t energyMeV = energykeV * 0.001;
 570   if (energyMeV >= 0.001) {
 571     return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
 572            GetMuAi(energyMeV) * fGapDens  * fGapThick  * GetTemp());
 573   }
 574   else {
 575     return 1.0e6;
 576   }
 577
 578 }
 579
 580 //_____________________________________________________________________________
 581 Double_t AliTRDsimTR::GetMuPo(Double_t energyMeV)
 582 {
 583   //
 584   // Returns the photon absorbtion cross section for polypropylene
 585   //
 586
 587   const Int_t kN = 36;
 588
 589   Double_t mu[kN] = { 1.894E+03, 5.999E+02, 2.593E+02
 590                     , 7.743E+01, 3.242E+01, 1.643E+01
 591                     , 9.432E+00, 3.975E+00, 2.088E+00
 592                     , 7.452E-01, 4.315E-01, 2.706E-01
 593                     , 2.275E-01, 2.084E-01, 1.970E-01
 594                     , 1.823E-01, 1.719E-01, 1.534E-01
 595                     , 1.402E-01, 1.217E-01, 1.089E-01
 596                     , 9.947E-02, 9.198E-02, 8.078E-02
 597                     , 7.262E-02, 6.495E-02, 5.910E-02
 598                     , 5.064E-02, 4.045E-02, 3.444E-02
 599                     , 3.045E-02, 2.760E-02, 2.383E-02
 600                     , 2.145E-02, 1.819E-02, 1.658E-02 };
 601
 602   Double_t en[kN] = { 1.000E-03, 1.500E-03, 2.000E-03
 603                     , 3.000E-03, 4.000E-03, 5.000E-03
 604                     , 6.000E-03, 8.000E-03, 1.000E-02
 605                     , 1.500E-02, 2.000E-02, 3.000E-02
 606                     , 4.000E-02, 5.000E-02, 6.000E-02
 607                     , 8.000E-02, 1.000E-01, 1.500E-01
 608                     , 2.000E-01, 3.000E-01, 4.000E-01
 609                     , 5.000E-01, 6.000E-01, 8.000E-01
 610                     , 1.000E+00, 1.250E+00, 1.500E+00
 611                     , 2.000E+00, 3.000E+00, 4.000E+00
 612                     , 5.000E+00, 6.000E+00, 8.000E+00
 613                     , 1.000E+01, 1.500E+01, 2.000E+01 };
 614
 615   return Interpolate(energyMeV,en,mu,kN);
 616
 617 }
 618
 619 //_____________________________________________________________________________
 620 Double_t AliTRDsimTR::GetMuCO(Double_t energyMeV)
 621 {
 622   //
 623   // Returns the photon absorbtion cross section for CO2
 624   //
 625
 626   const Int_t kN = 36;
 627
 628   Double_t mu[kN] = { 0.39383E+04, 0.13166E+04, 0.58750E+03
 629                     , 0.18240E+03, 0.77996E+02, 0.40024E+02
 630                     , 0.23116E+02, 0.96997E+01, 0.49726E+01
 631                     , 0.15543E+01, 0.74915E+00, 0.34442E+00
 632                     , 0.24440E+00, 0.20589E+00, 0.18632E+00
 633                     , 0.16578E+00, 0.15394E+00, 0.13558E+00
 634                     , 0.12336E+00, 0.10678E+00, 0.95510E-01
 635                     , 0.87165E-01, 0.80587E-01, 0.70769E-01
 636                     , 0.63626E-01, 0.56894E-01, 0.51782E-01
 637                     , 0.44499E-01, 0.35839E-01, 0.30825E-01
 638                     , 0.27555E-01, 0.25269E-01, 0.22311E-01
 639                     , 0.20516E-01, 0.18184E-01, 0.17152E-01 };
 640
 641   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
 642                     , 0.30000E-02, 0.40000E-02, 0.50000E-02
 643                     , 0.60000E-02, 0.80000E-02, 0.10000E-01
 644                     , 0.15000E-01, 0.20000E-01, 0.30000E-01
 645                     , 0.40000E-01, 0.50000E-01, 0.60000E-01
 646                     , 0.80000E-01, 0.10000E+00, 0.15000E+00
 647                     , 0.20000E+00, 0.30000E+00, 0.40000E+00
 648                     , 0.50000E+00, 0.60000E+00, 0.80000E+00
 649                     , 0.10000E+01, 0.12500E+01, 0.15000E+01
 650                     , 0.20000E+01, 0.30000E+01, 0.40000E+01
 651                     , 0.50000E+01, 0.60000E+01, 0.80000E+01
 652                     , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
 653
 654   return Interpolate(energyMeV,en,mu,kN);
 655
 656 }
 657
 658 //_____________________________________________________________________________
 659 Double_t AliTRDsimTR::GetMuXe(Double_t energyMeV)
 660 {
 661   //
 662   // Returns the photon absorbtion cross section for xenon
 663   //
 664
 665   const Int_t kN = 48;
 666
 667   Double_t mu[kN] = { 9.413E+03, 8.151E+03, 7.035E+03
 668                     , 7.338E+03, 4.085E+03, 2.088E+03
 669                     , 7.780E+02, 3.787E+02, 2.408E+02
 670                     , 6.941E+02, 6.392E+02, 6.044E+02
 671                     , 8.181E+02, 7.579E+02, 6.991E+02
 672                     , 8.064E+02, 6.376E+02, 3.032E+02
 673                     , 1.690E+02, 5.743E+01, 2.652E+01
 674                     , 8.930E+00, 6.129E+00, 3.316E+01
 675                     , 2.270E+01, 1.272E+01, 7.825E+00
 676                     , 3.633E+00, 2.011E+00, 7.202E-01
 677                     , 3.760E-01, 1.797E-01, 1.223E-01
 678                     , 9.699E-02, 8.281E-02, 6.696E-02
 679                     , 5.785E-02, 5.054E-02, 4.594E-02
 680                     , 4.078E-02, 3.681E-02, 3.577E-02
 681                     , 3.583E-02, 3.634E-02, 3.797E-02
 682                     , 3.987E-02, 4.445E-02, 4.815E-02 };
 683
 684   Double_t en[kN] = { 1.00000E-03, 1.07191E-03, 1.14900E-03
 685                     , 1.14900E-03, 1.50000E-03, 2.00000E-03
 686                     , 3.00000E-03, 4.00000E-03, 4.78220E-03
 687                     , 4.78220E-03, 5.00000E-03, 5.10370E-03
 688                     , 5.10370E-03, 5.27536E-03, 5.45280E-03
 689                     , 5.45280E-03, 6.00000E-03, 8.00000E-03
 690                     , 1.00000E-02, 1.50000E-02, 2.00000E-02
 691                     , 3.00000E-02, 3.45614E-02, 3.45614E-02
 692                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 693                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 694                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 695                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 696                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 697                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 698                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 699                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 700
 701   return Interpolate(energyMeV,en,mu,kN);
 702
 703 }
 704
 705 //_____________________________________________________________________________
 706 Double_t AliTRDsimTR::GetMuBu(Double_t energyMeV)
 707 {
 708   //
 709   // Returns the photon absorbtion cross section for isobutane
 710   //
 711
 712   const Int_t kN = 36;
 713
 714   Double_t mu[kN] = { 0.38846E+03, 0.12291E+03, 0.53225E+02
 715                     , 0.16091E+02, 0.69114E+01, 0.36541E+01
 716                     , 0.22282E+01, 0.11149E+01, 0.72887E+00
 717                     , 0.45053E+00, 0.38167E+00, 0.33920E+00
 718                     , 0.32155E+00, 0.30949E+00, 0.29960E+00
 719                     , 0.28317E+00, 0.26937E+00, 0.24228E+00
 720                     , 0.22190E+00, 0.19289E+00, 0.17288E+00
 721                     , 0.15789E+00, 0.14602E+00, 0.12829E+00
 722                     , 0.11533E+00, 0.10310E+00, 0.93790E-01
 723                     , 0.80117E-01, 0.63330E-01, 0.53229E-01
 724                     , 0.46390E-01, 0.41425E-01, 0.34668E-01
 725                     , 0.30267E-01, 0.23910E-01, 0.20509E-01 };
 726
 727   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
 728                     , 0.30000E-02, 0.40000E-02, 0.50000E-02
 729                     , 0.60000E-02, 0.80000E-02, 0.10000E-01
 730                     , 0.15000E-01, 0.20000E-01, 0.30000E-01
 731                     , 0.40000E-01, 0.50000E-01, 0.60000E-01
 732                     , 0.80000E-01, 0.10000E+00, 0.15000E+00
 733                     , 0.20000E+00, 0.30000E+00, 0.40000E+00
 734                     , 0.50000E+00, 0.60000E+00, 0.80000E+00
 735                     , 0.10000E+01, 0.12500E+01, 0.15000E+01
 736                     , 0.20000E+01, 0.30000E+01, 0.40000E+01
 737                     , 0.50000E+01, 0.60000E+01, 0.80000E+01
 738                     , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
 739
 740   return Interpolate(energyMeV,en,mu,kN);
 741
 742 }
 743
 744 //_____________________________________________________________________________
 745 Double_t AliTRDsimTR::GetMuMy(Double_t energyMeV)
 746 {
 747   //
 748   // Returns the photon absorbtion cross section for mylar
 749   //
 750
 751   const Int_t kN = 36;
 752
 753   Double_t mu[kN] = { 2.911E+03, 9.536E+02, 4.206E+02
 754                     , 1.288E+02, 5.466E+01, 2.792E+01
 755                     , 1.608E+01, 6.750E+00, 3.481E+00
 756                     , 1.132E+00, 5.798E-01, 3.009E-01
 757                     , 2.304E-01, 2.020E-01, 1.868E-01
 758                     , 1.695E-01, 1.586E-01, 1.406E-01
 759                     , 1.282E-01, 1.111E-01, 9.947E-02
 760                     , 9.079E-02, 8.395E-02, 7.372E-02
 761                     , 6.628E-02, 5.927E-02, 5.395E-02
 762                     , 4.630E-02, 3.715E-02, 3.181E-02
 763                     , 2.829E-02, 2.582E-02, 2.257E-02
 764                     , 2.057E-02, 1.789E-02, 1.664E-02 };
 765
 766   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 767                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 768                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 769                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 770                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 771                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 772                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 773                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 774                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 775                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 776                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 777                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 778
 779   return Interpolate(energyMeV,en,mu,kN);
 780
 781 }
 782
 783 //_____________________________________________________________________________
 784 Double_t AliTRDsimTR::GetMuN2(Double_t energyMeV)
 785 {
 786   //
 787   // Returns the photon absorbtion cross section for nitrogen
 788   //
 789
 790   const Int_t kN = 36;
 791
 792   Double_t mu[kN] = { 3.311E+03, 1.083E+03, 4.769E+02
 793                     , 1.456E+02, 6.166E+01, 3.144E+01
 794                     , 1.809E+01, 7.562E+00, 3.879E+00
 795                     , 1.236E+00, 6.178E-01, 3.066E-01
 796                     , 2.288E-01, 1.980E-01, 1.817E-01
 797                     , 1.639E-01, 1.529E-01, 1.353E-01
 798                     , 1.233E-01, 1.068E-01, 9.557E-02
 799                     , 8.719E-02, 8.063E-02, 7.081E-02
 800                     , 6.364E-02, 5.693E-02, 5.180E-02
 801                     , 4.450E-02, 3.579E-02, 3.073E-02
 802                     , 2.742E-02, 2.511E-02, 2.209E-02
 803                     , 2.024E-02, 1.782E-02, 1.673E-02 };
 804
 805   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 806                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 807                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 808                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 809                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 810                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 811                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 812                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 813                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 814                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 815                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 816                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 817
 818   return Interpolate(energyMeV,en,mu,kN);
 819
 820 }
 821
 822 //_____________________________________________________________________________
 823 Double_t AliTRDsimTR::GetMuO2(Double_t energyMeV)
 824 {
 825   //
 826   // Returns the photon absorbtion cross section for oxygen
 827   //
 828
 829   const Int_t kN = 36;
 830
 831   Double_t mu[kN] = { 4.590E+03, 1.549E+03, 6.949E+02
 832                     , 2.171E+02, 9.315E+01, 4.790E+01
 833                     , 2.770E+01, 1.163E+01, 5.952E+00
 834                     , 1.836E+00, 8.651E-01, 3.779E-01
 835                     , 2.585E-01, 2.132E-01, 1.907E-01
 836                     , 1.678E-01, 1.551E-01, 1.361E-01
 837                     , 1.237E-01, 1.070E-01, 9.566E-02
 838                     , 8.729E-02, 8.070E-02, 7.087E-02
 839                     , 6.372E-02, 5.697E-02, 5.185E-02
 840                     , 4.459E-02, 3.597E-02, 3.100E-02
 841                     , 2.777E-02, 2.552E-02, 2.263E-02
 842                     , 2.089E-02, 1.866E-02, 1.770E-02 };
 843
 844   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 845                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 846                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 847                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 848                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 849                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 850                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 851                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 852                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 853                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 854                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 855                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 856
 857   return Interpolate(energyMeV,en,mu,kN);
 858
 859 }
 860
 861 //_____________________________________________________________________________
 862 Double_t AliTRDsimTR::GetMuHe(Double_t energyMeV)
 863 {
 864   //
 865   // Returns the photon absorbtion cross section for helium
 866   //
 867
 868   const Int_t kN = 36;
 869
 870   Double_t mu[kN] = { 6.084E+01, 1.676E+01, 6.863E+00
 871                     , 2.007E+00, 9.329E-01, 5.766E-01
 872                     , 4.195E-01, 2.933E-01, 2.476E-01
 873                     , 2.092E-01, 1.960E-01, 1.838E-01
 874                     , 1.763E-01, 1.703E-01, 1.651E-01
 875                     , 1.562E-01, 1.486E-01, 1.336E-01
 876                     , 1.224E-01, 1.064E-01, 9.535E-02
 877                     , 8.707E-02, 8.054E-02, 7.076E-02
 878                     , 6.362E-02, 5.688E-02, 5.173E-02
 879                     , 4.422E-02, 3.503E-02, 2.949E-02
 880                     , 2.577E-02, 2.307E-02, 1.940E-02
 881                     , 1.703E-02, 1.363E-02, 1.183E-02 };
 882
 883   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 884                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 885                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 886                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 887                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 888                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 889                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 890                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 891                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 892                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 893                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 894                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 895
 896   return Interpolate(energyMeV,en,mu,kN);
 897
 898 }
 899
 900 //_____________________________________________________________________________
 901 Double_t AliTRDsimTR::GetMuAi(Double_t energyMeV)
 902 {
 903   //
 904   // Returns the photon absorbtion cross section for air
 905   // Implemented by Oliver Busch
 906   //
 907
 908   const Int_t kN = 38;
 909
 910   Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
 911                       0.16143E+03, 0.14250E+03, 0.15722E+03,
 912                       0.77538E+02, 0.40099E+02, 0.23313E+02,
 913                       0.98816E+01, 0.51000E+01, 0.16079E+01,
 914                       0.77536E+00, 0.35282E+00, 0.24790E+00,
 915                       0.20750E+00, 0.18703E+00, 0.16589E+00,
 916                       0.15375E+00, 0.13530E+00, 0.12311E+00,
 917                       0.10654E+00, 0.95297E-01, 0.86939E-01,
 918                       0.80390E-01, 0.70596E-01, 0.63452E-01,
 919                       0.56754E-01, 0.51644E-01, 0.44382E-01,
 920                       0.35733E-01, 0.30721E-01, 0.27450E-01,
 921                       0.25171E-01, 0.22205E-01, 0.20399E-01,
 922                       0.18053E-01, 0.18057E-01 };
 923
 924
 925
 926   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
 927                       0.30000E-02, 0.32029E-02, 0.32029E-02,
 928                       0.40000E-02, 0.50000E-02, 0.60000E-02,
 929                       0.80000E-02, 0.10000E-01, 0.15000E-01,
 930                       0.20000E-01, 0.30000E-01, 0.40000E-01,
 931                       0.50000E-01, 0.60000E-01, 0.80000E-01,
 932                       0.10000E+00, 0.15000E+00, 0.20000E+00,
 933                       0.30000E+00, 0.40000E+00, 0.50000E+00,
 934                       0.60000E+00, 0.80000E+00, 0.10000E+01,
 935                       0.12500E+01, 0.15000E+01, 0.20000E+01,
 936                       0.30000E+01, 0.40000E+01, 0.50000E+01,
 937                       0.60000E+01, 0.80000E+01, 0.10000E+02,
 938                       0.15000E+02, 0.20000E+02 };
 939
 940   return Interpolate(energyMeV,en,mu,kN);
 941
 942 }
 943
 944 //_____________________________________________________________________________
 945 Double_t AliTRDsimTR::Interpolate(Double_t energyMeV
 946                               , Double_t *en, Double_t *mu, Int_t n)
 947 {
 948   //
 949   // Interpolates the photon absorbtion cross section
 950   // for a given energy <energyMeV>.
 951   //
 952
 953   Double_t de    = 0;
 954   Int_t    index = 0;
 955   Int_t    istat = Locate(en,n,energyMeV,index,de);
 956   if (istat == 0) {
 957     return (mu[index] - de * (mu[index]   - mu[index+1])
 958                            / (en[index+1] - en[index]  ));
 959   }
 960   else {
 961     return 0.0;
 962   }
 963
 964 }
 965
 966 //_____________________________________________________________________________
 967 Int_t AliTRDsimTR::Locate(Double_t *xv, Int_t n, Double_t xval
 968                       , Int_t &kl, Double_t &dx)
 969 {
 970   //
 971   // Locates a point (xval) in a 1-dim grid (xv(n))
 972   //
 973
 974   if (xval >= xv[n-1]) {
 975     return  1;
 976   }
 977   if (xval <  xv[0]) {
 978     return -1;
 979   }
 980
 981   Int_t km;
 982   Int_t kh = n - 1;
 983
 984   kl = 0;
 985   while (kh - kl > 1) {
 986     if (xval < xv[km = (kl+kh)/2]) {
 987       kh = km;
 988     }
 989     else {
 990       kl = km;
 991     }
 992   }
 993   if ((xval <  xv[kl])   ||
 994       (xval >  xv[kl+1]) ||
 995       (kl   >= n-1)) {
 996     AliFatal(Form("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
 997                  ,kl,xv[kl],xval,kl+1,xv[kl+1]));
 998     exit(1);
 999   }
1000
1001   dx = xval - xv[kl];
1002
1003   return 0;
1004
1005 }
1006
1007 //_____________________________________________________________________________
1008 Int_t AliTRDsimTR::SelectNFoils(Float_t p) const
1009 {
1010   //
1011   // Selects the number of foils corresponding to the momentum
1012   //
1013
1014   Int_t foils = fNFoils[fNFoilsDim-1];
1015
1016   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
1017     if (p < fNFoilsUp[iFoil]) {
1018       foils = fNFoils[iFoil];
1019       break;
1020     }
1021   }
1022
1023   return foils;
1024
1025 }