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