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