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