Skip to content

Commit cefb94d

Browse files
committed
Clamp NIEL damage weights to the tabulated energy range
This fixes unbounded 1 MeV neutron-equivalent weights outside the RD50 tables. - GetWeight evaluated the spline beyond the last table point, giving a pion weight of 541 at 10 GeV and 5e11 at 1 TeV. - The energy is now clamped to the table range, so weights beyond it use the value at the table edge. - Graphs read from CSV are sorted, since the clamp takes the range from the first and last point. - An empty graph returns 0, as before.
1 parent d7aa492 commit cefb94d

1 file changed

Lines changed: 20 additions & 3 deletions

File tree

‎Common/SimConfig/src/FluenceWeightCalculator.cxx‎

Lines changed: 20 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,7 @@
1111

1212
#include "SimConfig/FluenceWeightCalculator.h"
1313
#include <TFile.h>
14+
#include <algorithm>
1415
#include <fstream>
1516
#include <sstream>
1617
#include <iostream>
@@ -19,6 +20,19 @@ std::unique_ptr<TGraph> FluenceWeightCalculator::neutronG;
1920
std::unique_ptr<TGraph> FluenceWeightCalculator::protonG;
2021
std::unique_ptr<TGraph> FluenceWeightCalculator::pionG;
2122

23+
namespace
24+
{
25+
// Damage weight at an energy clamped to the tabulated range
26+
double evalClamped(const TGraph& g, double kineticEnergy)
27+
{
28+
if (g.GetN() == 0) {
29+
return 0.;
30+
}
31+
const double e = std::clamp(kineticEnergy, g.GetX()[0], g.GetX()[g.GetN() - 1]);
32+
return g.Eval(e, nullptr, "S");
33+
}
34+
} // namespace
35+
2236
double FluenceWeightCalculator::GetWeight(const int pdg, const double kineticEnergy)
2337
{
2438
//
@@ -29,13 +43,13 @@ double FluenceWeightCalculator::GetWeight(const int pdg, const double kineticEne
2943
}
3044
switch (std::abs(pdg)) {
3145
case 2112: {
32-
return neutronG->Eval(kineticEnergy, nullptr, "S");
46+
return evalClamped(*neutronG, kineticEnergy);
3347
}
3448
case 2212: {
35-
return ((kineticEnergy > 1e-3) ? protonG->Eval(kineticEnergy, nullptr, "S") : 0.);
49+
return ((kineticEnergy > 1e-3) ? evalClamped(*protonG, kineticEnergy) : 0.);
3650
}
3751
case 211: {
38-
return ((kineticEnergy > 10.) ? pionG->Eval(kineticEnergy, nullptr, "S") : 0.);
52+
return ((kineticEnergy > 10.) ? evalClamped(*pionG, kineticEnergy) : 0.);
3953
}
4054
default:
4155
return 0.0;
@@ -130,6 +144,9 @@ void FluenceWeightCalculator::InitWeightsFromCSV(const std::string& filename)
130144
default:;
131145
}
132146
}
147+
neutronG->Sort();
148+
protonG->Sort();
149+
pionG->Sort();
133150
auto fout = new TFile("rd50_niel.root", "recreate");
134151
neutronG->Write();
135152
protonG->Write();

0 commit comments

Comments
 (0)