ShifterAndWeighter.cxx

Go to the documentation of this file.

//
//  ShifterAndWeighter.cxx
//  Implementation of singleton to calculate weights and keep track of shifts
//  for CovarianceMatrixFit fitting
//

#include "TH2.h"
#include "TFile.h"
#include "TKey.h"
#include <algorithm>

#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"

#include "CovarianceMatrixFit/core/ShifterAndWeighter.h"
#include "CovarianceMatrixFit/core/VarVals.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.h"
#include "CovarianceMatrixFit/utilities/SelectionUtility.h"

#include "OscLib/OscCalcPMNS_NSI.h"
#include "OscLib/OscCalcSterile.h"
#include "OscLib/OscCalcPMNSOptEigen.h"
#include "Utilities/func/MathUtil.h"


namespace cmf {

  static ShifterAndWeighter* gShifterAndWeighter = nullptr;

  //----------------------------------------------------------------------------
  ShifterAndWeighter* ShifterAndWeighter::Instance()
  {
    if(!gShifterAndWeighter) gShifterAndWeighter = new ShifterAndWeighter();

    return gShifterAndWeighter;
  }

  //----------------------------------------------------------------------------
  ShifterAndWeighter::ShifterAndWeighter()
  : fOscCalc (nullptr)
  , fCalcType(cmf::kPMNS)
  {}

  //---------------------------------------------------------------------------
  void ShifterAndWeighter::LoadRatiosFromFile(TFile* file, std::string systName){

    TH1D* ratioHist = nullptr;
    std::string ratioName = "";
     auto const& dbs = cmf::SelectionUtility::Instance()->SelectionsToUse();
     for(auto const& itr : dbs){
       for(auto const& sel : itr.Selections()){

         std::string strippedName = systName;
         if ((systName.substr(systName.length()-2, systName.length()) == "ND") |
             (systName.substr(systName.length()-2, systName.length()) == "FD"))
           strippedName = systName.substr(0, systName.length()-2);

         // TODO: get the biggest excursion of +/- 1 sigma
         ratioName = strippedName
                   + "Total"
                   + cmf::cBeamType_Strings[itr.BeamType()]
                   + cmf::cDetType_Strings[itr.Detector()]
                   + cmf::cSelectionType_Strings[sel]
                   + "Plus1Sigma";

         ratioHist = (TH1D*)file->Get(ratioName.c_str());

         if (ratioHist == nullptr)
           throw cet::exception("HistLookupError")
             << "Histogram: " 
             << ratioName 
             << " not found in file "
             << file->GetName()
             << ". Bailing";

         for(int i = 1; i < ratioHist->GetNbinsX() + 1; ++i){
           fileSystVec[ratioName].push_back(ratioHist->GetBinContent(i));

         }
      }
    }
  }

  //---------------------------------------------------------------------------
  void ShifterAndWeighter::InitCalibSystFile(){
    std::string filePath;
    cet::search_path sp("FW_SEARCH_PATH");
    if(!sp.find_file(calibFileName, filePath)){
      throw cet::exception("FileLookupError")
        << "Calibration systematics file: "
        << calibFileName
        << " cannot be found in FW_SEARCH_PATH. Bailing.";
    }

    calibSystFile = new TFile(filePath.c_str(), "READ");

    if (calibSystFile->IsZombie()){
      throw cet::exception("FileOpenError")
        << "Cannot open file: "
        << filePath
        << ". Bailing";
    }

  }

  //---------------------------------------------------------------------------
  void ShifterAndWeighter::InitEShiftSystFile(){

    std::string filePath;
    cet::search_path sp("FW_SEARCH_PATH");
    if(!sp.find_file(eshiftFileName, filePath)){
      throw cet::exception("FileLookupError")
        << "EShift systematics file: "
        << eshiftFileName
        << " cannot be found in FW_SEARCH_PATH. Bailing.";
    }

    eshiftSystFile = new TFile(filePath.c_str(), "READ");

    if (eshiftSystFile->IsZombie()){
      throw cet::exception("FileOpenError")
        << "Cannot open file: "
        << filePath
        << ". Bailing";
    }

  }


