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