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