  //----------------------------------------------------------------------------
  void ShifterAndWeighter::InitShiftsAndWeightsForDetector(cmf::DetType_t                 const& det,
                                                           cmf::ParameterSpaceLoc         const& detLoc,
                                                           std::vector<WeighterInfo>           & weighters,
                                                           std::map<std::string, uint8_t>      & shifters)
  {
    // loop over the systematic parameters from the input point and evaluate the
    // name associated with each one, then add the appropriate function to the
    // list.
    weighters.clear();
    shifters .clear();

    std::string systName;
    WeighterInfo weighter;
    cmf::ParType_t parType = cmf::kUnknownPar;
    for(auto const& itr : detLoc){
      if(itr.second.IsOscPar()) continue;
      systName = itr.first;
      parType  = itr.second.ParType();

      if(systName.find("Shifter") != std::string::npos){

        LOG_DEBUG("ShifterAndWeighter")
            << "Adding shifter for "
            << systName
            << " to detector "
            << det;

        this->FillShifterInfo(systName, det);
      } // end if this is a shifter      
      else{
        weighter.wgtName = systName;
        // check all the systematic parameter names and add the necessary functions
        // for weighting if there is a function for that name.
        if(parType == cmf::kCalibSystPar || parType == cmf::kEShiftSystPar){
          if (eshiftSystFile==nullptr && parType == cmf::kEShiftSystPar) {
            this->InitEShiftSystFile();
            this->LoadRatiosFromFile(eshiftSystFile, systName);
          }
          if (calibSystFile==nullptr && parType == cmf::kCalibSystPar){  
            this->InitCalibSystFile();
            this->LoadRatiosFromFile(calibSystFile, systName);
          }
          weighter.wgtFn = std::bind(&ShifterAndWeighter::FileSystWeight, this, systName);
        }
        else if(parType == cmf::kCMFNormSystPar)
          weighter.wgtFn = std::bind(&ShifterAndWeighter::NormSystWeight, this, systName);
        else if(parType == cmf::kCMFSystPar)
          weighter.wgtFn = std::bind(&ShifterAndWeighter::CMFSystVarWeight, this, systName);
        else{
          LOG_WARNING("ShifterAndWeighter")
              << "Unable to determine weighting function for "
              << systName;
          continue;
        }
        LOG_DEBUG("ShifterAndWeighter")
            << "Adding weighter "
            << weighter.wgtName
            << " to detector "
            << cmf::cDetType_Strings[det]
            << "\nthe current weight type is "
            << fWeightType;

        weighters.push_back(weighter);

      } // the parameter was for a weight
    } // end loop over systematic parameters

  }

  //----------------------------------------------------------------------------
  void ShifterAndWeighter::InitShiftsAndWeightsToUse(cmf::Location       const& loc,
                                                                                         fhicl::ParameterSet const& parset)
  {
    // this method sets which systematics are to be used in calculating weights
    // by creating a set of std::functions corresponding to the weighting functions
    // that are to be turned on
    // the code is a bit tedious, but it gets the job done

    calibFileName           = parset.get<std::string>("CalibSystFile"    , "CovarianceMatrixFit/data/calibration/CalibrationSystematics2020.root");
    eshiftFileName          = parset.get<std::string>("EShiftSystFile"   , "CovarianceMatrixFit/data/eshift/EShiftSystematics2020.root");
    auto calcString         = parset.get<std::string>("OscCalculatorType", "PMNS");
    fTrueEOscillationWeight = parset.get<bool       >("TrueEOscCalc"     , true);

    fWeightType = cmf::kOscWgtOnly;
    auto weightString = parset.get<std::string>("WeightType", "OscWeightsOnly");
    if     (weightString == "SystAndOscWeights") fWeightType = cmf::kSystAndOscWgt;
    else if(weightString == "SystWeightsOnly"  ) fWeightType = cmf::kSystWgtOnly;

    cmf::CalcType_t calcType = cmf::kPMNS;
    if     (calcString == "PMNS_NSI"    ) calcType = cmf::kNSI;
    else if(calcString == "PMNS_Sterile") calcType = cmf::kSterile;

    // setup the oscillation calculator to use
    this->OscillationCalculatorToUse(calcType);

    // go ahead and set the current values
    this->SetCurrentVals(loc);

    this->InitShiftsAndWeightsForDetector(cmf::kNEARDET,
                                          loc.NDLocation(),
                                          fWeightersToUseND,
                                          fShiftersToUseND);
    this->InitShiftsAndWeightsForDetector(cmf::kFARDET,
                                          loc.FDLocation(),
                                          fWeightersToUseFD,
                                          fShiftersToUseFD);

    /*
    for(auto const& itr : fShiftersToUseND){
      LOG_VERBATIM("ShifterAndWeighter")
      << "ND shifter param: "
      << itr.first;
      for(auto const& si : itr.second)
        LOG_VERBATIM("ShifterAndWeighter")
        << " associated shifter "
        << si.fShifterName
        << " "
        << si.fDetector;
    }
    for(auto const& itr : fShiftersToUseFD){
      LOG_VERBATIM("ShifterAndWeighter")
      << "FD shifter param: "
      << itr.first;
      for(auto const& si : itr.second)
        LOG_VERBATIM("ShifterAndWeighter")
        << " associated shifter "
        << si.fShifterName
        << " "
        << si.fDetector;
    }
    */

    LOG_DEBUG("ShifterAndWeighter")
      << "Done initialising shifters and weighters";

  }

