test_predictionscalecomp.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
4 
8 
12 
16 
17 #include "TH1.h"
18 #include "THStack.h"
19 #include "TCanvas.h"
20 #include "TLegend.h"
22 
23 #include <vector>
24 #include <string>
25 
26 
27 using namespace ana;
28 
29 const Cut kIsNDNumuCC([](const caf::SRProxy* sr)
30  {
31  if(sr->mc.nnu == 0) return false;
32  if(!sr->mc.nu[0].iscc) return false;
33  return (sr->mc.nu[0].pdg==14);
34  });
35 
36 
37 // This macro aims to demonstrate the PredictionScaleComp class, namely how
38 // it can be used to fit different components to a distribution in data.
39 // There's a constructor in CAFAna's MinuitFitter doesn't need an osccalc. We can
40 // use that to hook into Minuit so that analyzers won't have to spend time
41 // re-inventing the wheel. We're going to fit the SA ND numi spectrum to data
42 // in HadE just for example with MC separated by GENIE mode (from Kirk's
43 // initial studies. That's the first type of this anlaysis I remember.
44 
45 ///////////////////////////////////////////////////////////////////////////////
46 void MakeStackPlot(std::vector<TH1*> hists, TH1* hda, std::string title)
47 {
48  // Make a pretty stacked plot for the MC with data overlain
49 
50  // We have 6 different MC spectra that we need to assign fill colors to.
51  assert(hists.size()==6 && "Check the number of MC hists you're fitting");
52  hists[0]->SetFillColor(kGray+2);
53  hists[1]->SetFillColor(kRed);
54  hists[2]->SetFillColor(kYellow);
55  hists[3]->SetFillColor(kBlue);
56  hists[4]->SetFillColor(kMagenta);
57  hists[5]->SetFillColor(kGreen+2);
58 
59  THStack *hs = new THStack();
60  for (int i = 0; i < int(hists.size()); i++)
61  hs->Add(hists[i]);
62 
63  double max = std::max(hda->GetMaximum(), hs->GetMaximum());
64 
65  hda->SetMaximum(1.25*max);
66  hda->GetXaxis()->CenterTitle();
67  hda->GetYaxis()->CenterTitle();
68  hda->SetTitle(title.c_str());
69 
70  TCanvas *can = new TCanvas(UniqueName().c_str());
71  hda->Draw();
72  hs->Draw("hist,same");
73  hda->Draw("same");
74 
75  TLegend *leg = new TLegend(0.6,0.4,0.8,0.8);
76  leg->AddEntry(hists[0],"#nu_{#mu} CC Bkgd","f");
77  leg->AddEntry(hists[1],"#nu_{#mu} CC QE", "f");
78  leg->AddEntry(hists[2],"#nu_{#mu} CC MEC", "f");
79  leg->AddEntry(hists[3],"#nu_{#mu} CC Res", "f");
80  leg->AddEntry(hists[4],"#nu_{#mu} CC DIS", "f");
81  leg->AddEntry(hists[5],"#nu_{#mu} CC Coh", "f");
82  leg->AddEntry(hda, "Data", "le");
83  leg->Draw();
84 }
85 
86 ///////////////////////////////////////////////////////////////////////////////
87 void MakePlots(std::vector<Spectrum*> spects, Spectrum data,
88  std::vector<const ISyst*> systs, SystShifts shift)
89 {
90  double pot = data.POT();
91  // Translate MC spectra into TH1*, and scale by the fit output
92  std::vector<TH1*> nomhists, fithists;
93  for (int i = 0; i < int(spects.size()); i++){
94  nomhists.push_back(spects[i]->ToTH1(pot));
95  fithists.push_back(spects[i]->ToTH1(pot));
96  // Remember we didn't fit the first component of this syst vector
97  if (i!=0){
98  double val = shift.GetShift(systs[i-1]);
99  fithists[i]->Scale(1+val);
100  }
101  }
102  TH1* hda1 = data.ToTH1(pot);
103  TH1* hda2 = data.ToTH1(pot);
104 
105  MakeStackPlot(nomhists, hda1, "Nominal MC");
106  MakeStackPlot(fithists, hda2, "Fit MC");
107 }
108 
109 ///////////////////////////////////////////////////////////////////////////////
111 {
112  // Set up an experiment - most basic tool for setting up a CAFAna fit.
113  // It knows about the data that you observe, and a MC prediction
114  SingleSampleExperiment expt(&pred, data);
115 
116  // MinuitFitter will fit out the scale factor by treating the normalization as
117  // ISysts. We need to get the vector, with entries corresponding to each
118  // entry in your truthcuts vector to preform the fit.
119  std::vector<const ISyst*> systs = pred.GetSysts();
120 
121  MinuitFitter fit(&expt, {}, systs);
122  SystShifts seed; // fit needs this object to store the answer
123  fit.Fit(seed); // Run the fit
124 
125  // Make before / after data / mc plots for the fit
126  MakePlots(pred.GetSpectra(), data, systs, seed);
127 }
128 
129 ///////////////////////////////////////////////////////////////////////////////
131 {
132 
133  // Set up loaders and apply spill cuts. Get the files from prod's concats.
134  const std::string nameMC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
135  SpectrumLoader lNDMC(nameMC);
137 
138  const std::string nameDA = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
139  SpectrumLoader lNDDA(nameDA);
141 
142  // Set up your prediction and its inputs
143  HistAxis axis("Had E (GeV)",Binning::Simple(30,0,3),kHadE);
144  // We need to make sure that every event belongs to some category, otherwise
145  // PredictionScaleComp will give you an assert. Throw all NC / other bkgd
146  // into comp.
147  Cut comp = !kIsNDNumuCC||(!kIsQE&&!kIsDytmanMEC&&!kIsRes&&!kIsDIS&&!kIsCoh);
148  // The `bool' part of this pair dictates whether the component is floated
149  // in the fit. So, we hold the non-numu CC background fixed, along with
150  // the coherent contribution - they're small and we don't want to let the
151  // fitter confuse their shape with syst variations in the big components
152  std::vector< std::pair<Cut,bool> > truthcuts;
153  truthcuts.push_back(std::pair<Cut,bool>(!kIsNDNumuCC,false));
154  truthcuts.push_back(std::pair<Cut,bool>(kIsNDNumuCC&&kIsQE,true));
155  truthcuts.push_back(std::pair<Cut,bool>(kIsNDNumuCC&&kIsDytmanMEC,true));
156  truthcuts.push_back(std::pair<Cut,bool>(kIsNDNumuCC&&kIsRes,true));
157  truthcuts.push_back(std::pair<Cut,bool>(kIsNDNumuCC&&kIsDIS,true));
158  truthcuts.push_back(std::pair<Cut,bool>(kIsNDNumuCC&&kIsCoh,false));
159 
160  Cut sel = kNumuND; // This is your analysis cut. Here, the SA numu CC sel
162  PredictionScaleComp pred(lNDMC, axis, sel, truthcuts,
163  kNoShift, wei);
164 
165  // Now make your data spectrum
166  Spectrum data(lNDMC, axis, sel);
167 
168  // Once we've set up our prediction and our Data spectrum, press Go
169  lNDMC.Go();
170  lNDDA.Go();
171 
172  // This function sets up experiment + runs fitter, then makes some plots
173  RunFitter(pred, data);
174 
175  TFile *out = new TFile("out_predictionscalecomp.root","recreate");
176  pred.SaveTo(out, "pred");
177  data.SaveTo(out, "data");
178  out->Close();
179 }
const Var kHadE
Definition: NumuVars.h:23
const Cut kIsQE
Definition: TruthCuts.cxx:104
T max(const caf::Proxy< T > &a, T b)
const Cut kIsNDNumuCC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(!sr->mc.nu[0].iscc) return false;return(sr->mc.nu[0].pdg==14);})
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
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
const Cut kIsRes
Definition: TruthCuts.cxx:111
Proxy for caf::StandardRecord.
Definition: SRProxy.h:1993
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:570
caf::Proxy< short int > nnu
Definition: SRProxy.h:569
const Cut kIsDIS
Definition: TruthCuts.cxx:118
void test_predictionscalecomp()
void SetSpillCut(const SpillCut &cut)
TString hists[nhists]
Definition: bdt_com.C:3
virtual void SaveTo(TDirectory *dir, const std::string &name) const override
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
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
expt
Definition: demo5.py:34
unsigned int seed
Definition: runWimpSim.h:102
void MakeStackPlot(std::vector< TH1 * > hists, TH1 *hda, std::string title)
const Cut kNumuND
Definition: NumuCuts.h:55
#define pot
T GetShift(const ISyst *syst) const
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
double POT() const
Definition: Spectrum.h:231
const SystShifts kNoShift
Definition: SystShifts.h:115
std::vector< const ISyst * > GetSysts() const
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2005
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void RunFitter(PredictionScaleComp pred, Spectrum data)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Spectrum * > GetSpectra() const
assert(nhit_max >=nhit_nbins)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Cut kIsCoh
Definition: TruthCuts.cxx:133
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
static constexpr Double_t sr
Definition: Munits.h:164
Compare a single data spectrum to the MC + cosmics expectation.
void MakePlots(std::vector< Spectrum * > spects, Spectrum data, std::vector< const ISyst * > systs, SystShifts shift)
Prediction broken down into arbitrary components whose scales can be varied independently.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17