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