Parameter.cxx

Go to the documentation of this file.

//
// Created by Brian Rebel on 8/17/20.
//

#include "CovarianceMatrixFit/dataProducts/Constants.h"
#include "CovarianceMatrixFit/dataProducts/Parameter.h"
#include "CovarianceMatrixFit/dataProducts/StaticFuncs.h"

#include "cetlib_except/exception.h"

namespace cmf{

  //----------------------------------------------------------------------------
  size_t ParameterSpaceLoc::count(std::string const& parName) const
  {
    size_t cnt = 0;
    for(auto const& itr : fPars) if(itr.Name() == parName) ++cnt;

    if(cnt > 1)
      throw cet::exception("ParameterSpaceLoc") 
      << "Parameter "
      << parName
      << " shows up multiple times for detector "
      << cmf::cDetType_Strings[fDetector];

    return cnt;
  }

  //----------------------------------------------------------------------------
  cmf::Parameter const& ParameterSpaceLoc::Parameter(std::string const& parName) const
  {
    for(auto const& itr : fPars){
      if(itr.Name() == parName){
        return itr;
      }
    }

    // if we got here, we didn't find the parameter
    throw cet::exception("ParameterSpaceLoc")
    << "no parameter corresponding to "
    << parName
    << " found for detector "
    << cmf::cDetType_Strings[fDetector];

    // make the compiler happy, although we will never get here
    return fPars.front();
  }

  //----------------------------------------------------------------------------
  void ParameterSpaceLoc::SetParameterValue(std::string const& parName, 
                                            double      const& val)
  {
    bool foundPar = false;
    for(auto & itr : fPars){
      if(itr.Name() == parName){
        itr.SetValue(val);
        foundPar = true;
      }
    }

    if(!foundPar)
      throw cet::exception("ParameterSpaceLoc")
      << "no parameter corresponding to "
      << parName
      << " found for detector "
      << cmf::cDetType_Strings[fDetector]
      << " to set to value "
      << val;

    return;
  }

  //----------------------------------------------------------------------------
  void ParameterSpaceLoc::InsertParameter(cmf::Parameter const& par)
  {
    // first check that we are looking at the correct detector and 
    // if the parameter is already in there
    if(par.Detector() != fDetector) return;

    for(auto const& itr : fPars){
      if(itr.Name() == par.Name()) return;
    }

    fPars.emplace_back(par);
  }

  //------------------------------------------------------------
  std::ostream& operator<<(std::ostream& o, cmf::OscPar const& op) {
    o << cmf::cOscParams_Strings[op.Param()]
      << " " 
      << op.Value();
      
    return o;
  }

  //----------------------------------------------------------------------------
  Parameter::Parameter()
  : fDetector     (cmf::kFARDET       )
  , fVal          (cmf::kGarbageDouble)
  , fCentralVal   (cmf::kGarbageFloat )
  , fSigmaSize    (1.e-5              )
  , fVariance     (1.e-5              )
  , fIsFixed      (false              )
  , fIsNuisance   (false              )
  , fIsConstrained(false              )
  , fLowerBound   (cmf::kGarbageFloat )
  , fUpperBound   (cmf::kGarbageFloat )
  , fParType      (cmf::kUnknownPar   )
  {
  }

