CovarianceMatrixMaker_plugin.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \brief   Create a set of Covariance matrices to be used in fits
/// \author  brebel@fnal.gov
///////////////////////////////////////////////////////////////////////
#include <cmath>
#include <string>
#include <vector>
#include <sstream>

// Framework includes
#include "art/Framework/Core/ResultsProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Results.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvA includes
#include "NovaDAQConventions/DAQConventions.h"
#include "CovarianceMatrixFit/core/Event.h"
#include "CovarianceMatrixFit/core/VarVals.h"
#include "CovarianceMatrixFit/core/ShifterAndWeighter.h"
#include "CovarianceMatrixFit/core/SpectrumTools.h"
#include "CovarianceMatrixFit/modules/EventListManipulator.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/SelectionUtility.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.h"

#include "TH1.h"
#include "TH2.h"
#include "TRandom3.h"

namespace cmf {
  
  class CovarianceMatrixMaker : public art::ResultsProducer {
  public:
    explicit CovarianceMatrixMaker(fhicl::ParameterSet const& pset);
     ~CovarianceMatrixMaker() = default;
    
    // Plugins should not be copied or assigned.
    CovarianceMatrixMaker(CovarianceMatrixMaker const &)               = delete;
    CovarianceMatrixMaker(CovarianceMatrixMaker &&)                    = delete;
    CovarianceMatrixMaker & operator = (CovarianceMatrixMaker const &) = delete;
    CovarianceMatrixMaker & operator = (CovarianceMatrixMaker &&)      = delete;
    
    void readResults (art::Results const& r) override;
    void writeResults(art::Results      & r) override;
    void clear()                             override;
    void beginJob()                          override;

    void reconfigure (fhicl::ParameterSet const&pset);

  private:
    
    void SetSystematicParameter(cmf::Location & loc);
    void FillSpectrum          (cmf::EventListColl const& eventLists,
                                cmf::ExposureMap   const& exposureMap,
                                cmf::Location      const& loc,
                                TH1D*                     spectrumBin,
                                std::map<long, TH1D*>   & spectraEnergy);

    art::ServiceHandle<art::TFileService>      fTFS;                       ///< TFileService
    int                                        fIterations;                ///< number of iterations for picking sigmas
    TRandom3                                   fRand;                      ///< RNG for picking the variation in parameters
    fhicl::ParameterSet                        fManipulatorParameters;     ///< the fake point parameters
    TH2D*                                      fSystCovarianceMatrix;      ///< Syst Covariance matrix
    TH1D*                                      fNominalSpectrum;           ///< Nominal spectrum for each detector/beam/selection
    std::map<long, TH1D*>                      fNominalEnergySpectra;      ///< Nominal energy spectrum for each detector/file/beam/selection
    TH1D*                                      fPlus1SigmaSpectrum;        ///< Plus 1 Sigma shifted spectrum for each detector/beam/selection
    std::map<long, TH1D*>                      fPlus1SigmaEnergySpectra;   ///< Plus 1 Sigma shifted energy spectrum for each detector/file/beam/selection
    TH1D*                                      fMinus1SigmaSpectrum;       ///< Minus 1 Sigma shifted spectrum for each detector/beam/selection
    std::map<long, TH1D*>                      fMinus1SigmaEnergySpectra;  ///< Minus 1 Sigma shifted energy spectrum for each detector/file/beam/selection
    TH1D*                                      fShiftedSpectrum;           ///< Shifted spectrum for each detector/beam/selection
    TH2D*                                      fShiftedSpectrumByUniverse; ///< Shifted spectrum as a function of universe
    std::map<long, TH1D*>                      fShiftedEnergySpectra;      ///< Shifted energy spectrum for each detector/file/beam/selection
    std::map<std::string, TH1D*>               fShifts;                    ///< Shifts selected for each parameter adjusted
    TH2D*                                      fShiftWeights;              ///< Weights from shifts selected for each parameter adjusted
    fhicl::ParameterSet                        fSAWParameters;             ///< shifter and weighter parameters
  };
  
  //......................................................................
  CovarianceMatrixMaker::CovarianceMatrixMaker(fhicl::ParameterSet const& pset)
  : fRand(TRandom3(0))
  , fSystCovarianceMatrix(nullptr)
  , fNominalSpectrum(nullptr)
  , fPlus1SigmaSpectrum(nullptr)
  , fMinus1SigmaSpectrum(nullptr)
  , fShiftedSpectrum(nullptr)
  , fShiftedSpectrumByUniverse(nullptr)
  , fShiftWeights(nullptr)
  {
    this->reconfigure(pset);
  }
  
