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