  //----------------------------------------------------------------------------
  Parameter::Parameter(std::string    const& name,
                       cmf::DetType_t const& det,
                       double         const& val,
                       bool           const& fixed,
                       float          const& sigmaSize,
                       float          const& cv,
                       float          const& variance,
                       bool           const& nuisance,
                       bool           const& constrained,
                       float          const& lowerBound,
                       float          const& upperBound)
  : fDetector     (det             )
  , fVal          (val             )
  , fCentralVal   (cv              )
  , fSigmaSize    (sigmaSize       )
  , fVariance     (variance        )
  , fIsFixed      (fixed           )
  , fIsNuisance   (nuisance        )
  , fIsConstrained(constrained     )
  , fLowerBound   (lowerBound      )
  , fUpperBound   (upperBound      )
  , fParType      (cmf::kUnknownPar)
  {
    if(fVariance < 1.e-5) fVariance = 1.e-5;

    if(cmf::IsOscillationParameter(name)){
      fParType = cmf::kOscPar;
      fKey     = cmf::StringToOscParmType(name);
    }
    else{
      fKey = cmf::VarNameToKey(name);
      if(fKey >= cmf::kFirstSigmaVar && fKey <= cmf::kLastSigmaVar) 
        fParType = cmf::kCMFSystPar;
      else if (fKey >= kFirstNormalizationVar && fKey <= cmf::kLastNormalizationVar)
        fParType = cmf::kCMFNormSystPar;
      else if (fKey >= kFirstCalibrationVar && fKey <= cmf::kLastCalibrationVar) 
        fParType = cmf::kCalibSystPar;
      else if (fKey >= kFirstEShiftVar && fKey <= cmf::kLastEShiftVar)
        fParType = cmf::kEShiftSystPar;
      else if (fKey == kXSecCV2020Wgt || fKey == kXSecCVPPFX_Weight || fKey == kPPFXFluxCVWgt || fKey == kRPARES_Weight)
        fParType = cmf::kCMFWeightPar;
      else
        fParType = cmf::kOtherSystPar;
    }
      

    // sanity check warning - if the parameter is fixed and labeled as nuisance,
    // tell the user that the fixed part wins
    if(fIsFixed && fIsNuisance){
      MF_LOG_WARNING("Parameter")
          << name
          << " is set to be fixed but is also set to be a nuisance parameter, "
          << "assume that leaving it fixed is correct and set IsNuisance to "
          << " false";

      fIsNuisance = false;
    }

  }

  //----------------------------------------------------------------------------
  bool Parameter::operator<(cmf::Parameter const& other) const {
    return fKey < other.Key();
  }

  //----------------------------------------------------------------------------
  bool Parameter::operator==(cmf::Parameter const& other) const {
    return fKey == other.Key();
  }

  //------------------------------------------------------------
  std::ostream& operator<<(std::ostream& o, cmf::Parameter const& par) {

    o  << "-- parameter: "
       << std::setw(30)     << par.Name()  << std::setw(10)
       << " val: "          << par.Value()
       << " fixed: "        << par.IsFixed()
       << " oscillation "   << par.IsOscPar()
       << " nuisance "      << par.IsNuisance()
       << " constrained "   << par.IsConstrained()
       << " lower: "        << par.LowerBound()
       << " upper: "        << par.UpperBound()
       << std::endl;

    return o;
  }

  //------------------------------------------------------------
  std::ostream& operator<<(std::ostream& o, cmf::Location const& loc) 
  {
    auto const& detLocs = loc.DetectorLocations();

    std::set<std::string> parsDone;

    o << std::endl;
    for(auto const& itr : detLocs){
      for(auto const& parItr : itr.Parameters()){

        if(parsDone.find(parItr.Name()) != parsDone.end()) continue;

        parsDone.insert(parItr.Name());

        o << " -- parameter: "  
          << std::setw(30)
          << parItr.Name()
          << std::setw(15)
          << " val "
          << cmf::cDetType_Strings[itr.Detector()]
          << ": "
          << parItr.Value()
          << " fixed: "
          << std::setw(5)
          << parItr.IsFixed();

        // see if this parameter is set for other detectors too
        for(auto const& itr2 : detLocs){
          if(itr2.Detector() == itr.Detector()) continue;
          if(itr2.count(parItr.Name()) > 0){
            auto const& par2 = itr2.Parameter(parItr.Name());
            o << std::setw(15)
              << " val "
              << cmf::cDetType_Strings[itr2.Detector()]
              << ": "
              << std::setw(15)
              << " "
              << par2.Value()
              << " fixed: "
              << std::setw(5)
              << par2.IsFixed();
          }
        } // end loop over other detectors
              o << std::endl;
      } // end loop over parameters 
      o << std::endl;
    } // end loop over the current detector

    return o;
  }

  //----------------------------------------------------------------------------
  Location::Location(std::set<cmf::DetType_t> const& dets)
  {
    for(auto const& det : dets) fDetLocs.emplace_back(det, std::vector<cmf::Parameter>());
  }