  //----------------------------------------------------------------------------
  void ShifterAndWeighter::FillShifterInfo(std::string    const& name,
                                           cmf::DetType_t const& det)
  {
    uint8_t param;
    if     (name.find("HadE") != std::string::npos) param = cmf::kHad_RecoE;
    else if(name.find("LepE") != std::string::npos) param = cmf::kLep_RecoE;
    else if(name.find("NuE" ) != std::string::npos) param = cmf::kNu_RecoE;

    if     (det == cmf::kNEARDET) fShiftersToUseND.emplace(name, param);
    else if(det == cmf::kFARDET ) fShiftersToUseFD.emplace(name, param);
  }

  //----------------------------------------------------------------------------
  void ShifterAndWeighter::ReportCurrentVals()
  {
    LOG_VERBATIM("ShifterAndWeighter")
        << "Current shifter/weighter point:\n"
        << fCurrentLoc;
  }

  //----------------------------------------------------------------------------
  // Use this version if there are any parameters with different values in the ND
  // and FD
  void ShifterAndWeighter::SetCurrentVals(cmf::Location const& curLoc)
  {
    fCurrentLoc = cmf::Location(curLoc.NDLocation(),
                                curLoc.FDLocation());
    //this->ReportCurrentVals();

    this->SetOscillationVals();
  }

  //----------------------------------------------------------------------------
  // Only use this version if you are sure the ND and FD parameter values are the
  // same when the same parameter is present and changing during a fit in both
  // detectors
  void ShifterAndWeighter::SetCurrentVals(std::vector<std::string> const& pars,
                                          const double*                   vals)
  {

    for(size_t pv = 0; pv < pars.size(); ++pv){
      fCurrentLoc.SetParameterValue(pars[pv], vals[pv], cmf::kFARDET);
      fCurrentLoc.SetParameterValue(pars[pv], vals[pv], cmf::kNEARDET);
    }
    this->ReportCurrentVals();

    this->SetOscillationVals();
  }

