Weight of each event with wigner distribution corrected to account fo…

…r sign of maxweight
AchillesGen · Nov 22, 2024 · 55900ff · 55900ff
1 parent 0332034
commit 55900ff
Show file tree

Hide file tree

Showing 3 changed files with 11 additions and 10 deletions.
diff --git a/include/Achilles/WignerDistribution.hh b/include/Achilles/WignerDistribution.hh
@@ -20,7 +20,7 @@ class WignerDistribution {
     std::vector<double> &Radius() { return radius; }
     double MinMomentum() const { return mom.front(); }
     double MaxMomentum() const { return mom.back(); }
-    double MaxWeight(double r) const;
+    double MaxAbsWeightSign(double r) const;
     double MaxAbsWeight(double r) const;
     double MinRadius() const { return radius.front(); }
     double MaxRadius() const { return radius.back(); }

diff --git a/src/Achilles/Nucleus.cc b/src/Achilles/Nucleus.cc
@@ -220,8 +220,8 @@ double Nucleus::FermiGasWeight(const Particle &p) const {
     switch(fermi_gas.type) {
     case FermiGasType::Wigner: {
         auto position = p.Position().Magnitude();
-        auto max_wigner_value = wigner_d.MaxWeight(position);
-        return std::copysign(1.0, wigner_d(position, p.Momentum().P()) / max_wigner_value);
+        auto max_wigner_sign= wigner_d.MaxAbsWeightSign(position);
+        return std::copysign(1.0, wigner_d(position, p.Momentum().P())/max_wigner_sign);
     }
     case FermiGasType::Global:
     case FermiGasType::Local:

diff --git a/src/Achilles/WignerDistribution.cc b/src/Achilles/WignerDistribution.cc
@@ -85,16 +85,17 @@ double WignerDistribution::operator()(double r) const {
     return result > 0 ? result : 0;
 }
 
-double WignerDistribution::MaxWeight(double r) const {
-    auto wgt_func = [&](double p) { return -p * p * this->operator()(r, p); };
-    Brent brent(wgt_func);
-    auto maxmomweight = brent.Minimize(mom.front(), 200, mom.back());
-    return -wgt_func(maxmomweight);
+double WignerDistribution::MaxAbsWeightSign(double r) const {
+    auto absfunc = [&](double p) { return -p * p * std::abs(this->operator()(r, p)); };
+    Brent brent(absfunc);
+    auto maxweight_mom = brent.Minimize(mom.front(), 200, mom.back());
+    return std::copysign(1.0, this->operator()(r, maxweight_mom));
 }
 
 double WignerDistribution::MaxAbsWeight(double r) const {
     auto absfunc = [&](double p) { return -p * p * std::abs(this->operator()(r, p)); };
     Brent brent(absfunc);
-    auto maxmomweight = brent.Minimize(mom.front(), 200, mom.back());
-    return -absfunc(maxmomweight);
+    auto maxweight_mom = brent.Minimize(mom.front(), 200, mom.back());
+
+    return -absfunc(maxweight_mom);
 }