NumuCC2p2hBkgdEstimator.cxx

Go to the documentation of this file.

#include "CAFAna/Core/LoadFromFile.h"
#include "NDAna/numucc_2p2h/NumuCC2p2hBkgdEstimator.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/MinuitFitter.h"

#include "TH1.h"
#include "TObjString.h"
#include "TVectorD.h"

namespace ana
{
  Spectrum NumuCC2p2hBkgdEstimator::Background() const
  {
    assert(fCompPreds.size()>0 && "Must give at least one MC int type");
    assert(fCompPreds.size()==fDatas.size());

    // Set pot = 1 so that if there's a POT accounting error, it super obvious
    double pot = 1;

    // We're going to sum up the background templates
    // Grab the TH1* from fData to get a hist with the right binning
    TH1 *hbkgd = fData.ToTH1(pot);
    for (int i = 0; i <= hbkgd->GetNbinsX(); i++) hbkgd->SetBinContent(i,0);
    assert(hbkgd->GetNbinsX()==int(fCompPreds.size())
           && "Number of kinematic bins given to NumuCC2p2hBkgdEstimator"
              " does not match number of reconstructed bins in analysis");

    // Loop over Preds and fit, get the scale for background, and store in hbkgd
    for (int kinIdx = 0; kinIdx < int(fCompPreds.size()); kinIdx++){
      SingleSampleExperiment expt(fCompPreds[kinIdx],*fDatas[kinIdx]);
      std::vector<const ISyst*> systs = fCompPreds[kinIdx]->GetSysts();
      // Now we can make a fitter with the vector of systs
      MinuitFitter fit(&expt, {}, systs);
      SystShifts seed;
      // Do a fit
      fit.Fit(seed,MinuitFitter::kQuiet);
      double bkgdsum=0;
      for(int i=0; i<int(fCompPreds[kinIdx]->GetSpectra().size()-1);i++)//loop over non signal preds signal is last
        {
          double bkgdscale = seed.GetShift(systs[i]); 
          double bkgdnom   = fCompPreds[kinIdx]->GetSpectra()[i]->Integral(pot);
          bkgdsum+=(1+bkgdscale)*bkgdnom;
        }
      hbkgd->SetBinContent(kinIdx+1, bkgdsum);
    }
    Spectrum bkgdest(hbkgd, fData.GetLabels(), fData.GetBinnings(), pot, fData.Livetime());
    delete hbkgd; // ret copies hbkgd, prevent memory leak here
    return bkgdest;
  }

  NumuCC2p2hBkgdEstimator::~NumuCC2p2hBkgdEstimator()
  {}

  void NumuCC2p2hBkgdEstimator::SaveTo(TDirectory *dir) const
  {
    TDirectory *tmp = gDirectory;
    dir->cd();

    TObjString("NumuCC2p2hBkgdEstimator").Write("type");

    assert(fCompPreds.size()==fDatas.size());
    // Save number of analysis bins as a TVectorD for when we load this
    TVectorD nanabins(1);
    nanabins[0] = fCompPreds.size();
    nanabins.Write("nanabins");
    for (int i = 0; i < int(fCompPreds.size()); i++){
      std::string idxstr = std::to_string(i);
      fCompPreds[i]->SaveTo(dir->mkdir(("fCompPreds_"+idxstr).c_str()));
      fDatas[i]    ->SaveTo(dir->mkdir(("fDatas_"    +idxstr).c_str()));
    }
    fData.SaveTo(dir->mkdir("fData"));

    tmp->cd();
  }

  std::unique_ptr<NumuCC2p2hBkgdEstimator>
    NumuCC2p2hBkgdEstimator::LoadFrom(TDirectory *dir)
  {
    Spectrum* data = ana::LoadFrom<Spectrum>
      (dir->GetDirectory("fData")).release();

    TVectorD* nanabins = (TVectorD*)dir->Get("nanabins");
    assert(nanabins->GetNoElements()==1 && "N MCBkgds is not 1.");
    std::vector<Spectrum*> datas;
    std::vector<PredictionScaleComp*> ps;
    dir->cd();
    for (int i = 0; i < (*nanabins)[0]; i++){
      std::string idxstr = std::to_string(i);
      dir->cd(("fCompPreds_"+idxstr).c_str());
      ps.push_back(ana::LoadFrom<PredictionScaleComp>(gDirectory).release());
      dir->cd(("fDatas_"+idxstr).c_str());
      datas.push_back(ana::LoadFrom<Spectrum>(gDirectory).release());
    }
    return std::unique_ptr<NumuCC2p2hBkgdEstimator>(
      new NumuCC2p2hBkgdEstimator(ps, datas, *data));
  }

}