]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDv1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Using look-up table to calculate the integrated spectrum (Valgrind)
[u/mrichter/AliRoot.git] / TRD / AliTRDv1.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
88cb7938 16/* $Id$ */
4c039060 17
fe4da5cc 18///////////////////////////////////////////////////////////////////////////////
19// //
769257f4 20// Transition Radiation Detector version 1 -- slow simulator //
fe4da5cc 21// //
22//Begin_Html
23/*
5c7f4665 24<img src="picts/AliTRDfullClass.gif">
fe4da5cc 25*/
26//End_Html
27// //
28// //
29///////////////////////////////////////////////////////////////////////////////
30
769257f4 31#include <stdlib.h>
32
793ff80c 33#include <TF1.h>
1819f4bb 34#include <TLorentzVector.h>
88cb7938 35#include <TMath.h>
36#include <TRandom.h>
37#include <TVector.h>
38#include <TVirtualMC.h>
fe4da5cc 39
d3f347ff 40#include "AliConst.h"
45160b1f 41#include "AliLog.h"
42#include "AliMC.h"
88cb7938 43#include "AliRun.h"
44#include "AliTRDgeometry.h"
793ff80c 45#include "AliTRDhit.h"
793ff80c 46#include "AliTRDsim.h"
88cb7938 47#include "AliTRDv1.h"
851d3db9 48
fe4da5cc 49ClassImp(AliTRDv1)
8230f242 50
51//_____________________________________________________________________________
52AliTRDv1::AliTRDv1():AliTRD()
53{
54 //
55 // Default constructor
56 //
57
a328fff9 58 fSensSelect = 0;
59 fSensPlane = -1;
60 fSensChamber = -1;
61 fSensSector = -1;
62 fSensSectorRange = 0;
8230f242 63
a328fff9 64 fDeltaE = NULL;
65 fDeltaG = NULL;
66 fTR = NULL;
67
68 fStepSize = 0.1;
c4214bc0 69 fTypeOfStepManager = 1;
8230f242 70
71}
72
fe4da5cc 73//_____________________________________________________________________________
74AliTRDv1::AliTRDv1(const char *name, const char *title)
75 :AliTRD(name, title)
76{
77 //
851d3db9 78 // Standard constructor for Transition Radiation Detector version 1
fe4da5cc 79 //
82bbf98a 80
a328fff9 81 fSensSelect = 0;
82 fSensPlane = -1;
83 fSensChamber = -1;
84 fSensSector = -1;
85 fSensSectorRange = 0;
5c7f4665 86
a328fff9 87 fDeltaE = NULL;
88 fDeltaG = NULL;
89 fTR = NULL;
90 fStepSize = 0.1;
c4214bc0 91 fTypeOfStepManager = 1;
5c7f4665 92
93 SetBufferSize(128000);
94
95}
96
8230f242 97//_____________________________________________________________________________
73ae7b59 98AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
8230f242 99{
100 //
101 // Copy constructor
102 //
103
dd9a6ee3 104 ((AliTRDv1 &) trd).Copy(*this);
8230f242 105
106}
107
5c7f4665 108//_____________________________________________________________________________
109AliTRDv1::~AliTRDv1()
110{
dd9a6ee3 111 //
112 // AliTRDv1 destructor
113 //
82bbf98a 114
5c7f4665 115 if (fDeltaE) delete fDeltaE;
a328fff9 116 if (fDeltaG) delete fDeltaG;
793ff80c 117 if (fTR) delete fTR;
82bbf98a 118
fe4da5cc 119}
120
dd9a6ee3 121//_____________________________________________________________________________
122AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
123{
124 //
125 // Assignment operator
126 //
127
128 if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
129 return *this;
130
131}
8230f242 132
133//_____________________________________________________________________________
e0d47c25 134void AliTRDv1::Copy(TObject &trd) const
8230f242 135{
c4214bc0 136 printf("void AliTRDv1::Copy(TObject &trd) const\n");
8230f242 137 //
138 // Copy function
139 //
140
a328fff9 141 ((AliTRDv1 &) trd).fSensSelect = fSensSelect;
142 ((AliTRDv1 &) trd).fSensPlane = fSensPlane;
143 ((AliTRDv1 &) trd).fSensChamber = fSensChamber;
144 ((AliTRDv1 &) trd).fSensSector = fSensSector;
145 ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
146
147 ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
148 ((AliTRDv1 &) trd).fStepSize = fStepSize;
8230f242 149
793ff80c 150 fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
a328fff9 151 fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
793ff80c 152 fTR->Copy(*((AliTRDv1 &) trd).fTR);
8230f242 153
154}
155
fe4da5cc 156//_____________________________________________________________________________
157void AliTRDv1::CreateGeometry()
158{
159 //
851d3db9 160 // Create the GEANT geometry for the Transition Radiation Detector - Version 1
5c7f4665 161 // This version covers the full azimuth.
d3f347ff 162 //
163
82bbf98a 164 // Check that FRAME is there otherwise we have no place where to put the TRD
8230f242 165 AliModule* frame = gAlice->GetModule("FRAME");
166 if (!frame) return;
d3f347ff 167
82bbf98a 168 // Define the chambers
169 AliTRD::CreateGeometry();
d3f347ff 170
fe4da5cc 171}
172
173//_____________________________________________________________________________
174void AliTRDv1::CreateMaterials()
175{
176 //
851d3db9 177 // Create materials for the Transition Radiation Detector version 1
fe4da5cc 178 //
82bbf98a 179
d3f347ff 180 AliTRD::CreateMaterials();
82bbf98a 181
fe4da5cc 182}
183
793ff80c 184//_____________________________________________________________________________
185void AliTRDv1::CreateTRhit(Int_t det)
186{
187 //
188 // Creates an electron cluster from a TR photon.
189 // The photon is assumed to be created a the end of the radiator. The
190 // distance after which it deposits its energy takes into account the
191 // absorbtion of the entrance window and of the gas mixture in drift
192 // volume.
193 //
194
195 // PDG code electron
196 const Int_t kPdgElectron = 11;
197
198 // Ionization energy
bc327ce2 199 const Float_t kWion = 23.53;
793ff80c 200
201 // Maximum number of TR photons per track
202 const Int_t kNTR = 50;
203
204 TLorentzVector mom, pos;
793ff80c 205
793ff80c 206 // Create TR at the entrance of the chamber
207 if (gMC->IsTrackEntering()) {
208
f73816f5 209 // Create TR only for electrons
210 Int_t iPdg = gMC->TrackPid();
211 if (TMath::Abs(iPdg) != kPdgElectron) return;
212
793ff80c 213 Float_t eTR[kNTR];
214 Int_t nTR;
215
216 // Create TR photons
217 gMC->TrackMomentum(mom);
218 Float_t pTot = mom.Rho();
219 fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
220 if (nTR > kNTR) {
45160b1f 221 AliFatal(Form("Boundary error: nTR = %d, kNTR = %d",nTR,kNTR));
793ff80c 222 }
223
224 // Loop through the TR photons
225 for (Int_t iTR = 0; iTR < nTR; iTR++) {
226
227 Float_t energyMeV = eTR[iTR] * 0.001;
228 Float_t energyeV = eTR[iTR] * 1000.0;
229 Float_t absLength = 0;
230 Float_t sigma = 0;
231
232 // Take the absorbtion in the entrance window into account
233 Double_t muMy = fTR->GetMuMy(energyMeV);
234 sigma = muMy * fFoilDensity;
842287f2 235 if (sigma > 0.0) {
236 absLength = gRandom->Exp(1.0/sigma);
237 if (absLength < AliTRDgeometry::MyThick()) continue;
238 }
239 else {
240 continue;
241 }
793ff80c 242
243 // The absorbtion cross sections in the drift gas
3dac2b2d 244 // Gas-mixture (Xe/CO2)
245 Double_t muXe = fTR->GetMuXe(energyMeV);
246 Double_t muCO = fTR->GetMuCO(energyMeV);
247 sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
793ff80c 248
249 // The distance after which the energy of the TR photon
250 // is deposited.
842287f2 251 if (sigma > 0.0) {
252 absLength = gRandom->Exp(1.0/sigma);
a328fff9 253 if (absLength > (AliTRDgeometry::DrThick()
254 + AliTRDgeometry::AmThick())) {
255 continue;
256 }
842287f2 257 }
258 else {
259 continue;
260 }
793ff80c 261
262 // The position of the absorbtion
263 Float_t posHit[3];
264 gMC->TrackPosition(pos);
265 posHit[0] = pos[0] + mom[0] / pTot * absLength;
266 posHit[1] = pos[1] + mom[1] / pTot * absLength;
c4214bc0 267 posHit[2] = pos[2] + mom[2] / pTot * absLength;
793ff80c 268
269 // Create the charge
270 Int_t q = ((Int_t) (energyeV / kWion));
271
272 // Add the hit to the array. TR photon hits are marked
273 // by negative charge
c4214bc0 274 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
793ff80c 275
276 }
277
278 }
279
280}
281
5c7f4665 282//_____________________________________________________________________________
283void AliTRDv1::Init()
284{
285 //
286 // Initialise Transition Radiation Detector after geometry has been built.
5c7f4665 287 //
288
289 AliTRD::Init();
290
45160b1f 291 AliDebug(1,"Slow simulator\n");
851d3db9 292 if (fSensSelect) {
293 if (fSensPlane >= 0)
45160b1f 294 AliInfo(Form("Only plane %d is sensitive"));
851d3db9 295 if (fSensChamber >= 0)
45160b1f 296 AliInfo(Form("Only chamber %d is sensitive",fSensChamber));
9d0b222b 297 if (fSensSector >= 0) {
298 Int_t sens1 = fSensSector;
299 Int_t sens2 = fSensSector + fSensSectorRange;
793ff80c 300 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
301 * AliTRDgeometry::Nsect();
45160b1f 302 AliInfo(Form("Only sectors %d - %d are sensitive\n",sens1,sens2-1));
9d0b222b 303 }
851d3db9 304 }
793ff80c 305 if (fTR)
45160b1f 306 AliInfo("TR simulation on")
793ff80c 307 else
45160b1f 308 AliInfo("TR simulation off");
5c7f4665 309
310 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
311 const Float_t kPoti = 12.1;
312 // Maximum energy (50 keV);
313 const Float_t kEend = 50000.0;
314 // Ermilova distribution for the delta-ray spectrum
8230f242 315 Float_t poti = TMath::Log(kPoti);
316 Float_t eEnd = TMath::Log(kEend);
a328fff9 317
318 // Ermilova distribution for the delta-ray spectrum
c4214bc0 319 fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
a328fff9 320
321 // Geant3 distribution for the delta-ray spectrum
c4214bc0 322 fDeltaG = new TF1("deltag",IntSpecGeant,2.421257,28.536469,0);
5c7f4665 323
45160b1f 324 AliDebug(1,"+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
5c7f4665 325
fe4da5cc 326}
327
793ff80c 328//_____________________________________________________________________________
329AliTRDsim *AliTRDv1::CreateTR()
330{
331 //
332 // Enables the simulation of TR
333 //
334
335 fTR = new AliTRDsim();
336 return fTR;
337
338}
339
5c7f4665 340//_____________________________________________________________________________
341void AliTRDv1::SetSensPlane(Int_t iplane)
342{
343 //
851d3db9 344 // Defines the hit-sensitive plane (0-5)
5c7f4665 345 //
82bbf98a 346
851d3db9 347 if ((iplane < 0) || (iplane > 5)) {
45160b1f 348 AliWarning(Form("Wrong input value:%d",iplane));
349 AliWarning("Use standard setting");
851d3db9 350 fSensPlane = -1;
351 fSensSelect = 0;
5c7f4665 352 return;
353 }
82bbf98a 354
5c7f4665 355 fSensSelect = 1;
356 fSensPlane = iplane;
82bbf98a 357
5c7f4665 358}
359
360//_____________________________________________________________________________
361void AliTRDv1::SetSensChamber(Int_t ichamber)
362{
363 //
851d3db9 364 // Defines the hit-sensitive chamber (0-4)
5c7f4665 365 //
366
851d3db9 367 if ((ichamber < 0) || (ichamber > 4)) {
45160b1f 368 AliWarning(Form("Wrong input value: %d",ichamber));
369 AliWarning("Use standard setting");
851d3db9 370 fSensChamber = -1;
371 fSensSelect = 0;
5c7f4665 372 return;
373 }
374
375 fSensSelect = 1;
376 fSensChamber = ichamber;
377
378}
379
380//_____________________________________________________________________________
381void AliTRDv1::SetSensSector(Int_t isector)
382{
383 //
851d3db9 384 // Defines the hit-sensitive sector (0-17)
5c7f4665 385 //
386
9d0b222b 387 SetSensSector(isector,1);
388
389}
390
391//_____________________________________________________________________________
392void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
393{
394 //
395 // Defines a range of hit-sensitive sectors. The range is defined by
396 // <isector> (0-17) as the starting point and <nsector> as the number
397 // of sectors to be included.
398 //
399
851d3db9 400 if ((isector < 0) || (isector > 17)) {
45160b1f 401 AliWarning(Form("Wrong input value <isector>: %d",isector));
402 AliWarning("Use standard setting");
9d0b222b 403 fSensSector = -1;
404 fSensSectorRange = 0;
405 fSensSelect = 0;
5c7f4665 406 return;
407 }
408
9d0b222b 409 if ((nsector < 1) || (nsector > 18)) {
45160b1f 410 AliWarning(Form("Wrong input value <nsector>: %d",nsector));
411 AliWarning("Use standard setting");
9d0b222b 412 fSensSector = -1;
413 fSensSectorRange = 0;
414 fSensSelect = 0;
415 return;
416 }
417
418 fSensSelect = 1;
419 fSensSector = isector;
420 fSensSectorRange = nsector;
5c7f4665 421
422}
423
424//_____________________________________________________________________________
425void AliTRDv1::StepManager()
a328fff9 426{
427 //
c4214bc0 428 // Slow simulator. Every charged track produces electron cluster as hits
a328fff9 429 // along its path across the drift volume.
430 //
431
432 switch (fTypeOfStepManager) {
433 case 0 : StepManagerErmilova(); break; // 0 is Ermilova
434 case 1 : StepManagerGeant(); break; // 1 is Geant
435 case 2 : StepManagerFixedStep(); break; // 2 is fixed step
45160b1f 436 default : AliWarning("Not a valid Step Manager.");
a328fff9 437 }
438
439}
440
441//_____________________________________________________________________________
442void AliTRDv1::SelectStepManager(Int_t t)
443{
444 //
445 // Selects a step manager type:
446 // 0 - Ermilova
447 // 1 - Geant3
448 // 2 - Fixed step size
449 //
450
c4214bc0 451/* if (t == 1) {
45160b1f 452 AliWarning("Sorry, Geant parametrization step manager is not implemented yet. Please ask K.Oyama for detail.");
a328fff9 453 }
c4214bc0 454*/
a328fff9 455 fTypeOfStepManager = t;
45160b1f 456 AliInfo(Form("Step Manager type %d was selected",fTypeOfStepManager));
a328fff9 457
458}
459
460//_____________________________________________________________________________
461void AliTRDv1::StepManagerGeant()
462{
463 //
c4214bc0 464 // Slow simulator. Every charged track produces electron cluster as hits
a328fff9 465 // along its path across the drift volume. The step size is set acording
466 // to Bethe-Bloch. The energy distribution of the delta electrons follows
467 // a spectrum taken from Geant3.
468 //
c4214bc0 469 Int_t pla = 0;
470 Int_t cha = 0;
471 Int_t sec = 0;
472 Int_t det = 0;
473 Int_t iPdg;
474 Int_t qTot;
475
476 Float_t hits[3];
477 Float_t charge;
478 Float_t aMass;
479
480 Double_t pTot = 0;
481 Double_t eDelta;
482 Double_t betaGamma, pp;
483 Double_t stepSize=0;
484
485 Bool_t drRegion = kFALSE;
486 Bool_t amRegion = kFALSE;
487
488 TString cIdCurrent;
489 TString cIdSensDr = "J";
490 TString cIdSensAm = "K";
491 Char_t cIdChamber[3];
492 cIdChamber[2] = 0;
493
494 TLorentzVector pos, mom;
495
496 const Int_t kNplan = AliTRDgeometry::Nplan();
497 const Int_t kNcham = AliTRDgeometry::Ncham();
498 const Int_t kNdetsec = kNplan * kNcham;
499
500 const Double_t kBig = 1.0E+12; // Infinitely big
bc327ce2 501 const Float_t kWion = 23.53; // Ionization energy
c4214bc0 502 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
503
504 // Minimum energy for the step size adjustment
505 const Float_t kEkinMinStep = 1.0e-5;
506 // energy threshold for production of delta electrons
507 const Float_t kECut = 1.0e4;
508 // Parameters entering the parametrized range for delta electrons
509 const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
510 // Gas density -> To be made user adjustable !
511 const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
a328fff9 512
c4214bc0 513 // Plateau value of the energy-loss for electron in xenon
514 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
515 //const Double_t kPlateau = 1.70;
516 // the averaged value (26/3/99)
517 const Float_t kPlateau = 1.55;
518
519 const Float_t kPrim = 19.34; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
520 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
521 const Float_t kPoti = 12.1;
522
523 const Int_t kPdgElectron = 11; // PDG code electron
524
525 // Set the maximum step size to a very large number for all
526 // neutral particles and those outside the driftvolume
527 gMC->SetMaxStep(kBig);
528
529 // Use only charged tracks
530 if (( gMC->TrackCharge() ) &&
531 (!gMC->IsTrackStop() ) &&
532 (!gMC->IsTrackDisappeared())) {
533
534 // Inside a sensitive volume?
535 drRegion = kFALSE;
536 amRegion = kFALSE;
537 cIdCurrent = gMC->CurrentVolName();
538 if (cIdSensDr == cIdCurrent[1]) {
539 drRegion = kTRUE;
540 }
541 if (cIdSensAm == cIdCurrent[1]) {
542 amRegion = kTRUE;
543 }
544 if (drRegion || amRegion) {
a328fff9 545
c4214bc0 546 // The hit coordinates and charge
547 gMC->TrackPosition(pos);
548 hits[0] = pos[0];
549 hits[1] = pos[1];
550 hits[2] = pos[2];
551
552 // The sector number (0 - 17)
553 // The numbering goes clockwise and starts at y = 0
554 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
555 if (phi < 90.)
556 phi = phi + 270.;
557 else
558 phi = phi - 90.;
559 sec = ((Int_t) (phi / 20));
560
561 // The plane and chamber number
562 cIdChamber[0] = cIdCurrent[2];
563 cIdChamber[1] = cIdCurrent[3];
564 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
565 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
566 pla = ((Int_t) idChamber % kNplan);
567
568 // Check on selected volumes
569 Int_t addthishit = 1;
570 if (fSensSelect) {
571 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
572 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
573 if (fSensSector >= 0) {
574 Int_t sens1 = fSensSector;
575 Int_t sens2 = fSensSector + fSensSectorRange;
576 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
577 * AliTRDgeometry::Nsect();
578 if (sens1 < sens2) {
579 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
580 }
581 else {
582 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
583 }
584 }
585 }
586
587 // Add this hit
588 if (addthishit) {
589
590 // The detector number
591 det = fGeometry->GetDetector(pla,cha,sec);
592
593 // Special hits only in the drift region
594 if (drRegion) {
595 // Create a track reference at the entrance and
596 // exit of each chamber that contain the
597 // momentum components of the particle
598 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
599 gMC->TrackMomentum(mom);
600 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
601 }
602
603 if (gMC->IsTrackEntering() && !gMC->IsNewTrack()) {
604 // determine if hit belong to primary track
605 fPrimaryTrackPid=gAlice->GetMCApp()->GetCurrentTrackNumber();
606 //determine track length when entering the detector
607 fTrackLength0=gMC->TrackLength();
608 }
609
610 // Create the hits from TR photons
611 if (fTR) CreateTRhit(det);
612 }
613
614 // Calculate the energy of the delta-electrons
615 // modified by Alex Bercuci (A.Bercuci@gsi.de) on 26.01.06
616 // take into account correlation with the underlying GEANT tracking
617 // mechanism. see
618 // http://www-linux.gsi.de/~abercuci/Contributions/TRD/index.html
619
620 // determine the most significant process (last on the processes list)
621 // which caused this hit
622 TArrayI processes;
623 gMC->StepProcesses(processes);
624 int nofprocesses=processes.GetSize(), pid;
625 if(!nofprocesses) pid=0;
626 else pid= processes[nofprocesses-1];
627
628 // generate Edep according to GEANT parametrisation
629 eDelta =TMath::Exp(fDeltaG->GetRandom()) - kPoti;
630 eDelta=TMath::Max(eDelta,0.0);
631 float pr_range=0.;
632 float range=gMC->TrackLength()-fTrackLength0;
633 // merge GEANT tracker information with localy cooked one
634 if(gAlice->GetMCApp()->GetCurrentTrackNumber()==fPrimaryTrackPid) {
635// printf("primary pid=%d eDelta=%f\n",pid,eDelta);
636 if(pid==27){
637 if(eDelta>=kECut){
638 pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
639 if(pr_range>=(3.7-range)) eDelta*=.1;
640 }
641 } else if(pid==1){
642 if(eDelta<kECut) eDelta*=.5;
643 else {
644 pr_range=ra*eDelta*.001*(1.-rb/(1.+rc*eDelta*0.001))/rho;
645 if(pr_range>=((AliTRDgeometry::DrThick()
646 + AliTRDgeometry::AmThick())-range)) eDelta*=.05;
647 else eDelta*=.5;
648 }
649 } else eDelta=0.;
650 } else eDelta=0.;
651
652 // Generate the electron cluster size
653 if(eDelta==0.) qTot=0;
654 else qTot = ((Int_t) (eDelta / kWion) + 1);
655 // Create a new dEdx hit
656 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot, drRegion);
657
658 // Calculate the maximum step size for the next tracking step
659 // Produce only one hit if Ekin is below cutoff
660 aMass = gMC->TrackMass();
661 if ((gMC->Etot() - aMass) > kEkinMinStep) {
662
663 // The energy loss according to Bethe Bloch
664 iPdg = TMath::Abs(gMC->TrackPid());
665 if ( (iPdg != kPdgElectron) ||
666 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
667 gMC->TrackMomentum(mom);
668 pTot = mom.Rho();
669 betaGamma = pTot / aMass;
670 pp = BetheBlochGeant(betaGamma);
671 // Take charge > 1 into account
672 charge = gMC->TrackCharge();
673 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
674 } else { // Electrons above 20 Mev/c are at the plateau
675 pp = kPrim * kPlateau;
676 }
677
1b775a44 678 Int_t nsteps = 0;
679 do {nsteps = gRandom->Poisson(pp);} while(!nsteps);
680 stepSize = 1./nsteps;
681 gMC->SetMaxStep(stepSize);
682 }
c4214bc0 683 }
684 }
685 }
a328fff9 686}
687
688//_____________________________________________________________________________
689void AliTRDv1::StepManagerErmilova()
5c7f4665 690{
691 //
5c7f4665 692 // Slow simulator. Every charged track produces electron cluster as hits
693 // along its path across the drift volume. The step size is set acording
694 // to Bethe-Bloch. The energy distribution of the delta electrons follows
695 // a spectrum taken from Ermilova et al.
696 //
697
851d3db9 698 Int_t pla = 0;
699 Int_t cha = 0;
700 Int_t sec = 0;
793ff80c 701 Int_t det = 0;
851d3db9 702 Int_t iPdg;
793ff80c 703 Int_t qTot;
5c7f4665 704
793ff80c 705 Float_t hits[3];
a5cadd36 706 Double_t random[1];
5c7f4665 707 Float_t charge;
708 Float_t aMass;
709
f73816f5 710 Double_t pTot = 0;
5c7f4665 711 Double_t eDelta;
712 Double_t betaGamma, pp;
f73816f5 713 Double_t stepSize;
5c7f4665 714
332e9569 715 Bool_t drRegion = kFALSE;
716 Bool_t amRegion = kFALSE;
717
718 TString cIdCurrent;
719 TString cIdSensDr = "J";
720 TString cIdSensAm = "K";
593a9fc3 721 Char_t cIdChamber[3];
722 cIdChamber[2] = 0;
332e9569 723
5c7f4665 724 TLorentzVector pos, mom;
82bbf98a 725
332e9569 726 const Int_t kNplan = AliTRDgeometry::Nplan();
e644678a 727 const Int_t kNcham = AliTRDgeometry::Ncham();
728 const Int_t kNdetsec = kNplan * kNcham;
729
a328fff9 730 const Double_t kBig = 1.0E+12; // Infinitely big
bc327ce2 731 const Float_t kWion = 23.53; // Ionization energy
a328fff9 732 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
5c7f4665 733
c4214bc0 734 // energy threshold for production of delta electrons
735 //const Float_t kECut = 1.0e4;
736 // Parameters entering the parametrized range for delta electrons
737 //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
738 // Gas density -> To be made user adjustable !
739 //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
740
f73816f5 741 // Minimum energy for the step size adjustment
742 const Float_t kEkinMinStep = 1.0e-5;
a328fff9 743
5c7f4665 744 // Plateau value of the energy-loss for electron in xenon
745 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
746 //const Double_t kPlateau = 1.70;
747 // the averaged value (26/3/99)
a3c76cdc 748 const Float_t kPlateau = 1.55;
a328fff9 749
750 const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
5c7f4665 751 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
a3c76cdc 752 const Float_t kPoti = 12.1;
851d3db9 753
a328fff9 754 const Int_t kPdgElectron = 11; // PDG code electron
5c7f4665 755
756 // Set the maximum step size to a very large number for all
757 // neutral particles and those outside the driftvolume
758 gMC->SetMaxStep(kBig);
759
760 // Use only charged tracks
761 if (( gMC->TrackCharge() ) &&
762 (!gMC->IsTrackStop() ) &&
763 (!gMC->IsTrackDisappeared())) {
fe4da5cc 764
5c7f4665 765 // Inside a sensitive volume?
332e9569 766 drRegion = kFALSE;
767 amRegion = kFALSE;
768 cIdCurrent = gMC->CurrentVolName();
e6674585 769 if (cIdSensDr == cIdCurrent[1]) {
332e9569 770 drRegion = kTRUE;
771 }
e6674585 772 if (cIdSensAm == cIdCurrent[1]) {
332e9569 773 amRegion = kTRUE;
774 }
775 if (drRegion || amRegion) {
fe4da5cc 776
5c7f4665 777 // The hit coordinates and charge
778 gMC->TrackPosition(pos);
779 hits[0] = pos[0];
780 hits[1] = pos[1];
781 hits[2] = pos[2];
5c7f4665 782
851d3db9 783 // The sector number (0 - 17)
784 // The numbering goes clockwise and starts at y = 0
e15eb584 785 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
851d3db9 786 if (phi < 90.)
787 phi = phi + 270.;
788 else
789 phi = phi - 90.;
790 sec = ((Int_t) (phi / 20));
82bbf98a 791
332e9569 792 // The plane and chamber number
793 cIdChamber[0] = cIdCurrent[2];
794 cIdChamber[1] = cIdCurrent[3];
e644678a 795 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
a5cadd36 796 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
332e9569 797 pla = ((Int_t) idChamber % kNplan);
82bbf98a 798
5c7f4665 799 // Check on selected volumes
800 Int_t addthishit = 1;
801 if (fSensSelect) {
6f1e466d 802 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
803 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
9d0b222b 804 if (fSensSector >= 0) {
805 Int_t sens1 = fSensSector;
806 Int_t sens2 = fSensSector + fSensSectorRange;
793ff80c 807 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
808 * AliTRDgeometry::Nsect();
9d0b222b 809 if (sens1 < sens2) {
810 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
c4214bc0 811 }
9d0b222b 812 else {
813 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
c4214bc0 814 }
815 }
5c7f4665 816 }
817
818 // Add this hit
819 if (addthishit) {
820
f73816f5 821 // The detector number
793ff80c 822 det = fGeometry->GetDetector(pla,cha,sec);
823
a328fff9 824 // Special hits only in the drift region
332e9569 825 if (drRegion) {
f73816f5 826
c61f1a66 827 // Create a track reference at the entrance and
828 // exit of each chamber that contain the
829 // momentum components of the particle
f73816f5 830 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
831 gMC->TrackMomentum(mom);
5d12ce38 832 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
f73816f5 833 }
f73816f5 834 // Create the hits from TR photons
835 if (fTR) CreateTRhit(det);
c4214bc0 836 }
f73816f5 837
f73816f5 838
839 // Calculate the energy of the delta-electrons
840 eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
841 eDelta = TMath::Max(eDelta,0.0);
c4214bc0 842 // Generate the electron cluster size
843 if(eDelta==0.) qTot=0;
844 else qTot = ((Int_t) (eDelta / kWion) + 1);
f73816f5 845
c4214bc0 846 // Create a new dEdx hit
332e9569 847 if (drRegion) {
a328fff9 848 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
c4214bc0 849 ,det,hits,qTot, kTRUE);
850 }
5c7f4665 851 else {
a328fff9 852 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
c4214bc0 853 ,det,hits,qTot,kFALSE);
854 }
f73816f5 855
5c7f4665 856 // Calculate the maximum step size for the next tracking step
f73816f5 857 // Produce only one hit if Ekin is below cutoff
858 aMass = gMC->TrackMass();
859 if ((gMC->Etot() - aMass) > kEkinMinStep) {
860
861 // The energy loss according to Bethe Bloch
862 iPdg = TMath::Abs(gMC->TrackPid());
863 if ( (iPdg != kPdgElectron) ||
c4214bc0 864 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
f73816f5 865 gMC->TrackMomentum(mom);
866 pTot = mom.Rho();
867 betaGamma = pTot / aMass;
868 pp = kPrim * BetheBloch(betaGamma);
869 // Take charge > 1 into account
870 charge = gMC->TrackCharge();
871 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
c4214bc0 872 } else { // Electrons above 20 Mev/c are at the plateau
873 pp = kPrim * kPlateau;
f73816f5 874 }
875
876 if (pp > 0) {
877 do
b9d0a01d 878 gMC->GetRandom()->RndmArray(1, random);
f73816f5 879 while ((random[0] == 1.) || (random[0] == 0.));
880 stepSize = - TMath::Log(random[0]) / pp;
881 gMC->SetMaxStep(stepSize);
c4214bc0 882 }
883 }
5c7f4665 884 }
d3f347ff 885 }
5c7f4665 886 }
5c7f4665 887}
888
a328fff9 889//_____________________________________________________________________________
890void AliTRDv1::StepManagerFixedStep()
891{
892 //
893 // Slow simulator. Every charged track produces electron cluster as hits
894 // along its path across the drift volume. The step size is fixed in
895 // this version of the step manager.
896 //
897
898 Int_t pla = 0;
899 Int_t cha = 0;
900 Int_t sec = 0;
901 Int_t det = 0;
902 Int_t qTot;
903
904 Float_t hits[3];
905 Double_t eDep;
906
907 Bool_t drRegion = kFALSE;
908 Bool_t amRegion = kFALSE;
909
910 TString cIdCurrent;
911 TString cIdSensDr = "J";
912 TString cIdSensAm = "K";
913 Char_t cIdChamber[3];
914 cIdChamber[2] = 0;
915
916 TLorentzVector pos, mom;
917
918 const Int_t kNplan = AliTRDgeometry::Nplan();
919 const Int_t kNcham = AliTRDgeometry::Ncham();
920 const Int_t kNdetsec = kNplan * kNcham;
921
922 const Double_t kBig = 1.0E+12;
923
bc327ce2 924 const Float_t kWion = 23.53; // Ionization energy
a328fff9 925 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
926
927 // Set the maximum step size to a very large number for all
928 // neutral particles and those outside the driftvolume
929 gMC->SetMaxStep(kBig);
930
931 // If not charged track or already stopped or disappeared, just return.
932 if ((!gMC->TrackCharge()) ||
933 gMC->IsTrackStop() ||
934 gMC->IsTrackDisappeared()) return;
935
936 // Inside a sensitive volume?
937 cIdCurrent = gMC->CurrentVolName();
938
939 if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
940 if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
941
942 if ((!drRegion) && (!amRegion)) return;
943
944 // The hit coordinates and charge
945 gMC->TrackPosition(pos);
946 hits[0] = pos[0];
947 hits[1] = pos[1];
948 hits[2] = pos[2];
949
950 // The sector number (0 - 17)
951 // The numbering goes clockwise and starts at y = 0
952 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
953 if (phi < 90.) phi += 270.;
954 else phi -= 90.;
955 sec = ((Int_t) (phi / 20.));
956
957 // The plane and chamber number
958 cIdChamber[0] = cIdCurrent[2];
959 cIdChamber[1] = cIdCurrent[3];
960 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
a5cadd36 961 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
a328fff9 962 pla = ((Int_t) idChamber % kNplan);
e0d47c25 963
a328fff9 964 // Check on selected volumes
965 Int_t addthishit = 1;
966 if(fSensSelect) {
967 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
968 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
969 if (fSensSector >= 0) {
970 Int_t sens1 = fSensSector;
971 Int_t sens2 = fSensSector + fSensSectorRange;
972 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
973 if (sens1 < sens2) {
974 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
975 }
976 else {
977 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
978 }
979 }
980 }
981
982 if (!addthishit) return;
983
984 det = fGeometry->GetDetector(pla,cha,sec); // The detector number
985
986 Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting
987
988 // Special hits only in the drift region
989 if (drRegion) {
990
991 // Create a track reference at the entrance and exit of each
992 // chamber that contain the momentum components of the particle
993
994 if (gMC->IsTrackEntering()) {
995 gMC->TrackMomentum(mom);
996 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
997 trkStat = 1;
998 }
999 if (gMC->IsTrackExiting()) {
1000 gMC->TrackMomentum(mom);
1001 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1002 trkStat = 2;
1003 }
1004
1005 // Create the hits from TR photons
1006 if (fTR) CreateTRhit(det);
1007
1008 }
1009
1010 // Calculate the charge according to GEANT Edep
1011 // Create a new dEdx hit
1012 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1013 qTot = (Int_t) (eDep / kWion);
c4214bc0 1014 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1015 ,det,hits,qTot,drRegion);
a328fff9 1016
1017 // Set Maximum Step Size
1018 // Produce only one hit if Ekin is below cutoff
1019 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
1020 gMC->SetMaxStep(fStepSize);
1021
1022}
1023
5c7f4665 1024//_____________________________________________________________________________
1025Double_t AliTRDv1::BetheBloch(Double_t bg)
1026{
1027 //
1028 // Parametrization of the Bethe-Bloch-curve
1029 // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1030 //
1031
1032 // This parameters have been adjusted to averaged values from GEANT
1033 const Double_t kP1 = 7.17960e-02;
1034 const Double_t kP2 = 8.54196;
1035 const Double_t kP3 = 1.38065e-06;
1036 const Double_t kP4 = 5.30972;
1037 const Double_t kP5 = 2.83798;
1038
1039 // This parameters have been adjusted to Xe-data found in:
1040 // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
1041 //const Double_t kP1 = 0.76176E-1;
1042 //const Double_t kP2 = 10.632;
1043 //const Double_t kP3 = 3.17983E-6;
1044 //const Double_t kP4 = 1.8631;
1045 //const Double_t kP5 = 1.9479;
1046
f73816f5 1047 // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1048 const Double_t kBgMin = 0.8;
1049 const Double_t kBBMax = 6.83298;
1050 //const Double_t kBgMin = 0.6;
1051 //const Double_t kBBMax = 17.2809;
1052 //const Double_t kBgMin = 0.4;
1053 //const Double_t kBBMax = 82.0;
1054
1055 if (bg > kBgMin) {
5c7f4665 1056 Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
1057 Double_t aa = TMath::Power(yy,kP4);
1058 Double_t bb = TMath::Power((1./bg),kP5);
1059 bb = TMath::Log(kP3 + bb);
1060 return ((kP2 - aa - bb)*kP1 / aa);
1061 }
f73816f5 1062 else {
1063 return kBBMax;
1064 }
d3f347ff 1065
fe4da5cc 1066}
5c7f4665 1067
a328fff9 1068//_____________________________________________________________________________
c4214bc0 1069Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
a328fff9 1070{
1071 //
1072 // Return dN/dx (number of primary collisions per centimeter)
1073 // for given beta*gamma factor.
1074 //
1075 // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1076 // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1077 // This must be used as a set with IntSpecGeant.
1078 //
1079
1080 Double_t arr_g[20] = {
1081 1.100000, 1.200000, 1.300000, 1.500000,
1082 1.800000, 2.000000, 2.500000, 3.000000,
1083 4.000000, 7.000000, 10.000000, 20.000000,
1084 40.000000, 70.000000, 100.000000, 300.000000,
1085 600.000000, 1000.000000, 3000.000000, 10000.000000 };
1086
1087 Double_t arr_nc[20] = {
1088 75.009056, 45.508083, 35.299252, 27.116327,
1089 22.734999, 21.411915, 19.934095, 19.449375,
1090 19.344431, 20.185553, 21.027925, 22.912676,
1091 24.933352, 26.504053, 27.387468, 29.566597,
1092 30.353779, 30.787134, 31.129285, 31.157350 };
1093
1094 // betagamma to gamma
1095 Double_t g = TMath::Sqrt( 1. + bg*bg );
1096
1097 // Find the index just before the point we need.
1098 int i;
1099 for( i = 0 ; i < 18 ; i++ )
1100 if( arr_g[i] < g && arr_g[i+1] > g )
1101 break;
1102
1103 // Simple interpolation.
1104 Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
1105 (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
1106
c4214bc0 1107 return pp; //arr_nc[8];
a328fff9 1108
1109}
1110
c4214bc0 1111//_____________________________________________________________________________
1112void AliTRDv1::Stepping()
1113{
1114// Stepping info
1115// ---
1116
1117 cout << "X(cm) "
1118 << "Y(cm) "
1119 << "Z(cm) "
1120 << "KinE(MeV) "
1121 << "dE(MeV) "
1122 << "Step(cm) "
1123 << "TrackL(cm) "
1124 << "Volume "
1125 << "Process "
1126 << endl;
1127
1128 // Position
1129 //
1130 Double_t x, y, z;
1131 gMC->TrackPosition(x, y, z);
1132 cout << setw(8) << setprecision(3) << x << " "
1133 << setw(8) << setprecision(3) << y << " "
1134 << setw(8) << setprecision(3) << z << " ";
1135
1136 // Kinetic energy
1137 //
1138 Double_t px, py, pz, etot;
1139 gMC->TrackMomentum(px, py, pz, etot);
1140 Double_t ekin = etot - gMC->TrackMass();
1141 cout << setw(9) << setprecision(4) << ekin*1e03 << " ";
1142
1143 // Energy deposit
1144 //
1145 cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " ";
1146
1147 // Step length
1148 //
1149 cout << setw(8) << setprecision(3) << gMC->TrackStep() << " ";
1150
1151 // Track length
1152 //
1153 cout << setw(8) << setprecision(3) << gMC->TrackLength() << " ";
1154
1155 // Volume
1156 //
1157 if (gMC->CurrentVolName() != 0)
1158 cout << setw(4) << gMC->CurrentVolName() << " ";
1159 else
1160 cout << setw(4) << "None" << " ";
1161
1162 // Process
1163 //
1164 TArrayI processes;
1165 Int_t nofProcesses = gMC->StepProcesses(processes);
1166 for(int ip=0;ip<nofProcesses; ip++)
1167 cout << TMCProcessName[processes[ip]]<<" / ";
1168
1169 cout << endl;
1170}
1171
1172
5c7f4665 1173//_____________________________________________________________________________
1174Double_t Ermilova(Double_t *x, Double_t *)
1175{
1176 //
1177 // Calculates the delta-ray energy distribution according to Ermilova.
1178 // Logarithmic scale !
1179 //
1180
1181 Double_t energy;
1182 Double_t dpos;
1183 Double_t dnde;
1184
1185 Int_t pos1, pos2;
1186
8230f242 1187 const Int_t kNv = 31;
5c7f4665 1188
8230f242 1189 Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
1190 , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
1191 , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
1192 , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
1193 , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
1194 , 9.4727, 9.9035,10.3735,10.5966,10.8198
1195 ,11.5129 };
5c7f4665 1196
8230f242 1197 Float_t vye[kNv] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0
1198 , 20.9 , 20.8 , 20.0 , 16.0 , 11.0
1199 , 8.0 , 6.0 , 5.2 , 4.6 , 4.0
1200 , 3.5 , 3.0 , 1.4 , 0.67 , 0.44
1201 , 0.3 , 0.18 , 0.12 , 0.08 , 0.056
1202 , 0.04 , 0.023, 0.015, 0.011, 0.01
1203 , 0.004 };
5c7f4665 1204
1205 energy = x[0];
1206
1207 // Find the position
1208 pos1 = pos2 = 0;
1209 dpos = 0;
1210 do {
1211 dpos = energy - vxe[pos2++];
1212 }
1213 while (dpos > 0);
1214 pos2--;
f73816f5 1215 if (pos2 > kNv) pos2 = kNv - 1;
5c7f4665 1216 pos1 = pos2 - 1;
1217
1218 // Differentiate between the sampling points
1219 dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1220
1221 return dnde;
1222
1223}
a328fff9 1224
1225//_____________________________________________________________________________
1226Double_t IntSpecGeant(Double_t *x, Double_t *)
1227{
1228 //
1229 // Integrated spectrum from Geant3
1230 //
1231
96efaf83 1232 const Int_t npts = 83;
1233 Double_t arre[npts] = {
a328fff9 1234 2.421257, 2.483278, 2.534301, 2.592230,
1235 2.672067, 2.813299, 3.015059, 3.216819,
1236 3.418579, 3.620338, 3.868209, 3.920198,
1237 3.978284, 4.063923, 4.186264, 4.308605,
1238 4.430946, 4.553288, 4.724261, 4.837736,
1239 4.999842, 5.161949, 5.324056, 5.486163,
1240 5.679688, 5.752998, 5.857728, 5.962457,
1241 6.067185, 6.171914, 6.315653, 6.393674,
1242 6.471694, 6.539689, 6.597658, 6.655627,
1243 6.710957, 6.763648, 6.816338, 6.876198,
1244 6.943227, 7.010257, 7.106285, 7.252151,
1245 7.460531, 7.668911, 7.877290, 8.085670,
1246 8.302979, 8.353585, 8.413120, 8.483500,
1247 8.541030, 8.592857, 8.668865, 8.820485,
1248 9.037086, 9.253686, 9.470286, 9.686887,
1249 9.930838, 9.994655, 10.085822, 10.176990,
1250 10.268158, 10.359325, 10.503614, 10.627565,
1251 10.804637, 10.981709, 11.158781, 11.335854,
1252 11.593397, 11.781165, 12.049404, 12.317644,
1253 12.585884, 12.854123, 14.278421, 16.975889,
1254 20.829416, 24.682943, 28.536469
1255 };
96efaf83 1256 /*
1257 Double_t arrdndx[npts] = {
a328fff9 1258 19.344431, 18.664679, 18.136106, 17.567745,
1259 16.836426, 15.677382, 14.281277, 13.140237,
1260 12.207677, 11.445510, 10.697049, 10.562296,
1261 10.414673, 10.182341, 9.775256, 9.172330,
1262 8.240271, 6.898587, 4.808303, 3.889751,
1263 3.345288, 3.093431, 2.897347, 2.692470,
1264 2.436222, 2.340029, 2.208579, 2.086489,
1265 1.975535, 1.876519, 1.759626, 1.705024,
1266 1.656374, 1.502638, 1.330566, 1.200697,
1267 1.101168, 1.019323, 0.943867, 0.851951,
1268 0.755229, 0.671576, 0.570675, 0.449672,
1269 0.326722, 0.244225, 0.188225, 0.149608,
1270 0.121529, 0.116289, 0.110636, 0.103490,
1271 0.096147, 0.089191, 0.079780, 0.063927,
1272 0.047642, 0.036341, 0.028250, 0.022285,
1273 0.017291, 0.016211, 0.014802, 0.013533,
1274 0.012388, 0.011352, 0.009803, 0.008537,
1275 0.007039, 0.005829, 0.004843, 0.004034,
1276 0.003101, 0.002564, 0.001956, 0.001494,
1277 0.001142, 0.000873, 0.000210, 0.000014,
1278 0.000000, 0.000000, 0.000000
1279 };
96efaf83 1280 */
1281 // Differentiate
1282 // dnde = (arrdndx[i-1] - arrdndx[i]) / (arre[i] - arre[i-1]);
1283
1284 Double_t arrdnde[npts] = {
1285 10.960000, 10.960000, 10.359500, 9.811340,
1286 9.1601500, 8.206670, 6.919630, 5.655430,
1287 4.6221300, 3.777610, 3.019560, 2.591950,
1288 2.5414600, 2.712920, 3.327460, 4.928240,
1289 7.6185300, 10.966700, 12.225800, 8.094750,
1290 3.3586900, 1.553650, 1.209600, 1.263840,
1291 1.3241100, 1.312140, 1.255130, 1.165770,
1292 1.0594500, 0.945450, 0.813231, 0.699837,
1293 0.6235580, 2.260990, 2.968350, 2.240320,
1294 1.7988300, 1.553300, 1.432070, 1.535520,
1295 1.4429900, 1.247990, 1.050750, 0.829549,
1296 0.5900280, 0.395897, 0.268741, 0.185320,
1297 0.1292120, 0.103545, 0.0949525, 0.101535,
1298 0.1276380, 0.134216, 0.123816, 0.104557,
1299 0.0751843, 0.0521745, 0.0373546, 0.0275391,
1300 0.0204713, 0.0169234, 0.0154552, 0.0139194,
1301 0.0125592, 0.0113638, 0.0107354, 0.0102137,
1302 0.00845984, 0.00683338, 0.00556836, 0.00456874,
1303 0.0036227, 0.00285991, 0.00226664, 0.00172234,
1304 0.00131226, 0.00100284, 0.000465492, 7.26607e-05,
1305 3.63304e-06, 0.0000000, 0.0000000
1306 };
a328fff9 1307
1308 Int_t i;
1309 Double_t energy = x[0];
a328fff9 1310
96efaf83 1311 for( i = 0 ; i < npts ; i++ )
1312 if( energy < arre[i] ) break;
a328fff9 1313
1314 if( i == 0 )
45160b1f 1315 AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");
a328fff9 1316
96efaf83 1317 return arrdnde[i];
a328fff9 1318
1319}