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/modules/EventListManipulator.h"
#include "CovarianceMatrixFit/dataProducts/InputPoint.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.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::InputPoint & point);
    void FillSpectrum          (cmf::EventListMap const & eventLists,
                                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*                                      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                        fInputPointParameters;      ///< input parameters
    std::vector<std::string>                   fSystParameterNames;        ///< names of parameters we are floating from the input point
  };
  
  //......................................................................
  CovarianceMatrixMaker::CovarianceMatrixMaker(fhicl::ParameterSet const& pset)
  : fRand(TRandom3(0))
  , fSystCovarianceMatrix(nullptr)
  , fNominalSpectrum(nullptr)
  , fShiftedSpectrum(nullptr)
  , fShiftedSpectrumByUniverse(nullptr)
  , fShiftWeights(nullptr)
  {
    this->reconfigure(pset);
  }
  
  //......................................................................
  void CovarianceMatrixMaker::reconfigure(fhicl::ParameterSet const& pset)
  {
    // create the manipulator
    fManipulatorParameters = pset.get< fhicl::ParameterSet >("EventListManipulator");
    fIterations            = pset.get< int                 >("Iterations", 1000);
    fInputPointParameters  = pset.get< fhicl::ParameterSet >("InitialGuessInputPoint");
    auto randSeed          = pset.get< unsigned long int   >("RandomSeed",    0);

    if(randSeed > 0) fRand.SetSeed(randSeed);

    auto systematicParameters  = fInputPointParameters.get< fhicl::ParameterSet >("SystematicPullInitial");
    auto parNames = systematicParameters.get_names();

    fhicl::ParameterSet systPars;
    for(const auto &name : parNames){
      LOG_VERBATIM("cmf::CovarianceMatrixMaker")
        << "Running with systematic parameter: " << name;
      systPars = systematicParameters.get<fhicl::ParameterSet>(name);
      
      // don't bother with any systematics parameters that have been set to fixed
      if(systPars.get<bool>("Fixed")) continue;
      
      fSystParameterNames.push_back(systPars.get<std::string>("Parameter"));
    }
    LOG_VERBATIM("cmf::CovarianceMatrixMaker")
      << "Running with " << fSystParameterNames.size() << " systematic parameters";
  }
  
  //......................................................................
  // 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);

    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& name : fSystParameterNames){
      fShifts[name] = fTFS->make<TH1D>((name+"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;
    std::set<cmf::BeamType_t> beamTypes({cmf::kFHC,
                                         cmf::kRHC});
    std::set<novadaq::cnv::DetId> detectors({novadaq::cnv::kNEARDET,
                                             novadaq::cnv::kFARDET});
    std::set<cmf::FileType_t> fileTypes({cmf::kBeam,
                                         cmf::kSwap,
                                         cmf::kTauSwap,
                                         cmf::kRockNonSwap,
                                         cmf::kRockFluxSwap,
                                         cmf::kCosmicBackgroundFile});
    std::set<cmf::SelectionType_t> selectionTypes({kNuESelectionLowPID,
                                                   kNuESelectionHighPID,
                                                   kNuESelectionPeripheral,
                                                   kNuMuSelectionQ1,
                                                   kNuMuSelectionQ2,
                                                   kNuMuSelectionQ3,
                                                   kNuMuSelectionQ4,
                                                   kNCSelection});

    cmf::MetaData md;
    std::vector<double> bins;
    
    // loop over beams
    for(auto const& bt : beamTypes){
      md.epoch = (bt == cmf::kFHC) ? 1 : 4;

      // have to have the same bins in the ND and FD
      for(auto const& det : detectors){
        md.detector = det;

        title = (det == novadaq::cnv::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 : selectionTypes){

          // no peripheral nue selection in the ND
          if(det == novadaq::cnv::kNEARDET &&
             sel == cmf::kNuESelectionPeripheral) continue;

          md.selectionType = sel;

          auto const& highEdges = cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(md);
          bins.clear();
          if(md.IsNuESelected() &&
             md.selectionType != cmf::kNuESelectionPeripheral &&
             md.detector == novadaq::cnv::kFARDET)
            bins.insert(bins.begin(), 1.);
          else
            bins.insert(bins.begin(), 0.);

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

          // loop over file types - the ND only cares about
          // the nonswap files, FD has nonswap, fluxswap, tauswap,
          // rocknonswap, rockfluxswap and cosmic background
          for(auto fitr : fileTypes){

            // in the ND, we only care about the nonswap file type
            if(det == novadaq::cnv::kNEARDET && fitr != cmf::kBeam) continue;

            key  = md.DetectorBeamSelectionKey();
            key += cmf::FileTypeKey(fitr);

            name = (cmf::cBeamType_Strings[bt]          +
                    novadaq::cnv::DetInfo::GetName(det) +
                    cmf::cFileTypeStrings[fitr]         +
                    cmf::cSelectionType_Strings[sel]);

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

          } // end loop over file types
        } // end loop over selections
      } // end loop over detectors
    } // end loop over beam types

    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::EventListMap const& eventLists,
                                           TH1D*                    spectrumBin,
                                           std::map<long, TH1D*>  & spectraEnergy)
   {
    // first reset the spectra
    spectrumBin->Reset();
    for(auto hitr : spectraEnergy)
     hitr.second->Reset();
    
    // declare some variables we will be using a lot in the for loops below
    double weight;
    double energy;
    double normalization;
    long   key        = 0;
    int    bin        = 0;
    int    offset     = 0;
    double potForNorm = 1.e20; // spill summary POT values are in units of 10^12 POT,
                              // so this normalization would be to 10^20 POT

    // loop over each list and event in the list
    for(auto const& itr : eventLists){
      auto const& md = itr.first;
      
      normalization = potForNorm/itr.second.ListSpillSummary().goodPOT;
      
      // normalize the numu selected ND down to a reasonable number of events
      // to have about the same peak as the highest fd bins.  In the end we
      // are just looking to see relative changes in each bin anyway.
      if(md.detector == novadaq::cnv::kNEARDET){
        if     (md.IsNuMuSelected()) normalization *= 1.e-3;
        else if(md.IsNuESelected() ) normalization *= 1.e-2;
        else if(md.IsNCSelected()  ) normalization *= 1.e-3;
      }

      key = md.DetectorBeamSelectionKey();

      // figure out the offsets
      offset = cmf::CovarianceBinUtility::Instance()->KeyToOffset(key);
      
      LOG_DEBUG("CovarianceMatrixMaker")
      << md.ToString()
      << " key: "
      << key
      << " bin offset: "
      << offset;

      // add the file type to the key for the energy spectra
      key += cmf::FileTypeKey(md.fileType);

      for(auto const& evt : itr.second){
        
        energy = evt->DataVals().val_at(cmf::kNu_RecoE, md);
        bin    = cmf::CovarianceBinUtility::Instance()->EnergyToHistogramBin(energy, md);

        // if bin is negative then the energy is out of range for
        // the selection type and we'll ignore it.
        if(bin < 0) continue;

        weight = cmf::ShifterAndWeighter::Instance()->Weight(evt, md, cmf::kSystAndOscWgt);

        // now fill the spectrum
        spectrumBin->Fill(bin + offset,  weight * normalization);

        if(spectraEnergy.count(key) > 0)
          spectraEnergy[key]->Fill(energy, weight * normalization);

        // fill the histogram of weights as a function of energy bin
        fShiftWeights->Fill(bin + offset, weight * normalization);
      } // end loop over events
    } // end loop over lists of events
  }
  
  // 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::InputPoint & point)
  {
    
    // 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
    << " original point "
    << point;

    auto systematicPulls = point.SystematicPulls();

    for(auto const& name : fSystParameterNames) {
      ParameterDet tempPairFar  = std::make_pair(name, novadaq::cnv::kFARDET);
      ParameterDet tempPairNear = std::make_pair(name, novadaq::cnv::kNEARDET);

      // 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;
      
      if(systematicPulls.count(tempPairFar)  != 0){
        systematicPulls[tempPairFar] = sigFrac;
        
        LOG_DEBUG("CovarianceMatrixMaker")
        << "FD name: "
        << name
        << " "
        << systematicPulls[tempPairFar];
      }
      if(systematicPulls.count(tempPairNear) != 0){
        systematicPulls[tempPairNear] = sigFrac;
        
        LOG_DEBUG("CovarianceMatrixMaker")
        << "ND name: "
        << name
        << " "
        << systematicPulls[tempPairFar];
      }
    } // end loop over systematic parameters to change
    
    // Make the new input point
    //
    point = cmf::InputPoint(point.ParameterSpaceLocation(),
                            systematicPulls,
                            point.FixedParameters(),
                            point.FixedSystematics() );
    
    LOG_DEBUG("CovarianceMatrixMaker")
    << " new point "
    << point;
    
    cmf::ShifterAndWeighter::Instance()->SetCurrentVals(point);
  }
  
  //......................................................................
  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;

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


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

    // The shifter and weighter was initialized in the EventListManipulator::Deserialize() method
    // fill the nominal histograms first so we don't do this for every time we
    // change the systematic parameters
    this->FillSpectrum(eventLists,
                       fNominalSpectrum,
                       fNominalEnergySpectra);

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

    // 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(point);
      
      // this is a new shift for the parameter, so reset and refill the spectrum
      this->FillSpectrum(eventLists, 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
    if (fSystParameterNames.empty() || (fSystParameterNames.size() == 1 && fSystParameterNames[0] == "xseccvwgt")){
      for (int i_bin = 0; i_bin < fNominalSpectrum->GetNbinsX(); ++i_bin){

        float 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,
               binContent/(binContent * 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){
        double 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