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