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