Update HiggsMakeDecay_tool.cc

Fixed a bug with setting the boolean for the new form factor

Author: jberger7

sbncode/EventGenerator/MeVPrtl/Tools/Higgs/HiggsMakeDecay_tool.cc

@@ -73,7 +73,6 @@ class HiggsMakeDecay : public IMeVPrtlDecay {
   double fReferenceRayDistance;
   double fReferenceHiggsMass;
   double fReferenceHiggsMixing;
-  bool fNewFormFactor;
   double fReferenceHiggsEnergy;
   double fReferenceHiggsKaonEnergy;
 
@@ -83,6 +82,7 @@ class HiggsMakeDecay : public IMeVPrtlDecay {
   bool fAllowMuonDecay;
   bool fAllowPionDecay;
   bool fAllowPi0Decay;
+  bool fNewFormFactor;
   bool fAddTimeOfFlight;
 
 };
@@ -128,34 +128,34 @@ double MuonPartialWidth(double higs_mass, double mixing) {
   return LeptonPartialWidth(Constants::Instance().muon_mass, higs_mass, mixing);
 }
 
-double PionPartialWidth(double pion_mass, double higs_mass, double mixing) {
+double PionPartialWidth(double pion_mass, double higs_mass, double mixing, bool use_new_ff) {
   if (pion_mass * 2 >= higs_mass) return 0.;
 
   double higgs_vev = Constants::Instance().higgs_vev;
 
   double form_factor = 0.;
-  if (fNewFormFactor) {
+  if (use_new_ff) {
     // New, improved form factor based on the plot in arXiv:1909.11670v4
     // This form factor is significnatly different than in older versions
     // The expression is a fit to the Fig. 1 left panel in arXiv:1909.11670v4
     form_factor = pow(0.289 * pow(higs_mass - 2 * Constants::Instance().pizero_mass,1.5),0.5);
   }
   else {
     // Old form factor
-    (2. / 9.) * higs_mass * higs_mass + (11. / 9.) * pion_mass * pion_mass;
+    form_factor = (2. / 9.) * higs_mass * higs_mass + (11. / 9.) * pion_mass * pion_mass;
   }
 
   double width = (mixing * mixing * 3 * form_factor * form_factor / (32 * M_PI * higgs_vev * higgs_vev * higs_mass)) * pow(1- 4. * pion_mass * pion_mass / (higs_mass * higs_mass), 1. / 2.);
 
   return width;
 }
 
-double PiPlusPartialWidth(double higs_mass, double mixing) {
-  return 2*PionPartialWidth(Constants::Instance().piplus_mass, higs_mass, mixing);
+double PiPlusPartialWidth(double higs_mass, double mixing, bool use_new_ff) {
+  return 2*PionPartialWidth(Constants::Instance().piplus_mass, higs_mass, mixing, use_new_ff);
 }
 
-double PiZeroPartialWidth(double higs_mass, double mixing) {
-  return PionPartialWidth(Constants::Instance().pizero_mass, higs_mass, mixing);
+double PiZeroPartialWidth(double higs_mass, double mixing, bool use_new_ff) {
+  return PionPartialWidth(Constants::Instance().pizero_mass, higs_mass, mixing, use_new_ff);
 }
 
 HiggsMakeDecay::HiggsMakeDecay(fhicl::ParameterSet const &pset):
@@ -177,14 +177,14 @@ void HiggsMakeDecay::configure(fhicl::ParameterSet const &pset)
   fReferenceRayDistance = pset.get<double>("ReferenceRayDistance", 0.);
   fReferenceHiggsMass = pset.get<double>("ReferenceHiggsMass", -1);
   fReferenceHiggsMixing = pset.get<double>("ReferenceHiggsMixing", -1);
-  fNewFormFactor = pset.get<bool>("NewFormFactor", -1);
   fReferenceHiggsEnergy = pset.get<double>("ReferenceHiggsEnergy", -1);
   fReferenceHiggsKaonEnergy = pset.get<double>("ReferenceHiggsEnergyFromKaonEnergy", -1.);
 
   fAllowElectronDecay = pset.get<bool>("AllowElectronDecay", true);
   fAllowMuonDecay = pset.get<bool>("AllowMuonDecay", true);
   fAllowPionDecay = pset.get<bool>("AllowPionDecay", true);
   fAllowPi0Decay = pset.get<bool>("AllowPi0Decay", true);
+  fNewFormFactor = pset.get<bool>("NewFormFactor", true);
   fAddTimeOfFlight = pset.get<bool>("AddTimeOfFlight", true);
 
   if (fReferenceHiggsEnergy < 0. && fReferenceHiggsKaonEnergy > 0.) {
@@ -197,8 +197,8 @@ void HiggsMakeDecay::configure(fhicl::ParameterSet const &pset)
     // Get each partial width
     double width_elec = ElectronPartialWidth(fReferenceHiggsMass, fReferenceHiggsMixing);
     double width_muon = MuonPartialWidth(fReferenceHiggsMass, fReferenceHiggsMixing);
-    double width_piplus = PiPlusPartialWidth(fReferenceHiggsMass, fReferenceHiggsMixing);
-    double width_pizero = PiZeroPartialWidth(fReferenceHiggsMass, fReferenceHiggsMixing);
+    double width_piplus = PiPlusPartialWidth(fReferenceHiggsMass, fReferenceHiggsMixing, fNewFormFactor);
+    double width_pizero = PiZeroPartialWidth(fReferenceHiggsMass, fReferenceHiggsMixing, fNewFormFactor);
 
     // total lifetime
     double lifetime_ns = Constants::Instance().hbar / (width_elec + width_muon + width_piplus + width_pizero);
@@ -253,8 +253,8 @@ bool HiggsMakeDecay::Decay(const MeVPrtlFlux &flux, const TVector3 &in, const TV
   // Get each partial width
   double width_elec = ElectronPartialWidth(flux.mass, mixing);
   double width_muon = MuonPartialWidth(flux.mass, mixing);
-  double width_piplus = PiPlusPartialWidth(flux.mass, mixing);
-  double width_pizero = PiZeroPartialWidth(flux.mass, mixing);
+  double width_piplus = PiPlusPartialWidth(flux.mass, mixing, fNewFormFactor);
+  double width_pizero = PiZeroPartialWidth(flux.mass, mixing, fNewFormFactor);
 
   // total lifetime
   double lifetime_ns = Constants::Instance().hbar / (width_elec + width_muon + width_piplus + width_pizero);