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