  //......................................................................
  void CovarianceMatrixMaker::reconfigure(fhicl::ParameterSet const& pset)
  {
    cmf::SelectionUtility::Instance()    ->Initialize(pset.get< fhicl::ParameterSet >("SelectionsToUse", fhicl::ParameterSet()));
    cmf::CovarianceBinUtility::Instance()->Initialize(pset.get< fhicl::ParameterSet >("BinningToUse"                          ));
    cmf::ParameterUtility::Instance()    ->Initialize(pset.get< fhicl::ParameterSet >("ParametersToUse"                       ));

    fManipulatorParameters = pset.get< fhicl::ParameterSet >("EventListManipulator");
    fIterations            = pset.get< int                 >("Iterations", 1000);
    fSAWParameters         = pset.get< fhicl::ParameterSet >("ShifterAndWeighterParameters");
    auto randSeed          = pset.get< unsigned long int   >("RandomSeed",    0);

    cmf::ShifterAndWeighter::Instance()->InitShiftsAndWeightsToUse(cmf::ParameterUtility::Instance()->ParametersAsLocation(),
                                                                   fSAWParameters);

    if(randSeed > 0) fRand.SetSeed(randSeed);
  }
  
  //......................................................................
  // Method to clear out the collections of data products after the
  // writeResults method is called.
  void CovarianceMatrixMaker::clear()
  {
  }
  

  //......................................................................
  void CovarianceMatrixMaker::beginJob()
  {
    // Create histograms for the covariance matrix, the nominal spectrum and
    // the shifted values
    // Each axis accounts for both detectors, both beam types and all desired selections
    // Don't break out the interaction types or file types because they aren't known in
    // the data
    int numBins = cmf::CovarianceBinUtility::Instance()->TotalBins();

    fSystCovarianceMatrix = fTFS->make<TH2D>("CovarianceMatrix",
                                             ";Energy Bin;Energy Bin",
                                             numBins,
                                             1,
                                             numBins + 1,
                                             numBins,
                                             1,
                                             numBins + 1);

    fNominalSpectrum = fTFS->make<TH1D>("NominalBin",
                                        ";Energy Bin;Events",
                                        numBins,
                                        1,
                                        numBins + 1);
    
    fShiftedSpectrum = fTFS->make<TH1D>("ShiftedBin",
                                        ";Energy Bin;Events",
                                        numBins,
                                        1,
                                        numBins + 1);

    fPlus1SigmaSpectrum = fTFS->make<TH1D>("Plus1SigmaSpectrum",
                                           ";Energy Bin;Events",
                                           numBins,
                                           1,
                                           numBins + 1);

    fMinus1SigmaSpectrum = fTFS->make<TH1D>("Minus1SigmaSpectrum",
                                            ";Energy Bin;Events",
                                            numBins,
                                            1,
                                            numBins + 1);

    fShiftedSpectrumByUniverse = fTFS->make<TH2D>("ShiftedBinByUniverse",
                                                  ";Energy Bin;Universe",
                                                  numBins,
                                                  1,
                                                  numBins + 1,
                                                  fIterations,
                                                  0,
                                                  fIterations);

    std::stringstream shifttitles;
    shifttitles << "Seed: "
                << fRand.GetSeed()
                << ";Shifts;Iterations";
    for(auto const& itr : cmf::ParameterUtility::Instance()->SysParNames()){
      fShifts[itr] = fTFS->make<TH1D>((itr+"Shifts").c_str(),
                                       shifttitles.str().c_str(),
                                       500,
                                       -5.,
                                       5);

    }

    // histogram of the weights as a function of energy bin
    std::stringstream shiftwgttitles;
    shiftwgttitles << "Seed: "
                   << fRand.GetSeed()
                   << ";Energy Bin;Total Shift Weights";
    fShiftWeights = fTFS->make<TH2D>("TotalShiftWeights",
                                     shiftwgttitles.str().c_str(),
                                     numBins,
                                     1,
                                     numBins + 1,
                                     500,
                                     -1.,
                                     1);


    long        key  = 0;
    std::string name;
    std::string title;

    cmf::MetaData md;
    std::vector<double> bins;
    
    // loop over the selections to use
    auto const& dbs = cmf::SelectionUtility::Instance()->SelectionsToUse();

    for(auto const& itr : dbs){
      md.period   = (itr.BeamType() == cmf::kFHC) ? 1 : 4;
      md.detector = itr.Detector();

      title = (md.detector == cmf::kNEARDET) ? ";Energy (GeV);10^{3} Events/10^{20} POT" : ";Energy (GeV);Events/10^{20} POT";

      // loop over the selection types in each detector
      for(auto const& sel : itr.Selections()){

          // no peripheral nue selection in the ND
          if(md.detector == cmf::kNEARDET &&
             sel         == cmf::kNuESelectionPeripheral) continue;

          md.selectionType = sel;

          auto const& highEdges = cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(md);
          bins.clear();
          bins.insert(bins.begin(), cmf::CovarianceBinUtility::Instance()->SelectionLowEdge(md));

          for(auto const& itr : highEdges ) bins.push_back(itr.first);

          key  = md.DetectorBeamSelectionKey();

          name = (cmf::cBeamType_Strings[itr.BeamType()] +
                  cmf::cDetType_Strings[itr.Detector()]  +
                  cmf::cSelectionType_Strings[sel]);

//        LOG_VERBATIM("tpm")
//          << "key: " << cmf::KeyToString(key) << " bins: ";
//        for (auto& b : bins)
//          LOG_VERBATIM("tmp") << b;

          fNominalEnergySpectra[key] = fTFS->make<TH1D>(name.c_str(),
                                                        title.c_str(),
                                                        bins.size() - 1,
                                                        &bins[0]);

      } // end loop over selections
    } // end loop over dbs

    LOG_DEBUG("CovarianceMatrixMaker")
        << "histograms made";

  }
  
