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