FitInAnaBinsBkgdEstimator.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 subtract off the signal prediction from the data
21  // Grab the TH1* from fData to get a hist with the right binning
22  Eigen::ArrayXd asig = fData.GetEigen(pot);
23  asig.setZero();
24 
25  assert(asig.size()-2==int(fCompPreds.size())
26  && "Number of kinematic bins given to FitInAnaBinsBkgdEstimator"
27  " does not match number of reconstructed bins in analysis");
28 
29  // Loop over Preds and fit, get the scale for signal, and store in hsig
30  for (int kinIdx = 0; kinIdx < int(fCompPreds.size()); kinIdx++){
31  SingleSampleExperiment expt(fCompPreds[kinIdx],*fDatas[kinIdx]);
32  std::vector<const ISyst*> systs = fCompPreds[kinIdx]->GetSysts();
33  // Now we can make a fitter with the vector of systs
34  MinuitFitter fit(&expt, {}, systs);
36  // Do a fit
37  fit.Fit(seed,MinuitFitter::kQuiet);
38  double sigscale = seed.GetShift(systs.back()); // Signal scale is last
39  double signom = fCompPreds[kinIdx]->GetSpectra().back()->Integral(pot);
40  asig[kinIdx+1] = (1+sigscale)*signom;
41  }
42  Spectrum sigest(std::move(asig),
44  pot, fData.Livetime());
45  return fData - sigest;
46  }
47 
48  void FitInAnaBinsBkgdEstimator::SaveTo(TDirectory *dir, const std::string& name) const
49  {
50  TDirectory *tmp = gDirectory;
51 
52  dir = dir->mkdir(name.c_str()); // switch to subdir
53  dir->cd();
54 
55  TObjString("FitInAnaBinsBkgdEstimator").Write("type");
56 
57  assert(fCompPreds.size()==fDatas.size());
58  // Save number of analysis bins as a TVectorD for when we load this
59  TVectorD nanabins(1);
60  nanabins[0] = fCompPreds.size();
61  nanabins.Write("nanabins");
62  for (int i = 0; i < int(fCompPreds.size()); i++){
63  std::string idxstr = std::to_string(i);
64  fCompPreds[i]->SaveTo(dir, "fCompPreds_"+idxstr);
65  fDatas[i] ->SaveTo(dir, "fDatas_" +idxstr);
66  }
67  fData.SaveTo(dir, "fData");
68 
69  dir->Write();
70  delete dir;
71 
72  tmp->cd();
73  }
74 
75  std::unique_ptr<FitInAnaBinsBkgdEstimator>
77  {
78  dir = dir->GetDirectory(name.c_str()); // switch to subdir
79  assert(dir);
80 
81  Spectrum* data = ana::LoadFrom<Spectrum>(dir, "fData").release();
82 
83  TVectorD* nanabins = (TVectorD*)dir->Get("nanabins");
84  assert(nanabins->GetNoElements()==1 && "N MCBkgds is not 1.");
85  std::vector<Spectrum*> datas;
86  std::vector<PredictionScaleComp*> ps;
87  for (int i = 0; i < (*nanabins)[0]; i++){
88  std::string idxstr = std::to_string(i);
89  ps.push_back(ana::LoadFrom<PredictionScaleComp>(dir, "fCompPreds_"+idxstr).release());
90  datas.push_back(ana::LoadFrom<Spectrum>(dir, "fDatas_"+idxstr).release());
91  }
92 
93  delete dir;
94 
95  return std::unique_ptr<FitInAnaBinsBkgdEstimator>(
96  new FitInAnaBinsBkgdEstimator(ps, datas, *data));
97  }
98 
99 }
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:268
const XML_Char * name
Definition: expat.h:151
std::vector< PredictionScaleComp * > fCompPreds
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
FitInAnaBinsBkgdEstimator(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)
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:200
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
#define pot
void SaveTo(TDirectory *dir, const std::string &name) const override
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
TDirectory * dir
Definition: macro.C:5
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
static std::unique_ptr< FitInAnaBinsBkgdEstimator > LoadFrom(TDirectory *dir, const std::string &name)
assert(nhit_max >=nhit_nbins)
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