diff --git a/include/Achilles/WignerDistribution.hh b/include/Achilles/WignerDistribution.hh index 6e6f7ee1..1e0e9f83 100644 --- a/include/Achilles/WignerDistribution.hh +++ b/include/Achilles/WignerDistribution.hh @@ -20,7 +20,7 @@ class WignerDistribution { std::vector &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 index 1a7d6fc4..08dd353f 100644 --- 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 index ed9b7a5e..1e16581c 100644 --- 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); }