  //......................................................................
  void CovarianceMatrixMaker::readResults(art::Results const& /*r*/)
  {
  }
  
  // Function to fill the spectra for each meta data type
  //......................................................................
  void CovarianceMatrixMaker::FillSpectrum(cmf::EventListColl    const& eventLists,
                                           cmf::ExposureMap      const& exposureMap,
                                           cmf::Location         const& loc,
                                           TH1D*                        spectrumBin,
                                           std::map<long, TH1D*>      & spectraEnergy)
   {
     // first reset the spectra
     spectrumBin->Reset();
     for(auto hitr : spectraEnergy) hitr.second->Reset();

     // now call the SpectrumTools.
     std::vector<float> spectrum;
     std::vector<float> count;

     cmf::ShifterAndWeighter::Instance()->SetCurrentVals(loc);
     cmf::FillSpectrum(eventLists,
                       exposureMap,
                       spectrum,
                       count);

     // now fill the spectrum histogram
     long key;
     for(size_t b = 0; b < spectrum.size(); ++b){
       spectrumBin->Fill(b,  spectrum[b]);

       key = cmf::CovarianceBinUtility::Instance()->BinToKey(b);
       if(spectraEnergy.count(key) > 0)
         spectraEnergy[key]->Fill(cmf::CovarianceBinUtility::Instance()->BinToEnergy(b),
                                  spectrum[b]);

       // fill the histogram of weights as a function of energy bin
       if(count[b] > 0) fShiftWeights->Fill(b, spectrum[b]/count[b]);
     } // end loop over spectrum bins
   }
  