  //----------------------------------------------------------------------------
  Location::Location(cmf::OscillationParameterMap const& oscPars,
                     std::set<cmf::DetType_t>     const& dets)
  : Location(dets)
  {
    std::map<cmf::DetType_t, double> detBaseline({{cmf::kNEARDET,   1.  },
                                                  {cmf::kFARDET,    810.},
                                                  {cmf::kMINOSNEAR, 1.02},
                                                  {cmf::kMINOSFAR,  735.}});

    float parVal;

    for(auto const& itr : oscPars){
      parVal = itr.second;
      for(auto & locItr : fDetLocs){
        if(itr.first == cmf::kL) parVal = detBaseline.find(locItr.Detector())->second;

        locItr.InsertParameter(cmf::Parameter(cmf::cOscParams_Strings[itr.first],
                                              locItr.Detector(),
                                              parVal,
                                              true,
                                              0,
                                              parVal,
                                              0,
                                              false,
                                              false,
                                              0.,
                                              1.));
      } // end loop over locations
    } // end loop over oscillation parameters
  }

  //----------------------------------------------------------------------------
  void Location::InsertParameter(cmf::Parameter const& par)
  {
    for(auto & itr : fDetLocs){
      if(itr.Detector() == par.Detector()){
        if(itr.count(par.Name()) < 1) itr.InsertParameter(par);
        else itr.SetParameterValue(par.Name(), par.Value());
      }
    }
  }

  //----------------------------------------------------------------------------
  void Location::SetParameterValue(std::string    const& parName,
                                   double         const& value,
                                   cmf::DetType_t const& det)
  {
    for(auto & itr : fDetLocs){
      if(itr.Detector() == det || 
         cmf::kALLDET   == det ){
        if(itr.count(parName) > 0) 
          itr.SetParameterValue(parName, value);
        else
          MF_LOG_WARNING("Location")
          << parName
          << " is not a recognized parameter for the "
          << cmf::cDetType_Strings[det];
      }
    } // end loop over the ParameterSpaceLocs
  }

  //----------------------------------------------------------------------------
  cmf::OscillationParameterMap Location::AsOscillationParameterMap(cmf::DetType_t const& det) const
  {
    cmf::OscillationParameterMap oscMap;

    // The oscillation parameters should be the same for all detectors, 
    // except of course the baseline. Make the map with the NOvA far detector
    // by default
    for(auto const& itr : this->DetectorLocation(det).Parameters()){
      if(!itr.IsOscPar()) continue;

      oscMap.emplace((cmf::OscParm_t)itr.Key(),
                     itr.Value());
    }

    return std::move(oscMap);
  }

  //----------------------------------------------------------------------------
  cmf::ParameterSpaceLoc const& Location::DetectorLocation(cmf::DetType_t const& det) const
  {
    for(auto const& itr : fDetLocs){
      if(itr.Detector() == det) return itr;
    }
    
    // if we got here there is a problem because the detector is not in the 
    // collection
    throw cet::exception("Location")
    << "no cmf::ParameterSpaceLoc for "
    << cmf::cDetType_Strings[det];

    // the compiler will complain if we don't pretend to return what we 
    // say we will
    return fDetLocs.front();
  }

  //----------------------------------------------------------------------------
  std::vector<Parameter> Location::OscillationParameters(cmf::DetType_t const& det) const
  {
    std::vector<Parameter> oscPars;

    MF_LOG_DEBUG("Location")
      << "Get oscillation parameters from collection of "
      << fDetLocs.size()
      << " detector locations";
    
    for(auto const& itr : fDetLocs){
      if(itr.Detector() == det){
            MF_LOG_DEBUG("Location")
              << "Oscillation Parameters for "
              << cmf::cDetType_Strings[itr.Detector()];
        for(auto const& parItr : itr.Parameters()){
                if(parItr.IsOscPar()){
                  oscPars.emplace_back(parItr);
                  MF_LOG_DEBUG("Location")
                    << parItr;
                }
        }
      } 
    }

    return std::move(oscPars);
  }
  
  //----------------------------------------------------------------------------
  std::vector<Parameter> Location::SystematicParameters(cmf::DetType_t const& det) const
  {
    std::vector<Parameter> systPars;

    for(auto const& itr : fDetLocs){
      if(itr.Detector() == det){
        for(auto const& parItr : itr.Parameters()){
          if(!parItr.IsOscPar()) systPars.emplace_back(parItr);
        }
      } 
    }

    return std::move(systPars);
  }


}// namespace