RandomUniverseUtility.cxx

Go to the documentation of this file.

//
// Created by Brian Rebel on 4/20/20.
//

// Framework includes
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "art_root_io/TFileDirectory.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib_except/exception.h"

#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/RandomUniverseUtility.h"
#include "CovarianceMatrixFit/utilities/PlotUtilities.h"

#include "TH1.h"

namespace cmf{

  static RandomUniverseUtility *gRUU = nullptr;

  //----------------------------------------------------------------------------
  RandomUniverseUtility* RandomUniverseUtility::Instance()
  {
    if(gRUU == nullptr) gRUU = new RandomUniverseUtility();

    return gRUU;
  }

  //----------------------------------------------------------------------------
  RandomUniverseUtility::RandomUniverseUtility()
  {
  }

  //----------------------------------------------------------------------------
  RandomUniverseUtility::~RandomUniverseUtility()
  {
  }

  //----------------------------------------------------------------------------
  void RandomUniverseUtility::Initialize(fhicl::ParameterSet const& pset)
  {
    // set the seed for the random number to be what is configured, plus the
    // start universe value in case lots are jobs are submitted with the same seed
    // by accident
    fRand.SetSeed(pset.get<unsigned long int>("RandomSeed",    0) +
                  pset.get<unsigned int     >("StartUniverse", 0));

    fNumUniverses     = pset.get<unsigned int>("NumUniverses",                 1);
    fPoissonFluctuate = pset.get<bool        >("PoissonFluctuate",         false);
    fPoissonByHist    = pset.get<bool        >("PoissonDrawWithHistogram", false);
    fPoissonByBin     = pset.get<bool        >("PoissonDrawWithBin",       true );

    // fill histogram with systematic names and which values were used
    art::ServiceHandle<art::TFileService> tfs;

    fSystDrawHist = tfs->make<TH2D>("SystematicParameterDraws",
                                    ";Universe;Uncertainty Draw;",
                                    fNumUniverses,
                                    0,
                                    fNumUniverses,
                                    2,
                                    0,
                                    2);
    fSystDrawHist->SetCanExtend(TH1::kAllAxes);
    fSystDrawHist->SetStats(false);

    fOscParamDrawHist = tfs->make<TH2D>("OscillationParameterDraws",
                                        ";Unvierse;Parameter Draw",
                                        fNumUniverses,
                                        0,
                                        fNumUniverses,
                                        2,
                                        0,
                                        2);
    fOscParamDrawHist->SetCanExtend(TH1::kAllAxes);
    fOscParamDrawHist->SetStats(0);
  }

  //......................................................................
  void RandomUniverseUtility::CreateLocationsWithVariedSystematics(cmf::Location              const& nominalLoc,
                                                                   std::vector<cmf::Location>      & vecLoc)
  {
    // check that the vector is the correct size
    vecLoc.resize(fNumUniverses, nominalLoc);

    for(unsigned int u = 0; u < fNumUniverses; ++u){

      this->CreateLocationWithVariedSystematics(nominalLoc, vecLoc[u], u);
    } // end loop over the number of universes

    // fSystDrawHist->LabelsDeflate("Y");
    // fSystDrawHist->LabelsOption("h", "Y");

  }