  //----------------------------------------------------------------------------
  void ShifterAndWeighter::SetOscillationVals()
  {
    // get the oscillation parameters from the current point, we haven't set
    // the detector yet, but the only thing that changes is the
    // baseline which we change with the SetCurrentEvent method
    for(auto const& itr : fCurrentLoc.FDLocation()){

      if(fCalcType == cmf::kPMNS || fCalcType == cmf::kNSI){
        if     (itr.first == "L"     ) fOscCalc->SetL     (itr.second.Value());
        else if(itr.first == "Rho"   ) fOscCalc->SetRho   (itr.second.Value());
        else if(itr.first == "Dmsq21") fOscCalc->SetDmsq21(itr.second.Value());
        else if(itr.first == "Dmsq32") fOscCalc->SetDmsq32(itr.second.Value());
        else if(itr.first == "Th12"  ) fOscCalc->SetTh12  (itr.second.Value());
        else if(itr.first == "Th13"  ) fOscCalc->SetTh13  (itr.second.Value());
        else if(itr.first == "Th23"  ) fOscCalc->SetTh23  (itr.second.Value());
        else if(itr.first == "dCP"   ) fOscCalc->SetdCP   (itr.second.Value());
        if(fCalcType == cmf::kNSI){
          //Eps = NSI epsilon parameters
          if     (itr.first == "Eps_ee"     ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_ee     (itr.second.Value());
          else if(itr.first == "Eps_emu"    ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_emu    (itr.second.Value());
          else if(itr.first == "Eps_etau"   ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_etau   (itr.second.Value());
          else if(itr.first == "Eps_mumu"   ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mumu   (itr.second.Value());
          else if(itr.first == "Eps_mutau"  ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_mutau  (itr.second.Value());
          else if(itr.first == "Eps_tautau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetEps_tautau (itr.second.Value());
          else if(itr.first == "Delta_emu"  ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_emu  (itr.second.Value());
          else if(itr.first == "Delta_etau" ) dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_etau (itr.second.Value());
          else if(itr.first == "Delta_mutau") dynamic_cast<osc::OscCalcPMNS_NSI*>(fOscCalc)->SetDelta_mutau(itr.second.Value());
        } // end if NSI oscillations
      } // end if PMNS or NSI oscillations
      else if(fCalcType == cmf::kSterile){
        if     (itr.first == "Th12"  ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(1, 2, itr.second.Value());
        else if(itr.first == "Th23"  ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(2, 3, itr.second.Value());
        else if(itr.first == "Th24"  ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(2, 4, itr.second.Value());
        else if(itr.first == "Th34"  ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetAngle(3, 4, itr.second.Value());
        else if(itr.first == "dCP"   ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDelta(1, 3, itr.second.Value());
        else if(itr.first == "d24"   ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDelta(2, 4, itr.second.Value());
        else if(itr.first == "Dmsq21") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm   (2,    itr.second.Value());
        else if(itr.first == "Dmsq32") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm   (3,    itr.second.Value());
        else if(itr.first == "Dmsq41") dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetDm   (4,    itr.second.Value());
        else if(itr.first == "L"     ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetL    (      itr.second.Value());
        else if(itr.first == "Rho"   ) dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetRho  (      itr.second.Value());
      } // end if sterile oscillations
    } // end loop over FD location
  }

  //----------------------------------------------------------------------------
  void ShifterAndWeighter::SetCurrentEvent(cmf::EventPtr const& curEvent,
                                           cmf::MetaData const& md)
  {
    fCurrentEvent = curEvent.get();
    fCurrentMD    = md;

    LOG_DEBUG("SetCurrentEvent")
      << " Metadata is "
      << md.ToString();

    // for(auto const& itr : fCurrentEvent->MCVals().SystematicShifts()){
    //     LOG_VERBATIM("EventListManipulator")
    //       << "ShifterAndWeighter::SetCurrentEvent() systematic: " 
    //       << itr; 
    // }
    
    // check the validity of the mcVals before trying to get weights, if it is
    // a nullptr, we have a problem.
    if( !fCurrentEvent->MCValsFilled() && md.isMC)
      throw cet::exception("ShifterAndWeighter")
      << "Event does not have a valid MCVarVals unique_ptr, bail";

    // the data values better always be valid
    if( !fCurrentEvent->DataValsFilled() )
      throw cet::exception("ShifterAndWeighter")
      << "Event does not have a valid DataVarVals unique_ptr, bail";

    // reset the oscillation baseline based on the detector
    fOscCalc->SetL(this->CurrentParameterValue("L"));

  }

  //----------------------------------------------------------------------------
  float ShifterAndWeighter::FractionalShift(uint8_t const& param)
  {
    return (fCurrentMD.detector == cmf::kNEARDET) ?
    this->ShiftForParameter(fShiftersToUseND, fCurrentLoc.NDLocation(), param) : this->ShiftForParameter(fShiftersToUseFD, fCurrentLoc.FDLocation(), param);
  }

  //----------------------------------------------------------------------------
  float ShifterAndWeighter::ShiftForParameter(std::map<std::string, uint8_t> const& shifterMap,
                                              cmf::ParameterSpaceLoc         const& loc,
                                              uint8_t                        const& param)
  {

    float totalShift = 1.;
    float shiftScale;
    float curShift;

    for(auto const& itr : shifterMap){

      if(itr.second != param) continue;

      auto const& sigma = loc.find(itr.first)->second.Sigma();

      curShift = this->CurrentParameterValue(itr.first);

      shiftScale = 1. + sigma * curShift;

      // LOG_VERBATIM("ShifterAndWeighterTotalShift")
      // << " associated shifter "
      // << infoItr.fShifterName
      // << " "
      // << infoItr.fDetector
      // << " "
      // << fCurrentMD.detector
      // << " "
      // << sigma
      // << " "
      // << this->CurrentParameterValue(infoItr.fShifterName)
      // << " "
      // << totalShift;
      totalShift *= shiftScale;
    }

    return totalShift;
  }

  //----------------------------------------------------------------------------
  double ShifterAndWeighter::Weight(cmf::EventPtr     const& curEvent,
                                    cmf::MetaData     const& md)
  {
    // get ready to use the passed in event
    this->SetCurrentEvent(curEvent, md);

    // set the weight to be the fake weight value, which is 1 for
    // data.
    double wgt = curEvent->DataVals().val_at(cmf::kFake_Weight, md);


    // if the mcVals is null or the metadata claims it isn't MC,
    // return the fake weight for the event.  If the event is data
    // then the value is always 1, if it is fake data, then the
    // value is the proper weighting
    if(!curEvent->MCValsFilled() ||
       !md.isMC                  ||
       md.IsCosmicMuon()         )
      return wgt;

    //if(md.selectionType == cmf::kNuESelectionPeripheral)
    //  LOG_VERBATIM("ShifterAndWeighter")
    //    << "weight from systematics is "
    //    << wgt;
    
    // check if we only want weights from the systematic parameters or only
    // from the oscillation parameters
    if(fWeightType == cmf::kSystWgtOnly){
      wgt *= this->TotalWeightFromWeighters();
    }
    else if(fWeightType == cmf::kOscWgtOnly){
      wgt *= this->OscillationWeight();
    }
    else
      wgt *= this->OscillationWeight() * this->TotalWeightFromWeighters();

    //if(md.selectionType == cmf::kNuESelectionPeripheral)
    // LOG_VERBATIM("ShifterAndWeighter")
    //   << "weight from oscillations is "
    //   << this->OscillationWeight();

    LOG_DEBUG("ShifterAndWeighter")
    << "Oscillation weight is "
    << wgt
    << " energy: "
    << fCurrentEvent->MCVals().val_at(cmf::kTrueE)
    << " pdg: "
    << fCurrentEvent->MCVals().val_at(cmf::kTruePDG)
    << " orig pdg: "
    << fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig)
    << " detector: "
    << fCurrentMD.detector;

    //if(md.selectionType == cmf::kNuESelectionPeripheral)
    //  LOG_VERBATIM("ShifterAndWeighter")
    //    << "total weight is "
    //    << wgt;
   
    return wgt;
  }

  //----------------------------------------------------------------------------
  double ShifterAndWeighter::TotalWeightFromWeighters()
  {
    double wgt = 1.;

    auto & weighters = (fCurrentMD.detector == cmf::kNEARDET) ? fWeightersToUseND : fWeightersToUseFD;

    //std::string ss;
    for(auto itr : weighters){
      // ss += "\t" + itr.wgtName                       +
      //            " "  + std::to_string(itr.wgtFn())       +
      //            " "  + std::to_string(wgt)               +
      //            " "  + std::to_string(wgt * itr.wgtFn()) +
      //            "\n";
      
      wgt *= itr.wgtFn();
    }
    // LOG_VERBATIM("SAWWeighterInfo")
    // << ss;

    return wgt;
  }

  //----------------------------------------------------------------------------
  void ShifterAndWeighter::OscillationCalculatorToUse(cmf::CalcType_t calcType)
  {
    // nothing to change
    if(fOscCalc != nullptr && fCalcType == calcType) return;
    else{
      LOG_WARNING("ShifterAndWeighter")
          << "resetting oscillation calculator to type "
          << calcType
          << " from "
          << fCalcType
          << " are you sure you meant to change types in this job?";

      delete fOscCalc;
      fOscCalc = nullptr;
    }

    fCalcType = calcType;

    if     (fCalcType == cmf::kPMNS   ) fOscCalc = new osc::OscCalcPMNSOptEigen();
    else if(fCalcType == cmf::kNSI    ) fOscCalc = new osc::OscCalcPMNS_NSI();
    else if(fCalcType == cmf::kSterile){
      fOscCalc = new osc::OscCalcSterile();
      dynamic_cast<osc::OscCalcSterile*>(fOscCalc)->SetNFlavors(4);
    }
  }

  //----------------------------------------------------------------------------
  double ShifterAndWeighter::OscillationWeight()
  {
    // don't bother doing anything if we are looking at the near detector
    // and standard oscillations, it is just a waste of time.
    if(fCalcType           == cmf::kPMNS &&
       fCurrentMD.detector == cmf::kNEARDET) return 1.;

    if(fTrueEOscillationWeight) return this->OscillationWeightTrueE();

    return this->OscillationWeightBinCenter();
  }
  //----------------------------------------------------------------------------
  double ShifterAndWeighter::OscillationWeightTrueE()
  {
    double E       = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
    int    pdg     = fCurrentEvent->MCVals().val_at(cmf::kTruePDG);
    int    pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);

    // only attempt to weight neutrinos for oscillations
    if(std::abs(pdg) != 12 &&
       std::abs(pdg) != 14 &&
       std::abs(pdg) != 16) return 1.;

    return fOscCalc->P(pdgOrig, pdg, E);

  }

  //----------------------------------------------------------------------------
  double ShifterAndWeighter::OscillationWeightBinCenter()
  {
    if(fTrueEnergyBins.empty()){
      const int numTrueEnergyBins = 100;

      // How many edges to generate. Allow room for 0-Emin bin
      const double N = numTrueEnergyBins-1;
      const double A = N * 0.5; // 500 MeV: there's really no events below there

      fTrueEnergyBins.insert(0);
      for(int i = 1; i < N - 1; ++i) fTrueEnergyBins.insert(A/i);
      fTrueEnergyBins.insert(120.);
    }

    double E       = fCurrentEvent->MCVals().val_at(cmf::kTrueE);
    int    pdg     = fCurrentEvent->MCVals().val_at(cmf::kTruePDG);
    int    pdgOrig = fCurrentEvent->MCVals().val_at(cmf::kTruePDGOrig);

    // get the relevant bin for the current neutrino, then get the
    // energy at the center of that bin
    double eUpper = *(fTrueEnergyBins.rbegin());
    double eLower = *(fTrueEnergyBins.begin());

    // find the bin edge that goes after the value of E
    auto itr = fTrueEnergyBins.upper_bound(E);
    if(itr != fTrueEnergyBins.end()) eUpper = *itr;
    else --itr;

    // now find the edge before by just decrementing the iterator
    --itr;
    if(itr != fTrueEnergyBins.end()) eLower = *itr;

    // get the midpoint between the upper and lower values,
    // ie take their mean
    E = 0.5 * (eUpper + eLower);

    LOG_DEBUG("ShifterAndWeighter")
    << "true energy: "
    << fCurrentEvent->MCVals().val_at(cmf::kTrueE)
    << " central bin energy "
    << E;

    // only attempt to weight neutrinos for oscillations
    if(std::abs(pdg) != 12 &&
       std::abs(pdg) != 14 &&
       std::abs(pdg) != 16) return 1.;

    return fOscCalc->P(pdgOrig, pdg, E);

  }

  //----------------------------------------------------------------------------
  double ShifterAndWeighter::CalcLinearInterpWgt(double sigma,
                                                 std::map<float, float> const& wgtsBySigma)
  {
    auto   mapEnd    = wgtsBySigma.end();

    // no way this is a valid number and is obviously stupid
    // not using a std::numeric_limits value to avoid rounding to 0
    double badVal    = -99999999999999.; 
    double minus2Val = (wgtsBySigma.find(-2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-2.)->second) : badVal;
    double minus1Val = (wgtsBySigma.find(-1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find(-1.)->second) : badVal;
    //    double zeroVal   = 1.;
    double zeroVal   = (wgtsBySigma.find( 0.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 0.)->second) : 1.;
    double plus1Val  = (wgtsBySigma.find( 1.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 1.)->second) : badVal;
    double plus2Val  = (wgtsBySigma.find( 2.) != mapEnd) ? static_cast<double>(wgtsBySigma.find( 2.)->second) : badVal;
    double slope     = 0.;
    double intercept = 0.;

    //Debugging RPA weights
    LOG_DEBUG("ShifterAndWeighter")
      << "\tCalcLinearInterpWgt: "
      <<  " -2:" << minus2Val
      << ", -1:" << minus1Val
      <<  ", 0:" << zeroVal
      <<  ", 1:" << plus1Val
      <<  ", 2:" << plus2Val;

    /*
    if(wgtsBySigma.find(-1.) != mapEnd)
      LOG_VERBATIM("ShifterAndWeighter")
       << "\t\t 0 val from "
       << wgtsBySigma.find(0.)->second;
    */

    // The following code assumes that at least the -1 sigma, +1 sigma
    // values are in the map.  If they are not, we have a problem
    if(minus1Val == badVal ||
       plus1Val  == badVal){
      LOG_WARNING("ShifterAndWeighterCalcLinearInterpWgt")
      << "No values associated with either -1 sigma or +1 sigma "
      << "so we cannot calculate a weight, return 1.";
      return 1.;
    }

    if(sigma <= -1){
      // It is possible there is no -2 sigma value from the map,
      // in that case just interpolate from 0 to -1 sigma
      if(minus2Val != badVal){
        slope     = (minus1Val - minus2Val);
        intercept = minus2Val - slope * (-2.);
      }
      else{
        slope     = (zeroVal - minus1Val);
        intercept = minus1Val - slope * (-1.);
      }
    }
    else if(sigma > -1. && sigma <= 0.){
      slope     = (zeroVal - minus1Val);
      intercept = minus1Val - slope * (-1.);
    }
    else if(sigma > 0. && sigma <= 1.){
      slope     = (plus1Val - zeroVal);
      intercept = plus1Val - slope * 1.;
    }
    else if(sigma > 1.){
      // It is possible there is no +2 sigma value from the map,
      // in that case just interpolate from 0 to +1 sigma
      if(plus2Val != badVal){
        slope     = (plus2Val - plus1Val);
        intercept = plus2Val - slope * 2.;
      }
      else{
        slope     = (plus1Val - zeroVal);
        intercept = plus1Val - slope * 1.;
      }
    }

    double wgt = sigma * slope + intercept;
   
    // Don't let weights be negative
    if(wgt < 0) wgt = 0.;
 
    return wgt;
  }

  //-------------------------------------------------------------------------
  double ShifterAndWeighter::NormSystWeight(std::string const& systName){
    // we have a separate handling for normalisation systematics because 
    // there's a bunch of numbers it's just not worth storing in the event 
    // list
    double oneSigmaWgt = 1.0;
    
    double normPOT              = 0.0055;
    double normNDMass           = 0.003;
    double normFDMass           = 0.003;
    double normML               = 0.008;
    double normNDNumuContext = fCurrentMD.BeamType() == cmf::kFHC ? 0.0021 : 0.0048;
    double normNDNueContext  = fCurrentMD.BeamType() == cmf::kFHC ? 0.0041 : 0.0053;
    double normNDNCContext   = fCurrentMD.BeamType() == cmf::kFHC ? 0.01   : 0.00;

    bool isND   = fCurrentMD.detector == cmf::kNEARDET ? true : false;
    bool isFD   = fCurrentMD.detector == cmf::kFARDET  ? true : false; 
    bool isFHC  = fCurrentMD.BeamType() == cmf::kFHC   ? true : false;
    bool isRHC  = fCurrentMD.BeamType() == cmf::kRHC   ? true : false;
    bool isNuMu = false; 
    bool isNuE  = false; 
    bool isNC   = false;

    if (fCurrentMD.selectionType == cmf::kNuMuSelection  ||
       fCurrentMD.selectionType == cmf::kNuMuSelectionQ1 ||
       fCurrentMD.selectionType == cmf::kNuMuSelectionQ2 ||
       fCurrentMD.selectionType == cmf::kNuMuSelectionQ3 ||
       fCurrentMD.selectionType == cmf::kNuMuSelectionQ4) 
      isNuMu = true;
    if (fCurrentMD.selectionType == cmf::kNuESelectionHighPID ||
       fCurrentMD.selectionType == cmf::kNuESelectionLowPID   ||
       fCurrentMD.selectionType == cmf::kNuESelectionPeripheral)
      isNuE = true;
    if (fCurrentMD.selectionType == cmf::kNCSelection)
      isNC = true;

    if      (systName == "NormalizationPOT")                                       oneSigmaWgt += normPOT;
    else if (systName == "NormalizationNDMass" && isND)                            oneSigmaWgt += normNDMass;
    else if (systName == "NormalizationFDMass" && isFD)                            oneSigmaWgt += normFDMass;
    else if (systName == "NormalizationML" && isFD)                                oneSigmaWgt += normML;
    else if (systName == "NormalizationFHCNuMuContext" && isND && isNuMu && isFHC) oneSigmaWgt += normNDNumuContext;
    else if (systName == "NormalizationRHCNuMuContext" && isND && isNuMu && isRHC) oneSigmaWgt += normNDNumuContext;
    else if (systName == "NormalizationFHCNuEContext" && isND && isNuE && isFHC)   oneSigmaWgt += normNDNueContext;
    else if (systName == "NormalizationRHCNuEContext" && isND && isNuE && isRHC)   oneSigmaWgt += normNDNueContext;
    else if (systName == "NormalizationFHCNCContext" && isND && isNC && isFHC)     oneSigmaWgt += normNDNCContext;
    else if (systName == "NormalizationRHCNCContext" && isND && isNC && isRHC)     oneSigmaWgt += normNDNCContext;

    float sigShift = this->CurrentParameterValue(systName);
    return 1. - (1. - oneSigmaWgt) * sigShift;

  }

  //-------------------------------------------------------------------------
  double ShifterAndWeighter::FileSystWeight(std::string const& systName)
  {
    bool isWrongDet = false;
    if (systName.find("ND") != std::string::npos && fCurrentMD.detector == cmf::kFARDET)  isWrongDet = true;
    if (systName.find("FD") != std::string::npos && fCurrentMD.detector == cmf::kNEARDET) isWrongDet = true;

    std::string strippedName = systName;
    if ((systName.substr(systName.length()-2, systName.length()) == "ND") |
        (systName.substr(systName.length()-2, systName.length()) == "FD"))
      strippedName = systName.substr(0, systName.length()-2);

    std::string ratioName = strippedName
                          + "Total"
                          + cmf::cBeamType_Strings[fCurrentMD.BeamType()]
                          + cmf::cDetType_Strings[fCurrentMD.detector]
                          + cmf::cSelectionType_Strings[fCurrentMD.selectionType]
                          + "Plus1Sigma";

    float energy = fCurrentEvent->DataVals().val_at(cmf::kNu_RecoE, fCurrentMD);
    int   bin    = cmf::CovarianceBinUtility::Instance()->EnergyToBin(energy, fCurrentMD);

    float oneSigmaWgt = isWrongDet ? 1 : fileSystVec[ratioName][bin];
    float sigShift    = this->CurrentParameterValue(systName);

    return 1. - (1. - oneSigmaWgt) * sigShift;
  }
  
  //----------------------------------------------------------------------------
  double ShifterAndWeighter::CMFSystVarWeight(std::string const& systName)
  {
    // return 1 if this is a cosmic muon
    if(fCurrentMD.interactionType == cmf::kCosmicMuon) return 1.;

    // if this weight has NC or CC in the name and this interaction type
    // is not that, return 1
    if((systName.find("CC") != std::string::npos && fCurrentMD.interactionType == cmf::kNC) ||
       (systName.find("NC") == std::string::npos && fCurrentMD.interactionType == cmf::kNC)  ) return 1.;

    // may need to come up with default for cosmic muons, however I think the fact that
    // the default weights for every GENIE variable is a collection of 1's should
    // mean this is OK
   
    std::map<float, float> wgtBySigma = fCurrentEvent->MCVals().SystVarWgts(cmf::VarNameToKey(systName));

    //for(auto itr : wgtBySigma)
    //  LOG_VERBATIM("CMFSystVarWeight")
    //    << systName
    //    << " weights "
    //    << itr.first
    //    << " "
    //    << itr.second;

    return this->CalcLinearInterpWgt(this->CurrentParameterValue(systName),
                                     wgtBySigma);
  }

  //----------------------------------------------------------------------------
  double ShifterAndWeighter::XSecCVPPFXWeight()
  {
    return fCurrentEvent->MCVals().val_at(cmf::kXSecCV2020Wgt) * fCurrentEvent->MCVals().val_at(cmf::kPPFXFluxCVWgt);
  }

} // cmf