FitInAnaBinsBkgdEstimator.cxx

Go to the documentation of this file.

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/XSec/FitInAnaBinsBkgdEstimator.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/MinuitFitter.h"

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

namespace ana
{
  Spectrum FitInAnaBinsBkgdEstimator::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 subtract off the signal prediction from the data
    // Grab the TH1* from fData to get a hist with the right binning
    Eigen::ArrayXd asig = fData.GetEigen(pot);
    asig.setZero();

    assert(asig.size()-2==int(fCompPreds.size())
           && "Number of kinematic bins given to FitInAnaBinsBkgdEstimator"
              " does not match number of reconstructed bins in analysis");

    // Loop over Preds and fit, get the scale for signal, and store in hsig
    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 sigscale = seed.GetShift(systs.back()); // Signal scale is last
      double signom   = fCompPreds[kinIdx]->GetSpectra().back()->Integral(pot);
      asig[kinIdx+1] = (1+sigscale)*signom;
    }
    Spectrum sigest(std::move(asig),
                    HistAxis(fData.GetLabels(), fData.GetBinnings()),
                    pot, fData.Livetime());
    return fData - sigest;
  }

  void FitInAnaBinsBkgdEstimator::SaveTo(TDirectory *dir, const std::string& name) const
  {
    TDirectory *tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("FitInAnaBinsBkgdEstimator").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, "fCompPreds_"+idxstr);
      fDatas[i]    ->SaveTo(dir, "fDatas_"    +idxstr);
    }
    fData.SaveTo(dir, "fData");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  std::unique_ptr<FitInAnaBinsBkgdEstimator>
  FitInAnaBinsBkgdEstimator::LoadFrom(TDirectory *dir, const std::string& name)
  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    Spectrum* data = ana::LoadFrom<Spectrum>(dir, "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;
    for (int i = 0; i < (*nanabins)[0]; i++){
      std::string idxstr = std::to_string(i);
      ps.push_back(ana::LoadFrom<PredictionScaleComp>(dir, "fCompPreds_"+idxstr).release());
      datas.push_back(ana::LoadFrom<Spectrum>(dir, "fDatas_"+idxstr).release());
    }

    delete dir;

    return std::unique_ptr<FitInAnaBinsBkgdEstimator>(
      new FitInAnaBinsBkgdEstimator(ps, datas, *data));
  }

}