diff --git a/examples/WIMP.rml b/examples/WIMP.rml
index 03a0bc7..c944a6a 100644
--- a/examples/WIMP.rml
+++ b/examples/WIMP.rml
@@ -9,7 +9,7 @@
-
+
diff --git a/inc/TRestWimpSensitivity.h b/inc/TRestWimpSensitivity.h
index e7e7479..a800f0e 100644
--- a/inc/TRestWimpSensitivity.h
+++ b/inc/TRestWimpSensitivity.h
@@ -70,11 +70,35 @@ class TRestWimpSensitivity : public TRestMetadata {
void ReadNuclei();
const Double_t GetSensitivity(const double wimpMass);
void CalculateQuenchingFactor();
+ bool isEnergySpectraWideEnough();
const std::string BuildOutputFileName(const std::string& extension = ".txt");
std::map GetRecoilSpectra(const double wimpMass, const double crossSection);
std::map GetFormFactor();
inline auto GetQuenchingFactor() { return quenchingFactor; }
+ inline std::vector& GetNuclei() { return fNuclei; };
+ inline auto GetEnergySpectra() { return fEnergySpectra; }
+ inline auto GetEnergySpectraStep() { return fEnergySpectraStep; }
+ inline auto GetEnergyRange() { return fEnergyRange; }
+ inline auto GetWimpDensity() { return fWimpDensity; }
+ inline auto GetLabVelocity() { return fLabVelocity; }
+ inline auto GetEscapeVelocity() { return fEscapeVelocity; }
+ inline auto GetRmsVelocity() { return fRmsVelocity; }
+ inline auto GetExposure() { return fExposure; }
+ inline auto GetBackground() { return fBackground; }
+ inline auto GetUseQuenchingFactor() { return fUseQuenchingFactor; }
+
+ void SetNuclei(const std::vector& nuclei) { fNuclei = nuclei; }
+ void SetEnergySpectra(const TVector2& energySpectra) { fEnergySpectra = energySpectra; }
+ void SetEnergySpectraStep(const Double_t energySpectraStep) { fEnergySpectraStep = energySpectraStep; }
+ void SetEnergyRange(const TVector2& energyRange) { fEnergyRange = energyRange; }
+ void SetWimpDensity(const Double_t wimpDensity) { fWimpDensity = wimpDensity; }
+ void SetLabVelocity(const Double_t labVelocity) { fLabVelocity = labVelocity; }
+ void SetEscapeVelocity(const Double_t escapeVelocity) { fEscapeVelocity = escapeVelocity; }
+ void SetRmsVelocity(const Double_t rmsVelocity) { fRmsVelocity = rmsVelocity; }
+ void SetExposure(const Double_t exposure) { fExposure = exposure; }
+ void SetBackground(const Double_t background) { fBackground = background; }
+ void SetUseQuenchingFactor(const Bool_t useQuenchingFactor) { fUseQuenchingFactor = useQuenchingFactor; }
ClassDefOverride(TRestWimpSensitivity, 1);
};
diff --git a/inc/TRestWimpUtils.h b/inc/TRestWimpUtils.h
index baa5949..38f8030 100644
--- a/inc/TRestWimpUtils.h
+++ b/inc/TRestWimpUtils.h
@@ -31,7 +31,7 @@ namespace TRestWimpUtils {
/// Physics constants
constexpr double GEV_PER_UMA = 0.93149432;
constexpr double HC_KEV_FM = 197327.053;
-constexpr double LIGHT_SPEED = 300000.0; // km/s
+constexpr double LIGHT_SPEED = 299792.458; // km/s
constexpr double SECONDS_PER_DAY = 86400;
constexpr double N_AVOGADRO = 6.0221367E23;
constexpr double MBARN_PER_GEVM2 = 0.38937966;
@@ -46,6 +46,8 @@ const double Bessel(const double x);
const double GetVMin(const double wimpMass, const double Anum, const double recoilEnergy);
const double GetVelocityDistribution(const double v, const double vLab, const double vRMS,
const double vEscape);
+const double GetDifferentialCrossSectionNoHelmFormFactor(const double wimpMass, const double crossSection,
+ const double velocity, const double Anum);
const double GetDifferentialCrossSection(const double wimpMass, const double crossSection,
const double velocity, const double recoilEnergy, const double Anum);
const double GetRecoilRate(const double wimpMass, const double crossSection, const double recoilEnergy,
diff --git a/src/TRestWimpSensitivity.cxx b/src/TRestWimpSensitivity.cxx
index d60cad9..7987e6c 100644
--- a/src/TRestWimpSensitivity.cxx
+++ b/src/TRestWimpSensitivity.cxx
@@ -61,7 +61,7 @@
///
///
///
-///
+///
///
///
///
@@ -177,6 +177,7 @@ void TRestWimpSensitivity::ReadNuclei() {
///////////////////////////////////////////////
/// \brief Get recoil spectra for a given WIMP
/// mass and cross section
+/// Better performance version
///
std::map TRestWimpSensitivity::GetRecoilSpectra(const double wimpMass,
const double crossSection) {
@@ -188,16 +189,76 @@ std::map TRestWimpSensitivity::GetRecoilSpectra(const double
std::string histName = "RecoilSpc_" + std::string(nucl.fNucleusName);
TH1D* recoilSpc =
new TH1D(histName.c_str(), histName.c_str(), nBins, fEnergySpectra.X(), fEnergySpectra.Y());
- for (int i = 1; i < recoilSpc->GetNbinsX(); i++) {
+
+ // Build vector of tuples=(recoilEnergy, minimum velocity, rate) used in further calculations
+ std::vector> tEnergyVminRate;
+ for (int i = 0; i < recoilSpc->GetNbinsX(); i++) {
+ double E = recoilSpc->GetBinCenter(i);
+ if (E <= 0) continue;
+ tEnergyVminRate.push_back(
+ std::make_tuple(E, TRestWimpUtils::GetVMin(wimpMass, nucl.fAnum, E), 0));
+ }
+
+ const double nNuclei =
+ nucl.fAbundance * TRestWimpUtils::N_AVOGADRO * 1E3 / nucl.fAnum; // Number of atoms
+ const double vMin = std::get<1>(
+ tEnergyVminRate.at(0)); // element 0 should be the lowest (positive) energy -> lowest vMin
+ const double vMax = fEscapeVelocity + fLabVelocity;
+
+ // calculation of the rate for each recoil energy
+ double rate{0}; // will contain integral from minimun vMin to vMax, idem integral_min(vMin)^vMax
+ const double velStep = 0.1; // km/s
+ int j = 0;
+ double flux = 0, diffRate = 0, v = 0;
+ // vMax+velStep to save the rate when v is in interval (vMax-velStep, vMax]
+ for (v = vMin; v < vMax + velStep; v += velStep) {
+ // save (in 3rd element of tEnergyVminRate tuples) the integral from minimun vMin to each vMin,
+ // idem integral_min(vMin)^vMin
+ while (j < (int)tEnergyVminRate.size()) {
+ const double vmin = std::get<1>(tEnergyVminRate.at(j));
+ if (vmin < v) {
+ // std::get<2>(tEnergyVminRate.at(j)) = rate; //les precise
+ std::get<2>(tEnergyVminRate.at(j)) = rate - diffRate * (v - vmin); // more precise
+ j++;
+ } else
+ break;
+ }
+ flux = 1E5 * v * fWimpDensity / wimpMass;
+ diffRate =
+ flux *
+ TRestWimpUtils::GetDifferentialCrossSectionNoHelmFormFactor(wimpMass, crossSection, v,
+ nucl.fAnum) *
+ TRestWimpUtils::GetVelocityDistribution(v, fLabVelocity, fRmsVelocity, fEscapeVelocity);
+ rate += diffRate * velStep;
+ }
+ rate -=
+ diffRate * (v - vMax); // substract last diffRate*(v - vMax) to obtain the rate from vMin to vMax
+
+ /*obtain the rate (integral from each vMin to vMax) by substracting integral from minimun vMin to each
+ vMin to the integral from minimun vMin to vMax
+ idem: integral_vMin^vMax = integral_min(vMin)^vMax - integral_min(vMin)^vMin */
+ for (auto& [E, vmin, r] : tEnergyVminRate) {
+ if (vmin > vMax) continue; // r=0
+ const double formFactor = TRestWimpUtils::GetHelmFormFactor(E, nucl.fAnum);
+ r = (rate - r) * formFactor * formFactor * TRestWimpUtils::SECONDS_PER_DAY * nNuclei;
+ }
+
+ // copy results to recoilMap
+ j = 0;
+ for (int i = 0; i < recoilSpc->GetNbinsX(); i++) {
const double recoilEnergy = recoilSpc->GetBinCenter(i);
- const double recoilRate =
- TRestWimpUtils::GetRecoilRate(wimpMass, crossSection, recoilEnergy, nucl.fAnum, fLabVelocity,
- fRmsVelocity, fEscapeVelocity, fWimpDensity, nucl.fAbundance);
- recoilSpc->SetBinContent(i, recoilRate);
+ // const double recoilRate = std::get<2> (tEnergyVminRate.at(i));
+ while (j < (int)tEnergyVminRate.size()) {
+ if (recoilEnergy == std::get<0>(tEnergyVminRate.at(j))) {
+ recoilSpc->SetBinContent(i, std::get<2>(tEnergyVminRate.at(j)));
+ j++;
+ } else
+ break;
+ }
}
+
recoilMap[std::string(nucl.fNucleusName)] = recoilSpc;
}
-
return recoilMap;
}
@@ -211,19 +272,48 @@ const Double_t TRestWimpSensitivity::GetSensitivity(const double wimpMass) {
if (fUseQuenchingFactor) CalculateQuenchingFactor();
- const int nBins = (fEnergySpectra.Y() - fEnergySpectra.X()) / fEnergySpectraStep;
+ if (!isEnergySpectraWideEnough()) {
+ RESTError << "Energy spectra range is not wide enough to match the energy range given." << RESTendl;
+ // return 0;
+ }
+
+ double nMeas = 0;
- TH1D recoilSpc("recoilSpc", "recoilSpc", nBins, fEnergySpectra.X(), fEnergySpectra.Y());
+ const double crossSection = 1E-45;
+ auto rSpc = GetRecoilSpectra(wimpMass, crossSection);
+
+ for (auto& nucl : fNuclei) {
+ auto recoilSpc = rSpc[std::string(nucl.fNucleusName)];
+
+ for (int i = 1; i < recoilSpc->GetNbinsX(); i++) {
+ double recoilEnergy = recoilSpc->GetBinCenter(i);
+ const double recoilRate = recoilSpc->GetBinContent(i);
+
+ if (fUseQuenchingFactor)
+ recoilEnergy *= quenchingFactor[std::string(nucl.fNucleusName)]->GetBinContent(i);
+
+ if (recoilEnergy < fEnergyRange.X() || recoilEnergy > fEnergyRange.Y()) continue;
+ nMeas += recoilRate * fEnergySpectraStep;
+ }
+ }
double bckCounts = 0;
- for (int i = 1; i < recoilSpc.GetNbinsX(); i++) {
- const double en = recoilSpc.GetBinCenter(i);
+ auto recSpc = rSpc[std::string(fNuclei.front().fNucleusName)];
+ for (int i = 1; i < recSpc->GetNbinsX(); i++) {
+ const double en = recSpc->GetBinCenter(i);
if (en < fEnergyRange.X() || en > fEnergyRange.Y()) continue;
bckCounts += fBackground * fEnergySpectraStep;
}
bckCounts *= fExposure;
+ for (auto& [name, histo] : rSpc) delete histo;
+ rSpc.clear();
+
+ RESTExtreme << "nMeas = " << nMeas << " c/kg/day" << RESTendl;
+ RESTExtreme << "bckCounts = " << bckCounts << RESTendl;
+ if (nMeas == 0) return 0;
+
double signalCounts = 0, prob = 0;
do {
@@ -237,35 +327,10 @@ const Double_t TRestWimpSensitivity::GetSensitivity(const double wimpMass) {
signalCounts++;
} while (fabs(prob - 0.1) > 0.01 && signalCounts < 1E6);
- double nMeas = 0;
- const double crossSection = 1E-45;
-
- for (auto& nucl : fNuclei) {
- recoilSpc.Reset();
-
- for (int i = 1; i < recoilSpc.GetNbinsX(); i++) {
- const double recoilEnergy = recoilSpc.GetBinCenter(i);
- const double recoilRate =
- TRestWimpUtils::GetRecoilRate(wimpMass, crossSection, recoilEnergy, nucl.fAnum, fLabVelocity,
- fRmsVelocity, fEscapeVelocity, fWimpDensity, nucl.fAbundance);
- recoilSpc.SetBinContent(i, recoilRate);
- }
-
- for (int i = 1; i < recoilSpc.GetNbinsX(); i++) {
- double recoilEnergy = recoilSpc.GetBinCenter(i);
- // const double recoilRate = recoilSpc.GetBinContent(i);
- if (fUseQuenchingFactor)
- recoilEnergy *= quenchingFactor[std::string(nucl.fNucleusName)]->GetBinContent(i);
-
- if (recoilEnergy < fEnergyRange.X() || recoilEnergy > fEnergyRange.Y()) continue;
- nMeas += recoilSpc.GetBinContent(i) * fEnergySpectraStep;
- }
- }
-
- if (nMeas == 0) return 0;
-
const double sensitivity = signalCounts * 1E-45 / (nMeas * fExposure);
+ RESTExtreme << "sigCounts = " << signalCounts << RESTendl;
+
return sensitivity;
}
@@ -274,7 +339,17 @@ const Double_t TRestWimpSensitivity::GetSensitivity(const double wimpMass) {
/// stores in a map
///
void TRestWimpSensitivity::CalculateQuenchingFactor() {
- if (!quenchingFactor.empty()) return;
+ // do not calculate if already calculated (with same energy spectra limits)
+ if (!quenchingFactor.empty()) {
+ bool same = true;
+ for (auto& [name, histo] : quenchingFactor)
+ if (histo->GetXaxis()->GetXmin() != fEnergySpectra.X() ||
+ histo->GetXaxis()->GetXmax() != fEnergySpectra.Y()) {
+ same = false;
+ break;
+ }
+ if (same) return;
+ }
std::cout << "Calculating quenching factor " << std::endl;
@@ -296,6 +371,21 @@ void TRestWimpSensitivity::CalculateQuenchingFactor() {
}
}
+bool TRestWimpSensitivity::isEnergySpectraWideEnough() {
+ if (!fUseQuenchingFactor)
+ return fEnergySpectra.X() <= fEnergyRange.X() && fEnergySpectra.Y() >= fEnergyRange.Y();
+
+ CalculateQuenchingFactor();
+ for (auto& nucl : fNuclei) {
+ auto qf = quenchingFactor[std::string(nucl.fNucleusName)];
+ // assuming that Energy_nr * QF(Energy_nr) is a monotonically increasing function
+ if (qf->GetBinContent(1) * qf->GetBinCenter(1) > fEnergyRange.X() ||
+ qf->GetBinContent(qf->GetNbinsX() - 1) * qf->GetBinCenter(qf->GetNbinsX() - 1) < fEnergyRange.Y())
+ return false;
+ }
+ return true;
+}
+
///////////////////////////////////////////////
/// \brief Return output file format with different
/// parameters used in the calculation.
@@ -309,6 +399,7 @@ const std::string TRestWimpSensitivity::BuildOutputFileName(const std::string& e
ss << "WD_" << fWimpDensity << "_";
ss << "Vel_" << fLabVelocity << "_" << fRmsVelocity << "_" << fEscapeVelocity << "_";
ss << "Bck_" << fBackground << "_";
+ ss << "Exp_" << fExposure << "_";
ss << "RecEn_" << fEnergySpectra.X() << "_" << fEnergySpectra.Y() << "_" << fEnergySpectraStep << "_";
ss << "EnRange_" << fEnergyRange.X() << "_" << fEnergyRange.Y() << "_";
@@ -366,6 +457,6 @@ void TRestWimpSensitivity::PrintMetadata() {
<< ") Step: " << fEnergySpectraStep << " keV" << RESTendl;
RESTMetadata << "Sensitivity energy range: (" << fEnergyRange.X() << ", " << fEnergyRange.Y() << ") keV"
<< RESTendl;
- RESTMetadata << "Use quenching factor: " << fUseQuenchingFactor << RESTendl;
+ RESTMetadata << "Use quenching factor: " << (fUseQuenchingFactor ? "true" : "false") << RESTendl;
RESTMetadata << "+++++" << RESTendl;
}
diff --git a/src/TRestWimpUtils.cxx b/src/TRestWimpUtils.cxx
index 68447d4..816c95b 100644
--- a/src/TRestWimpUtils.cxx
+++ b/src/TRestWimpUtils.cxx
@@ -23,7 +23,7 @@
//////////////////////////////////////////////////////////////////////////
/// TRestWimpUtils defines several functions to calculate different WIMP
/// parameters, based on WimpRoot code written by I.G. Irastorza, some
-/// functions are derived from a D. Diez-Ibañez phython code.
+/// functions are derived from a D. Diez-Ibañez python code.
///
/// Note units are keV for the recoil energy, GeV for the WIMP mass,
/// km/s for the WIMP velocity and cm-2 for the cross section.
@@ -53,7 +53,7 @@
const double TRestWimpUtils::GetRelativeNuclearCS(const double wimpMass, const double Anum) {
const double reducedMass = GetReducedMass(wimpMass, Anum);
const double reducedMassSingle = GetReducedMass(wimpMass, 1.);
- return pow(Anum * GEV_PER_UMA * reducedMass / reducedMassSingle, 2.);
+ return (Anum * reducedMass / reducedMassSingle) * (Anum * reducedMass / reducedMassSingle);
}
//////////////////////////////////////////////////
@@ -99,35 +99,61 @@ const double TRestWimpUtils::GetVMin(const double wimpMass, const double Anum, c
//////////////////////////////////////////////////
/// \brief Get velocity distribution for a given
-/// WIMP velocity
+/// Velocity distribution in Earth frame, f(velocity) in units of 1/(km/s), velocity in km/s
+/// Result is the integral for all solid angle of the Boltzmann velocity distribution in Earth frame*/
///
const double TRestWimpUtils::GetVelocityDistribution(const double v, const double vLab, const double vRMS,
const double vEscape) {
- const double vAdim = vRMS / vLab;
- const double Nesc = erf(vAdim) - (2. / sqrt(TMath::Pi())) * (vAdim)*exp(-vAdim * vAdim);
- const double xMax = std::min(1., (vEscape * vEscape - vLab * vLab - v * v) / (2. * vLab * v));
- return v * Nesc / (vLab * vRMS * sqrt(TMath::Pi())) *
+ if (v > vLab + vEscape) return 0;
+
+ const double vAdim = vEscape / vRMS;
+ const double Nesc =
+ erf(vAdim) - 2. / sqrt(TMath::Pi()) * vAdim *
+ exp(-vAdim * vAdim); // Nesc=1 for vEscape=infinity (see Lewin&Smith appendix 1a)
+ // xMax = max(cosTheta) to meet [vec(v) + vec(vLab)]^2 < vEscape^2 boundary condition (also xMax in
+ // [-1,+1])
+ const double xMax =
+ std::max(-1., std::min(1., (vEscape * vEscape - vLab * vLab - v * v) / (2. * vLab * v)));
+
+ return v / Nesc / (vLab * vRMS * sqrt(TMath::Pi())) *
(exp(-(v - vLab) * (v - vLab) / (vRMS * vRMS)) -
exp(-(v * v + vLab * vLab + 2 * v * vLab * xMax) / (vRMS * vRMS)));
}
//////////////////////////////////////////////////
-/// \brief Differential cross section in energy,
+/// \brief Differential cross section without Helm form factor in energy,
/// in [cm/keV]. E in keV, velocity in km/s,
/// wimp mass in Gev/c^2, cross section per
/// nucleon in cm^2, Anum in atomic units (amu)
/// (or atomic weight)
+/// Useful function for performance enhancement
///
-const double TRestWimpUtils::GetDifferentialCrossSection(const double wimpMass, const double crossSection,
- const double velocity, const double recoilEnergy,
- const double Anum) {
+const double TRestWimpUtils::GetDifferentialCrossSectionNoHelmFormFactor(const double wimpMass,
+ const double crossSection,
+ const double velocity,
+ const double Anum) {
const double cs = GetRelativeNuclearCS(wimpMass, Anum) * crossSection;
const double reducedMass = GetReducedMass(wimpMass, Anum);
const double Emax = 1E6 / LIGHT_SPEED / LIGHT_SPEED * 2. * reducedMass * reducedMass * velocity *
velocity / (Anum * GEV_PER_UMA);
+
+ return cs / Emax;
+}
+
+//////////////////////////////////////////////////
+/// \brief Differential cross section in energy,
+/// in [cm/keV]. E in keV, velocity in km/s,
+/// wimp mass in Gev/c^2, cross section per
+/// nucleon in cm^2, Anum in atomic units (amu)
+/// (or atomic weight)
+///
+const double TRestWimpUtils::GetDifferentialCrossSection(const double wimpMass, const double crossSection,
+ const double velocity, const double recoilEnergy,
+ const double Anum) {
const double formFactor = GetHelmFormFactor(recoilEnergy, Anum);
- return cs * formFactor * formFactor / Emax;
+ return GetDifferentialCrossSectionNoHelmFormFactor(wimpMass, crossSection, velocity, Anum) * formFactor *
+ formFactor;
}
//////////////////////////////////////////////////
diff --git a/test/files/wimp.rml b/test/files/wimp.rml
index 03a0bc7..abdf50e 100644
--- a/test/files/wimp.rml
+++ b/test/files/wimp.rml
@@ -9,8 +9,8 @@
-
-
+
+
diff --git a/test/src/TRestWimpSensitivity.cxx b/test/src/TRestWimpSensitivity.cxx
index a49642f..aa3d305 100644
--- a/test/src/TRestWimpSensitivity.cxx
+++ b/test/src/TRestWimpSensitivity.cxx
@@ -30,5 +30,6 @@ TEST(TRestWimpSensitivity, FromRml) {
WS.PrintMetadata();
- EXPECT_TRUE(abs(WS.GetSensitivity(1) - 2.9732007e-39) < 1E-45);
+ const double val = WS.GetSensitivity(1);
+ EXPECT_NEAR(val, 9.764469e-40, 1E-46);
}