CMFSpectraMaker_module.cc

Go to the documentation of this file.

// Spectra making module

#include <cmath>
#include <string>
#include <vector>

#include "TFile.h"
#include "TTree.h"
#include "TH1.h"

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

// CMF includes
#include "CovarianceMatrixFit/core/VarVals.h"
#include "CovarianceMatrixFit/core/ShifterAndWeighter.h"
#include "CovarianceMatrixFit/core/SpectrumTools.cxx"
#include "CovarianceMatrixFit/dataProducts/Structs.h"
#include "CovarianceMatrixFit/modules/EventListManipulator.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.h"
#include "CovarianceMatrixFit/utilities/PlotUtilities.h"

namespace cmf {
  class CMFSpectraMaker : public art::EDAnalyzer {
    public:
      explicit CMFSpectraMaker(fhicl::ParameterSet const& pset);
      ~CMFSpectraMaker() = default;

      // Required fucntions
      void analyze(art::Event const& evt) override;
      void beginJob()                     override;
      void endJob()                       override;

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

    private:

      void Initialize();
      void FillHistograms(art::TFileDirectory& td,
                          int                  sigmaVal,
                          cmf::EventListColl & eventLists,
                          cmf::ExposureMap   & exposureMap);

      // Member data
      fhicl::ParameterSet fManipulatorParameters; ///< event list manipulator parameters
      fhicl::ParameterSet fSAWParameters;         ///< shifter and weighter parameters

      art::ServiceHandle<art::TFileService> fTFS;

      cmf::Location                                   loc;
      std::vector<std::string>                        fITypeStrings;
      std::map<std::string, cmf::InteractionType_t>   ITypes;
      std::vector<int>                                fSigmaValues;
      bool                                            fApplyExposureNorm;
  };

  //----------------------------------------------------------------//
  CMFSpectraMaker::CMFSpectraMaker(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
  {
    this->reconfigure(pset);
  }

  //----------------------------------------------------------------//
  void CMFSpectraMaker::reconfigure(fhicl::ParameterSet const& pset)
  {
    fManipulatorParameters = pset.get< fhicl::ParameterSet >("EventListManipulator"    );
    fSAWParameters         = pset.get< fhicl::ParameterSet >("ShifterAndWeighterConfig");
    fSigmaValues           = pset.get< std::vector<int>    >("SigmaValues", {0}        );
    fApplyExposureNorm     = pset.get< bool                >("ApplyExposureNorm", true );

    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"                       ));
    cmf::PlotUtilities::Instance()       ->Initialize(pset.get< fhicl::ParameterSet >("PlotConfig"                            ));
    loc = cmf::ParameterUtility::Instance()->ParametersAsLocation();

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

    fITypeStrings = pset.get< std::vector<std::string>> ("InteractionTypes");

    ITypes["Total"] = cmf::kUnknownInteraction;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NuMuCC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNuMuCC]] = cmf::kNuMuCC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NuMuBarCC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNuMuBarCC]] = cmf::kNuMuBarCC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NuECC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNuECC]] = cmf::kNuECC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NuEBarCC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNuEBarCC]] = cmf::kNuEBarCC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NuTauCC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNuTauCC]] = cmf::kNuTauCC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NuTauBarCC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNuTauBarCC]] = cmf::kNuTauBarCC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "NC") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kNC]] = cmf::kNC;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "CosmicMuon") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kCosmicMuon]] = cmf::kCosmicMuon;
    if (std::find(fITypeStrings.begin(), fITypeStrings.end(), "RockMuon") != fITypeStrings.end())
      ITypes[cInteractionType_Strings[cmf::kRockMuon]] = cmf::kRockMuon;

  }

  //----------------------------------------------------------------//
  void CMFSpectraMaker::beginJob()
  {
  }

  //----------------------------------------------------------------//
  void CMFSpectraMaker::analyze(art::Event const& evt)
  {
  }

  //----------------------------------------------------------------//
  void CMFSpectraMaker::FillHistograms(art::TFileDirectory& td,
                                       int                  sigmaValue,
                                       cmf::EventListColl & eventLists,
                                       cmf::ExposureMap   & exposureMap)
  {
    for (auto& iType : ITypes){

      std::vector<float> spectrum;
      std::vector<float> count;

      if (iType.second == cmf::kUnknownInteraction)
        cmf::FillSpectrum(eventLists,
                          exposureMap,
                          spectrum,
                          count,
                          fApplyExposureNorm);
      else
        cmf::FillSpectrum(eventLists,
                          exposureMap,
                          spectrum,
                          count,
                          iType.second,
                          fApplyExposureNorm);

      std::map<long, TH1D*> energySpecMap;

      std::string preName = "";
      if (sigmaValue != 0)
        preName = std::to_string(sigmaValue)+"Sigma";

      cmf::PlotUtilities::Instance()->MakeEnergySpectraFromBins(spectrum,
                                                                count,
                                                                energySpecMap,
                                                                iType.first,
                                                                td,
                                                                preName);
    }
  }

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

    cmf::EventListColl eventLists;
    manipulator.Deserialize(eventLists, cmf::kMC);

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

    // this allows us to get systematically spectra by setting
    // the systematic sigma
    auto parNames = cmf::ParameterUtility::Instance()->SysParNames();
    if (parNames.size() > 0){
      for (auto& parItr : parNames) {

        art::TFileDirectory td = fTFS->mkdir(parItr, parItr);

        for (size_t is = 0; is < fSigmaValues.size(); ++is){

          loc.SetParameterValue(parItr, fSigmaValues.at(is), cmf::kNEARDET);
          loc.SetParameterValue(parItr, fSigmaValues.at(is), cmf::kFARDET );

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

          this->FillHistograms(td, fSigmaValues.at(is), eventLists, exposureMap);
        }
      }
    }
    else{
      art::TFileDirectory td = fTFS->mkdir("nosyst", "Not systematically shifted");
      this->FillHistograms(td, 0, eventLists, exposureMap);
    }

  }

}

DEFINE_ART_MODULE(cmf::CMFSpectraMaker)