  // This function will pick a random value for the specified parameter
  // from a guassian distribution.  The random value will be within +/- 1 sigma
  // of the nominal value
  //......................................................................
  void CovarianceMatrixMaker::SetSystematicParameter(cmf::Location & loc)
  {
    
    // we will pick a random value from a unit gaussian, ie gaussian with
    // sigma = 1.  That will be the fractional change on the sigma we want
    // to apply to the new input point
    double sigFrac = 0.;

    LOG_DEBUG("CovarianceMatrixMaker")
    << "sigfrac = "
    << sigFrac;

    for(auto const& name : cmf::ParameterUtility::Instance()->SysParNames()){

      // each systematic parameter should have a different random number
      // change in sigma so that we don't artificially make them change
      // in a correlated way if doing more than one parameter in a job
      sigFrac = fRand.Gaus();

      // The calibration and light level systematic uncertainties are
      // really only defined for a 1 sigma shift because their excursions
      // are calculated from special data sets with the shift applied
      // Also pull them from a uniform distribution as they are definitely
      // not Gaussian.
      if(name.find("lightmodel")  != std::string::npos ||
         name.find("ckv-proton" ) != std::string::npos ||
         name.find("Calibration") != std::string::npos ){
        sigFrac = fRand.Rndm();

        LOG_DEBUG("CovarianceMatrixMaker")
            << "pulling from uniform distribution "
            << sigFrac;
      }
      
      fShifts[name]->Fill(sigFrac);

      LOG_DEBUG("CovarianceMatrixMaker")
      << " parameter: "
      << name
      << " sigFrac "
      << sigFrac;

      loc.SetParameterValue(name, sigFrac, cmf::kFARDET);
      LOG_DEBUG("CovarianceMatrixMaker")
          << "FD name: "
          << name
          << " "
          << loc.FDLocation().find(name)->second.Value();

      loc.SetParameterValue(name, sigFrac, cmf::kNEARDET);
      LOG_DEBUG("CovarianceMatrixMaker")
          << "ND name: "
          << name
          << " "
          << loc.NDLocation().find(name)->second.Value();

    } // end loop over systematic parameters to change
  }
  
