]> 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
Copy arrays in assignment instead of the pointer; avoid double delete.
[u/mrichter/AliRoot.git] / TRD / AliTRDsimTR.cxx
CommitLineData
99336540 1
46d29e70 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
88cb7938 17/* $Id$ */
46d29e70 18
c8ab4518 19////////////////////////////////////////////////////////////////////////////
20// //
21// TRD simulation - multimodule (regular rad.) //
22// after: M. CASTELLANO et al., COMP. PHYS. COMM. 51 (1988) 431 //
23// + COMP. PHYS. COMM. 61 (1990) 395 //
24// //
25// 17.07.1998 - A.Andronic //
26// 08.12.1998 - simplified version //
27// 11.07.2000 - Adapted code to aliroot environment (C.Blume) //
28// 04.06.2004 - Momentum dependent parameters implemented (CBL) //
29// //
30////////////////////////////////////////////////////////////////////////////
46d29e70 31
32#include <stdlib.h>
33
0e9c2ad5 34#include <TH1.h>
35#include <TRandom.h>
36#include <TMath.h>
37#include <TParticle.h>
a16be6c3 38#include <TVirtualMC.h>
39#include <TVirtualMCStack.h>
46d29e70 40
46d29e70 41#include "AliModule.h"
3bc9d03e 42#include "AliLog.h"
a16be6c3 43#include "AliMC.h"
46d29e70 44
cb2f9e9b 45#include "AliTRDsimTR.h"
0e9c2ad5 46
cb2f9e9b 47ClassImp(AliTRDsimTR)
46d29e70 48
49//_____________________________________________________________________________
cb2f9e9b 50AliTRDsimTR::AliTRDsimTR()
3bc9d03e 51 :TObject()
52 ,fNFoilsDim(0)
53 ,fNFoils(0)
54 ,fNFoilsUp(0)
55 ,fFoilThick(0)
56 ,fGapThick(0)
57 ,fFoilDens(0)
58 ,fGapDens(0)
59 ,fFoilOmega(0)
60 ,fGapOmega()
61 ,fFoilZ(0)
62 ,fGapZ(0)
63 ,fFoilA(0)
64 ,fGapA(0)
65 ,fTemp(0)
66 ,fSpNBins(0)
67 ,fSpRange(0)
68 ,fSpBinWidth(0)
69 ,fSpLower(0)
70 ,fSpUpper(0)
71 ,fSigma(0)
72 ,fSpectrum(0)
46d29e70 73{
74 //
cb2f9e9b 75 // AliTRDsimTR default constructor
46d29e70 76 //
77
78 Init();
79
80}
81
82//_____________________________________________________________________________
cb2f9e9b 83AliTRDsimTR::AliTRDsimTR(AliModule *mod, Int_t foil, Int_t gap)
3bc9d03e 84 :TObject()
85 ,fNFoilsDim(0)
86 ,fNFoils(0)
87 ,fNFoilsUp(0)
88 ,fFoilThick(0)
89 ,fGapThick(0)
90 ,fFoilDens(0)
91 ,fGapDens(0)
92 ,fFoilOmega(0)
93 ,fGapOmega()
94 ,fFoilZ(0)
95 ,fGapZ(0)
96 ,fFoilA(0)
97 ,fGapA(0)
98 ,fTemp(0)
99 ,fSpNBins(0)
100 ,fSpRange(0)
101 ,fSpBinWidth(0)
102 ,fSpLower(0)
103 ,fSpUpper(0)
104 ,fSigma(0)
105 ,fSpectrum(0)
46d29e70 106{
107 //
cb2f9e9b 108 // AliTRDsimTR constructor. Takes the material properties of the radiator
46d29e70 109 // foils and the gas in the gaps from AliModule <mod>.
110 // The default number of foils is 100 with a thickness of 20 mu. The
111 // thickness of the gaps is 500 mu.
112 //
113
3bc9d03e 114 Float_t aFoil;
115 Float_t zFoil;
116 Float_t rhoFoil;
117
118 Float_t aGap;
119 Float_t zGap;
120 Float_t rhoGap;
46d29e70 121
3bc9d03e 122 Float_t rad;
123 Float_t abs;
124
125 Char_t name[21];
fa5e892a 126
46d29e70 127 Init();
128
129 mod->AliGetMaterial(foil,name,aFoil,zFoil,rhoFoil,rad,abs);
130 mod->AliGetMaterial(gap ,name,aGap ,zGap ,rhoGap ,rad,abs);
131
132 fFoilDens = rhoFoil;
133 fFoilA = aFoil;
134 fFoilZ = zFoil;
135 fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA);
136
137 fGapDens = rhoGap;
138 fGapA = aGap;
139 fGapZ = zGap;
140 fGapOmega = Omega(fGapDens ,fGapZ ,fGapA );
141
142}
143
144//_____________________________________________________________________________
cb2f9e9b 145AliTRDsimTR::AliTRDsimTR(const AliTRDsimTR &s)
3bc9d03e 146 :TObject(s)
147 ,fNFoilsDim(s.fNFoilsDim)
148 ,fNFoils(0)
149 ,fNFoilsUp(0)
150 ,fFoilThick(s.fFoilThick)
151 ,fGapThick(s.fGapThick)
152 ,fFoilDens(s.fFoilDens)
153 ,fGapDens(s.fGapDens)
154 ,fFoilOmega(s.fFoilOmega)
155 ,fGapOmega(s.fGapOmega)
156 ,fFoilZ(s.fFoilZ)
157 ,fGapZ(s.fGapZ)
158 ,fFoilA(s.fFoilA)
159 ,fGapA(s.fGapA)
160 ,fTemp(s.fTemp)
161 ,fSpNBins(s.fSpNBins)
162 ,fSpRange(s.fSpRange)
163 ,fSpBinWidth(s.fSpBinWidth)
164 ,fSpLower(s.fSpLower)
165 ,fSpUpper(s.fSpUpper)
166 ,fSigma(0)
167 ,fSpectrum(0)
46d29e70 168{
169 //
cb2f9e9b 170 // AliTRDsimTR copy constructor
46d29e70 171 //
172
cb2f9e9b 173 if (((AliTRDsimTR &) s).fNFoils) {
174 delete [] ((AliTRDsimTR &) s).fNFoils;
3bc9d03e 175 }
cb2f9e9b 176 ((AliTRDsimTR &) s).fNFoils = new Int_t[fNFoilsDim];
3bc9d03e 177 for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
cb2f9e9b 178 ((AliTRDsimTR &) s).fNFoils[iFoil] = fNFoils[iFoil];
3bc9d03e 179 }
180
cb2f9e9b 181 if (((AliTRDsimTR &) s).fNFoilsUp) {
182 delete [] ((AliTRDsimTR &) s).fNFoilsUp;
3bc9d03e 183 }
cb2f9e9b 184 ((AliTRDsimTR &) s).fNFoilsUp = new Double_t[fNFoilsDim];
3bc9d03e 185 for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
cb2f9e9b 186 ((AliTRDsimTR &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
3bc9d03e 187 }
188
cb2f9e9b 189 if (((AliTRDsimTR &) s).fSigma) {
190 delete [] ((AliTRDsimTR &) s).fSigma;
3bc9d03e 191 }
cb2f9e9b 192 ((AliTRDsimTR &) s).fSigma = new Double_t[fSpNBins];
3bc9d03e 193 for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
cb2f9e9b 194 ((AliTRDsimTR &) s).fSigma[iBin] = fSigma[iBin];
3bc9d03e 195 }
196
cb2f9e9b 197 fSpectrum->Copy(*((AliTRDsimTR &) s).fSpectrum);
46d29e70 198
199}
200
201//_____________________________________________________________________________
cb2f9e9b 202AliTRDsimTR::~AliTRDsimTR()
46d29e70 203{
204 //
cb2f9e9b 205 // AliTRDsimTR destructor
46d29e70 206 //
207
3bc9d03e 208 if (fSigma) {
209 delete [] fSigma;
210 fSigma = 0;
211 }
212
213 if (fNFoils) {
214 delete [] fNFoils;
215 fNFoils = 0;
216 }
217
218 if (fNFoilsUp) {
219 delete [] fNFoilsUp;
220 fNFoilsUp = 0;
221 }
46d29e70 222
120272f4 223 if (fSpectrum) {
224 delete fSpectrum;
225 fSpectrum = 0;
226 }
227
46d29e70 228}
229
230//_____________________________________________________________________________
cb2f9e9b 231AliTRDsimTR &AliTRDsimTR::operator=(const AliTRDsimTR &s)
46d29e70 232{
233 //
234 // Assignment operator
235 //
236
cb2f9e9b 237 if (this != &s) ((AliTRDsimTR &) s).Copy(*this);
3bc9d03e 238
46d29e70 239 return *this;
240
241}
242
243//_____________________________________________________________________________
cb2f9e9b 244void AliTRDsimTR::Copy(TObject &s) const
46d29e70 245{
246 //
247 // Copy function
248 //
249
cb2f9e9b 250 ((AliTRDsimTR &) s).fFoilThick = fFoilThick;
251 ((AliTRDsimTR &) s).fFoilDens = fFoilDens;
252 ((AliTRDsimTR &) s).fFoilOmega = fFoilOmega;
253 ((AliTRDsimTR &) s).fFoilZ = fFoilZ;
254 ((AliTRDsimTR &) s).fFoilA = fFoilA;
255 ((AliTRDsimTR &) s).fGapThick = fGapThick;
256 ((AliTRDsimTR &) s).fGapDens = fGapDens;
257 ((AliTRDsimTR &) s).fGapOmega = fGapOmega;
258 ((AliTRDsimTR &) s).fGapZ = fGapZ;
259 ((AliTRDsimTR &) s).fGapA = fGapA;
260 ((AliTRDsimTR &) s).fTemp = fTemp;
261 ((AliTRDsimTR &) s).fSpNBins = fSpNBins;
262 ((AliTRDsimTR &) s).fSpRange = fSpRange;
263 ((AliTRDsimTR &) s).fSpBinWidth = fSpBinWidth;
264 ((AliTRDsimTR &) s).fSpLower = fSpLower;
265 ((AliTRDsimTR &) s).fSpUpper = fSpUpper;
266
267 if (((AliTRDsimTR &) s).fNFoils) {
268 delete [] ((AliTRDsimTR &) s).fNFoils;
3bc9d03e 269 }
cb2f9e9b 270 ((AliTRDsimTR &) s).fNFoils = new Int_t[fNFoilsDim];
0142cb22 271 for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
cb2f9e9b 272 ((AliTRDsimTR &) s).fNFoils[iFoil] = fNFoils[iFoil];
0142cb22 273 }
274
cb2f9e9b 275 if (((AliTRDsimTR &) s).fNFoilsUp) {
276 delete [] ((AliTRDsimTR &) s).fNFoilsUp;
3bc9d03e 277 }
cb2f9e9b 278 ((AliTRDsimTR &) s).fNFoilsUp = new Double_t[fNFoilsDim];
0142cb22 279 for (Int_t iFoil = 0; iFoil < fNFoilsDim; iFoil++) {
cb2f9e9b 280 ((AliTRDsimTR &) s).fNFoilsUp[iFoil] = fNFoilsUp[iFoil];
0142cb22 281 }
282
cb2f9e9b 283 if (((AliTRDsimTR &) s).fSigma) {
284 delete [] ((AliTRDsimTR &) s).fSigma;
3bc9d03e 285 }
cb2f9e9b 286 ((AliTRDsimTR &) s).fSigma = new Double_t[fSpNBins];
46d29e70 287 for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
cb2f9e9b 288 ((AliTRDsimTR &) s).fSigma[iBin] = fSigma[iBin];
46d29e70 289 }
290
cb2f9e9b 291 fSpectrum->Copy(*((AliTRDsimTR &) s).fSpectrum);
46d29e70 292
293}
294
295//_____________________________________________________________________________
cb2f9e9b 296void AliTRDsimTR::Init()
46d29e70 297{
298 //
299 // Initialization
0142cb22 300 // The default radiator are prolypropilene foils of 10 mu thickness
301 // with gaps of 80 mu filled with N2.
46d29e70 302 //
303
0142cb22 304 fNFoilsDim = 7;
305
3bc9d03e 306 if (fNFoils) {
307 delete [] fNFoils;
308 }
0142cb22 309 fNFoils = new Int_t[fNFoilsDim];
310 fNFoils[0] = 170;
3bc9d03e 311 fNFoils[1] = 225;
312 fNFoils[2] = 275;
313 fNFoils[3] = 305;
314 fNFoils[4] = 325;
315 fNFoils[5] = 340;
316 fNFoils[6] = 350;
317
318 if (fNFoilsUp) {
319 delete [] fNFoilsUp;
320 }
0142cb22 321 fNFoilsUp = new Double_t[fNFoilsDim];
322 fNFoilsUp[0] = 1.25;
323 fNFoilsUp[1] = 1.75;
324 fNFoilsUp[2] = 2.50;
325 fNFoilsUp[3] = 3.50;
326 fNFoilsUp[4] = 4.50;
327 fNFoilsUp[5] = 5.50;
328 fNFoilsUp[6] = 10000.0;
46d29e70 329
db30bf0f 330 fFoilThick = 0.0013;
46d29e70 331 fFoilDens = 0.92;
332 fFoilZ = 5.28571;
333 fFoilA = 10.4286;
334 fFoilOmega = Omega(fFoilDens,fFoilZ,fFoilA);
335
db30bf0f 336 fGapThick = 0.0060;
0142cb22 337 fGapDens = 0.00125;
338 fGapZ = 7.0;
339 fGapA = 14.00674;
46d29e70 340 fGapOmega = Omega(fGapDens ,fGapZ ,fGapA );
341
342 fTemp = 293.16;
343
344 fSpNBins = 200;
345 fSpRange = 100;
346 fSpBinWidth = fSpRange / fSpNBins;
347 fSpLower = 1.0 - 0.5 * fSpBinWidth;
348 fSpUpper = fSpLower + fSpRange;
349
350 if (fSpectrum) delete fSpectrum;
351 fSpectrum = new TH1D("TRspectrum","TR spectrum",fSpNBins,fSpLower,fSpUpper);
abaf1f1d 352 fSpectrum->SetDirectory(0);
46d29e70 353
354 // Set the sigma values
355 SetSigma();
356
357}
358
359//_____________________________________________________________________________
cb2f9e9b 360Int_t AliTRDsimTR::CreatePhotons(Int_t pdg, Float_t p
46d29e70 361 , Int_t &nPhoton, Float_t *ePhoton)
362{
363 //
364 // Create TRD photons for a charged particle of type <pdg> with the total
365 // momentum <p>.
366 // Number of produced TR photons: <nPhoton>
367 // Energies of the produced TR photons: <ePhoton>
368 //
369
370 // PDG codes
371 const Int_t kPdgEle = 11;
372 const Int_t kPdgMuon = 13;
373 const Int_t kPdgPion = 211;
374 const Int_t kPdgKaon = 321;
375
376 Float_t mass = 0;
377 switch (TMath::Abs(pdg)) {
378 case kPdgEle:
379 mass = 5.11e-4;
380 break;
381 case kPdgMuon:
382 mass = 0.10566;
383 break;
384 case kPdgPion:
385 mass = 0.13957;
386 break;
387 case kPdgKaon:
388 mass = 0.4937;
389 break;
390 default:
391 return 0;
392 break;
393 };
394
46d29e70 395 // Calculate the TR photons
0142cb22 396 return TrPhotons(p, mass, nPhoton, ePhoton);
46d29e70 397
398}
399
400//_____________________________________________________________________________
cb2f9e9b 401Int_t AliTRDsimTR::TrPhotons(Float_t p, Float_t mass
0142cb22 402 , Int_t &nPhoton, Float_t *ePhoton)
46d29e70 403{
404 //
a16be6c3 405 // Produces TR photons using a parametric model for regular radiator. Photons
406 // with energy larger than 15 keV are included in the MC stack and tracked by VMC
407 // machinary.
46d29e70 408 //
a16be6c3 409 // Input parameters:
410 // p - parent momentum [GeV/c]
411 // mass - parent mass
412 //
413 // Output :
414 // nPhoton - number of photons which have to be processed by custom code
415 // ePhoton - energy of this photons in keV.
416 //
46d29e70 417
418 const Double_t kAlpha = 0.0072973;
99336540 419 const Int_t kSumMax = 30;
420
3bc9d03e 421 Double_t tau = fGapThick / fFoilThick;
46d29e70 422
0142cb22 423 // Calculate gamma
424 Double_t gamma = TMath::Sqrt(p*p + mass*mass) / mass;
425
426 // Select the number of foils corresponding to momentum
427 Int_t foils = SelectNFoils(p);
428
46d29e70 429 fSpectrum->Reset();
430
431 // The TR spectrum
3bc9d03e 432 Double_t csi1;
433 Double_t csi2;
434 Double_t rho1;
435 Double_t rho2;
ad4aeaf4 436 Double_t sigma;
3bc9d03e 437 Double_t sum;
438 Double_t nEqu;
439 Double_t thetaN;
440 Double_t aux;
441 Double_t energyeV;
442 Double_t energykeV;
3bc9d03e 443 for (Int_t iBin = 1; iBin <= fSpNBins; iBin++) {
46d29e70 444
3bc9d03e 445 energykeV = fSpectrum->GetBinCenter(iBin);
a16be6c3 446 energyeV = energykeV * 1.0e3;
447
ad4aeaf4 448 sigma = Sigma(energykeV);
3bc9d03e 449
450 csi1 = fFoilOmega / energyeV;
451 csi2 = fGapOmega / energyeV;
452
453 rho1 = 2.5 * energyeV * fFoilThick * 1.0e4
a16be6c3 454 * (1.0 / (gamma*gamma) + csi1*csi1);
3bc9d03e 455 rho2 = 2.5 * energyeV * fFoilThick * 1.0e4
a16be6c3 456 * (1.0 / (gamma*gamma) + csi2 *csi2);
46d29e70 457
458 // Calculate the sum
a16be6c3 459 sum = 0.0;
99336540 460 for (Int_t n = 1; n <= kSumMax; n++) {
3bc9d03e 461 thetaN = (TMath::Pi() * 2.0 * n - (rho1 + tau * rho2)) / (1.0 + tau);
462 if (thetaN < 0.0) {
463 thetaN = 0.0;
464 }
a16be6c3 465 aux = 1.0 / (rho1 + thetaN) - 1.0 / (rho2 + thetaN);
3bc9d03e 466 sum += thetaN * (aux*aux) * (1.0 - TMath::Cos(rho1 + thetaN));
46d29e70 467 }
468
99336540 469 // Equivalent number of foils
ad4aeaf4 470 nEqu = (1.0 - TMath::Exp(-foils * sigma)) / (1.0 - TMath::Exp(-sigma));
46d29e70 471
472 // dN / domega
a16be6c3 473 fSpectrum->SetBinContent(iBin,4.0 * kAlpha * nEqu * sum / (energykeV * (1.0 + tau)));
3bc9d03e 474
46d29e70 475 }
476
477 // <nTR> (binsize corr.)
a16be6c3 478 Float_t nTr = fSpBinWidth * fSpectrum->Integral();
479 // Number of TR photons from Poisson distribution with mean <nTr>
480 Int_t nPhCand = gRandom->Poisson(nTr);
481
482 // Link the MC stack and get info about parent electron
483 TVirtualMCStack *stack = gMC->GetStack();
a16be6c3 484 Int_t track = stack->GetCurrentTrackNumber();
3bf98338 485 Double_t px,py,pz,ptot;
486 gMC->TrackMomentum(px,py,pz,ptot);
487 ptot = TMath::Sqrt(px*px+py*py+pz*pz);
488 px /= ptot;
489 py /= ptot;
490 pz /= ptot;
a16be6c3 491
492 // Current position of electron
493 Double_t x;
494 Double_t y;
495 Double_t z;
496 gMC->TrackPosition(x,y,z);
497
498 // Counter for TR analysed in custom code (e < 15keV)
499 nPhoton = 0;
500
501 for (Int_t iPhoton = 0; iPhoton < nPhCand; iPhoton++) {
502
503 // Energy of the TR photon
504 Double_t e = fSpectrum->GetRandom();
505
506 // Put TR photon on particle stack
507 if (e > 15.0 ) {
508
509 e *= 1.0e-6; // Convert it to GeV
510
511 Int_t phtrack;
512 stack-> PushTrack(1 // Must be 1
513 ,track // Identifier of the parent track, -1 for a primary
514 ,22 // Particle code.
515 ,px*e // 4 momentum (The photon is generated on the same
516 ,py*e // direction as the parent. For irregular radiator one
517 ,pz*e // can calculate also the angle but this is a secondary
518 ,e // order effect)
519 ,x,y,z,0.0 // 4 vertex
520 ,0.0,0.0,0.0 // Polarisation
521 ,kPFeedBackPhoton // Production mechanism (there is no TR in G3 so one
522 // has to make some convention)
523 ,phtrack // On output the number of the track stored
524 ,1.0
525 ,1);
526
527 }
528 // Custom treatment of TR photons
529 else {
530
531 ePhoton[nPhoton++] = e;
532
533 }
534
46d29e70 535 }
536
537 return 1;
538
539}
540
541//_____________________________________________________________________________
cb2f9e9b 542void AliTRDsimTR::SetSigma()
46d29e70 543{
544 //
545 // Sets the absorbtion crosssection for the energies of the TR spectrum
546 //
547
3bc9d03e 548 if (fSigma) {
549 delete [] fSigma;
550 }
46d29e70 551 fSigma = new Double_t[fSpNBins];
3bc9d03e 552
46d29e70 553 for (Int_t iBin = 0; iBin < fSpNBins; iBin++) {
554 Double_t energykeV = iBin * fSpBinWidth + 1.0;
555 fSigma[iBin] = Sigma(energykeV);
46d29e70 556 }
557
558}
559
560//_____________________________________________________________________________
cb2f9e9b 561Double_t AliTRDsimTR::Sigma(Double_t energykeV)
46d29e70 562{
563 //
564 // Calculates the absorbtion crosssection for a one-foil-one-gap-radiator
565 //
566
46d29e70 567 // keV -> MeV
568 Double_t energyMeV = energykeV * 0.001;
569 if (energyMeV >= 0.001) {
842287f2 570 return(GetMuPo(energyMeV) * fFoilDens * fFoilThick +
571 GetMuAi(energyMeV) * fGapDens * fGapThick * GetTemp());
46d29e70 572 }
573 else {
3bc9d03e 574 return 1.0e6;
46d29e70 575 }
576
577}
578
579//_____________________________________________________________________________
cb2f9e9b 580Double_t AliTRDsimTR::GetMuPo(Double_t energyMeV)
46d29e70 581{
582 //
583 // Returns the photon absorbtion cross section for polypropylene
584 //
585
586 const Int_t kN = 36;
587
588 Double_t mu[kN] = { 1.894E+03, 5.999E+02, 2.593E+02
589 , 7.743E+01, 3.242E+01, 1.643E+01
590 , 9.432E+00, 3.975E+00, 2.088E+00
591 , 7.452E-01, 4.315E-01, 2.706E-01
592 , 2.275E-01, 2.084E-01, 1.970E-01
593 , 1.823E-01, 1.719E-01, 1.534E-01
594 , 1.402E-01, 1.217E-01, 1.089E-01
595 , 9.947E-02, 9.198E-02, 8.078E-02
596 , 7.262E-02, 6.495E-02, 5.910E-02
597 , 5.064E-02, 4.045E-02, 3.444E-02
598 , 3.045E-02, 2.760E-02, 2.383E-02
599 , 2.145E-02, 1.819E-02, 1.658E-02 };
600
601 Double_t en[kN] = { 1.000E-03, 1.500E-03, 2.000E-03
602 , 3.000E-03, 4.000E-03, 5.000E-03
603 , 6.000E-03, 8.000E-03, 1.000E-02
604 , 1.500E-02, 2.000E-02, 3.000E-02
605 , 4.000E-02, 5.000E-02, 6.000E-02
606 , 8.000E-02, 1.000E-01, 1.500E-01
607 , 2.000E-01, 3.000E-01, 4.000E-01
608 , 5.000E-01, 6.000E-01, 8.000E-01
609 , 1.000E+00, 1.250E+00, 1.500E+00
610 , 2.000E+00, 3.000E+00, 4.000E+00
611 , 5.000E+00, 6.000E+00, 8.000E+00
612 , 1.000E+01, 1.500E+01, 2.000E+01 };
613
614 return Interpolate(energyMeV,en,mu,kN);
615
616}
617
618//_____________________________________________________________________________
cb2f9e9b 619Double_t AliTRDsimTR::GetMuCO(Double_t energyMeV)
46d29e70 620{
621 //
622 // Returns the photon absorbtion cross section for CO2
623 //
624
625 const Int_t kN = 36;
626
627 Double_t mu[kN] = { 0.39383E+04, 0.13166E+04, 0.58750E+03
628 , 0.18240E+03, 0.77996E+02, 0.40024E+02
629 , 0.23116E+02, 0.96997E+01, 0.49726E+01
630 , 0.15543E+01, 0.74915E+00, 0.34442E+00
631 , 0.24440E+00, 0.20589E+00, 0.18632E+00
632 , 0.16578E+00, 0.15394E+00, 0.13558E+00
633 , 0.12336E+00, 0.10678E+00, 0.95510E-01
634 , 0.87165E-01, 0.80587E-01, 0.70769E-01
635 , 0.63626E-01, 0.56894E-01, 0.51782E-01
636 , 0.44499E-01, 0.35839E-01, 0.30825E-01
637 , 0.27555E-01, 0.25269E-01, 0.22311E-01
638 , 0.20516E-01, 0.18184E-01, 0.17152E-01 };
639
640 Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
641 , 0.30000E-02, 0.40000E-02, 0.50000E-02
642 , 0.60000E-02, 0.80000E-02, 0.10000E-01
643 , 0.15000E-01, 0.20000E-01, 0.30000E-01
644 , 0.40000E-01, 0.50000E-01, 0.60000E-01
645 , 0.80000E-01, 0.10000E+00, 0.15000E+00
646 , 0.20000E+00, 0.30000E+00, 0.40000E+00
647 , 0.50000E+00, 0.60000E+00, 0.80000E+00
648 , 0.10000E+01, 0.12500E+01, 0.15000E+01
649 , 0.20000E+01, 0.30000E+01, 0.40000E+01
650 , 0.50000E+01, 0.60000E+01, 0.80000E+01
651 , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
652
653 return Interpolate(energyMeV,en,mu,kN);
654
655}
656
657//_____________________________________________________________________________
cb2f9e9b 658Double_t AliTRDsimTR::GetMuXe(Double_t energyMeV)
46d29e70 659{
660 //
661 // Returns the photon absorbtion cross section for xenon
662 //
663
664 const Int_t kN = 48;
665
666 Double_t mu[kN] = { 9.413E+03, 8.151E+03, 7.035E+03
667 , 7.338E+03, 4.085E+03, 2.088E+03
668 , 7.780E+02, 3.787E+02, 2.408E+02
669 , 6.941E+02, 6.392E+02, 6.044E+02
670 , 8.181E+02, 7.579E+02, 6.991E+02
671 , 8.064E+02, 6.376E+02, 3.032E+02
672 , 1.690E+02, 5.743E+01, 2.652E+01
673 , 8.930E+00, 6.129E+00, 3.316E+01
674 , 2.270E+01, 1.272E+01, 7.825E+00
675 , 3.633E+00, 2.011E+00, 7.202E-01
676 , 3.760E-01, 1.797E-01, 1.223E-01
677 , 9.699E-02, 8.281E-02, 6.696E-02
678 , 5.785E-02, 5.054E-02, 4.594E-02
679 , 4.078E-02, 3.681E-02, 3.577E-02
680 , 3.583E-02, 3.634E-02, 3.797E-02
681 , 3.987E-02, 4.445E-02, 4.815E-02 };
682
683 Double_t en[kN] = { 1.00000E-03, 1.07191E-03, 1.14900E-03
684 , 1.14900E-03, 1.50000E-03, 2.00000E-03
685 , 3.00000E-03, 4.00000E-03, 4.78220E-03
686 , 4.78220E-03, 5.00000E-03, 5.10370E-03
687 , 5.10370E-03, 5.27536E-03, 5.45280E-03
688 , 5.45280E-03, 6.00000E-03, 8.00000E-03
689 , 1.00000E-02, 1.50000E-02, 2.00000E-02
690 , 3.00000E-02, 3.45614E-02, 3.45614E-02
691 , 4.00000E-02, 5.00000E-02, 6.00000E-02
692 , 8.00000E-02, 1.00000E-01, 1.50000E-01
693 , 2.00000E-01, 3.00000E-01, 4.00000E-01
694 , 5.00000E-01, 6.00000E-01, 8.00000E-01
695 , 1.00000E+00, 1.25000E+00, 1.50000E+00
696 , 2.00000E+00, 3.00000E+00, 4.00000E+00
697 , 5.00000E+00, 6.00000E+00, 8.00000E+00
698 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
699
700 return Interpolate(energyMeV,en,mu,kN);
701
702}
703
704//_____________________________________________________________________________
cb2f9e9b 705Double_t AliTRDsimTR::GetMuBu(Double_t energyMeV)
46d29e70 706{
707 //
708 // Returns the photon absorbtion cross section for isobutane
709 //
710
711 const Int_t kN = 36;
712
713 Double_t mu[kN] = { 0.38846E+03, 0.12291E+03, 0.53225E+02
714 , 0.16091E+02, 0.69114E+01, 0.36541E+01
715 , 0.22282E+01, 0.11149E+01, 0.72887E+00
716 , 0.45053E+00, 0.38167E+00, 0.33920E+00
717 , 0.32155E+00, 0.30949E+00, 0.29960E+00
718 , 0.28317E+00, 0.26937E+00, 0.24228E+00
719 , 0.22190E+00, 0.19289E+00, 0.17288E+00
720 , 0.15789E+00, 0.14602E+00, 0.12829E+00
721 , 0.11533E+00, 0.10310E+00, 0.93790E-01
722 , 0.80117E-01, 0.63330E-01, 0.53229E-01
723 , 0.46390E-01, 0.41425E-01, 0.34668E-01
724 , 0.30267E-01, 0.23910E-01, 0.20509E-01 };
725
726 Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02
727 , 0.30000E-02, 0.40000E-02, 0.50000E-02
728 , 0.60000E-02, 0.80000E-02, 0.10000E-01
729 , 0.15000E-01, 0.20000E-01, 0.30000E-01
730 , 0.40000E-01, 0.50000E-01, 0.60000E-01
731 , 0.80000E-01, 0.10000E+00, 0.15000E+00
732 , 0.20000E+00, 0.30000E+00, 0.40000E+00
733 , 0.50000E+00, 0.60000E+00, 0.80000E+00
734 , 0.10000E+01, 0.12500E+01, 0.15000E+01
735 , 0.20000E+01, 0.30000E+01, 0.40000E+01
736 , 0.50000E+01, 0.60000E+01, 0.80000E+01
737 , 0.10000E+02, 0.15000E+02, 0.20000E+02 };
738
739 return Interpolate(energyMeV,en,mu,kN);
740
741}
742
743//_____________________________________________________________________________
cb2f9e9b 744Double_t AliTRDsimTR::GetMuMy(Double_t energyMeV)
46d29e70 745{
746 //
747 // Returns the photon absorbtion cross section for mylar
748 //
749
750 const Int_t kN = 36;
751
752 Double_t mu[kN] = { 2.911E+03, 9.536E+02, 4.206E+02
753 , 1.288E+02, 5.466E+01, 2.792E+01
754 , 1.608E+01, 6.750E+00, 3.481E+00
755 , 1.132E+00, 5.798E-01, 3.009E-01
756 , 2.304E-01, 2.020E-01, 1.868E-01
757 , 1.695E-01, 1.586E-01, 1.406E-01
758 , 1.282E-01, 1.111E-01, 9.947E-02
759 , 9.079E-02, 8.395E-02, 7.372E-02
760 , 6.628E-02, 5.927E-02, 5.395E-02
761 , 4.630E-02, 3.715E-02, 3.181E-02
762 , 2.829E-02, 2.582E-02, 2.257E-02
763 , 2.057E-02, 1.789E-02, 1.664E-02 };
764
765 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
766 , 3.00000E-03, 4.00000E-03, 5.00000E-03
767 , 6.00000E-03, 8.00000E-03, 1.00000E-02
768 , 1.50000E-02, 2.00000E-02, 3.00000E-02
769 , 4.00000E-02, 5.00000E-02, 6.00000E-02
770 , 8.00000E-02, 1.00000E-01, 1.50000E-01
771 , 2.00000E-01, 3.00000E-01, 4.00000E-01
772 , 5.00000E-01, 6.00000E-01, 8.00000E-01
773 , 1.00000E+00, 1.25000E+00, 1.50000E+00
774 , 2.00000E+00, 3.00000E+00, 4.00000E+00
775 , 5.00000E+00, 6.00000E+00, 8.00000E+00
776 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
777
778 return Interpolate(energyMeV,en,mu,kN);
779
780}
781
782//_____________________________________________________________________________
cb2f9e9b 783Double_t AliTRDsimTR::GetMuN2(Double_t energyMeV)
46d29e70 784{
785 //
786 // Returns the photon absorbtion cross section for nitrogen
787 //
788
789 const Int_t kN = 36;
790
791 Double_t mu[kN] = { 3.311E+03, 1.083E+03, 4.769E+02
792 , 1.456E+02, 6.166E+01, 3.144E+01
793 , 1.809E+01, 7.562E+00, 3.879E+00
794 , 1.236E+00, 6.178E-01, 3.066E-01
795 , 2.288E-01, 1.980E-01, 1.817E-01
796 , 1.639E-01, 1.529E-01, 1.353E-01
797 , 1.233E-01, 1.068E-01, 9.557E-02
798 , 8.719E-02, 8.063E-02, 7.081E-02
799 , 6.364E-02, 5.693E-02, 5.180E-02
800 , 4.450E-02, 3.579E-02, 3.073E-02
801 , 2.742E-02, 2.511E-02, 2.209E-02
802 , 2.024E-02, 1.782E-02, 1.673E-02 };
803
804 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
805 , 3.00000E-03, 4.00000E-03, 5.00000E-03
806 , 6.00000E-03, 8.00000E-03, 1.00000E-02
807 , 1.50000E-02, 2.00000E-02, 3.00000E-02
808 , 4.00000E-02, 5.00000E-02, 6.00000E-02
809 , 8.00000E-02, 1.00000E-01, 1.50000E-01
810 , 2.00000E-01, 3.00000E-01, 4.00000E-01
811 , 5.00000E-01, 6.00000E-01, 8.00000E-01
812 , 1.00000E+00, 1.25000E+00, 1.50000E+00
813 , 2.00000E+00, 3.00000E+00, 4.00000E+00
814 , 5.00000E+00, 6.00000E+00, 8.00000E+00
815 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
816
817 return Interpolate(energyMeV,en,mu,kN);
818
819}
820
821//_____________________________________________________________________________
cb2f9e9b 822Double_t AliTRDsimTR::GetMuO2(Double_t energyMeV)
46d29e70 823{
824 //
825 // Returns the photon absorbtion cross section for oxygen
826 //
827
828 const Int_t kN = 36;
829
830 Double_t mu[kN] = { 4.590E+03, 1.549E+03, 6.949E+02
831 , 2.171E+02, 9.315E+01, 4.790E+01
832 , 2.770E+01, 1.163E+01, 5.952E+00
833 , 1.836E+00, 8.651E-01, 3.779E-01
834 , 2.585E-01, 2.132E-01, 1.907E-01
835 , 1.678E-01, 1.551E-01, 1.361E-01
836 , 1.237E-01, 1.070E-01, 9.566E-02
837 , 8.729E-02, 8.070E-02, 7.087E-02
838 , 6.372E-02, 5.697E-02, 5.185E-02
839 , 4.459E-02, 3.597E-02, 3.100E-02
840 , 2.777E-02, 2.552E-02, 2.263E-02
841 , 2.089E-02, 1.866E-02, 1.770E-02 };
842
843 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
844 , 3.00000E-03, 4.00000E-03, 5.00000E-03
845 , 6.00000E-03, 8.00000E-03, 1.00000E-02
846 , 1.50000E-02, 2.00000E-02, 3.00000E-02
847 , 4.00000E-02, 5.00000E-02, 6.00000E-02
848 , 8.00000E-02, 1.00000E-01, 1.50000E-01
849 , 2.00000E-01, 3.00000E-01, 4.00000E-01
850 , 5.00000E-01, 6.00000E-01, 8.00000E-01
851 , 1.00000E+00, 1.25000E+00, 1.50000E+00
852 , 2.00000E+00, 3.00000E+00, 4.00000E+00
853 , 5.00000E+00, 6.00000E+00, 8.00000E+00
854 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
855
856 return Interpolate(energyMeV,en,mu,kN);
857
858}
859
860//_____________________________________________________________________________
cb2f9e9b 861Double_t AliTRDsimTR::GetMuHe(Double_t energyMeV)
46d29e70 862{
863 //
864 // Returns the photon absorbtion cross section for helium
865 //
866
867 const Int_t kN = 36;
868
869 Double_t mu[kN] = { 6.084E+01, 1.676E+01, 6.863E+00
870 , 2.007E+00, 9.329E-01, 5.766E-01
871 , 4.195E-01, 2.933E-01, 2.476E-01
872 , 2.092E-01, 1.960E-01, 1.838E-01
873 , 1.763E-01, 1.703E-01, 1.651E-01
874 , 1.562E-01, 1.486E-01, 1.336E-01
875 , 1.224E-01, 1.064E-01, 9.535E-02
876 , 8.707E-02, 8.054E-02, 7.076E-02
877 , 6.362E-02, 5.688E-02, 5.173E-02
878 , 4.422E-02, 3.503E-02, 2.949E-02
879 , 2.577E-02, 2.307E-02, 1.940E-02
880 , 1.703E-02, 1.363E-02, 1.183E-02 };
881
882 Double_t en[kN] = { 1.00000E-03, 1.50000E-03, 2.00000E-03
883 , 3.00000E-03, 4.00000E-03, 5.00000E-03
884 , 6.00000E-03, 8.00000E-03, 1.00000E-02
885 , 1.50000E-02, 2.00000E-02, 3.00000E-02
886 , 4.00000E-02, 5.00000E-02, 6.00000E-02
887 , 8.00000E-02, 1.00000E-01, 1.50000E-01
888 , 2.00000E-01, 3.00000E-01, 4.00000E-01
889 , 5.00000E-01, 6.00000E-01, 8.00000E-01
890 , 1.00000E+00, 1.25000E+00, 1.50000E+00
891 , 2.00000E+00, 3.00000E+00, 4.00000E+00
892 , 5.00000E+00, 6.00000E+00, 8.00000E+00
893 , 1.00000E+01, 1.50000E+01, 2.00000E+01 };
894
895 return Interpolate(energyMeV,en,mu,kN);
896
897}
898
842287f2 899//_____________________________________________________________________________
cb2f9e9b 900Double_t AliTRDsimTR::GetMuAi(Double_t energyMeV)
842287f2 901{
902 //
903 // Returns the photon absorbtion cross section for air
904 // Implemented by Oliver Busch
905 //
906
907 const Int_t kN = 38;
908
909 Double_t mu[kN] = { 0.35854E+04, 0.11841E+04, 0.52458E+03,
910 0.16143E+03, 0.14250E+03, 0.15722E+03,
911 0.77538E+02, 0.40099E+02, 0.23313E+02,
912 0.98816E+01, 0.51000E+01, 0.16079E+01,
913 0.77536E+00, 0.35282E+00, 0.24790E+00,
914 0.20750E+00, 0.18703E+00, 0.16589E+00,
915 0.15375E+00, 0.13530E+00, 0.12311E+00,
916 0.10654E+00, 0.95297E-01, 0.86939E-01,
917 0.80390E-01, 0.70596E-01, 0.63452E-01,
918 0.56754E-01, 0.51644E-01, 0.44382E-01,
919 0.35733E-01, 0.30721E-01, 0.27450E-01,
920 0.25171E-01, 0.22205E-01, 0.20399E-01,
921 0.18053E-01, 0.18057E-01 };
922
923
924
925 Double_t en[kN] = { 0.10000E-02, 0.15000E-02, 0.20000E-02,
926 0.30000E-02, 0.32029E-02, 0.32029E-02,
927 0.40000E-02, 0.50000E-02, 0.60000E-02,
928 0.80000E-02, 0.10000E-01, 0.15000E-01,
929 0.20000E-01, 0.30000E-01, 0.40000E-01,
930 0.50000E-01, 0.60000E-01, 0.80000E-01,
931 0.10000E+00, 0.15000E+00, 0.20000E+00,
932 0.30000E+00, 0.40000E+00, 0.50000E+00,
933 0.60000E+00, 0.80000E+00, 0.10000E+01,
934 0.12500E+01, 0.15000E+01, 0.20000E+01,
935 0.30000E+01, 0.40000E+01, 0.50000E+01,
936 0.60000E+01, 0.80000E+01, 0.10000E+02,
937 0.15000E+02, 0.20000E+02 };
938
939 return Interpolate(energyMeV,en,mu,kN);
940
941}
942
46d29e70 943//_____________________________________________________________________________
cb2f9e9b 944Double_t AliTRDsimTR::Interpolate(Double_t energyMeV
46d29e70 945 , Double_t *en, Double_t *mu, Int_t n)
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
46d29e70 967 , Int_t &kl, Double_t &dx)
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}