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