]> 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
Change of C02 density
[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
c4214bc0 199 const Float_t kWion = 22.57;
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
501 const Float_t kWion = 22.57; // Ionization energy
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
678 stepSize = 1./gRandom->Poisson(pp);
679 gMC->SetMaxStep(stepSize);
680 }
681 }
682 }
683 }
a328fff9 684}
685
686//_____________________________________________________________________________
687void AliTRDv1::StepManagerErmilova()
5c7f4665 688{
689 //
5c7f4665 690 // Slow simulator. Every charged track produces electron cluster as hits
691 // along its path across the drift volume. The step size is set acording
692 // to Bethe-Bloch. The energy distribution of the delta electrons follows
693 // a spectrum taken from Ermilova et al.
694 //
695
851d3db9 696 Int_t pla = 0;
697 Int_t cha = 0;
698 Int_t sec = 0;
793ff80c 699 Int_t det = 0;
851d3db9 700 Int_t iPdg;
793ff80c 701 Int_t qTot;
5c7f4665 702
793ff80c 703 Float_t hits[3];
a5cadd36 704 Double_t random[1];
5c7f4665 705 Float_t charge;
706 Float_t aMass;
707
f73816f5 708 Double_t pTot = 0;
5c7f4665 709 Double_t eDelta;
710 Double_t betaGamma, pp;
f73816f5 711 Double_t stepSize;
5c7f4665 712
332e9569 713 Bool_t drRegion = kFALSE;
714 Bool_t amRegion = kFALSE;
715
716 TString cIdCurrent;
717 TString cIdSensDr = "J";
718 TString cIdSensAm = "K";
593a9fc3 719 Char_t cIdChamber[3];
720 cIdChamber[2] = 0;
332e9569 721
5c7f4665 722 TLorentzVector pos, mom;
82bbf98a 723
332e9569 724 const Int_t kNplan = AliTRDgeometry::Nplan();
e644678a 725 const Int_t kNcham = AliTRDgeometry::Ncham();
726 const Int_t kNdetsec = kNplan * kNcham;
727
a328fff9 728 const Double_t kBig = 1.0E+12; // Infinitely big
c4214bc0 729 const Float_t kWion = 22.57; // Ionization energy
a328fff9 730 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
5c7f4665 731
c4214bc0 732 // energy threshold for production of delta electrons
733 //const Float_t kECut = 1.0e4;
734 // Parameters entering the parametrized range for delta electrons
735 //const float ra=5.37E-4, rb=0.9815, rc=3.123E-3;
736 // Gas density -> To be made user adjustable !
737 //const float rho=0.004945 ; //[0.85*0.00549+0.15*0.00186 (Xe-CO2 85-15)]
738
f73816f5 739 // Minimum energy for the step size adjustment
740 const Float_t kEkinMinStep = 1.0e-5;
a328fff9 741
5c7f4665 742 // Plateau value of the energy-loss for electron in xenon
743 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
744 //const Double_t kPlateau = 1.70;
745 // the averaged value (26/3/99)
a3c76cdc 746 const Float_t kPlateau = 1.55;
a328fff9 747
748 const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
5c7f4665 749 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
a3c76cdc 750 const Float_t kPoti = 12.1;
851d3db9 751
a328fff9 752 const Int_t kPdgElectron = 11; // PDG code electron
5c7f4665 753
754 // Set the maximum step size to a very large number for all
755 // neutral particles and those outside the driftvolume
756 gMC->SetMaxStep(kBig);
757
758 // Use only charged tracks
759 if (( gMC->TrackCharge() ) &&
760 (!gMC->IsTrackStop() ) &&
761 (!gMC->IsTrackDisappeared())) {
fe4da5cc 762
5c7f4665 763 // Inside a sensitive volume?
332e9569 764 drRegion = kFALSE;
765 amRegion = kFALSE;
766 cIdCurrent = gMC->CurrentVolName();
e6674585 767 if (cIdSensDr == cIdCurrent[1]) {
332e9569 768 drRegion = kTRUE;
769 }
e6674585 770 if (cIdSensAm == cIdCurrent[1]) {
332e9569 771 amRegion = kTRUE;
772 }
773 if (drRegion || amRegion) {
fe4da5cc 774
5c7f4665 775 // The hit coordinates and charge
776 gMC->TrackPosition(pos);
777 hits[0] = pos[0];
778 hits[1] = pos[1];
779 hits[2] = pos[2];
5c7f4665 780
851d3db9 781 // The sector number (0 - 17)
782 // The numbering goes clockwise and starts at y = 0
e15eb584 783 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
851d3db9 784 if (phi < 90.)
785 phi = phi + 270.;
786 else
787 phi = phi - 90.;
788 sec = ((Int_t) (phi / 20));
82bbf98a 789
332e9569 790 // The plane and chamber number
791 cIdChamber[0] = cIdCurrent[2];
792 cIdChamber[1] = cIdCurrent[3];
e644678a 793 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
a5cadd36 794 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
332e9569 795 pla = ((Int_t) idChamber % kNplan);
82bbf98a 796
5c7f4665 797 // Check on selected volumes
798 Int_t addthishit = 1;
799 if (fSensSelect) {
6f1e466d 800 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
801 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
9d0b222b 802 if (fSensSector >= 0) {
803 Int_t sens1 = fSensSector;
804 Int_t sens2 = fSensSector + fSensSectorRange;
793ff80c 805 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
806 * AliTRDgeometry::Nsect();
9d0b222b 807 if (sens1 < sens2) {
808 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
c4214bc0 809 }
9d0b222b 810 else {
811 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
c4214bc0 812 }
813 }
5c7f4665 814 }
815
816 // Add this hit
817 if (addthishit) {
818
f73816f5 819 // The detector number
793ff80c 820 det = fGeometry->GetDetector(pla,cha,sec);
821
a328fff9 822 // Special hits only in the drift region
332e9569 823 if (drRegion) {
f73816f5 824
c61f1a66 825 // Create a track reference at the entrance and
826 // exit of each chamber that contain the
827 // momentum components of the particle
f73816f5 828 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
829 gMC->TrackMomentum(mom);
5d12ce38 830 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
f73816f5 831 }
f73816f5 832 // Create the hits from TR photons
833 if (fTR) CreateTRhit(det);
c4214bc0 834 }
f73816f5 835
f73816f5 836
837 // Calculate the energy of the delta-electrons
838 eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
839 eDelta = TMath::Max(eDelta,0.0);
c4214bc0 840 // Generate the electron cluster size
841 if(eDelta==0.) qTot=0;
842 else qTot = ((Int_t) (eDelta / kWion) + 1);
f73816f5 843
c4214bc0 844 // Create a new dEdx hit
332e9569 845 if (drRegion) {
a328fff9 846 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
c4214bc0 847 ,det,hits,qTot, kTRUE);
848 }
5c7f4665 849 else {
a328fff9 850 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
c4214bc0 851 ,det,hits,qTot,kFALSE);
852 }
f73816f5 853
5c7f4665 854 // Calculate the maximum step size for the next tracking step
f73816f5 855 // Produce only one hit if Ekin is below cutoff
856 aMass = gMC->TrackMass();
857 if ((gMC->Etot() - aMass) > kEkinMinStep) {
858
859 // The energy loss according to Bethe Bloch
860 iPdg = TMath::Abs(gMC->TrackPid());
861 if ( (iPdg != kPdgElectron) ||
c4214bc0 862 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
f73816f5 863 gMC->TrackMomentum(mom);
864 pTot = mom.Rho();
865 betaGamma = pTot / aMass;
866 pp = kPrim * BetheBloch(betaGamma);
867 // Take charge > 1 into account
868 charge = gMC->TrackCharge();
869 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
c4214bc0 870 } else { // Electrons above 20 Mev/c are at the plateau
871 pp = kPrim * kPlateau;
f73816f5 872 }
873
874 if (pp > 0) {
875 do
b9d0a01d 876 gMC->GetRandom()->RndmArray(1, random);
f73816f5 877 while ((random[0] == 1.) || (random[0] == 0.));
878 stepSize = - TMath::Log(random[0]) / pp;
879 gMC->SetMaxStep(stepSize);
c4214bc0 880 }
881 }
5c7f4665 882 }
d3f347ff 883 }
5c7f4665 884 }
5c7f4665 885}
886
a328fff9 887//_____________________________________________________________________________
888void AliTRDv1::StepManagerFixedStep()
889{
890 //
891 // Slow simulator. Every charged track produces electron cluster as hits
892 // along its path across the drift volume. The step size is fixed in
893 // this version of the step manager.
894 //
895
896 Int_t pla = 0;
897 Int_t cha = 0;
898 Int_t sec = 0;
899 Int_t det = 0;
900 Int_t qTot;
901
902 Float_t hits[3];
903 Double_t eDep;
904
905 Bool_t drRegion = kFALSE;
906 Bool_t amRegion = kFALSE;
907
908 TString cIdCurrent;
909 TString cIdSensDr = "J";
910 TString cIdSensAm = "K";
911 Char_t cIdChamber[3];
912 cIdChamber[2] = 0;
913
914 TLorentzVector pos, mom;
915
916 const Int_t kNplan = AliTRDgeometry::Nplan();
917 const Int_t kNcham = AliTRDgeometry::Ncham();
918 const Int_t kNdetsec = kNplan * kNcham;
919
920 const Double_t kBig = 1.0E+12;
921
c4214bc0 922 const Float_t kWion = 22.57; // Ionization energy
a328fff9 923 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
924
925 // Set the maximum step size to a very large number for all
926 // neutral particles and those outside the driftvolume
927 gMC->SetMaxStep(kBig);
928
929 // If not charged track or already stopped or disappeared, just return.
930 if ((!gMC->TrackCharge()) ||
931 gMC->IsTrackStop() ||
932 gMC->IsTrackDisappeared()) return;
933
934 // Inside a sensitive volume?
935 cIdCurrent = gMC->CurrentVolName();
936
937 if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
938 if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
939
940 if ((!drRegion) && (!amRegion)) return;
941
942 // The hit coordinates and charge
943 gMC->TrackPosition(pos);
944 hits[0] = pos[0];
945 hits[1] = pos[1];
946 hits[2] = pos[2];
947
948 // The sector number (0 - 17)
949 // The numbering goes clockwise and starts at y = 0
950 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
951 if (phi < 90.) phi += 270.;
952 else phi -= 90.;
953 sec = ((Int_t) (phi / 20.));
954
955 // The plane and chamber number
956 cIdChamber[0] = cIdCurrent[2];
957 cIdChamber[1] = cIdCurrent[3];
958 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
a5cadd36 959 cha = kNcham - ((Int_t) idChamber / kNplan) - 1;
a328fff9 960 pla = ((Int_t) idChamber % kNplan);
e0d47c25 961
a328fff9 962 // Check on selected volumes
963 Int_t addthishit = 1;
964 if(fSensSelect) {
965 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
966 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
967 if (fSensSector >= 0) {
968 Int_t sens1 = fSensSector;
969 Int_t sens2 = fSensSector + fSensSectorRange;
970 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
971 if (sens1 < sens2) {
972 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
973 }
974 else {
975 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
976 }
977 }
978 }
979
980 if (!addthishit) return;
981
982 det = fGeometry->GetDetector(pla,cha,sec); // The detector number
983
984 Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting
985
986 // Special hits only in the drift region
987 if (drRegion) {
988
989 // Create a track reference at the entrance and exit of each
990 // chamber that contain the momentum components of the particle
991
992 if (gMC->IsTrackEntering()) {
993 gMC->TrackMomentum(mom);
994 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
995 trkStat = 1;
996 }
997 if (gMC->IsTrackExiting()) {
998 gMC->TrackMomentum(mom);
999 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
1000 trkStat = 2;
1001 }
1002
1003 // Create the hits from TR photons
1004 if (fTR) CreateTRhit(det);
1005
1006 }
1007
1008 // Calculate the charge according to GEANT Edep
1009 // Create a new dEdx hit
1010 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
1011 qTot = (Int_t) (eDep / kWion);
c4214bc0 1012 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
1013 ,det,hits,qTot,drRegion);
a328fff9 1014
1015 // Set Maximum Step Size
1016 // Produce only one hit if Ekin is below cutoff
1017 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
1018 gMC->SetMaxStep(fStepSize);
1019
1020}
1021
5c7f4665 1022//_____________________________________________________________________________
1023Double_t AliTRDv1::BetheBloch(Double_t bg)
1024{
1025 //
1026 // Parametrization of the Bethe-Bloch-curve
1027 // The parametrization is the same as for the TPC and is taken from Lehrhaus.
1028 //
1029
1030 // This parameters have been adjusted to averaged values from GEANT
1031 const Double_t kP1 = 7.17960e-02;
1032 const Double_t kP2 = 8.54196;
1033 const Double_t kP3 = 1.38065e-06;
1034 const Double_t kP4 = 5.30972;
1035 const Double_t kP5 = 2.83798;
1036
1037 // This parameters have been adjusted to Xe-data found in:
1038 // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
1039 //const Double_t kP1 = 0.76176E-1;
1040 //const Double_t kP2 = 10.632;
1041 //const Double_t kP3 = 3.17983E-6;
1042 //const Double_t kP4 = 1.8631;
1043 //const Double_t kP5 = 1.9479;
1044
f73816f5 1045 // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
1046 const Double_t kBgMin = 0.8;
1047 const Double_t kBBMax = 6.83298;
1048 //const Double_t kBgMin = 0.6;
1049 //const Double_t kBBMax = 17.2809;
1050 //const Double_t kBgMin = 0.4;
1051 //const Double_t kBBMax = 82.0;
1052
1053 if (bg > kBgMin) {
5c7f4665 1054 Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
1055 Double_t aa = TMath::Power(yy,kP4);
1056 Double_t bb = TMath::Power((1./bg),kP5);
1057 bb = TMath::Log(kP3 + bb);
1058 return ((kP2 - aa - bb)*kP1 / aa);
1059 }
f73816f5 1060 else {
1061 return kBBMax;
1062 }
d3f347ff 1063
fe4da5cc 1064}
5c7f4665 1065
a328fff9 1066//_____________________________________________________________________________
c4214bc0 1067Double_t AliTRDv1::BetheBlochGeant(Double_t bg)
a328fff9 1068{
1069 //
1070 // Return dN/dx (number of primary collisions per centimeter)
1071 // for given beta*gamma factor.
1072 //
1073 // Implemented by K.Oyama according to GEANT 3 parametrization shown in
1074 // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
1075 // This must be used as a set with IntSpecGeant.
1076 //
1077
1078 Double_t arr_g[20] = {
1079 1.100000, 1.200000, 1.300000, 1.500000,
1080 1.800000, 2.000000, 2.500000, 3.000000,
1081 4.000000, 7.000000, 10.000000, 20.000000,
1082 40.000000, 70.000000, 100.000000, 300.000000,
1083 600.000000, 1000.000000, 3000.000000, 10000.000000 };
1084
1085 Double_t arr_nc[20] = {
1086 75.009056, 45.508083, 35.299252, 27.116327,
1087 22.734999, 21.411915, 19.934095, 19.449375,
1088 19.344431, 20.185553, 21.027925, 22.912676,
1089 24.933352, 26.504053, 27.387468, 29.566597,
1090 30.353779, 30.787134, 31.129285, 31.157350 };
1091
1092 // betagamma to gamma
1093 Double_t g = TMath::Sqrt( 1. + bg*bg );
1094
1095 // Find the index just before the point we need.
1096 int i;
1097 for( i = 0 ; i < 18 ; i++ )
1098 if( arr_g[i] < g && arr_g[i+1] > g )
1099 break;
1100
1101 // Simple interpolation.
1102 Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
1103 (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
1104
c4214bc0 1105 return pp; //arr_nc[8];
a328fff9 1106
1107}
1108
c4214bc0 1109//_____________________________________________________________________________
1110void AliTRDv1::Stepping()
1111{
1112// Stepping info
1113// ---
1114
1115 cout << "X(cm) "
1116 << "Y(cm) "
1117 << "Z(cm) "
1118 << "KinE(MeV) "
1119 << "dE(MeV) "
1120 << "Step(cm) "
1121 << "TrackL(cm) "
1122 << "Volume "
1123 << "Process "
1124 << endl;
1125
1126 // Position
1127 //
1128 Double_t x, y, z;
1129 gMC->TrackPosition(x, y, z);
1130 cout << setw(8) << setprecision(3) << x << " "
1131 << setw(8) << setprecision(3) << y << " "
1132 << setw(8) << setprecision(3) << z << " ";
1133
1134 // Kinetic energy
1135 //
1136 Double_t px, py, pz, etot;
1137 gMC->TrackMomentum(px, py, pz, etot);
1138 Double_t ekin = etot - gMC->TrackMass();
1139 cout << setw(9) << setprecision(4) << ekin*1e03 << " ";
1140
1141 // Energy deposit
1142 //
1143 cout << setw(9) << setprecision(4) << gMC->Edep()*1e03 << " ";
1144
1145 // Step length
1146 //
1147 cout << setw(8) << setprecision(3) << gMC->TrackStep() << " ";
1148
1149 // Track length
1150 //
1151 cout << setw(8) << setprecision(3) << gMC->TrackLength() << " ";
1152
1153 // Volume
1154 //
1155 if (gMC->CurrentVolName() != 0)
1156 cout << setw(4) << gMC->CurrentVolName() << " ";
1157 else
1158 cout << setw(4) << "None" << " ";
1159
1160 // Process
1161 //
1162 TArrayI processes;
1163 Int_t nofProcesses = gMC->StepProcesses(processes);
1164 for(int ip=0;ip<nofProcesses; ip++)
1165 cout << TMCProcessName[processes[ip]]<<" / ";
1166
1167 cout << endl;
1168}
1169
1170
5c7f4665 1171//_____________________________________________________________________________
1172Double_t Ermilova(Double_t *x, Double_t *)
1173{
1174 //
1175 // Calculates the delta-ray energy distribution according to Ermilova.
1176 // Logarithmic scale !
1177 //
1178
1179 Double_t energy;
1180 Double_t dpos;
1181 Double_t dnde;
1182
1183 Int_t pos1, pos2;
1184
8230f242 1185 const Int_t kNv = 31;
5c7f4665 1186
8230f242 1187 Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
1188 , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
1189 , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
1190 , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
1191 , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
1192 , 9.4727, 9.9035,10.3735,10.5966,10.8198
1193 ,11.5129 };
5c7f4665 1194
8230f242 1195 Float_t vye[kNv] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0
1196 , 20.9 , 20.8 , 20.0 , 16.0 , 11.0
1197 , 8.0 , 6.0 , 5.2 , 4.6 , 4.0
1198 , 3.5 , 3.0 , 1.4 , 0.67 , 0.44
1199 , 0.3 , 0.18 , 0.12 , 0.08 , 0.056
1200 , 0.04 , 0.023, 0.015, 0.011, 0.01
1201 , 0.004 };
5c7f4665 1202
1203 energy = x[0];
1204
1205 // Find the position
1206 pos1 = pos2 = 0;
1207 dpos = 0;
1208 do {
1209 dpos = energy - vxe[pos2++];
1210 }
1211 while (dpos > 0);
1212 pos2--;
f73816f5 1213 if (pos2 > kNv) pos2 = kNv - 1;
5c7f4665 1214 pos1 = pos2 - 1;
1215
1216 // Differentiate between the sampling points
1217 dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
1218
1219 return dnde;
1220
1221}
a328fff9 1222
1223//_____________________________________________________________________________
1224Double_t IntSpecGeant(Double_t *x, Double_t *)
1225{
1226 //
1227 // Integrated spectrum from Geant3
1228 //
1229
1230 const Int_t n_pts = 83;
1231 Double_t arr_e[n_pts] = {
1232 2.421257, 2.483278, 2.534301, 2.592230,
1233 2.672067, 2.813299, 3.015059, 3.216819,
1234 3.418579, 3.620338, 3.868209, 3.920198,
1235 3.978284, 4.063923, 4.186264, 4.308605,
1236 4.430946, 4.553288, 4.724261, 4.837736,
1237 4.999842, 5.161949, 5.324056, 5.486163,
1238 5.679688, 5.752998, 5.857728, 5.962457,
1239 6.067185, 6.171914, 6.315653, 6.393674,
1240 6.471694, 6.539689, 6.597658, 6.655627,
1241 6.710957, 6.763648, 6.816338, 6.876198,
1242 6.943227, 7.010257, 7.106285, 7.252151,
1243 7.460531, 7.668911, 7.877290, 8.085670,
1244 8.302979, 8.353585, 8.413120, 8.483500,
1245 8.541030, 8.592857, 8.668865, 8.820485,
1246 9.037086, 9.253686, 9.470286, 9.686887,
1247 9.930838, 9.994655, 10.085822, 10.176990,
1248 10.268158, 10.359325, 10.503614, 10.627565,
1249 10.804637, 10.981709, 11.158781, 11.335854,
1250 11.593397, 11.781165, 12.049404, 12.317644,
1251 12.585884, 12.854123, 14.278421, 16.975889,
1252 20.829416, 24.682943, 28.536469
1253 };
1254 Double_t arr_dndx[n_pts] = {
1255 19.344431, 18.664679, 18.136106, 17.567745,
1256 16.836426, 15.677382, 14.281277, 13.140237,
1257 12.207677, 11.445510, 10.697049, 10.562296,
1258 10.414673, 10.182341, 9.775256, 9.172330,
1259 8.240271, 6.898587, 4.808303, 3.889751,
1260 3.345288, 3.093431, 2.897347, 2.692470,
1261 2.436222, 2.340029, 2.208579, 2.086489,
1262 1.975535, 1.876519, 1.759626, 1.705024,
1263 1.656374, 1.502638, 1.330566, 1.200697,
1264 1.101168, 1.019323, 0.943867, 0.851951,
1265 0.755229, 0.671576, 0.570675, 0.449672,
1266 0.326722, 0.244225, 0.188225, 0.149608,
1267 0.121529, 0.116289, 0.110636, 0.103490,
1268 0.096147, 0.089191, 0.079780, 0.063927,
1269 0.047642, 0.036341, 0.028250, 0.022285,
1270 0.017291, 0.016211, 0.014802, 0.013533,
1271 0.012388, 0.011352, 0.009803, 0.008537,
1272 0.007039, 0.005829, 0.004843, 0.004034,
1273 0.003101, 0.002564, 0.001956, 0.001494,
1274 0.001142, 0.000873, 0.000210, 0.000014,
1275 0.000000, 0.000000, 0.000000
1276 };
1277
1278 Int_t i;
1279 Double_t energy = x[0];
1280 Double_t dnde;
1281
1282 for( i = 0 ; i < n_pts ; i++ )
1283 if( energy < arr_e[i] ) break;
1284
1285 if( i == 0 )
45160b1f 1286 AliErrorGeneral("AliTRDv1","Given energy value is too small or zero");
a328fff9 1287
c4214bc0 1288 // Differentiate
a328fff9 1289 dnde = (arr_dndx[i-1] - arr_dndx[i]) / (arr_e[i] - arr_e[i-1]);
1290
1291 return dnde;
1292
1293}