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