CMFRandomUniverses_plugin.cc

Go to the documentation of this file.

//
// Created by Brian Rebel on 10/29/19.
//

// Module to create fake data for multiple random universes
//

#include <memory>
#include <string>
#include <vector>
#include <map>
#include <set>

// 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 "CovarianceMatrixFit/core/ShifterAndWeighter.h"
#include "CovarianceMatrixFit/core/SpectrumTools.h"
#include "CovarianceMatrixFit/core/VarVals.h"
#include "CovarianceMatrixFit/dataProducts/FakeUniverse.h"
#include "CovarianceMatrixFit/dataProducts/Constants.h"
#include "CovarianceMatrixFit/dataProducts/StaticFuncs.h"
#include "CovarianceMatrixFit/modules/EventListManipulator.h"
#include "CovarianceMatrixFit/utilities/RandomUniverseUtility.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.h"
#include "CovarianceMatrixFit/utilities/SelectionUtility.h"

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

namespace cmf {

  class RandomUniverses : public art::ResultsProducer {
  public:
    explicit RandomUniverses(fhicl::ParameterSet const& pset);
    virtual ~RandomUniverses();

    // Plugins should not be copied or assigned.
    RandomUniverses(RandomUniverses const &)               = delete;
    RandomUniverses(RandomUniverses &&)                    = delete;
    RandomUniverses & operator = (RandomUniverses const &) = delete;
    RandomUniverses & operator = (RandomUniverses &&)      = 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& p);

  private:

