ParameterUtility.cxx

Go to the documentation of this file.

//
// Created by Brian Rebel on 8/17/20.
//
// NB: This code assumes the SelectionUtility is available to it

#include "CovarianceMatrixFit/utilities/ParameterUtility.h"
#include "CovarianceMatrixFit/utilities/SelectionUtility.h"

namespace cmf{

  //-----------------------------------------------------------------------------
  void ParameterSetToLocation(fhicl::ParameterSet const& pset,
                              cmf::DetType_t      const& det,
                              cmf::Location            & loc)
  {
    // nothing to be done with an empty parameter set
    if(pset.is_empty()) return;

    fhicl::ParameterSet param;
    bool        fixed;
    bool        nuisance;
    bool        constraint;
    float       centralVal;
    float       variance;
    float       sigma;
    double      val;
    float       lowerBound;
    float       upperBound;
    std::string parstr;
    std::string detstr;

    auto const& parNames = pset.get_names();

    for(auto const& itr : parNames){

      param = pset.get<fhicl::ParameterSet>(itr);

      MF_LOG_DEBUG("ParameterSetToLocation")
          << param.to_indented_string();

      detstr     = param.get<std::string>("Detector",     "BOTH");
      parstr     = param.get<std::string>("Parameter"           );
      fixed      = param.get<bool       >("Fixed",        true  );
      val        = param.get<double     >("Value",        0.    );
      nuisance   = param.get<bool       >("Nuisance",     false );
      sigma      = param.get<float      >("Sigma",        1.e6  );
      centralVal = param.get<float      >("CentralValue", 0.    );
      variance   = param.get<float      >("Variance",     1.e-5 );
      constraint = param.get<bool       >("Constraint",   false );
      lowerBound = param.get<float      >("LoVal",        0.    );
      upperBound = param.get<float      >("HiVal",        1.e6  );

      // remove any Wgt from the parameter name
      if(parstr.find("Wgt") != std::string::npos) parstr.erase(parstr.find("Wgt"), 3);

      if(detstr == "BOTH" || cmf::StringToDetectorType(detstr) == det){
        loc.InsertParameter(cmf::Parameter(parstr,
                                           det,
                                           val,
                                           fixed,
                                           sigma,
                                           centralVal,
                                           variance,
                                           nuisance,
                                           constraint,
                                           lowerBound,
                                           upperBound));

        MF_LOG_DEBUG("ParameterSpaceLoc")
            << "inserting parameter "
            << parstr;
      } // end if using the correct detector
    } // end loop over parameter names

    MF_LOG_DEBUG("ParameterSpaceLoc")
      << " full location is "
      << loc;

  }

  static ParameterUtility *gParUtil = nullptr;

  //-----------------------------------------------------------------------------
  ParameterUtility *ParameterUtility::Instance() {
    if(gParUtil == nullptr)
      gParUtil = new ParameterUtility();

    return gParUtil;
  }

  //-----------------------------------------------------------------------------
  ParameterUtility::ParameterUtility()
  {
  }

  //-----------------------------------------------------------------------------
  void ParameterUtility::Initialize(fhicl::ParameterSet const & pset)
  {
    // have default empty vectors for these which is what we expect if
    // making the event lists
    std::vector<std::string> psetNames({"OscillationParameters",
                                        "SystematicParameters"});

    auto const& dets = cmf::SelectionUtility::Instance()->DetectorsToUse();
    for(auto const& detItr : dets){
      for(auto const& nameItr : psetNames){
        cmf::Location loc(dets);
        cmf::ParameterSetToLocation(pset.get<fhicl::ParameterSet>(nameItr), detItr,  loc);

        if(nameItr.find("Oscillation") != std::string::npos)
          fDetOscParameters.emplace_back(detItr, loc.DetectorLocation(detItr).Parameters());
        else if(nameItr.find("Systematic") != std::string::npos)
          fDetSysParameters.emplace_back(detItr, loc.DetectorLocation(detItr).Parameters());

      } // end loop over parameter set names
    } // end loop over detectors

    MF_LOG_VERBATIM("cmf::CovarianceMatrixMaker")
      << "\n-------------------------------------------------------"
      << "\nPARAMETER SUMMARY"
      << "\n-------------------------------------------------------";

    for(auto const& itr : fDetSysParameters){
      MF_LOG_VERBATIM("ParameterUtility")
      << "\n=======================================================\n"
      << cmf::cDetType_Strings[itr.Detector()]
      << "\n======================================================="
      << "\n\nRunning with systematic parameters set to \n";

      for(auto const& parItr : itr.Parameters()){
        MF_LOG_VERBATIM("ParameterUtility")
        << " -- parameter: "
        << std::setw(30)
        << parItr.Name()
        << " val: "
        << std::setw(15)
        << parItr.Value()
        << " fixed: "
        << std::setw(3)
        << parItr.IsFixed();
      } // end loop over parameters for detector

      MF_LOG_VERBATIM("ParameterUtility")
      << "\nRunning with oscillation parameters set to \n";
      for(auto const& odItr : fDetOscParameters){
        if(odItr.Detector() != itr.Detector()) continue;
        for(auto const& oscItr : odItr.Parameters()){
          MF_LOG_VERBATIM("ParameterUtility")
          << " -- parameter: "
          << std::setw(30)
          << oscItr.Name()
          << " val: "
          << std::setw(15)
          << oscItr.Value()
          << " fixed: "
          << std::setw(3)
          << oscItr.IsFixed();
        } // end loop over oscillation parameters
      } // end loop over detectors for oscillation parameters
    } // end loop over systematic parameter detectors
  }