  //......................................................................
  void CovarianceMatrixMaker::writeResults(art::Results & /*r*/)
  {
    // declare some variables we will be using a lot in the for loops below
    int    nbins = fNominalSpectrum->GetNbinsX();
    double diffI;
    double diffJ;
    double iCenter;

    // containers for different oscillation and systematic parameter locations
    auto loc = cmf::ParameterUtility::Instance()->ParametersAsLocation();

    // get the lists of events from the TTree and put them into a map
    // of event lists
    cmf::EventListManipulator manipulator(fManipulatorParameters);

    // we only want the MC lists
    cmf::EventListColl eventLists;
    manipulator.Deserialize(eventLists, cmf::kMC);

    cmf::ExposureMap exposureMap = manipulator.ExposureMap();

    // fill the nominal histograms first so we don't do this for every time we
    // change the systematic parameters
    this->FillSpectrum(eventLists,
                       exposureMap,
                       loc,
                       fNominalSpectrum,
                       fNominalEnergySpectra);

    LOG_DEBUG("CovarianceMatrixMaker")
    << "nominal spectrum filled";

    // now fake input points for the +/- 1 sigma shifts

    auto parNames = cmf::ParameterUtility::Instance()->SysParNames();

    // plus 1 sigma
    for(auto const& itr : parNames){
      loc.SetParameterValue(itr, 1, cmf::kNEARDET);
      loc.SetParameterValue(itr, 1, cmf::kFARDET);
    }

    this->FillSpectrum(eventLists,
                       exposureMap,
                       loc,
                       fPlus1SigmaSpectrum,
                       fPlus1SigmaEnergySpectra);

    // minus 1 sigma
    for(auto const& itr : parNames){
      loc.SetParameterValue(itr, -1, cmf::kNEARDET);
      loc.SetParameterValue(itr, -1, cmf::kFARDET);
    }

    this->FillSpectrum(eventLists,
                       exposureMap,
                       loc,
                       fMinus1SigmaSpectrum,
                       fMinus1SigmaEnergySpectra);

    // loop over all the events and fill the matrices for many different
    // values of the systematic parameters
    for(int k = 0; k < fIterations; ++k){

      LOG_VERBATIM("CovarianceMatrixMaker")
      << "iteration: "
      << k;
      
      this->SetSystematicParameter(loc);

      cmf::ShifterAndWeighter::Instance()->SetCurrentVals(loc);

      // this is a new shift for the parameter, so reset and refill the spectrum
      this->FillSpectrum(eventLists,
                         exposureMap,
                         loc,
                         fShiftedSpectrum,
                         fShiftedEnergySpectra);
    
      // now to fill the covariance matrices
      // For ROOT histograms bin 0 is the underflow bin, so start at bin 1
      // if either diffI or diffJ are 0, then we don't need to fill
      // the bin with a 0 (ie the product) and can continue on
      for(int i = 1; i < nbins + 1; ++i){
        iCenter = fNominalSpectrum->GetBinCenter(i);
        diffI   = (fShiftedSpectrum->GetBinContent(i) - fNominalSpectrum->GetBinContent(i));
        if(diffI != 0.) diffI /= fNominalSpectrum->GetBinContent(i);
        else continue;
        
        for(int j = 1; j < nbins + 1; ++j){
          diffJ = (fShiftedSpectrum->GetBinContent(j) - fNominalSpectrum->GetBinContent(j));
          if (diffJ != 0.) diffJ /= fNominalSpectrum->GetBinContent(j);
          else continue;
          /*
           if(i == 1)
           LOG_VERBATIM("CovarianceMatrixMaker")
           << "i: "
           << i
           << " shifted i: "
           << fShiftedSpectrum->GetBinContent(i)
           << " nominal i: "
           << fNominalSpectrum->GetBinContent(i)
           << " j: "
           << j
           << " shifted j: "
           << fShiftedSpectrum->GetBinContent(j)
           << " nominal j: "
           << fNominalSpectrum->GetBinContent(j)
           << " diffI: "
           << diffI
           << " diffJ: "
           << diffJ
           << " product: "
           << diffI * diffJ;
           */
          fSystCovarianceMatrix->Fill(iCenter, fNominalSpectrum->GetBinCenter(j), diffI * diffJ);

        } // end loop over j bins of spectrum
      } // end loop over i bins of spectrum
    
      // loop over each bin and fill fShiftedSpectrumByUniverse
      for (int i_bin = 1; i_bin < fShiftedSpectrum->GetNbinsX(); ++i_bin){
        fShiftedSpectrumByUniverse->SetBinContent(i_bin, k + 1, fShiftedSpectrum->GetBinContent(i_bin));
      }
    } // end loop over the iterations on the systematic parameter

    // if we're not running any systematic uncertainties, 
    // then assume we want stat only
    // n.b. we're usually running the xsec weights 
    // so take that into account
    float binContent;
    if (parNames.empty() || (parNames.size() == 1 && *(parNames.begin()) == "XSecCVWgt")){
      for (int i_bin = 0; i_bin < fNominalSpectrum->GetNbinsX(); ++i_bin){

        binContent = fNominalSpectrum->GetBinContent(i_bin);

        if(binContent == 0){
          fSystCovarianceMatrix->SetBinContent(i_bin,
                                               i_bin,
                                               std::numeric_limits<float>::min());
        }
        else {
          fSystCovarianceMatrix->SetBinContent(i_bin,
                                               i_bin,
                                               1./binContent);
        }
      }
    }

    // protect against very low and very high covariances
    // these usually happen because nentries in these bins for nominal
    // and shifted spectra are identially 0
    for (int i = 1; i < nbins + 1; ++i){
      for (int j = 1; j < fSystCovarianceMatrix->GetNbinsY()+1; ++j){
        binContent = fSystCovarianceMatrix->GetBinContent(i, j);
        if (std::abs(binContent) < std::pow(10,-14) ||
            std::abs(binContent) > std::pow(10,14)  ||
            std::isnan(binContent)                  ){
          // LOG_VERBATIM("CovarianceMatrixMaker")
          // << "bin " << i 
          // << ", covariance: " << binContent
          // << "\n -- nominal i: " << fNominalSpectrum->GetBinContent(i)
          // << "\n -- shifted i: " << fShiftedSpectrum->GetBinContent(i)
          // << "\n -- nominal minus shifted i: " << fNominalSpectrum->GetBinContent(i) - fShiftedSpectrum->GetBinContent(i)
          // << "\n -- nominal j: " << fNominalSpectrum->GetBinContent(j)
          // << "\n -- shifted j: " << fShiftedSpectrum->GetBinContent(j)
          // << "\n -- nominal minus shifted j: " << fNominalSpectrum->GetBinContent(j) - fShiftedSpectrum->GetBinContent(j)
          // << "\n -- diffI*diffJ" << (fNominalSpectrum->GetBinContent(i)-fShiftedSpectrum->GetBinContent(i))*(fNominalSpectrum->GetBinContent(j)-fShiftedSpectrum->GetBinContent(j));

          fSystCovarianceMatrix->SetBinContent(i, j, 0);
          fSystCovarianceMatrix->SetBinContent(j, i, 0);
        } // end if the bin content is acceptable
      } // end loop over j bins
    } // end loop over i bins
  }
  
  DEFINE_ART_RESULTS_PLUGIN(CovarianceMatrixMaker)
  
} // end cmf namespace