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/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "cetlib_except/exception.h"

#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/RandomUniverseUtility.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& itr : nominalLoc.NDLocation()){
      if(itr.second.IsOscPar() ||
         itr.second.IsFixed()    ) 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.first)) ? fRand.Rndm() : fRand.Gaus();

      systParVals.emplace(itr.first, sigFrac);
    }
    for(auto const& itr : nominalLoc.FDLocation()){
      if(itr.second.IsOscPar()           ||
         itr.second.IsFixed()            ||
         systParVals.count(itr.first) > 0 ) continue;

      sigFrac = (this->UniformDrawParameters(itr.first)) ? fRand.Rndm() : fRand.Gaus();

      systParVals.emplace(itr.first, 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);

      vecLoc.SetParameterValue(itr.first, itr.second, cmf::kNEARDET);
      vecLoc.SetParameterValue(itr.first, itr.second, cmf::kFARDET );
    } // end loop over systematic parameters to change

    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. The parameters are the same
    // in both the ND and FD
    std::map<std::string, double> oscParVals;
    for(auto const& itr : nominalLoc.FDLocation()){

      name = itr.first;
      val  = itr.second.Value();

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

      // The distance the neutrino travels should never be adjusted, so
      // don't bother putting it into oscParVals
      if(itr.second.IsOscPar() && (cmf::OscParm_t)itr.second.Key() != cmf::kL){
        if(!itr.second.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.second.IsNuisance()){
            // treat the parameter as Gaussian distributed
            val = fRand.Gaus(itr.second.CentralValue(),
                             itr.second.Variance());

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

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

          oscStr += name + ": " + std::to_string(val) + "\t";
        }
        oscParVals.emplace(name, val);
      } // 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(!vecLoc.FDLocation().find(itr.first)->second.IsFixed()){
        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
      vecLoc.SetParameterValue(itr.first, itr.second, cmf::kNEARDET);
      vecLoc.SetParameterValue(itr.first, itr.second, cmf::kFARDET );
    } // end loop over oscillation parameters to change

    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);

      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(std::vector<float> & spectrum,
                                                  std::vector<float> & 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){
    //   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(std::vector<float> & spectrum,
                                                      std::vector<float> & poissonSpectrum)
  {
    for(size_t b = 0; b < spectrum.size(); ++b){
      poissonSpectrum[b] = fRand.Poisson(spectrum[b]);

      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";
    TH1D *poissonHist = dir.make<TH1D>(name.data(),
                                       title.data(),
                                       fakeU.PoissonSpectrum().size() - 1,
                                       0.,
                                       fakeU.PoissonSpectrum().size() - 1);
    poissonHist->Sumw2();

    // 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;
    TH1D *asimovHist = dir.make<TH1D>(name.data(),
                                      title.data(),
                                      fakeU.AsimovSpectrum().size() - 1,
                                      0.,
                                      fakeU.AsimovSpectrum().size() - 1);
    asimovHist->Sumw2();

    name = "RatioToAsimovUniverse_" + unumStr;
    TH1D *ratioHist = dir.make<TH1D>(name.data(),
                                     title.data(),
                                     fakeU.AsimovSpectrum().size() - 1,
                                     0.,
                                     fakeU.AsimovSpectrum().size() - 1);
    ratioHist->Sumw2();

    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(poissonHist, asimovHist);
  }

}