  //-----------------------------------------------------------------------------
  cmf::Parameter const& ParameterUtility::OscParameterInfo(std::string    const& parName,
                                                           cmf::DetType_t const& det) const
  {
    // TODO this method returns the first parameter from any detector matching the name
    // we should probably match the detector too, but as of March 2021, nothing calls it
    for(auto const& dopItr : fDetOscParameters){
      if(dopItr.Detector() != det) continue;

      if(dopItr.count(parName) < 1)
        throw cet::exception("ParameterUtility")
            << "oscillation parameter "
            << parName
            << " is not in the list for the "
            << cmf::cDetType_Strings[det]
            << " detector";

      return dopItr.Parameter(parName);
    } // end loop over detectors

    // we should never get here, but make the compiler happy
    return fDetOscParameters.front().Parameter(parName);
  }

  //-----------------------------------------------------------------------------
  cmf::Parameter const& ParameterUtility::SysParameterInfo(std::string const& parName) const
  {
    // TODO this method returns the first parameter from any detector matching the name
    // we should probably match the detector too, but as of March 2021, nothing calls it
    for(auto const& dspItr : fDetSysParameters){
      if(dspItr.count(parName) < 1)
        throw cet::exception("ParameterUtility")
            << "systematic parameter "
            << parName
            << " is not in the list";

      return dspItr.Parameter(parName);
    } // end loop over detectors

    // we should never get here, but make the compiler happy
    return fDetSysParameters.front().Parameter(parName);
  }

  //-----------------------------------------------------------------------------
  cmf::Parameter const& ParameterUtility::ParameterInfo(std::string const& parName) const
  {
    // the oscillation parameters are the same in all detectors, so the 
    // default of the NOvA FD is fine
    if(cmf::IsOscillationParameter(parName))
      return this->OscParameterInfo(parName);

    // if we don't have an oscillation parameter it must be a systematic parameter
    return this->SysParameterInfo(parName);
  }

  //-----------------------------------------------------------------------------
  bool ParameterUtility::UseSysPar(std::string const& parName)
  {
    // the systematic parameters could be different, so check all the lists
    for(auto const& dspItr : fDetSysParameters){
      if(dspItr.count(parName) > 0 ) return true;
    }

    return false;
  }

  //-----------------------------------------------------------------------------
  cmf::Location ParameterUtility::ParametersAsLocation()
  {

    // loop over the detectors we are using and combine the systematic
    // and oscillation parameters
    std::vector<cmf::ParameterSpaceLoc> pslVec;

    for(auto const& itr : fDetOscParameters){
      std::vector<cmf::Parameter> detPars;

      for(auto const& oscParItr : itr.Parameters()){
        detPars.emplace_back(oscParItr);
      }

      for(auto const& sysDetItr : fDetSysParameters){
        if(sysDetItr.Detector() == itr.Detector()){
          for(auto const& sysParItr : sysDetItr.Parameters()){
            detPars.emplace_back(sysParItr);
          }
        }
      } // end loop to get systematic parameters for the detector

      pslVec.emplace_back(itr.Detector(), detPars);
    }

    cmf::Location loc(pslVec);

    return std::move(loc);
  }

  //-----------------------------------------------------------------------------
  std::set<std::string> ParameterUtility::SysParNames() const
  {
    std::set<std::string> parNames;
    for(auto const& itr : fDetSysParameters){
      for(auto const& parItr : itr.Parameters()) 
        parNames.emplace(parItr.Name());
    } 

    return std::move(parNames);
  }

  //-----------------------------------------------------------------------------
  std::set<uint8_t> ParameterUtility::SysParKeys() const
  {
    std::set<uint8_t> parKeys;
    for(auto const& itr : fDetSysParameters){
      for(auto const& parItr : itr.Parameters()) 
        parKeys.emplace(parItr.Key());
    } 

    return std::move(parKeys);
  }

  //-----------------------------------------------------------------------------
  std::vector<cmf::Parameter> const& ParameterUtility::SysParameters(cmf::DetType_t const& det) const
  {
    for(auto const& itr : fDetSysParameters){
      if(itr.Detector() == det) return itr.Parameters();
    } 

    // if there are no systematic parameters for this detector return an
    // empty vector.
    std::vector<cmf::Parameter> emptyParVec;
    return std::move(emptyParVec);
  }

} // end namespace