  // 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 RandomUniverseUtility::CreateLocationWithVariedSystematics(cmf::Location const& nominalLoc,
                                                                  cmf::Location      & vecLoc,
                                                                  size_t               universeNum)
  {
    // 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.;

    // get the names of the systematic parameters in the ND and FD
    std::map<std::string, double> systParVals;
    for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
      for(auto const& itr : nominalLoc.DetectorLocation(detItr).Parameters()){
        if(itr.IsOscPar() ||
           itr.IsFixed()  ||
           systParVals.count(itr.Name()) ) continue;
  
        // 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 = (this->UniformDrawParameters(itr.Name())) ? fRand.Rndm() : fRand.Gaus();
  
        systParVals.emplace(itr.Name(), sigFrac);
      }
    }

    // prepare some output for the random values in this universe
    std::string systStr("Systematic Parameters " + std::to_string(universeNum));

    for(auto const& itr : systParVals){

      systStr += "\n\t" + itr.first + ": " + std::to_string(itr.second);
      
      fSystDrawHist->Fill(universeNum, itr.first.c_str(), itr.second);

      for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
        vecLoc.SetParameterValue(itr.first, itr.second, detItr);
      }
    } // end loop over systematic parameters to change

    MF_LOG_VERBATIM("RandomUniverseUtility")
        << systStr;
  }

  //......................................................................
  void RandomUniverseUtility::CreateLocationWithVariedOscParameters(cmf::Location const& nominalLoc,
                                                                    cmf::Location      & vecLoc,
                                                                    size_t               universeNum)
  {
    // loop over all the parameters and vary the ones we want
    double      val;
    std::string name;

    // set the location to the nominal location
    vecLoc = nominalLoc;

    // prepare some output for the random values in this universe
    std::string oscStr("Oscillation Parameters " + std::to_string(universeNum) + ": ");

    // get the names of the oscillation parameters. Only the baseline parameter
    // varies from one detector to the next.  By default we get the NOvA FD
    std::map<std::string, double> oscParVals;
    std::map<std::string, bool>   oscParFixed;
    for(auto const& itr : nominalLoc.OscillationParameters()){

      name = itr.Name();
      val  = itr.Value();

      MF_LOG_DEBUG("RandomUniverseUtility")
          << name
          << " is osc par: "
          << itr.IsOscPar()
          << " is fixed: "
          << itr.IsFixed();

      // The distance the neutrino travels should never be adjusted, so
      // don't bother putting it into oscParVals
      if(itr.IsOscPar() && (cmf::OscParm_t)itr.Key() != cmf::kL){
        if(!itr.IsFixed()){
          // we have a parameter that is not fixed, so let's figure out how to
          // vary it
          // first check if it is in the list of nuisance parameters, if not
          // look to see if it is constrained.
          if(itr.IsNuisance()){
            // treat the parameter as Gaussian distributed
            val = fRand.Gaus(itr.CentralValue(),
                             itr.Variance());

            MF_LOG_DEBUG("RandomUniverseUtility")
                << "nuisance parameter "
                << name
                << " val "
                << val
                << " "
                << itr.CentralValue()
                << " "
                << itr.Variance()
                << " "
                << itr.CentralValue() - itr.Variance();
          }
          else if(itr.IsConstrained()){
            val = fRand.Rndm() * (itr.UpperBound() - itr.LowerBound()) +
                  itr.LowerBound();

            MF_LOG_DEBUG("RandomUniverseUtility")
                << "non-nuisance parameter "
                << name
                << " val "
                << val
                << " "
                << itr.LowerBound()
                << " "
                << itr.UpperBound()
                << " "
                << itr.UpperBound() - itr.LowerBound();
          }

          oscStr += name + ": " + std::to_string(val) + "\t";
        }
        oscParVals.emplace(name, val);
        oscParFixed.emplace(name, itr.IsFixed());
      } // end if we have an oscillation parameter
    } // end loop over parameters

    //loop over oscParVals and set the oscillation parameter values that have
    // changed
    for(auto const& itr : oscParVals){

      // only add the parameter to the draw histogram if it isn't fixed
      if(!oscParFixed.find(itr.first)->second){
        fOscParamDrawHist->Fill(universeNum, itr.first.c_str(), itr.second);
      }

      // do this for both fixed and not fixed parameters to ensure we have
      // all oscillation parameters set
      for(auto const& detItr : cmf::SelectionUtility::Instance()->DetectorsToUse()){
        vecLoc.SetParameterValue(itr.first, itr.second, detItr);
      }
    } // end loop over oscillation parameters to change

    MF_LOG_DEBUG("RandomUniverseUtility")
      << vecLoc;
  }

  //......................................................................
  void RandomUniverseUtility::CreateLocationsWithVariedOscParameters(cmf::Location              const& nominalLoc,
                                                                     std::vector<cmf::Location>      & vecLocs)
  {
    // check that the vector is the correct size
    vecLocs.resize(fNumUniverses, nominalLoc);

    for(unsigned int u = 0; u < fNumUniverses; ++u){

      this->CreateLocationWithVariedOscParameters(nominalLoc, vecLocs[u], u);

      MF_LOG_DEBUG("RandomUniverseUtility")
          << " new point "
          << u;
    } // end loop over the number of universes

    // fOscParamDrawHist->LabelsDeflate("Y");
    // fOscParamDrawHist->LabelsOption("h", "Y");

  }


  // Function to fill the spectra for each meta data type
  //......................................................................
  void RandomUniverseUtility::FillPoissonSpectrum(cmf::Spectrum & spectrum,
                                                  cmf::Spectrum & poissonSpectrum)
  {
    // just copy the spectrum into the poissonSpectrum if we aren't configured for this
    if(!fPoissonFluctuate){
      poissonSpectrum = spectrum;
      return;
    }

    // first reset the poisson spectrum to be all 0's
    poissonSpectrum.resize(spectrum.size(), 0.);

    if(fPoissonByBin) this->PoissonSpectrumBinByBin(spectrum, poissonSpectrum);

    // for(size_t b = 0; b < poissonSpectrum.size(); ++b){
    //   MF_LOG_VERBATIM("RandomUniverseUtility")
    //  << "FillPoisonSpectrum bin "
    //  << b
    //  << " detector "
    //  << cmf::cDetType_Strings[cmf::CovarianceBinUtility::Instance()->BinToDetector(b)]
    //  << " without fluctuations "
    //  << spectrum[b]
    //  << " poisson expectation "
    //  << poissonSpectrum[b];
    // }

  }

  // poisson fluctuates each bin in the spectra
  //......................................................................
  void RandomUniverseUtility::PoissonSpectrumBinByBin(cmf::Spectrum & spectrum,
                                                      cmf::Spectrum & poissonSpectrum)
  {
    for(size_t b = 0; b < spectrum.size(); ++b){
      poissonSpectrum[b] = fRand.Poisson(spectrum[b]);

      MF_LOG_DEBUG("RandomUniverseUtility")
          << "PoissonSpectrumBinByBin, bin: "
          << b
          << " "
          << cmf::KeyToString(cmf::CovarianceBinUtility::Instance()->BinToKey(b))
          << " "
          << spectrum[b]
          << "->"
          << poissonSpectrum[b]
          << " ratio: "
          << 1. - poissonSpectrum[b] / spectrum[b];

    }
  }
        
  //......................................................................
  void RandomUniverseUtility::MakeUniverseHistograms(cmf::FakeUniverse & fakeU,
                                                     unsigned int        uNumber)
  {
    art::ServiceHandle<art::TFileService> tfs;

    auto unumStr = std::to_string(uNumber);
    // first off, let's make a new directory for this universe
    art::TFileDirectory dir = tfs->mkdir("Universe_" + unumStr);

    std::string title = "Universe " + unumStr;
    std::string name  = "PoissonUniverse_" + unumStr;

    title += ";Bin;Events";
    auto poissonHist = cmf::PlotUtilities::Instance()->MakeSpectrumHistogram(dir,
                                                                            name,
                                                                            title,
                                                                            fakeU.PoissonSpectrum());

    // comment out the next lines as FakeUniverse no longer carries the raw count
    // information.  Leave these here in case we want to pass in the vector from
    // elsewhere
//    name = "RawCount_" + unumStr;
//    TH1D *countHist = dir.make<TH1D>(name.data(),
//                                     title.data(),
//                                     fakeU.EventCounts().size() - 1,
//                                     0.,
//                                     fakeU.EventCounts().size() - 1);

    name = "AsimovUniverse_" + unumStr;
    auto asimovHist = cmf::PlotUtilities::Instance()->MakeSpectrumHistogram(dir, 
                                                                            name,
                                                                            title,
                                                                            fakeU.AsimovSpectrum());

    name = "RatioToAsimovUniverse_" + unumStr;
    auto ratioHist = cmf::PlotUtilities::Instance()->MakeSpectrumHistogram(dir,
                                                                           name,
                                                                           title,
                                                                           fakeU.PoissonSpectrum());

    // for(size_t b = 0; b < fakeU.PoissonSpectrum().size(); ++b){
    //   poissonHist->SetBinContent(b + 1, fakeU.PoissonSpectrum()[b]);
    //   asimovHist ->SetBinContent(b + 1, fakeU.AsimovSpectrum()[b] );
    //   //countHist  ->SetBinContent(b + 1, fakeU.EventCounts()[b]    );
    // }

    ratioHist->Divide(asimovHist);
  }

}