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