demo_trivial_xsec_analysis.C
Go to the documentation of this file.
1 #pragma once
2 
3 // ND cross section framework
7 #include "CAFAna/XSec/Plotting.h"
8 
9 // Borrow cuts and binning from another analysis
12 
13 // CAFAna includes
14 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/Vars/Vars.h"
17 
18 // Multiverse includes
20 #include "CAFAna/Vars/XsecTunes.h"
24 
25 // ROOT includes
26 #include "TColor.h"
27 
28 using namespace ana;
29 
30 // run with cafe --ndxsec --numuccinc -bq -l 1 demo_trivial_xsec_analysis.C
31 void demo_trivial_xsec_analysis(bool fill_spectra = false)
32 {
33  // nominal mc and fake data dataset definitions
34  std::string fdata = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_v1 with stride 10 offset 0";
36 
37  // systematic samples
38  std::string fcalibneg = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-neg-offset_v1_good_seventh with stride 10 offset 0";
39  std::string fcalibpos = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-xyview-pos-offset_v1_good_seventh with stride 10 offset 0";
40  std::string fcalibshape = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1_good_seventh with stride 10 offset 0";
41  std::string flightupcalibdown = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1_good_seventh with stride 10 offset 0";
42  std::string flightdowncalibup = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1_good_seventh with stride 10 offset 0";
43  std::string fcherenkov = "def_snapshot prod_caf_R17-11-14-prod4reco.remid-hotfix.b_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1_good_seventh with stride 10 offset 0";
44 // std::string testdir = "/nova/app/users/ddoyle/testfiles/";
45 // std::string fmc = testdir + "neardet_genie_nonswap_genierw_fhc_v08_2157_r00010565_s07_c003_R17-11-14-prod4reco.d_v1_20170322_204739_sim.caf.root";
46 // std::string fdata = testdir + "neardet_genie_nonswap_genierw_fhc_v08_2157_r00010565_s07_c003_R17-11-14-prod4reco.d_v1_20170322_204739_sim.caf.root";
47 // std::string fcalibneg = testdir + "neardet_genie_nonswap_genierw_fhc_v08_1151_r00012018_s07_c000_R17-11-14-prod4reco.g_v1_20170322_204739_sim_calib-shift-nd-xyview-neg-offset.caf.root";
48 // std::string fcalibpos = testdir + "neardet_genie_nonswap_genierw_fhc_v08_1558_r00011380_s14_c000_R17-11-14-prod4reco.g_v1_20170322_204739_sim_calib-shift-nd-xyview-pos-offset.caf.root";
49 // std::string fcalibshape = testdir + "neardet_genie_nonswap_genierw_fhc_v08_1566_r00010923_s17_c000_R17-11-14-prod4reco.g_v1_20170322_204739_sim_calib-shift-nd-func.caf.root";
50 // std::string flightupcalibdown = testdir + "neardet_genie_nonswap_genierw_fhc_v08_1207_r00011492_s09_c000_R17-11-14-prod4reco.g_v1_20170322_204739_sim_lightmodel-lightup-calibdown.caf.root";
51 // std::string flightdowncalibup = testdir + "neardet_genie_nonswap_genierw_fhc_v08_1309_r00011444_s17_c000_R17-11-14-prod4reco.g_v1_20170322_204739_sim_lightmodel-lightdown-calibup.caf.root";
52 // std::string fcherenkov = testdir + "neardet_genie_nonswap_genierw_fhc_v08_1549_r00011393_s04_c000_R17-11-14-prod4reco.g_v1_20170322_204739_sim_ckv-proton-shift-down.caf.root";
53 
54  std::vector<XSecSysts::Syst_t> systs = XSecSysts::GetAllSystematics();
55  // same order as GetAllSystematics()
56  std::vector<std::string> syst_samples = {fmc,
57  fcalibneg,
58  fcalibpos,
59  fcalibshape,
60  flightupcalibdown,
61  flightdowncalibup,
62  fcherenkov};
63  if(fill_spectra) {
64  std::vector<SpectrumLoader*> loaders;
65  std::vector<const Var *> weights;
66  std::vector<const NuTruthVar *> nt_weights;
67  std::vector<const SystShifts *> shifts;
68  for(int isyst = 0; isyst < (int) syst_samples.size(); isyst++) {
69  // Define the loader, weights, and syst shift for each systematic sample
70  loaders.push_back(new SpectrumLoader(syst_samples[isyst])); // loaders
71  weights.push_back(new Var(kPPFXFluxCVWgt * kXSecCVWgt2018)); // slice weights
72  nt_weights.push_back(new NuTruthVar(kPPFXFluxCVWgtST * kXSecCVWgt2018_smallerDISScale_NT)); // truth weights
73  shifts.push_back(&kNoShift); // syst shifts
74  loaders.back()->SetSpillCut(kStandardSpillCuts);
75  }
76 
77  // Loader for data
78  SpectrumLoader * data_loader = new SpectrumLoader(fdata);
79  data_loader->SetSpillCut(kStandardSpillCuts);
80 
81  // TrivialSignalEstimator uses a TrivialBkgdEstimator.
82  // TrivialBkgdEstimator subtracts a sum of MC distributions
83  // from the selected data.
84  std::vector<Cut> bkgd_components = {kIsNueCC,
87  kIsNCCut,
89  !kIsNueCC && !kIsAntiNueCC && !kIsAntiNumuCC &&
90  !kIsNCCut && !kAllNumuCCCutsTrueNotFiducial &&
92 
93  // Binning for flux and in this case, also the analysis
94  const Binning enu_binning = Binning::Simple(30, 0, 6);
95 
96  // CrossSectionAnalysis uses HistAxes
97  const HistAxis * reco_axis = new HistAxis("Reco E_{#nu} (GeV)",
98  enu_binning,
99  kNumuE2018);
100  const NuTruthHistAxis * truth_axis = new NuTruthHistAxis("True E_{#nu} (GeV)",
101  enu_binning,
102  kTrueE_NT);
103  // first do initial setup of the analysis object
104  int pdg = 14;
105  bool is_differential = false;
106  TrivialCrossSectionAnalysis xsec_ana(kAllNumuCCCuts, // selection cut
107  reco_axis, // reco axis
108  kIsTrueSigST, // signal cut
109  truth_axis, // truth axis
110  vtxmax, // fiducial max
111  vtxmin, // fiducial min
112  &enu_binning, // binning for flux
113  bkgd_components, // background components
114  pdg, // pdg for flux
115  is_differential); //
116 
117  // now setup all of the spectra we'll be using in the analysis
118  for(int isyst = 0; isyst < (int) systs.size(); isyst++) {
119  xsec_ana.InitSystematicSpectra(systs[isyst],
120  loaders[isyst],
121  weights[isyst],
122  nt_weights[isyst],
123  shifts[isyst]);
124  }
125 
126  // Setup the genie and ppfx multiverses
127  std::vector<Var> ppfx_weights = GetkPPFXFluxUnivWgt();
128  std::vector<NuTruthVar> ppfx_weights_nt = GetkPPFXFluxUnivWgtST();
129  std::vector<Var> ppfx_weights_sparse(&ppfx_weights[0], &ppfx_weights[50]);
130  std::vector<NuTruthVar> ppfx_weights_sparse_nt(&ppfx_weights_nt[0], &ppfx_weights_nt[50]);
131  GenieMultiverseParameters genie_multiverse(50, getAllXsecSysts_2018());
132  std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();
133 
134  // Give them to the CrossSectionAnalysis object
135  xsec_ana.InitPPFXUniverses(loaders[0],
136  ppfx_weights_sparse,
137  ppfx_weights_sparse_nt);
138  xsec_ana.InitGenieUniverses(loaders[0],
139  genie_shifts);
140 
141  // and finally we'll add the data
142  xsec_ana.InitDataSpectrum(data_loader);
143 
144  // analysis setup is done, now we can let the loaders go
145  for(int iload = 0; iload < (int) loaders.size(); iload++) {
146  loaders[iload]->Go();
147  }
148  data_loader->Go();
149 
150  // save the analysis once the spectra are done loading
151  TFile * out = new TFile("demo_trivial_xsec_analysis.root", "recreate");
152  xsec_ana.SaveTo(out->mkdir("trivial_xsec_ana"));
153  }
154  else {
155  TString in_dir = "/nova/ana/users/ddoyle/TrivialCrossSectionAnalysis/";
156  TString in_file = "demo_trivial_xsec_analysis.root";
157 
158  TFile * in = TFile::Open((in_dir + in_file).Data());
159  TrivialCrossSectionAnalysis * xsec_ana
160  = TrivialCrossSectionAnalysis::LoadFrom(in->GetDirectory("trivial_xsec_ana")).release();
161 
163  xsec_ana->SetUnfoldingReg(1);
164  xsec_ana->CalculateCrossSections();
165 
166  std::vector<XSecSysts::Syst_t> systs = XSecSysts::GetAllSystematics();
167  std::vector<TString> syst_labels;
168  for(int i = 0; i < (int) systs.size(); i++) {
169  syst_labels.push_back(XSecSysts::LabelFromSystematic(systs[i]));
170  }
171 
172  std::vector<TString> bkgd_labels = {"NueCC",
173  "AntiNueCC",
174  "AntiNumuCC",
175  "NC",
176  "Nonfiducial NumuCC",
177  "other"};
178  std::vector<Int_t> syst_colors =
179  {kRed,
180  kAzure - 3,
181  kAzure - 3,
182  kMagenta,
183  kGreen -3,
184  kGreen -3,
185  kOrange + 7
186  };
187  std::vector<Int_t> bkgd_colors =
188  {
189  kAzure + 10,
190  kOrange +10,
191  kMagenta,
192  kGreen,
193  kGray + 1,
194  kRed
195  };
196 
197  PlotAllCrossSectionOverENu("trivial_xsec_ana_xsec_over_enu.pdf",
198  xsec_ana,
199  systs,
200  syst_labels,
201  syst_colors,
202  true);
203  PlotAllCrossSectionResult("trivial_xsec_ana_result.pdf",
204  "Neutrino Energy (GeV)",
205  "#sigma (cm^{2})",
206  xsec_ana,
207  systs,
208  syst_labels,
209  syst_colors,
210  -1);
211  PlotAllSignalEstimates("trivial_xsec_ana_signalest.pdf",
212  xsec_ana,
213  systs,
214  syst_labels,
215  syst_colors);
216  PlotAllUnfoldedSignalEstimates("trivial_xsec_ana_unfolded_signalest.pdf",
217  xsec_ana,
218  systs,
219  syst_labels,
220  syst_colors);
221  PlotAllRecoTrue("trivial_xsec_ana_recotrue.pdf",
222  xsec_ana,
223  systs,
224  syst_labels);
225  PlotAllEfficiency("trivial_xsec_ana_eff.pdf",
226  xsec_ana,
227  systs,
228  syst_labels,
229  syst_colors);
230  PlotAllPurity("trivial_xsec_ana_purity.pdf",
231  xsec_ana,
232  systs,
233  syst_labels,
234  syst_colors);
235  PlotAllFlux("trivial_xsec_ana_flux.pdf",
236  xsec_ana,
237  systs,
238  syst_labels,
239  syst_colors);
240  PlotAllSelectionDecomposition("trivial_xsec_ana_sel_decomp.pdf",
241  xsec_ana,
242  systs,
243  syst_labels,
244  bkgd_labels,
245  bkgd_colors);
246 
247  }
248 }
static TString LabelFromSystematic(Syst_t)
Definition: XSecSysts.cxx:25
void PlotAllSignalEstimates(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > labels, std::vector< Int_t > colors)
Plot all of the signal estimates in a ICrossSectionAnalysis object.
Definition: Plotting.h:156
const NuTruthCut kIsTrueSigST
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NuTruthVar kTrueE_NT([](const caf::SRNeutrinoProxy *sr){return float(sr->E);})
True neutrino energy.
Definition: Vars.h:95
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
GenericVar< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:76
const TVector3 vtxmin(-130,-176, 225)
const Cut kAllNumuCCCutsTrueNotFiducial
const Var kNumuE2018
Definition: NumuEFxs.h:171
void PlotAllSelectionDecomposition(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > systs_labels, std::vector< TString > bkgd_labels, std::vector< Int_t > bkgd_colors)
Definition: Plotting.h:61
void PlotAllPurity(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > labels, std::vector< Int_t > colors)
Plot all purity distributions in an ICrossSectionAnalysis object.
Definition: Plotting.h:201
TFile * fmc
Definition: bdt_com.C:14
void SetSpillCut(const SpillCut &cut)
void PlotAllRecoTrue(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > labels)
Plot all reco to true matrices.
Definition: Plotting.h:135
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
void PlotAllCrossSectionOverENu(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > syst_labels, std::vector< Int_t > colors, bool logx=true)
Plot cross section over Enu for each systematic shift.
Definition: Plotting.h:222
void PlotAllEfficiency(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > labels, std::vector< Int_t > colors)
Plot all efficiencies.
Definition: Plotting.h:93
void InitDataSpectrum(SpectrumLoaderBase *)
void PlotAllCrossSectionResult(TString outfile, TString xtitle, TString ytitle, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > syst_labels, std::vector< Int_t > colors, double scale=1e-38)
Plot cross section result of each systematic shift.
Definition: Plotting.h:254
const Cut kIsNCCut
const Cut kIsAntiNueCC
void SaveTo(TDirectory *, const std::string &name) const override
void PlotAllFlux(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > labels, std::vector< Int_t > colors)
Plot all flux spectra.
Definition: Plotting.h:115
static std::vector< Syst_t > GetAllSystematics()
Definition: XSecSysts.cxx:9
void demo_trivial_xsec_analysis(bool fill_spectra=false)
static std::unique_ptr< TrivialCrossSectionAnalysis > LoadFrom(TDirectory *, const std::string &name)
Var weights
virtual void Go() override
Load all the registered spectra.
std::vector< NuTruthVar > GetkPPFXFluxUnivWgtST()
Definition: PPFXWeights.cxx:33
const NuTruthVar kXSecCVWgt2018_smallerDISScale_NT
Definition: XsecTunes.h:48
void InitPPFXUniverses(SpectrumLoaderBase *, std::vector< Var >, std::vector< NuTruthVar >)
const SystShifts kNoShift
Definition: SystShifts.h:115
void SetUnfoldingMethod(UnfoldMethod_t)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ifstream in
Definition: comparison.C:7
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
void PlotAllUnfoldedSignalEstimates(TString outfile, ICrossSectionAnalysis *xsec_ana, std::vector< XSecSysts::Syst_t > systs, std::vector< TString > labels, std::vector< Int_t > colors)
Plot all unfolded signal estimates in an ICrossSectionAnalysis object.
Definition: Plotting.h:180
void InitGenieUniverses(SpectrumLoaderBase *, std::vector< SystShifts >)
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::vector< SystShifts > GetSystShifts()
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
GenericVar< caf::SRNeutrinoProxy > NuTruthVar
Var designed to be used over the nuTree, ie all neutrinos, not just those that got slices...
Definition: Var.h:84
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
const TVector3 vtxmax(160, 160, 1000)
const Cut kIsAntiNumuCC
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Cut kIsNueCC
void InitSystematicSpectra(XSecSysts::Syst_t, SpectrumLoaderBase *)
const NuTruthVar kPPFXFluxCVWgtST([](const caf::SRNeutrinoProxy *nu){ if(nu->rwgt.ppfx.cv!=nu->rwgt.ppfx.cv){return 1.f;}if(nu->rwgt.ppfx.cv >90){return 1.f;}return float(nu->rwgt.ppfx.cv);})
weight events with the flux PPFX Central value correction.
Definition: PPFXWeights.h:12
const Cut kAllNumuCCCuts
Cut CutFromNuTruthCut(const NuTruthCut &stc)
Definition: Cut.cxx:128
GenericHistAxis< NuTruthVar > NuTruthHistAxis
Definition: HistAxis.h:114