NumuCC2p2hBkgdEstimator.cxx
Go to the documentation of this file.
5 
6 #include "TH1.h"
7 #include "TObjString.h"
8 #include "TVectorD.h"
9 
10 namespace ana
11 {
13  {
14  assert(fCompPreds.size()>0 && "Must give at least one MC int type");
15  assert(fCompPreds.size()==fDatas.size());
16 
17  // Set pot = 1 so that if there's a POT accounting error, it super obvious
18  double pot = 1;
19 
20  // We're going to sum up the background templates
21  // Grab the TH1* from fData to get a hist with the right binning
22  TH1 *hbkgd = fData.ToTH1(pot);
23  for (int i = 0; i <= hbkgd->GetNbinsX(); i++) hbkgd->SetBinContent(i,0);
24  assert(hbkgd->GetNbinsX()==int(fCompPreds.size())
25  && "Number of kinematic bins given to NumuCC2p2hBkgdEstimator"
26  " does not match number of reconstructed bins in analysis");
27 
28  // Loop over Preds and fit, get the scale for background, and store in hbkgd
29  for (int kinIdx = 0; kinIdx < int(fCompPreds.size()); kinIdx++){
30  SingleSampleExperiment expt(fCompPreds[kinIdx],*fDatas[kinIdx]);
31  std::vector<const ISyst*> systs = fCompPreds[kinIdx]->GetSysts();
32  // Now we can make a fitter with the vector of systs
33  MinuitFitter fit(&expt, {}, systs);
35  // Do a fit
36  fit.Fit(seed,MinuitFitter::kQuiet);
37  double bkgdsum=0;
38  for(int i=0; i<int(fCompPreds[kinIdx]->GetSpectra().size()-1);i++)//loop over non signal preds signal is last
39  {
40  double bkgdscale = seed.GetShift(systs[i]);
41  double bkgdnom = fCompPreds[kinIdx]->GetSpectra()[i]->Integral(pot);
42  bkgdsum+=(1+bkgdscale)*bkgdnom;
43  }
44  hbkgd->SetBinContent(kinIdx+1, bkgdsum);
45  }
46  Spectrum bkgdest(hbkgd, fData.GetLabels(), fData.GetBinnings(), pot, fData.Livetime());
47  delete hbkgd; // ret copies hbkgd, prevent memory leak here
48  return bkgdest;
49  }
50 
52  {}
53 
54  void NumuCC2p2hBkgdEstimator::SaveTo(TDirectory *dir) const
55  {
56  TDirectory *tmp = gDirectory;
57  dir->cd();
58 
59  TObjString("NumuCC2p2hBkgdEstimator").Write("type");
60 
61  assert(fCompPreds.size()==fDatas.size());
62  // Save number of analysis bins as a TVectorD for when we load this
63  TVectorD nanabins(1);
64  nanabins[0] = fCompPreds.size();
65  nanabins.Write("nanabins");
66  for (int i = 0; i < int(fCompPreds.size()); i++){
67  std::string idxstr = std::to_string(i);
68  fCompPreds[i]->SaveTo(dir->mkdir(("fCompPreds_"+idxstr).c_str()));
69  fDatas[i] ->SaveTo(dir->mkdir(("fDatas_" +idxstr).c_str()));
70  }
71  fData.SaveTo(dir->mkdir("fData"));
72 
73  tmp->cd();
74  }
75 
76  std::unique_ptr<NumuCC2p2hBkgdEstimator>
78  {
79  Spectrum* data = ana::LoadFrom<Spectrum>
80  (dir->GetDirectory("fData")).release();
81 
82  TVectorD* nanabins = (TVectorD*)dir->Get("nanabins");
83  assert(nanabins->GetNoElements()==1 && "N MCBkgds is not 1.");
84  std::vector<Spectrum*> datas;
85  std::vector<PredictionScaleComp*> ps;
86  dir->cd();
87  for (int i = 0; i < (*nanabins)[0]; i++){
88  std::string idxstr = std::to_string(i);
89  dir->cd(("fCompPreds_"+idxstr).c_str());
90  ps.push_back(ana::LoadFrom<PredictionScaleComp>(gDirectory).release());
91  dir->cd(("fDatas_"+idxstr).c_str());
92  datas.push_back(ana::LoadFrom<Spectrum>(gDirectory).release());
93  }
94  return std::unique_ptr<NumuCC2p2hBkgdEstimator>(
95  new NumuCC2p2hBkgdEstimator(ps, datas, *data));
96  }
97 
98 }
Spectrum Background() const override
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:268
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:209
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< PredictionScaleComp * > fCompPreds
Float_t tmp
Definition: plot.C:36
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const XML_Char const XML_Char * data
Definition: expat.h:268
expt
Definition: demo5.py:34
unsigned int seed
Definition: runWimpSim.h:102
void SaveTo(TDirectory *dir) const override
NumuCC2p2hBkgdEstimator(SpectrumLoaderBase &lMC, SpectrumLoaderBase &lDA, const HistAxis axisana, const HistAxis axisfit, const Cut sel, std::vector< Cut > kinsel, std::vector< std::pair< Cut, bool > > bkgdtypes, Cut isSig, const SystShifts &shiftmc=kNoShift, const SystShifts &shiftda=kNoShift, const Var &weimc=kUnweighted, const Var &weida=kUnweighted)
#define pot
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void GetSpectra(IPrediction *pred, osc::IOscCalc *calc, TH1 *&hSig, TH1 *&hNC, TH1 *&hCC, TH1 *&hBeam, TH1 *&hTotBkg)
Definition: nue_pid_effs.C:427
std::vector< Spectrum * > fDatas
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
static std::unique_ptr< NumuCC2p2hBkgdEstimator > LoadFrom(TDirectory *dir)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:267
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17