1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // Transition Radiation Detector version 1 -- slow simulator //
24 <img src="picts/AliTRDfullClass.gif">
29 ///////////////////////////////////////////////////////////////////////////////
34 #include <TLorentzVector.h>
38 #include <TVirtualMC.h>
42 #include "AliTRDgeometry.h"
43 #include "AliTRDhit.h"
44 #include "AliTRDsim.h"
50 //_____________________________________________________________________________
51 AliTRDv1::AliTRDv1():AliTRD()
54 // Default constructor
68 fTypeOfStepManager = 2;
72 //_____________________________________________________________________________
73 AliTRDv1::AliTRDv1(const char *name, const char *title)
77 // Standard constructor for Transition Radiation Detector version 1
90 fTypeOfStepManager = 2;
92 SetBufferSize(128000);
96 //_____________________________________________________________________________
97 AliTRDv1::AliTRDv1(const AliTRDv1 &trd):AliTRD(trd)
103 ((AliTRDv1 &) trd).Copy(*this);
107 //_____________________________________________________________________________
108 AliTRDv1::~AliTRDv1()
111 // AliTRDv1 destructor
114 if (fDeltaE) delete fDeltaE;
115 if (fDeltaG) delete fDeltaG;
120 //_____________________________________________________________________________
121 AliTRDv1 &AliTRDv1::operator=(const AliTRDv1 &trd)
124 // Assignment operator
127 if (this != &trd) ((AliTRDv1 &) trd).Copy(*this);
132 //_____________________________________________________________________________
133 void AliTRDv1::Copy(TObject &trd)
139 ((AliTRDv1 &) trd).fSensSelect = fSensSelect;
140 ((AliTRDv1 &) trd).fSensPlane = fSensPlane;
141 ((AliTRDv1 &) trd).fSensChamber = fSensChamber;
142 ((AliTRDv1 &) trd).fSensSector = fSensSector;
143 ((AliTRDv1 &) trd).fSensSectorRange = fSensSectorRange;
145 ((AliTRDv1 &) trd).fTypeOfStepManager = fTypeOfStepManager;
146 ((AliTRDv1 &) trd).fStepSize = fStepSize;
148 fDeltaE->Copy(*((AliTRDv1 &) trd).fDeltaE);
149 fDeltaG->Copy(*((AliTRDv1 &) trd).fDeltaG);
150 fTR->Copy(*((AliTRDv1 &) trd).fTR);
154 //_____________________________________________________________________________
155 void AliTRDv1::CreateGeometry()
158 // Create the GEANT geometry for the Transition Radiation Detector - Version 1
159 // This version covers the full azimuth.
162 // Check that FRAME is there otherwise we have no place where to put the TRD
163 AliModule* frame = gAlice->GetModule("FRAME");
166 // Define the chambers
167 AliTRD::CreateGeometry();
171 //_____________________________________________________________________________
172 void AliTRDv1::CreateMaterials()
175 // Create materials for the Transition Radiation Detector version 1
178 AliTRD::CreateMaterials();
182 //_____________________________________________________________________________
183 void AliTRDv1::CreateTRhit(Int_t det)
186 // Creates an electron cluster from a TR photon.
187 // The photon is assumed to be created a the end of the radiator. The
188 // distance after which it deposits its energy takes into account the
189 // absorbtion of the entrance window and of the gas mixture in drift
194 const Int_t kPdgElectron = 11;
197 const Float_t kWion = 22.04;
199 // Maximum number of TR photons per track
200 const Int_t kNTR = 50;
202 TLorentzVector mom, pos;
204 // Create TR at the entrance of the chamber
205 if (gMC->IsTrackEntering()) {
207 // Create TR only for electrons
208 Int_t iPdg = gMC->TrackPid();
209 if (TMath::Abs(iPdg) != kPdgElectron) return;
215 gMC->TrackMomentum(mom);
216 Float_t pTot = mom.Rho();
217 fTR->CreatePhotons(iPdg,pTot,nTR,eTR);
219 printf("AliTRDv1::CreateTRhit -- ");
220 printf("Boundary error: nTR = %d, kNTR = %d\n",nTR,kNTR);
224 // Loop through the TR photons
225 for (Int_t iTR = 0; iTR < nTR; iTR++) {
227 Float_t energyMeV = eTR[iTR] * 0.001;
228 Float_t energyeV = eTR[iTR] * 1000.0;
229 Float_t absLength = 0;
232 // Take the absorbtion in the entrance window into account
233 Double_t muMy = fTR->GetMuMy(energyMeV);
234 sigma = muMy * fFoilDensity;
236 absLength = gRandom->Exp(1.0/sigma);
237 if (absLength < AliTRDgeometry::MyThick()) continue;
243 // The absorbtion cross sections in the drift gas
245 // Gas-mixture (Xe/CO2)
246 Double_t muXe = fTR->GetMuXe(energyMeV);
247 Double_t muCO = fTR->GetMuCO(energyMeV);
248 sigma = (0.85 * muXe + 0.15 * muCO) * fGasDensity * fTR->GetTemp();
251 // Gas-mixture (Xe/Isobutane)
252 Double_t muXe = fTR->GetMuXe(energyMeV);
253 Double_t muBu = fTR->GetMuBu(energyMeV);
254 sigma = (0.97 * muXe + 0.03 * muBu) * fGasDensity * fTR->GetTemp();
257 // The distance after which the energy of the TR photon
260 absLength = gRandom->Exp(1.0/sigma);
261 if (absLength > (AliTRDgeometry::DrThick()
262 + AliTRDgeometry::AmThick())) {
270 // The position of the absorbtion
272 gMC->TrackPosition(pos);
273 posHit[0] = pos[0] + mom[0] / pTot * absLength;
274 posHit[1] = pos[1] + mom[1] / pTot * absLength;
275 posHit[2] = pos[2] + mom[2] / pTot * absLength;
278 Int_t q = ((Int_t) (energyeV / kWion));
280 // Add the hit to the array. TR photon hits are marked
281 // by negative charge
282 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,posHit,-q,kTRUE);
290 //_____________________________________________________________________________
291 void AliTRDv1::Init()
294 // Initialise Transition Radiation Detector after geometry has been built.
299 if(fDebug) printf("%s: Slow simulator\n",ClassName());
302 printf(" Only plane %d is sensitive\n",fSensPlane);
303 if (fSensChamber >= 0)
304 printf(" Only chamber %d is sensitive\n",fSensChamber);
305 if (fSensSector >= 0) {
306 Int_t sens1 = fSensSector;
307 Int_t sens2 = fSensSector + fSensSectorRange;
308 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
309 * AliTRDgeometry::Nsect();
310 printf(" Only sectors %d - %d are sensitive\n",sens1,sens2-1);
314 printf("%s: TR simulation on\n",ClassName());
316 printf("%s: TR simulation off\n",ClassName());
319 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
320 const Float_t kPoti = 12.1;
321 // Maximum energy (50 keV);
322 const Float_t kEend = 50000.0;
323 // Ermilova distribution for the delta-ray spectrum
324 Float_t poti = TMath::Log(kPoti);
325 Float_t eEnd = TMath::Log(kEend);
327 // Ermilova distribution for the delta-ray spectrum
328 fDeltaE = new TF1("deltae" ,Ermilova ,poti,eEnd,0);
330 // Geant3 distribution for the delta-ray spectrum
331 fDeltaG = new TF1("deltaeg",IntSpecGeant,poti,eEnd,0);
334 printf("%s: ",ClassName());
335 for (Int_t i = 0; i < 80; i++) printf("*");
341 //_____________________________________________________________________________
342 AliTRDsim *AliTRDv1::CreateTR()
345 // Enables the simulation of TR
348 fTR = new AliTRDsim();
353 //_____________________________________________________________________________
354 void AliTRDv1::SetSensPlane(Int_t iplane)
357 // Defines the hit-sensitive plane (0-5)
360 if ((iplane < 0) || (iplane > 5)) {
361 printf("Wrong input value: %d\n",iplane);
362 printf("Use standard setting\n");
373 //_____________________________________________________________________________
374 void AliTRDv1::SetSensChamber(Int_t ichamber)
377 // Defines the hit-sensitive chamber (0-4)
380 if ((ichamber < 0) || (ichamber > 4)) {
381 printf("Wrong input value: %d\n",ichamber);
382 printf("Use standard setting\n");
389 fSensChamber = ichamber;
393 //_____________________________________________________________________________
394 void AliTRDv1::SetSensSector(Int_t isector)
397 // Defines the hit-sensitive sector (0-17)
400 SetSensSector(isector,1);
404 //_____________________________________________________________________________
405 void AliTRDv1::SetSensSector(Int_t isector, Int_t nsector)
408 // Defines a range of hit-sensitive sectors. The range is defined by
409 // <isector> (0-17) as the starting point and <nsector> as the number
410 // of sectors to be included.
413 if ((isector < 0) || (isector > 17)) {
414 printf("Wrong input value <isector>: %d\n",isector);
415 printf("Use standard setting\n");
417 fSensSectorRange = 0;
422 if ((nsector < 1) || (nsector > 18)) {
423 printf("Wrong input value <nsector>: %d\n",nsector);
424 printf("Use standard setting\n");
426 fSensSectorRange = 0;
432 fSensSector = isector;
433 fSensSectorRange = nsector;
437 //_____________________________________________________________________________
438 void AliTRDv1::StepManager()
441 // Slow simulator. Every charged track produces electron cluster as hits
442 // along its path across the drift volume.
445 switch (fTypeOfStepManager) {
446 case 0 : StepManagerErmilova(); break; // 0 is Ermilova
447 case 1 : StepManagerGeant(); break; // 1 is Geant
448 case 2 : StepManagerFixedStep(); break; // 2 is fixed step
449 default : printf("<AliTRDv1::StepManager>: Not a valid Step Manager.\n");
454 //_____________________________________________________________________________
455 void AliTRDv1::SelectStepManager(Int_t t)
458 // Selects a step manager type:
461 // 2 - Fixed step size
465 printf("<AliTRDv1::SelectStepManager>: Sorry, Geant parametrization step"
466 "manager is not implemented yet. Please ask K.Oyama for detail.\n");
469 fTypeOfStepManager = t;
470 printf("<AliTRDv1::SelectStepManager>: Step Manager type %d was selected.\n"
471 , fTypeOfStepManager);
475 //_____________________________________________________________________________
476 void AliTRDv1::StepManagerGeant()
479 // Slow simulator. Every charged track produces electron cluster as hits
480 // along its path across the drift volume. The step size is set acording
481 // to Bethe-Bloch. The energy distribution of the delta electrons follows
482 // a spectrum taken from Geant3.
485 printf("AliTRDv1::StepManagerGeant: Not implemented yet.\n");
489 //_____________________________________________________________________________
490 void AliTRDv1::StepManagerErmilova()
493 // Slow simulator. Every charged track produces electron cluster as hits
494 // along its path across the drift volume. The step size is set acording
495 // to Bethe-Bloch. The energy distribution of the delta electrons follows
496 // a spectrum taken from Ermilova et al.
513 Double_t betaGamma, pp;
516 Bool_t drRegion = kFALSE;
517 Bool_t amRegion = kFALSE;
520 TString cIdSensDr = "J";
521 TString cIdSensAm = "K";
522 Char_t cIdChamber[3];
525 TLorentzVector pos, mom;
527 const Int_t kNplan = AliTRDgeometry::Nplan();
528 const Int_t kNcham = AliTRDgeometry::Ncham();
529 const Int_t kNdetsec = kNplan * kNcham;
531 const Double_t kBig = 1.0E+12; // Infinitely big
532 const Float_t kWion = 22.04; // Ionization energy
533 const Float_t kPTotMaxEl = 0.002; // Maximum momentum for e+ e- g
535 // Minimum energy for the step size adjustment
536 const Float_t kEkinMinStep = 1.0e-5;
538 // Plateau value of the energy-loss for electron in xenon
539 // taken from: Allison + Comb, Ann. Rev. Nucl. Sci. (1980), 30, 253
540 //const Double_t kPlateau = 1.70;
541 // the averaged value (26/3/99)
542 const Float_t kPlateau = 1.55;
544 const Float_t kPrim = 48.0; // dN1/dx|min for the gas mixture (90% Xe + 10% CO2)
545 // First ionization potential (eV) for the gas mixture (90% Xe + 10% CO2)
546 const Float_t kPoti = 12.1;
548 const Int_t kPdgElectron = 11; // PDG code electron
550 // Set the maximum step size to a very large number for all
551 // neutral particles and those outside the driftvolume
552 gMC->SetMaxStep(kBig);
554 // Use only charged tracks
555 if (( gMC->TrackCharge() ) &&
556 (!gMC->IsTrackStop() ) &&
557 (!gMC->IsTrackDisappeared())) {
559 // Inside a sensitive volume?
562 cIdCurrent = gMC->CurrentVolName();
563 if (cIdSensDr == cIdCurrent[1]) {
566 if (cIdSensAm == cIdCurrent[1]) {
569 if (drRegion || amRegion) {
571 // The hit coordinates and charge
572 gMC->TrackPosition(pos);
577 // The sector number (0 - 17)
578 // The numbering goes clockwise and starts at y = 0
579 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
584 sec = ((Int_t) (phi / 20));
586 // The plane and chamber number
587 cIdChamber[0] = cIdCurrent[2];
588 cIdChamber[1] = cIdCurrent[3];
589 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
590 cha = ((Int_t) idChamber / kNplan);
591 pla = ((Int_t) idChamber % kNplan);
593 // Check on selected volumes
594 Int_t addthishit = 1;
596 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
597 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
598 if (fSensSector >= 0) {
599 Int_t sens1 = fSensSector;
600 Int_t sens2 = fSensSector + fSensSectorRange;
601 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
602 * AliTRDgeometry::Nsect();
604 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
607 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
615 // The detector number
616 det = fGeometry->GetDetector(pla,cha,sec);
618 // Special hits only in the drift region
621 // Create a track reference at the entrance and
622 // exit of each chamber that contain the
623 // momentum components of the particle
624 if (gMC->IsTrackEntering() || gMC->IsTrackExiting()) {
625 gMC->TrackMomentum(mom);
626 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
629 // Create the hits from TR photons
630 if (fTR) CreateTRhit(det);
634 // Calculate the energy of the delta-electrons
635 eDelta = TMath::Exp(fDeltaE->GetRandom()) - kPoti;
636 eDelta = TMath::Max(eDelta,0.0);
638 // The number of secondary electrons created
639 qTot = ((Int_t) (eDelta / kWion) + 1);
641 // Create a new dEdx hit
643 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
644 ,det,hits,qTot,kTRUE);
647 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber()
648 ,det,hits,qTot,kFALSE);
651 // Calculate the maximum step size for the next tracking step
652 // Produce only one hit if Ekin is below cutoff
653 aMass = gMC->TrackMass();
654 if ((gMC->Etot() - aMass) > kEkinMinStep) {
656 // The energy loss according to Bethe Bloch
657 iPdg = TMath::Abs(gMC->TrackPid());
658 if ( (iPdg != kPdgElectron) ||
659 ((iPdg == kPdgElectron) && (pTot < kPTotMaxEl))) {
660 gMC->TrackMomentum(mom);
662 betaGamma = pTot / aMass;
663 pp = kPrim * BetheBloch(betaGamma);
664 // Take charge > 1 into account
665 charge = gMC->TrackCharge();
666 if (TMath::Abs(charge) > 1) pp = pp * charge*charge;
668 // Electrons above 20 Mev/c are at the plateau
670 pp = kPrim * kPlateau;
675 gMC->GetRandom()->RndmArray(1, random);
676 while ((random[0] == 1.) || (random[0] == 0.));
677 stepSize = - TMath::Log(random[0]) / pp;
678 gMC->SetMaxStep(stepSize);
691 //_____________________________________________________________________________
692 void AliTRDv1::StepManagerFixedStep()
695 // Slow simulator. Every charged track produces electron cluster as hits
696 // along its path across the drift volume. The step size is fixed in
697 // this version of the step manager.
709 Bool_t drRegion = kFALSE;
710 Bool_t amRegion = kFALSE;
713 TString cIdSensDr = "J";
714 TString cIdSensAm = "K";
715 Char_t cIdChamber[3];
718 TLorentzVector pos, mom;
720 const Int_t kNplan = AliTRDgeometry::Nplan();
721 const Int_t kNcham = AliTRDgeometry::Ncham();
722 const Int_t kNdetsec = kNplan * kNcham;
724 const Double_t kBig = 1.0E+12;
726 const Float_t kWion = 22.04; // Ionization energy
727 const Float_t kEkinMinStep = 1.0e-5; // Minimum energy for the step size adjustment
729 // Set the maximum step size to a very large number for all
730 // neutral particles and those outside the driftvolume
731 gMC->SetMaxStep(kBig);
733 // If not charged track or already stopped or disappeared, just return.
734 if ((!gMC->TrackCharge()) ||
735 gMC->IsTrackStop() ||
736 gMC->IsTrackDisappeared()) return;
738 // Inside a sensitive volume?
739 cIdCurrent = gMC->CurrentVolName();
741 if (cIdSensDr == cIdCurrent[1]) drRegion = kTRUE;
742 if (cIdSensAm == cIdCurrent[1]) amRegion = kTRUE;
744 if ((!drRegion) && (!amRegion)) return;
746 // The hit coordinates and charge
747 gMC->TrackPosition(pos);
752 // The sector number (0 - 17)
753 // The numbering goes clockwise and starts at y = 0
754 Float_t phi = kRaddeg*TMath::ATan2(pos[0],pos[1]);
755 if (phi < 90.) phi += 270.;
757 sec = ((Int_t) (phi / 20.));
759 // The plane and chamber number
760 cIdChamber[0] = cIdCurrent[2];
761 cIdChamber[1] = cIdCurrent[3];
762 Int_t idChamber = (atoi(cIdChamber) % kNdetsec);
763 cha = ((Int_t) idChamber / kNplan);
764 pla = ((Int_t) idChamber % kNplan);
766 // Check on selected volumes
767 Int_t addthishit = 1;
769 if ((fSensPlane >= 0) && (pla != fSensPlane )) addthishit = 0;
770 if ((fSensChamber >= 0) && (cha != fSensChamber)) addthishit = 0;
771 if (fSensSector >= 0) {
772 Int_t sens1 = fSensSector;
773 Int_t sens2 = fSensSector + fSensSectorRange;
774 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) * AliTRDgeometry::Nsect();
776 if ((sec < sens1) || (sec >= sens2)) addthishit = 0;
779 if ((sec < sens1) && (sec >= sens2)) addthishit = 0;
784 if (!addthishit) return;
786 det = fGeometry->GetDetector(pla,cha,sec); // The detector number
788 Int_t trkStat = 0; // 0: InFlight 1:Entering 2:Exiting
790 // Special hits only in the drift region
793 // Create a track reference at the entrance and exit of each
794 // chamber that contain the momentum components of the particle
796 if (gMC->IsTrackEntering()) {
797 gMC->TrackMomentum(mom);
798 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
801 if (gMC->IsTrackExiting()) {
802 gMC->TrackMomentum(mom);
803 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber());
807 // Create the hits from TR photons
808 if (fTR) CreateTRhit(det);
812 // Calculate the charge according to GEANT Edep
813 // Create a new dEdx hit
814 eDep = TMath::Max(gMC->Edep(),0.0) * 1.0e+09;
815 qTot = (Int_t) (eDep / kWion);
816 AddHit(gAlice->GetMCApp()->GetCurrentTrackNumber(),det,hits,qTot,drRegion);
818 // Set Maximum Step Size
819 // Produce only one hit if Ekin is below cutoff
820 if ((gMC->Etot() - gMC->TrackMass()) < kEkinMinStep) return;
821 gMC->SetMaxStep(fStepSize);
825 //_____________________________________________________________________________
826 Double_t AliTRDv1::BetheBloch(Double_t bg)
829 // Parametrization of the Bethe-Bloch-curve
830 // The parametrization is the same as for the TPC and is taken from Lehrhaus.
833 // This parameters have been adjusted to averaged values from GEANT
834 const Double_t kP1 = 7.17960e-02;
835 const Double_t kP2 = 8.54196;
836 const Double_t kP3 = 1.38065e-06;
837 const Double_t kP4 = 5.30972;
838 const Double_t kP5 = 2.83798;
840 // This parameters have been adjusted to Xe-data found in:
841 // Allison & Cobb, Ann. Rev. Nucl. Sci. (1980), 30, 253
842 //const Double_t kP1 = 0.76176E-1;
843 //const Double_t kP2 = 10.632;
844 //const Double_t kP3 = 3.17983E-6;
845 //const Double_t kP4 = 1.8631;
846 //const Double_t kP5 = 1.9479;
848 // Lower cutoff of the Bethe-Bloch-curve to limit step sizes
849 const Double_t kBgMin = 0.8;
850 const Double_t kBBMax = 6.83298;
851 //const Double_t kBgMin = 0.6;
852 //const Double_t kBBMax = 17.2809;
853 //const Double_t kBgMin = 0.4;
854 //const Double_t kBBMax = 82.0;
857 Double_t yy = bg / TMath::Sqrt(1. + bg*bg);
858 Double_t aa = TMath::Power(yy,kP4);
859 Double_t bb = TMath::Power((1./bg),kP5);
860 bb = TMath::Log(kP3 + bb);
861 return ((kP2 - aa - bb)*kP1 / aa);
869 //_____________________________________________________________________________
870 Double_t BetheBlochGeant(Double_t bg)
873 // Return dN/dx (number of primary collisions per centimeter)
874 // for given beta*gamma factor.
876 // Implemented by K.Oyama according to GEANT 3 parametrization shown in
877 // A.Andronic's webpage: http://www-alice.gsi.de/trd/papers/dedx/dedx.html
878 // This must be used as a set with IntSpecGeant.
881 Double_t arr_g[20] = {
882 1.100000, 1.200000, 1.300000, 1.500000,
883 1.800000, 2.000000, 2.500000, 3.000000,
884 4.000000, 7.000000, 10.000000, 20.000000,
885 40.000000, 70.000000, 100.000000, 300.000000,
886 600.000000, 1000.000000, 3000.000000, 10000.000000 };
888 Double_t arr_nc[20] = {
889 75.009056, 45.508083, 35.299252, 27.116327,
890 22.734999, 21.411915, 19.934095, 19.449375,
891 19.344431, 20.185553, 21.027925, 22.912676,
892 24.933352, 26.504053, 27.387468, 29.566597,
893 30.353779, 30.787134, 31.129285, 31.157350 };
895 // betagamma to gamma
896 Double_t g = TMath::Sqrt( 1. + bg*bg );
898 // Find the index just before the point we need.
900 for( i = 0 ; i < 18 ; i++ )
901 if( arr_g[i] < g && arr_g[i+1] > g )
904 // Simple interpolation.
905 Double_t pp = ((arr_nc[i+1] - arr_nc[i]) /
906 (arr_g[i+1]-arr_g[i])) * (g-arr_g[i]) + arr_nc[i];
912 //_____________________________________________________________________________
913 Double_t Ermilova(Double_t *x, Double_t *)
916 // Calculates the delta-ray energy distribution according to Ermilova.
917 // Logarithmic scale !
926 const Int_t kNv = 31;
928 Float_t vxe[kNv] = { 2.3026, 2.9957, 3.4012, 3.6889, 3.9120
929 , 4.0943, 4.2485, 4.3820, 4.4998, 4.6052
930 , 4.7005, 5.0752, 5.2983, 5.7038, 5.9915
931 , 6.2146, 6.5221, 6.9078, 7.3132, 7.6009
932 , 8.0064, 8.5172, 8.6995, 8.9872, 9.2103
933 , 9.4727, 9.9035,10.3735,10.5966,10.8198
936 Float_t vye[kNv] = { 80.0 , 31.0 , 23.3 , 21.1 , 21.0
937 , 20.9 , 20.8 , 20.0 , 16.0 , 11.0
938 , 8.0 , 6.0 , 5.2 , 4.6 , 4.0
939 , 3.5 , 3.0 , 1.4 , 0.67 , 0.44
940 , 0.3 , 0.18 , 0.12 , 0.08 , 0.056
941 , 0.04 , 0.023, 0.015, 0.011, 0.01
950 dpos = energy - vxe[pos2++];
954 if (pos2 > kNv) pos2 = kNv - 1;
957 // Differentiate between the sampling points
958 dnde = (vye[pos1] - vye[pos2]) / (vxe[pos2] - vxe[pos1]);
964 //_____________________________________________________________________________
965 Double_t IntSpecGeant(Double_t *x, Double_t *)
968 // Integrated spectrum from Geant3
971 const Int_t n_pts = 83;
972 Double_t arr_e[n_pts] = {
973 2.421257, 2.483278, 2.534301, 2.592230,
974 2.672067, 2.813299, 3.015059, 3.216819,
975 3.418579, 3.620338, 3.868209, 3.920198,
976 3.978284, 4.063923, 4.186264, 4.308605,
977 4.430946, 4.553288, 4.724261, 4.837736,
978 4.999842, 5.161949, 5.324056, 5.486163,
979 5.679688, 5.752998, 5.857728, 5.962457,
980 6.067185, 6.171914, 6.315653, 6.393674,
981 6.471694, 6.539689, 6.597658, 6.655627,
982 6.710957, 6.763648, 6.816338, 6.876198,
983 6.943227, 7.010257, 7.106285, 7.252151,
984 7.460531, 7.668911, 7.877290, 8.085670,
985 8.302979, 8.353585, 8.413120, 8.483500,
986 8.541030, 8.592857, 8.668865, 8.820485,
987 9.037086, 9.253686, 9.470286, 9.686887,
988 9.930838, 9.994655, 10.085822, 10.176990,
989 10.268158, 10.359325, 10.503614, 10.627565,
990 10.804637, 10.981709, 11.158781, 11.335854,
991 11.593397, 11.781165, 12.049404, 12.317644,
992 12.585884, 12.854123, 14.278421, 16.975889,
993 20.829416, 24.682943, 28.536469
995 Double_t arr_dndx[n_pts] = {
996 19.344431, 18.664679, 18.136106, 17.567745,
997 16.836426, 15.677382, 14.281277, 13.140237,
998 12.207677, 11.445510, 10.697049, 10.562296,
999 10.414673, 10.182341, 9.775256, 9.172330,
1000 8.240271, 6.898587, 4.808303, 3.889751,
1001 3.345288, 3.093431, 2.897347, 2.692470,
1002 2.436222, 2.340029, 2.208579, 2.086489,
1003 1.975535, 1.876519, 1.759626, 1.705024,
1004 1.656374, 1.502638, 1.330566, 1.200697,
1005 1.101168, 1.019323, 0.943867, 0.851951,
1006 0.755229, 0.671576, 0.570675, 0.449672,
1007 0.326722, 0.244225, 0.188225, 0.149608,
1008 0.121529, 0.116289, 0.110636, 0.103490,
1009 0.096147, 0.089191, 0.079780, 0.063927,
1010 0.047642, 0.036341, 0.028250, 0.022285,
1011 0.017291, 0.016211, 0.014802, 0.013533,
1012 0.012388, 0.011352, 0.009803, 0.008537,
1013 0.007039, 0.005829, 0.004843, 0.004034,
1014 0.003101, 0.002564, 0.001956, 0.001494,
1015 0.001142, 0.000873, 0.000210, 0.000014,
1016 0.000000, 0.000000, 0.000000
1020 Double_t energy = x[0];
1023 for( i = 0 ; i < n_pts ; i++ )
1024 if( energy < arr_e[i] ) break;
1027 printf("Error in AliTRDv1::IntSpecGeant: "
1028 "given energy value is too small or zero.\n");
1031 dnde = (arr_dndx[i-1] - arr_dndx[i]) / (arr_e[i] - arr_e[i-1]);