]> 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
Updated version: Performance and error calculation
[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 <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
47ClassImp(AliTRDsimTR)
48
49//_____________________________________________________________________________
50AliTRDsimTR::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//_____________________________________________________________________________
83AliTRDsimTR::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//_____________________________________________________________________________
145AliTRDsimTR::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//_____________________________________________________________________________
202AliTRDsimTR::~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//_____________________________________________________________________________
226AliTRDsimTR &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//_____________________________________________________________________________
239void 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//_____________________________________________________________________________
291void 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//_____________________________________________________________________________
355Int_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//_____________________________________________________________________________
396Int_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//_____________________________________________________________________________
537void 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//_____________________________________________________________________________
556Double_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//_____________________________________________________________________________
575Double_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//_____________________________________________________________________________
614Double_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//_____________________________________________________________________________
653Double_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//_____________________________________________________________________________
700Double_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//_____________________________________________________________________________
739Double_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//_____________________________________________________________________________
778Double_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//_____________________________________________________________________________
817Double_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//_____________________________________________________________________________
856Double_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//_____________________________________________________________________________
895Double_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//_____________________________________________________________________________
939Double_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//_____________________________________________________________________________
961Int_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//_____________________________________________________________________________
1002Int_t AliTRDsimTR::SelectNFoils(Float_t p) const
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}