    std::map<cmf::DetType_t, fhicl::ParameterSet> fManipulatorPars;   ///< parameters to make object for deserializing events
    fhicl::ParameterSet                                fSAWPars;           ///< parameters for the shifter and weighter
    bool                                               fIsFCJob;           ///< are we making Feldman Cousins points or not?
    unsigned int                                       fStartUniverse;     ///< which universe to start?
    TH2D*                                              fPoissonDraws;      ///< histogram recording the Poisson draws for each logical bin
    std::vector<cmf::Location>                         fUniverseLocs;      ///< vector of input points for the random universes
    std::vector<cmf::Location>                         fAsimovLocs;        ///< the Asimov input points, ie no systematic shifts
  };

  //......................................................................
  RandomUniverses::RandomUniverses(fhicl::ParameterSet const& pset)
   : fSAWPars     (pset.get<fhicl::ParameterSet>("ShifterAndWeighterParameters"))
   , fPoissonDraws(nullptr)
  {
    this->reconfigure(pset);

    LOG_DEBUG("RandomUniverses")
      << "ShifterAndWeighter parameters: "
      << fSAWPars.to_string();
    
    // produce a vector of FakeUniverses
    produces< std::vector<cmf::FakeUniverse> >();
    produces< std::vector<cmf::Location>     >();
  }

  //......................................................................
  RandomUniverses::~RandomUniverses()
  = default;

  //......................................................................
  void RandomUniverses::reconfigure(fhicl::ParameterSet const& p)
  {
    cmf::SelectionUtility::Instance()    ->Initialize(p.get< fhicl::ParameterSet >("SelectionsToUse", fhicl::ParameterSet() ));
    cmf::CovarianceBinUtility::Instance()->Initialize(p.get< fhicl::ParameterSet >("BinningToUse"                           ));
    cmf::ParameterUtility::Instance()    ->Initialize(p.get< fhicl::ParameterSet >("ParametersToUse"                        ));

    fManipulatorPars.emplace(cmf::kFARDET,  p.get<fhicl::ParameterSet>("FarDetEventListManipulator" ));
    fManipulatorPars.emplace(cmf::kNEARDET, p.get<fhicl::ParameterSet>("NearDetEventListManipulator"));
    fIsFCJob       = p.get<bool        >("IsFeldmanCousinsJob", false);
    fStartUniverse = p.get<unsigned int>("StartUniverse");

    // 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
    cmf::RandomUniverseUtility::Instance()->Initialize(p);
  }

  //......................................................................
  void RandomUniverses::beginJob()
  {
    cmf::Location nominalLoc = cmf::ParameterUtility::Instance()->ParametersAsLocation();

    LOG_DEBUG("RandomUniverse")
        << "nominal location is\n"
        << nominalLoc;

    // get a list of the oscillation parameter names
    std::set<std::string> oscParNames;
    for(auto const& itr : nominalLoc.FDLocation()){
      if(itr.second.IsOscPar())
        oscParNames.emplace(itr.first);
    }

    if(fIsFCJob){
      // use the RandomUniverseUtility to get our random values for the hidden
      // parameters in oscillation space first then use those oscillation
      // parameters to create systematically varied points as well
      cmf::RandomUniverseUtility::Instance()->CreateLocationsWithVariedOscParameters(nominalLoc,
                                                                                     fAsimovLocs);

      // now loop over the Asimov points and for each one create a template
      // to be used to get the systematic uncertainty parameter variations
      fUniverseLocs.resize(fAsimovLocs.size(), nominalLoc);

      for(size_t u = 0; u < fAsimovLocs.size(); ++u){
        // grab the oscillation parameters from this Asimov point
        // the oscillation parameter values are the same in both detectors
        for(auto const& itr : oscParNames){
          fUniverseLocs[u].SetParameterValue(itr, fAsimovLocs[u].NDLocation().find(itr)->second.Value(), cmf::kNEARDET);
          fUniverseLocs[u].SetParameterValue(itr, fAsimovLocs[u].FDLocation().find(itr)->second.Value(), cmf::kFARDET );
        }

        cmf::RandomUniverseUtility::Instance()->CreateLocationWithVariedSystematics(nominalLoc,
                                                                                    fUniverseLocs[u],
                                                                                    u);
      }
    }
    else{
      // In this mode we are only going to change the systematic parameters but not the
      // oscillation parameters so we only need 1 Asimov point
      fAsimovLocs.emplace_back(nominalLoc);

      // get the set of locations for our universes
      cmf::RandomUniverseUtility::Instance()->CreateLocationsWithVariedSystematics(nominalLoc,
                                                                                   fUniverseLocs);
    }

    for(size_t u = 0; u < fUniverseLocs.size(); ++u)
      LOG_VERBATIM("RandomUniverses")
          << " new location "
          << u
          << "\n"
          << fUniverseLocs[u];

    // now make some histograms
    art::ServiceHandle<art::TFileService> tfs;

    fPoissonDraws = tfs->make<TH2D>("PoissonDraws",
                                    ";Energy Bin;Poisson Draw",
                                    cmf::CovarianceBinUtility::Instance()->TotalBins(),
                                    0,
                                    cmf::CovarianceBinUtility::Instance()->TotalBins(),
                                    fUniverseLocs.size(),
                                    0.,
                                    fUniverseLocs.size());

  }

  //......................................................................
  // Method to clear out the collections of data products after the
  // writeResults method is called.
  void RandomUniverses::clear()
  {
  }

  //......................................................................
  void RandomUniverses::readResults(art::Results const& /*r*/)
  {
  }


  //......................................................................
  void RandomUniverses::writeResults(art::Results & r)
  {
    cmf::EventListColl mcLists;

    std::vector<cmf::Spectrum> spectrum       (fUniverseLocs.size());
    std::vector<cmf::Spectrum> poissonSpectrum(fUniverseLocs.size());
    std::vector<cmf::Spectrum> fakeCount      (fUniverseLocs.size());
    std::vector<cmf::Spectrum> count          (fUniverseLocs.size());
    std::vector<cmf::Spectrum> asimov         (fAsimovLocs.size()  );

    std::set<cmf::DetType_t> dets;
    for(auto const& mpItr : fManipulatorPars){

      mcLists.clear();

      dets.clear();
      dets.insert(mpItr.first);

      LOG_VERBATIM("RandomUniverses")
      << "Fill spectrum from "
      << cmf::cDetType_Strings[mpItr.first]
      << " starting universe "
      << fStartUniverse
      << " "
      << fUniverseLocs.size();

      // get the manipulator to deserialize the data and MC
      // we'll keep the MC lists around, but will drop the data lists as
      // soon as we leave this method
      cmf::EventListManipulator elm(mpItr.second);

      elm.Deserialize(mcLists, cmf::kMC, dets);
      auto const& exposureMap = elm.ExposureMap();

      // now make the spectra for all the universes
      unsigned int uni = 0;
      for(unsigned int u = fStartUniverse; u < fStartUniverse + fUniverseLocs.size(); ++u){

        uni = u - fStartUniverse;

        if(fIsFCJob){
          LOG_DEBUG("RandomUniverses")
              << "Filling FC Asimov Spectra";

          cmf::ShifterAndWeighter::Instance()->InitShiftsAndWeightsToUse(fAsimovLocs[uni],
                                                                         fSAWPars);

          // fill the asimov spectrum
          cmf::FillSpectrum(mcLists,
                            exposureMap,
                            asimov[uni],
                            count[uni]);
        }
        else if(!fIsFCJob && u == 0){
          LOG_DEBUG("RandomUniverses")
              << fIsFCJob
              << " "
              << u
              << " fill Asimov spectrum for "
              << fAsimovLocs.size()
              << " locations";

          cmf::ShifterAndWeighter::Instance()->InitShiftsAndWeightsToUse(fAsimovLocs[u],
                                                                         fSAWPars);
          // fill the asimov spectrum
          cmf::FillSpectrum(mcLists,
                            exposureMap,
                            asimov[u],
                            count[u]);

          LOG_DEBUG("RandomUniverses")
              << "Asimov spectrum filled "
              << asimov.size()
              << " "
              << count.size();

        }

        LOG_DEBUG("CMFRandomUniverses")
            << "Finished filling asimov spectrum";

        cmf::ShifterAndWeighter::Instance()->InitShiftsAndWeightsToUse(fUniverseLocs[uni],
                                                                       fSAWPars);

        LOG_DEBUG("CMFRandomUniverses")
            << "Filling universe spectrum";

        cmf::FillSpectrum(mcLists,
                          exposureMap,
                          spectrum[uni],
                          fakeCount[uni]);

        cmf::RandomUniverseUtility::Instance()->FillPoissonSpectrum(spectrum[uni], poissonSpectrum[uni]);

        for(size_t b = 0; b < poissonSpectrum[uni].size(); ++b){
          // the ND bins fill twice without the following check because
          // the ND portion of the spectrum is filled before the FD portion
          if(mcLists.begin()->ListMetaData().detector == cmf::CovarianceBinUtility::Instance()->BinToDetector(b))
            fPoissonDraws->Fill(b, uni, poissonSpectrum[uni][b]);
        }

        // for(size_t b = 0; b < spectrum[uni].size(); ++b){
        //   LOG_VERBATIM("RandomUniverses")
        //   << "universe "
            //   << u
            //   << " spectrum bin "
        //   << b
        //   << " fluctuated spectrum: "
        //   << poissonSpectrum[uni][b]
        //   << " unfluctuated/asimov: "
        //   << spectrum[uni][b]
        //   << " / "
        //   << asimov[0][b];
        // }

        LOG_VERBATIM("RandomUniverses")
            << " Finished universe "
            << uni
            << " for "
            << cmf::cDetType_Strings[mpItr.first];

      } // end loop to make fake universe spectra
    } // end loop over the manipulators

    // actually make the fake universes and their locations
    // we are storing these separately so that we don't have to
    // read in the locations (which use a lot of memory) unless we want to
    std::unique_ptr<std::vector<cmf::FakeUniverse>> fakeUniverses    = std::make_unique<std::vector<cmf::FakeUniverse>>();
    std::unique_ptr<std::vector<cmf::Location>>     fakeUniverseLocs = std::make_unique<std::vector<cmf::Location>>();

    for(size_t u = 0; u < fUniverseLocs.size(); ++u){
      cmf::Spectrum asimovSpec(asimov[0]);
      if(fIsFCJob) asimovSpec.swap(asimov[u]);
      fakeUniverses->emplace_back(fUniverseLocs[u],
                                  asimovSpec,
                                  poissonSpectrum[u]);
      fakeUniverseLocs->emplace_back(fUniverseLocs[u]);

      LOG_DEBUG("RandomUniverses")
          << "added fake universe "
          << u;
      
      // make histograms for the fake spectra and the total spectrum to compare
      // them against each other.
      cmf::RandomUniverseUtility::Instance()->MakeUniverseHistograms(fakeUniverses->back(),
                                                                     u + fStartUniverse);

      LOG_DEBUG("RandomUniverses")
          << fakeUniverses->back();
    }

    r.put(std::move(fakeUniverses   ));
    r.put(std::move(fakeUniverseLocs));
  } // end writeResults

  DEFINE_ART_RESULTS_PLUGIN(RandomUniverses)

} // end cmf namespace