]> 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
Moved from Reve to Alieve.
[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
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) //
0142cb22 28// 04.06.2004 - Momentum dependent parameters implemented (CBL) //
46d29e70 29// //
30///////////////////////////////////////////////////////////////////////////////
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
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
ad4aeaf4 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
478 TVirtualMCStack *stack = gMC->GetStack();
a16be6c3 479 Int_t track = stack->GetCurrentTrackNumber();
3bf98338 480 Double_t px,py,pz,ptot;
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);
492
493 // Counter for TR analysed in custom code (e < 15keV)
494 nPhoton = 0;
495
496 for (Int_t iPhoton = 0; iPhoton < nPhCand; iPhoton++) {
497
498 // Energy of the TR photon
499 Double_t e = fSpectrum->GetRandom();
500
501 // Put TR photon on particle stack
502 if (e > 15.0 ) {
503
504 e *= 1.0e-6; // Convert it to GeV
505
506 Int_t phtrack;
507 stack-> PushTrack(1 // Must be 1
508 ,track // Identifier of the parent track, -1 for a primary
509 ,22 // Particle code.
510 ,px*e // 4 momentum (The photon is generated on the same
511 ,py*e // direction as the parent. For irregular radiator one
512 ,pz*e // can calculate also the angle but this is a secondary
513 ,e // order effect)
514 ,x,y,z,0.0 // 4 vertex
515 ,0.0,0.0,0.0 // Polarisation
516 ,kPFeedBackPhoton // Production mechanism (there is no TR in G3 so one
517 // has to make some convention)
518 ,phtrack // On output the number of the track stored
519 ,1.0
520 ,1);
521
522 }
523 // Custom treatment of TR photons
524 else {
525
526 ePhoton[nPhoton++] = e;
527
528 }
529
46d29e70 530 }
531
532 return 1;
533
534}
535
536//_____________________________________________________________________________
cb2f9e9b 537void AliTRDsimTR::SetSigma()
46d29e70 538{
539 //
540 // Sets the absorbtion crosssection for the energies of the TR spectrum
541 //
542
3bc9d03e 543 if (fSigma) {
544 delete [] fSigma;
545 }
46d29e70 546 fSigma = new Double_t[fSpNBins];
3bc9d03e 547
46d29e70 548 for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
549 Double_t energykeV = iBin * fSpBinWidth + 1.0;
550 fSigma[iBin] = Sigma(energykeV);
46d29e70 551 }
552
553}
554
555//_____________________________________________________________________________
cb2f9e9b 556Double_t AliTRDsimTR::Sigma(Double_t energykeV)
46d29e70 557{
558 //
559 // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
560 //
561
46d29e70 562 // keV -> MeV
563 Double_t energyMeV = energykeV * 0.001;
564 if (energyMeV >= 0.001) {
842287f2 565 return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
566 GetMuAi(energyMeV) * fGapDens * fGapThick * GetTemp());
46d29e70 567 }
568 else {
3bc9d03e 569 return 1.0e6;
46d29e70 570 }
571
572}
573
574//_____________________________________________________________________________
cb2f9e9b 575Double_t AliTRDsimTR::GetMuPo(Double_t energyMeV)
46d29e70 576{
577 //
578 // Returns the photon absorbtion cross section for polypropylene
579 //
580
581 const Int_t kN = 36;
582
583 Double_t mu[kN] = { 1.894E+03, 5.999E+02, 2.593E+02
584 , 7.743E+01, 3.242E+01, 1.643E+01
585 , 9.432E+00, 3.975E+00, 2.088E+00
586 , 7.452E-01, 4.315E-01, 2.706E-01
587 , 2.275E-01, 2.084E-01, 1.970E-01
588 , 1.823E-01, 1.719E-01, 1.534E-01
589 , 1.402E-01, 1.217E-01, 1.089E-01
590 , 9.947E-02, 9.198E-02, 8.078E-02
591 , 7.262E-02, 6.495E-02, 5.910E-02
592 , 5.064E-02, 4.045E-02, 3.444E-02
593 , 3.045E-02, 2.760E-02, 2.383E-02
594 , 2.145E-02, 1.819E-02, 1.658E-02 };
595
596 Double_t en[kN] = { 1.000E-03, 1.500E-03, 2.000E-03
597 , 3.000E-03, 4.000E-03, 5.000E-03
598 , 6.000E-03, 8.000E-03, 1.000E-02
599 , 1.500E-02, 2.000E-02, 3.000E-02
600 , 4.000E-02, 5.000E-02, 6.000E-02
601 , 8.000E-02, 1.000E-01, 1.500E-01
602 , 2.000E-01, 3.000E-01, 4.000E-01
603 , 5.000E-01, 6.000E-01, 8.000E-01
604 , 1.000E+00, 1.250E+00, 1.500E+00
605 , 2.000E+00, 3.000E+00, 4.000E+00
606 , 5.000E+00, 6.000E+00, 8.000E+00
607 , 1.000E+01, 1.500E+01, 2.000E+01 };
608
609 return Interpolate(energyMeV,en,mu,kN);
610
611}
612
613//_____________________________________________________________________________
cb2f9e9b 614Double_t AliTRDsimTR::GetMuCO(Double_t energyMeV)
46d29e70 615{
616 //
617 // Returns the photon absorbtion cross section for CO2
618 //
619
620 const Int_t kN = 36;
621
622 Double_t mu[kN] = { 0.39383E+04, 0.13166E+04, 0.58750E+03
623 , 0.18240E+03, 0.77996E+02, 0.40024E+02
624 , 0.23116E+02, 0.96997E+01, 0.49726E+01
625 , 0.15543E+01, 0.74915E+00, 0.34442E+00
626 , 0.24440E+00, 0.20589E+00, 0.18632E+00
627 , 0.16578E+00, 0.15394E+00, 0.13558E+00
628 , 0.12336E+00, 0.10678E+00, 0.95510E-01
629 , 0.87165E-01, 0.80587E-01, 0.70769E-01
630 , 0.63626E-01, 0.56894E-01, 0.51782E-01
631 , 0.44499E-01, 0.35839E-01, 0.30825E-01
632 , 0.27555E-01, 0.25269E-01, 0.22311E-01
633 , 0.20516E-01, 0.18184E-01, 0.17152E-01 };
634
635 Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
636 , 0.30000E-02, 0.40000E-02, 0.50000E-02
637 , 0.60000E-02, 0.80000E-02, 0.10000E-01
638 , 0.15000E-01, 0.20000E-01, 0.30000E-01
639 , 0.40000E-01, 0.50000E-01, 0.60000E-01
640 , 0.80000E-01, 0.10000E+00, 0.15000E+00
641 , 0.20000E+00, 0.30000E+00, 0.40000E+00
642 , 0.50000E+00, 0.60000E+00, 0.80000E+00
643 , 0.10000E+01, 0.12500E+01, 0.15000E+01
644 , 0.20000E+01, 0.30000E+01, 0.40000E+01
645 , 0.50000E+01, 0.60000E+01, 0.80000E+01
646 , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
647
648 return Interpolate(energyMeV,en,mu,kN);
649
650}
651
652//_____________________________________________________________________________
cb2f9e9b 653Double_t AliTRDsimTR::GetMuXe(Double_t energyMeV)
46d29e70 654{
655 //
656 // Returns the photon absorbtion cross section for xenon
657 //
658
659 const Int_t kN = 48;
660
661 Double_t mu[kN] = { 9.413E+03, 8.151E+03, 7.035E+03
662 , 7.338E+03, 4.085E+03, 2.088E+03
663 , 7.780E+02, 3.787E+02, 2.408E+02
664 , 6.941E+02, 6.392E+02, 6.044E+02
665 , 8.181E+02, 7.579E+02, 6.991E+02
666 , 8.064E+02, 6.376E+02, 3.032E+02
667 , 1.690E+02, 5.743E+01, 2.652E+01
668 , 8.930E+00, 6.129E+00, 3.316E+01
669 , 2.270E+01, 1.272E+01, 7.825E+00
670 , 3.633E+00, 2.011E+00, 7.202E-01
671 , 3.760E-01, 1.797E-01, 1.223E-01
672 , 9.699E-02, 8.281E-02, 6.696E-02
673 , 5.785E-02, 5.054E-02, 4.594E-02
674 , 4.078E-02, 3.681E-02, 3.577E-02
675 , 3.583E-02, 3.634E-02, 3.797E-02
676 , 3.987E-02, 4.445E-02, 4.815E-02 };
677
678 Double_t en[kN] = { 1.00000E-03, 1.07191E-03, 1.14900E-03
679 , 1.14900E-03, 1.50000E-03, 2.00000E-03
680 , 3.00000E-03, 4.00000E-03, 4.78220E-03
681 , 4.78220E-03, 5.00000E-03, 5.10370E-03
682 , 5.10370E-03, 5.27536E-03, 5.45280E-03
683 , 5.45280E-03, 6.00000E-03, 8.00000E-03
684 , 1.00000E-02, 1.50000E-02, 2.00000E-02
685 , 3.00000E-02, 3.45614E-02, 3.45614E-02
686 , 4.00000E-02, 5.00000E-02, 6.00000E-02
687 , 8.00000E-02, 1.00000E-01, 1.50000E-01
688 , 2.00000E-01, 3.00000E-01, 4.00000E-01
689 , 5.00000E-01, 6.00000E-01, 8.00000E-01
690 , 1.00000E+00, 1.25000E+00, 1.50000E+00
691 , 2.00000E+00, 3.00000E+00, 4.00000E+00
692 , 5.00000E+00, 6.00000E+00, 8.00000E+00
693 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
694
695 return Interpolate(energyMeV,en,mu,kN);
696
697}
698
699//_____________________________________________________________________________
cb2f9e9b 700Double_t AliTRDsimTR::GetMuBu(Double_t energyMeV)
46d29e70 701{
702 //
703 // Returns the photon absorbtion cross section for isobutane
704 //
705
706 const Int_t kN = 36;
707
708 Double_t mu[kN] = { 0.38846E+03, 0.12291E+03, 0.53225E+02
709 , 0.16091E+02, 0.69114E+01, 0.36541E+01
710 , 0.22282E+01, 0.11149E+01, 0.72887E+00
711 , 0.45053E+00, 0.38167E+00, 0.33920E+00
712 , 0.32155E+00, 0.30949E+00, 0.29960E+00
713 , 0.28317E+00, 0.26937E+00, 0.24228E+00
714 , 0.22190E+00, 0.19289E+00, 0.17288E+00
715 , 0.15789E+00, 0.14602E+00, 0.12829E+00
716 , 0.11533E+00, 0.10310E+00, 0.93790E-01
717 , 0.80117E-01, 0.63330E-01, 0.53229E-01
718 , 0.46390E-01, 0.41425E-01, 0.34668E-01
719 , 0.30267E-01, 0.23910E-01, 0.20509E-01 };
720
721 Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
722 , 0.30000E-02, 0.40000E-02, 0.50000E-02
723 , 0.60000E-02, 0.80000E-02, 0.10000E-01
724 , 0.15000E-01, 0.20000E-01, 0.30000E-01
725 , 0.40000E-01, 0.50000E-01, 0.60000E-01
726 , 0.80000E-01, 0.10000E+00, 0.15000E+00
727 , 0.20000E+00, 0.30000E+00, 0.40000E+00
728 , 0.50000E+00, 0.60000E+00, 0.80000E+00
729 , 0.10000E+01, 0.12500E+01, 0.15000E+01
730 , 0.20000E+01, 0.30000E+01, 0.40000E+01
731 , 0.50000E+01, 0.60000E+01, 0.80000E+01
732 , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
733
734 return Interpolate(energyMeV,en,mu,kN);
735
736}
737
738//_____________________________________________________________________________
cb2f9e9b 739Double_t AliTRDsimTR::GetMuMy(Double_t energyMeV)
46d29e70 740{
741 //
742 // Returns the photon absorbtion cross section for mylar
743 //
744
745 const Int_t kN = 36;
746
747 Double_t mu[kN] = { 2.911E+03, 9.536E+02, 4.206E+02
748 , 1.288E+02, 5.466E+01, 2.792E+01
749 , 1.608E+01, 6.750E+00, 3.481E+00
750 , 1.132E+00, 5.798E-01, 3.009E-01
751 , 2.304E-01, 2.020E-01, 1.868E-01
752 , 1.695E-01, 1.586E-01, 1.406E-01
753 , 1.282E-01, 1.111E-01, 9.947E-02
754 , 9.079E-02, 8.395E-02, 7.372E-02
755 , 6.628E-02, 5.927E-02, 5.395E-02
756 , 4.630E-02, 3.715E-02, 3.181E-02
757 , 2.829E-02, 2.582E-02, 2.257E-02
758 , 2.057E-02, 1.789E-02, 1.664E-02 };
759
760 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
761 , 3.00000E-03, 4.00000E-03, 5.00000E-03
762 , 6.00000E-03, 8.00000E-03, 1.00000E-02
763 , 1.50000E-02, 2.00000E-02, 3.00000E-02
764 , 4.00000E-02, 5.00000E-02, 6.00000E-02
765 , 8.00000E-02, 1.00000E-01, 1.50000E-01
766 , 2.00000E-01, 3.00000E-01, 4.00000E-01
767 , 5.00000E-01, 6.00000E-01, 8.00000E-01
768 , 1.00000E+00, 1.25000E+00, 1.50000E+00
769 , 2.00000E+00, 3.00000E+00, 4.00000E+00
770 , 5.00000E+00, 6.00000E+00, 8.00000E+00
771 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
772
773 return Interpolate(energyMeV,en,mu,kN);
774
775}
776
777//_____________________________________________________________________________
cb2f9e9b 778Double_t AliTRDsimTR::GetMuN2(Double_t energyMeV)
46d29e70 779{
780 //
781 // Returns the photon absorbtion cross section for nitrogen
782 //
783
784 const Int_t kN = 36;
785
786 Double_t mu[kN] = { 3.311E+03, 1.083E+03, 4.769E+02
787 , 1.456E+02, 6.166E+01, 3.144E+01
788 , 1.809E+01, 7.562E+00, 3.879E+00
789 , 1.236E+00, 6.178E-01, 3.066E-01
790 , 2.288E-01, 1.980E-01, 1.817E-01
791 , 1.639E-01, 1.529E-01, 1.353E-01
792 , 1.233E-01, 1.068E-01, 9.557E-02
793 , 8.719E-02, 8.063E-02, 7.081E-02
794 , 6.364E-02, 5.693E-02, 5.180E-02
795 , 4.450E-02, 3.579E-02, 3.073E-02
796 , 2.742E-02, 2.511E-02, 2.209E-02
797 , 2.024E-02, 1.782E-02, 1.673E-02 };
798
799 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
800 , 3.00000E-03, 4.00000E-03, 5.00000E-03
801 , 6.00000E-03, 8.00000E-03, 1.00000E-02
802 , 1.50000E-02, 2.00000E-02, 3.00000E-02
803 , 4.00000E-02, 5.00000E-02, 6.00000E-02
804 , 8.00000E-02, 1.00000E-01, 1.50000E-01
805 , 2.00000E-01, 3.00000E-01, 4.00000E-01
806 , 5.00000E-01, 6.00000E-01, 8.00000E-01
807 , 1.00000E+00, 1.25000E+00, 1.50000E+00
808 , 2.00000E+00, 3.00000E+00, 4.00000E+00
809 , 5.00000E+00, 6.00000E+00, 8.00000E+00
810 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
811
812 return Interpolate(energyMeV,en,mu,kN);
813
814}
815
816//_____________________________________________________________________________
cb2f9e9b 817Double_t AliTRDsimTR::GetMuO2(Double_t energyMeV)
46d29e70 818{
819 //
820 // Returns the photon absorbtion cross section for oxygen
821 //
822
823 const Int_t kN = 36;
824
825 Double_t mu[kN] = { 4.590E+03, 1.549E+03, 6.949E+02
826 , 2.171E+02, 9.315E+01, 4.790E+01
827 , 2.770E+01, 1.163E+01, 5.952E+00
828 , 1.836E+00, 8.651E-01, 3.779E-01
829 , 2.585E-01, 2.132E-01, 1.907E-01
830 , 1.678E-01, 1.551E-01, 1.361E-01
831 , 1.237E-01, 1.070E-01, 9.566E-02
832 , 8.729E-02, 8.070E-02, 7.087E-02
833 , 6.372E-02, 5.697E-02, 5.185E-02
834 , 4.459E-02, 3.597E-02, 3.100E-02
835 , 2.777E-02, 2.552E-02, 2.263E-02
836 , 2.089E-02, 1.866E-02, 1.770E-02 };
837
838 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
839 , 3.00000E-03, 4.00000E-03, 5.00000E-03
840 , 6.00000E-03, 8.00000E-03, 1.00000E-02
841 , 1.50000E-02, 2.00000E-02, 3.00000E-02
842 , 4.00000E-02, 5.00000E-02, 6.00000E-02
843 , 8.00000E-02, 1.00000E-01, 1.50000E-01
844 , 2.00000E-01, 3.00000E-01, 4.00000E-01
845 , 5.00000E-01, 6.00000E-01, 8.00000E-01
846 , 1.00000E+00, 1.25000E+00, 1.50000E+00
847 , 2.00000E+00, 3.00000E+00, 4.00000E+00
848 , 5.00000E+00, 6.00000E+00, 8.00000E+00
849 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
850
851 return Interpolate(energyMeV,en,mu,kN);
852
853}
854
855//_____________________________________________________________________________
cb2f9e9b 856Double_t AliTRDsimTR::GetMuHe(Double_t energyMeV)
46d29e70 857{
858 //
859 // Returns the photon absorbtion cross section for helium
860 //
861
862 const Int_t kN = 36;
863
864 Double_t mu[kN] = { 6.084E+01, 1.676E+01, 6.863E+00
865 , 2.007E+00, 9.329E-01, 5.766E-01
866 , 4.195E-01, 2.933E-01, 2.476E-01
867 , 2.092E-01, 1.960E-01, 1.838E-01
868 , 1.763E-01, 1.703E-01, 1.651E-01
869 , 1.562E-01, 1.486E-01, 1.336E-01
870 , 1.224E-01, 1.064E-01, 9.535E-02
871 , 8.707E-02, 8.054E-02, 7.076E-02
872 , 6.362E-02, 5.688E-02, 5.173E-02
873 , 4.422E-02, 3.503E-02, 2.949E-02
874 , 2.577E-02, 2.307E-02, 1.940E-02
875 , 1.703E-02, 1.363E-02, 1.183E-02 };
876
877 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
878 , 3.00000E-03, 4.00000E-03, 5.00000E-03
879 , 6.00000E-03, 8.00000E-03, 1.00000E-02
880 , 1.50000E-02, 2.00000E-02, 3.00000E-02
881 , 4.00000E-02, 5.00000E-02, 6.00000E-02
882 , 8.00000E-02, 1.00000E-01, 1.50000E-01
883 , 2.00000E-01, 3.00000E-01, 4.00000E-01
884 , 5.00000E-01, 6.00000E-01, 8.00000E-01
885 , 1.00000E+00, 1.25000E+00, 1.50000E+00
886 , 2.00000E+00, 3.00000E+00, 4.00000E+00
887 , 5.00000E+00, 6.00000E+00, 8.00000E+00
888 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
889
890 return Interpolate(energyMeV,en,mu,kN);
891
892}
893
842287f2 894//_____________________________________________________________________________
cb2f9e9b 895Double_t AliTRDsimTR::GetMuAi(Double_t energyMeV)
842287f2 896{
897 //
898 // Returns the photon absorbtion cross section for air
899 // Implemented by Oliver Busch
900 //
901
902 const Int_t kN = 38;
903
904 Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
905 0.16143E+03, 0.14250E+03, 0.15722E+03,
906 0.77538E+02, 0.40099E+02, 0.23313E+02,
907 0.98816E+01, 0.51000E+01, 0.16079E+01,
908 0.77536E+00, 0.35282E+00, 0.24790E+00,
909 0.20750E+00, 0.18703E+00, 0.16589E+00,
910 0.15375E+00, 0.13530E+00, 0.12311E+00,
911 0.10654E+00, 0.95297E-01, 0.86939E-01,
912 0.80390E-01, 0.70596E-01, 0.63452E-01,
913 0.56754E-01, 0.51644E-01, 0.44382E-01,
914 0.35733E-01, 0.30721E-01, 0.27450E-01,
915 0.25171E-01, 0.22205E-01, 0.20399E-01,
916 0.18053E-01, 0.18057E-01 };
917
918
919
920 Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
921 0.30000E-02, 0.32029E-02, 0.32029E-02,
922 0.40000E-02, 0.50000E-02, 0.60000E-02,
923 0.80000E-02, 0.10000E-01, 0.15000E-01,
924 0.20000E-01, 0.30000E-01, 0.40000E-01,
925 0.50000E-01, 0.60000E-01, 0.80000E-01,
926 0.10000E+00, 0.15000E+00, 0.20000E+00,
927 0.30000E+00, 0.40000E+00, 0.50000E+00,
928 0.60000E+00, 0.80000E+00, 0.10000E+01,
929 0.12500E+01, 0.15000E+01, 0.20000E+01,
930 0.30000E+01, 0.40000E+01, 0.50000E+01,
931 0.60000E+01, 0.80000E+01, 0.10000E+02,
932 0.15000E+02, 0.20000E+02 };
933
934 return Interpolate(energyMeV,en,mu,kN);
935
936}
937
46d29e70 938//_____________________________________________________________________________
cb2f9e9b 939Double_t AliTRDsimTR::Interpolate(Double_t energyMeV
46d29e70 940 , Double_t *en, Double_t *mu, Int_t n)
941{
942 //
943 // Interpolates the photon absorbtion cross section
944 // for a given energy <energyMeV>.
945 //
946
947 Double_t de = 0;
948 Int_t index = 0;
949 Int_t istat = Locate(en,n,energyMeV,index,de);
950 if (istat == 0) {
951 return (mu[index] - de * (mu[index] - mu[index+1])
952 / (en[index+1] - en[index] ));
953 }
954 else {
955 return 0.0;
956 }
957
958}
959
960//_____________________________________________________________________________
cb2f9e9b 961Int_t AliTRDsimTR::Locate(Double_t *xv, Int_t n, Double_t xval
46d29e70 962 , Int_t &kl, Double_t &dx)
963{
964 //
965 // Locates a point (xval) in a 1-dim grid (xv(n))
966 //
967
3bc9d03e 968 if (xval >= xv[n-1]) {
969 return 1;
970 }
971 if (xval < xv[0]) {
972 return -1;
973 }
46d29e70 974
975 Int_t km;
976 Int_t kh = n - 1;
977
978 kl = 0;
979 while (kh - kl > 1) {
3bc9d03e 980 if (xval < xv[km = (kl+kh)/2]) {
981 kh = km;
982 }
983 else {
984 kl = km;
985 }
46d29e70 986 }
3bc9d03e 987 if ((xval < xv[kl]) ||
988 (xval > xv[kl+1]) ||
989 (kl >= n-1)) {
a16be6c3 990 AliFatal(Form("Locate failed xv[%d] %f xval %f xv[%d] %f!!!\n"
3bc9d03e 991 ,kl,xv[kl],xval,kl+1,xv[kl+1]));
46d29e70 992 exit(1);
993 }
994
995 dx = xval - xv[kl];
996
997 return 0;
998
999}
0142cb22 1000
1001//_____________________________________________________________________________
cb2f9e9b 1002Int_t AliTRDsimTR::SelectNFoils(Float_t p)
0142cb22 1003{
1004 //
1005 // Selects the number of foils corresponding to the momentum
1006 //
1007
1008 Int_t foils = fNFoils[fNFoilsDim-1];
1009
1010 for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
1011 if (p < fNFoilsUp[iFoil]) {
1012 foils = fNFoils[iFoil];
1013 break;
1014 }
1015 }
1016
1017 return foils;
1018
1019}