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::GetMuAr(Double_t energyMeV)
 707 {
 708   //
 709   // Returns the photon absorbtion cross section for argon
 710   //
 711
 712   const Int_t kN = 38;
 713
 714   Double_t mu[kN] = { 3.184E+03, 1.105E+03, 5.120E+02
 715                     , 1.703E+02, 1.424E+02, 1.275E+03
 716                     , 7.572E+02, 4.225E+02, 2.593E+02
 717                     , 1.180E+02, 6.316E+01, 1.983E+01
 718                     , 8.629E+00, 2.697E+00, 1.228E+00
 719                     , 7.012E-01, 4.664E-01, 2.760E-01
 720                     , 2.043E-01, 1.427E-01, 1.205E-01
 721                     , 9.953E-02, 8.776E-02, 7.958E-02
 722                     , 7.335E-02, 6.419E-02, 5.762E-02
 723                     , 5.150E-02, 4.695E-02, 4.074E-02
 724                     , 3.384E-02, 3.019E-02, 2.802E-02
 725                     , 2.667E-02, 2.517E-02, 2.451E-02
 726                     , 2.418E-02, 2.453E-02 };
 727
 728   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 729                     , 3.00000E-03, 3.20290E-03, 3.20290E-03
 730                     , 4.00000E-03, 5.00000E-03, 6.00000E-03
 731                     , 8.00000E-03, 1.00000E-02, 1.50000E-02
 732                     , 2.00000E-02, 3.00000E-02, 4.00000E-02
 733                     , 5.00000E-02, 6.00000E-02, 8.00000E-02
 734                     , 1.00000E-01, 1.50000E-01, 2.00000E-01
 735                     , 3.00000E-01, 4.00000E-01, 5.00000E-01
 736                     , 6.00000E-01, 8.00000E-01, 1.00000E+00
 737                     , 1.25000E+00, 1.50000E+00, 2.00000E+00
 738                     , 3.00000E+00, 4.00000E+00, 5.00000E+00
 739                     , 6.00000E+00, 8.00000E+00, 1.00000E+01
 740                     , 1.50000E+01, 2.00000E+01 };
 741
 742   return Interpolate(energyMeV,en,mu,kN);
 743
 744 }
 745
 746 //_____________________________________________________________________________
 747 Double_t AliTRDsimTR::GetMuMy(Double_t energyMeV)
 748 {
 749   //
 750   // Returns the photon absorbtion cross section for mylar
 751   //
 752
 753   const Int_t kN = 36;
 754
 755   Double_t mu[kN] = { 2.911E+03, 9.536E+02, 4.206E+02
 756                     , 1.288E+02, 5.466E+01, 2.792E+01
 757                     , 1.608E+01, 6.750E+00, 3.481E+00
 758                     , 1.132E+00, 5.798E-01, 3.009E-01
 759                     , 2.304E-01, 2.020E-01, 1.868E-01
 760                     , 1.695E-01, 1.586E-01, 1.406E-01
 761                     , 1.282E-01, 1.111E-01, 9.947E-02
 762                     , 9.079E-02, 8.395E-02, 7.372E-02
 763                     , 6.628E-02, 5.927E-02, 5.395E-02
 764                     , 4.630E-02, 3.715E-02, 3.181E-02
 765                     , 2.829E-02, 2.582E-02, 2.257E-02
 766                     , 2.057E-02, 1.789E-02, 1.664E-02 };
 767
 768   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 769                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 770                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 771                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 772                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 773                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 774                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 775                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 776                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 777                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 778                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 779                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 780
 781   return Interpolate(energyMeV,en,mu,kN);
 782
 783 }
 784
 785 //_____________________________________________________________________________
 786 Double_t AliTRDsimTR::GetMuN2(Double_t energyMeV)
 787 {
 788   //
 789   // Returns the photon absorbtion cross section for nitrogen
 790   //
 791
 792   const Int_t kN = 36;
 793
 794   Double_t mu[kN] = { 3.311E+03, 1.083E+03, 4.769E+02
 795                     , 1.456E+02, 6.166E+01, 3.144E+01
 796                     , 1.809E+01, 7.562E+00, 3.879E+00
 797                     , 1.236E+00, 6.178E-01, 3.066E-01
 798                     , 2.288E-01, 1.980E-01, 1.817E-01
 799                     , 1.639E-01, 1.529E-01, 1.353E-01
 800                     , 1.233E-01, 1.068E-01, 9.557E-02
 801                     , 8.719E-02, 8.063E-02, 7.081E-02
 802                     , 6.364E-02, 5.693E-02, 5.180E-02
 803                     , 4.450E-02, 3.579E-02, 3.073E-02
 804                     , 2.742E-02, 2.511E-02, 2.209E-02
 805                     , 2.024E-02, 1.782E-02, 1.673E-02 };
 806
 807   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 808                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 809                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 810                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 811                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 812                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 813                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 814                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 815                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 816                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 817                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 818                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 819
 820   return Interpolate(energyMeV,en,mu,kN);
 821
 822 }
 823
 824 //_____________________________________________________________________________
 825 Double_t AliTRDsimTR::GetMuO2(Double_t energyMeV)
 826 {
 827   //
 828   // Returns the photon absorbtion cross section for oxygen
 829   //
 830
 831   const Int_t kN = 36;
 832
 833   Double_t mu[kN] = { 4.590E+03, 1.549E+03, 6.949E+02
 834                     , 2.171E+02, 9.315E+01, 4.790E+01
 835                     , 2.770E+01, 1.163E+01, 5.952E+00
 836                     , 1.836E+00, 8.651E-01, 3.779E-01
 837                     , 2.585E-01, 2.132E-01, 1.907E-01
 838                     , 1.678E-01, 1.551E-01, 1.361E-01
 839                     , 1.237E-01, 1.070E-01, 9.566E-02
 840                     , 8.729E-02, 8.070E-02, 7.087E-02
 841                     , 6.372E-02, 5.697E-02, 5.185E-02
 842                     , 4.459E-02, 3.597E-02, 3.100E-02
 843                     , 2.777E-02, 2.552E-02, 2.263E-02
 844                     , 2.089E-02, 1.866E-02, 1.770E-02 };
 845
 846   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 847                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 848                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 849                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 850                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 851                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 852                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 853                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 854                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 855                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 856                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 857                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 858
 859   return Interpolate(energyMeV,en,mu,kN);
 860
 861 }
 862
 863 //_____________________________________________________________________________
 864 Double_t AliTRDsimTR::GetMuHe(Double_t energyMeV)
 865 {
 866   //
 867   // Returns the photon absorbtion cross section for helium
 868   //
 869
 870   const Int_t kN = 36;
 871
 872   Double_t mu[kN] = { 6.084E+01, 1.676E+01, 6.863E+00
 873                     , 2.007E+00, 9.329E-01, 5.766E-01
 874                     , 4.195E-01, 2.933E-01, 2.476E-01
 875                     , 2.092E-01, 1.960E-01, 1.838E-01
 876                     , 1.763E-01, 1.703E-01, 1.651E-01
 877                     , 1.562E-01, 1.486E-01, 1.336E-01
 878                     , 1.224E-01, 1.064E-01, 9.535E-02
 879                     , 8.707E-02, 8.054E-02, 7.076E-02
 880                     , 6.362E-02, 5.688E-02, 5.173E-02
 881                     , 4.422E-02, 3.503E-02, 2.949E-02
 882                     , 2.577E-02, 2.307E-02, 1.940E-02
 883                     , 1.703E-02, 1.363E-02, 1.183E-02 };
 884
 885   Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
 886                     , 3.00000E-03, 4.00000E-03, 5.00000E-03
 887                     , 6.00000E-03, 8.00000E-03, 1.00000E-02
 888                     , 1.50000E-02, 2.00000E-02, 3.00000E-02
 889                     , 4.00000E-02, 5.00000E-02, 6.00000E-02
 890                     , 8.00000E-02, 1.00000E-01, 1.50000E-01
 891                     , 2.00000E-01, 3.00000E-01, 4.00000E-01
 892                     , 5.00000E-01, 6.00000E-01, 8.00000E-01
 893                     , 1.00000E+00, 1.25000E+00, 1.50000E+00
 894                     , 2.00000E+00, 3.00000E+00, 4.00000E+00
 895                     , 5.00000E+00, 6.00000E+00, 8.00000E+00
 896                     , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
 897
 898   return Interpolate(energyMeV,en,mu,kN);
 899
 900 }
 901
 902 //_____________________________________________________________________________
 903 Double_t AliTRDsimTR::GetMuAi(Double_t energyMeV)
 904 {
 905   //
 906   // Returns the photon absorbtion cross section for air
 907   // Implemented by Oliver Busch
 908   //
 909
 910   const Int_t kN = 38;
 911
 912   Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
 913                       0.16143E+03, 0.14250E+03, 0.15722E+03,
 914                       0.77538E+02, 0.40099E+02, 0.23313E+02,
 915                       0.98816E+01, 0.51000E+01, 0.16079E+01,
 916                       0.77536E+00, 0.35282E+00, 0.24790E+00,
 917                       0.20750E+00, 0.18703E+00, 0.16589E+00,
 918                       0.15375E+00, 0.13530E+00, 0.12311E+00,
 919                       0.10654E+00, 0.95297E-01, 0.86939E-01,
 920                       0.80390E-01, 0.70596E-01, 0.63452E-01,
 921                       0.56754E-01, 0.51644E-01, 0.44382E-01,
 922                       0.35733E-01, 0.30721E-01, 0.27450E-01,
 923                       0.25171E-01, 0.22205E-01, 0.20399E-01,
 924                       0.18053E-01, 0.18057E-01 };
 925
 926
 927
 928   Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
 929                       0.30000E-02, 0.32029E-02, 0.32029E-02,
 930                       0.40000E-02, 0.50000E-02, 0.60000E-02,
 931                       0.80000E-02, 0.10000E-01, 0.15000E-01,
 932                       0.20000E-01, 0.30000E-01, 0.40000E-01,
 933                       0.50000E-01, 0.60000E-01, 0.80000E-01,
 934                       0.10000E+00, 0.15000E+00, 0.20000E+00,
 935                       0.30000E+00, 0.40000E+00, 0.50000E+00,
 936                       0.60000E+00, 0.80000E+00, 0.10000E+01,
 937                       0.12500E+01, 0.15000E+01, 0.20000E+01,
 938                       0.30000E+01, 0.40000E+01, 0.50000E+01,
 939                       0.60000E+01, 0.80000E+01, 0.10000E+02,
 940                       0.15000E+02, 0.20000E+02 };
 941
 942   return Interpolate(energyMeV,en,mu,kN);
 943
 944 }
 945
 946 //_____________________________________________________________________________
 947 Double_t AliTRDsimTR::Interpolate(Double_t energyMeV
 948                               , Double_t *en, Double_t *mu, Int_t n)
 949 {
 950   //
 951   // Interpolates the photon absorbtion cross section
 952   // for a given energy <energyMeV>.
 953   //
 954
 955   Double_t de    = 0;
 956   Int_t    index = 0;
 957   Int_t    istat = Locate(en,n,energyMeV,index,de);
 958   if (istat == 0) {
 959     return (mu[index] - de * (mu[index]   - mu[index+1])
 960                            / (en[index+1] - en[index]  ));
 961   }
 962   else {
 963     return 0.0;
 964   }
 965
 966 }
 967
 968 //_____________________________________________________________________________
 969 Int_t AliTRDsimTR::Locate(Double_t *xv, Int_t n, Double_t xval
 970                       , Int_t &kl, Double_t &dx)
 971 {
 972   //
 973   // Locates a point (xval) in a 1-dim grid (xv(n))
 974   //
 975
 976   if (xval >= xv[n-1]) {
 977     return  1;
 978   }
 979   if (xval <  xv[0]) {
 980     return -1;
 981   }
 982
 983   Int_t km;
 984   Int_t kh = n - 1;
 985
 986   kl = 0;
 987   while (kh - kl > 1) {
 988     if (xval < xv[km = (kl+kh)/2]) {
 989       kh = km;
 990     }
 991     else {
 992       kl = km;
 993     }
 994   }
 995   if ((xval <  xv[kl])   ||
 996       (xval >  xv[kl+1]) ||
 997       (kl   >= n-1)) {
 998     AliFatal(Form("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
 999                  ,kl,xv[kl],xval,kl+1,xv[kl+1]));
1000     exit(1);
1001   }
1002
1003   dx = xval - xv[kl];
1004
1005   return 0;
1006
1007 }
1008
1009 //_____________________________________________________________________________
1010 Int_t AliTRDsimTR::SelectNFoils(Float_t p) const
1011 {
1012   //
1013   // Selects the number of foils corresponding to the momentum
1014   //
1015
1016   Int_t foils = fNFoils[fNFoilsDim-1];
1017
1018   for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
1019     if (p < fNFoilsUp[iFoil]) {
1020       foils = fNFoils[iFoil];
1021       break;
1022     }
1023   }
1024
1025   return foils;
1026
1027 }