uncertainty1png.C
Go to the documentation of this file.
1 
2 
5 #include "CAFAna/Core/HistAxis.h"
6 #include "CAFAna/Vars/Vars.h"
7 #include "CAFAna/Core/Binning.h"
8 #include "CAFAna/Core/Cut.h"
9 #include "CAFAna/Core/Ratio.h"
10 #include "CAFAna/Cuts/Cuts.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
13 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "TFile.h"
18 #include "CAFAna/Core/Spectrum.h"
19 #include "TH1.h"
20 #include "TH2.h"
22 
25 
26 #include "CAFAna/Analysis/Plots.h"
27 #include "TLegend.h"
28 #include "TCanvas.h"
29 #include "TColor.h"
30 #include <stdio.h>
31 #include <string.h>
32 
33 using namespace ana;
34 
35 std::string fname = "BDTGOutput_optimizationTest.root";
36 
37 TFile fout(fname.c_str(),"RECREATE");
38 
39 /******************************************************************/
40 
41 //Declare vectors of all Spectra that we want to generate
42 std::vector<std::string> parse_as_vector(std::string s)
43 {
44  std::vector<std::string> ret;
45  char *pch;
46  char *cstr = &s[0u];
47  // strtok splits at delimiter, but input has to be char *
48  pch = strtok (cstr,",");
49 
50  while (pch != NULL)
51  {
52  // cout << pch << endl;
53  ret.push_back(pch);
54  pch = strtok (NULL, ",");
55  }
56  return ret;
57 }
58 
59 /****************************************************************/
60 
62 {
63  //Define the datasets and indices to be used.
64  const std::string& datasetString =
65  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1,"
66  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1,"
67  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1,"
68  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1,"
69  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-neg-offset_v1,"
70  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-pos-offset_v1,"
71  "dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.j_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1";
72 
73  const int dataIdx = 0;
74  const int nominalIdx = 1;
75 
76  /*****************************************************************/
77 
78  //Label and load in the spectrum objects for each dataset.
79  std::vector<std::string> dataset_labels = {"nominal", "lightdown", "lightup", "ckv", "calibneg", "calibpos", "calibshape"};
80  std::vector<ana::SpectrumLoaderBase*> loaders;
81  std::vector<std::string> datasets = parse_as_vector(datasetString);
82 
83  // Print all the inputs for easy checking.
84  std::cout<<"\nDataset names:\n";
85  for(auto &id: datasets)
86  std::cout<<id<<"\n";
87 
88  /*****************************************************************/
89 
90  /*Define vectors of spectra of different variables for different loaders.*/
91 
92  //Total truth information.
93  std::vector<std::vector<ana::Spectrum*> > spec_totmc;
94 
95  //True signal with fiducial cuts only.
96  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only;
97 
98  //True signal that passes preselection cuts.
99  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts;
100 
101  //Total background.
102  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts;
103 
104  /******************************************************************/
105 
106  /*Define vectors of spectra of different variables for
107  nominal loader but different flux universes.*/
108 
109  //Total flux truth information.
110  std::vector<std::vector<ana::Spectrum*> > spec_totmc_flux;
111 
112  //True signal with fiducial cuts only.
113  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only_flux;
114 
115  //True signal that passes preselection cuts.
116  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts_flux;
117 
118  //True signal that passes all cuts (Total Background)
119  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts_flux;
120 
121  /*******************************************************************/
122 
123  /*Define vectors of spectra of different variables
124  for nominal loader, but different xsec universes.*/
125 
126  //Total MC Preselection.
127  std::vector<std::vector<ana::Spectrum*> > spec_totmc_xsec;
128 
129  //True signal with fiducial cuts only.
130  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only_xsec;
131 
132  //Signal that passes preselection cuts.
133  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts_xsec;
134 
135  //True signal that passes all cuts (Total Background)
136  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts_xsec;
137 
138  /********************************************************************/
139 
140  //Set spill cuts on all spectra by looping over all loaders.
141  for(unsigned int idataset = 0; idataset < datasets.size(); ++idataset)
142  {
144  loader = new SpectrumLoader(datasets[idataset]);
146  loaders.push_back(loader);
147  }
148 
149  const int nvar = 1;
150 
151  const Binning BDTBins = Binning::Simple(100,-1,1);
152 
153  HistAxis taxes[nvar] =
154  {
155  HistAxis ("Single Prong BDTG Output", BDTBins, kNCPi0BDTID)
156  };
157 
158  /************************************************************************/
159 
160  //Begin filling spectra with multiverse parameters.
161  std::vector<std::string> taxes_labels = {"BDTG Output"};
162  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
163  GenieMultiverseParameters verse(100, "/nova/app/users/rbowles/tag_releaseS18-04-23/CAFAna/XSec/Utilities/knob_config.txt");
164  std::vector<std::map<const ISyst*, double> > genie_shifts = verse.ShiftTables();
165 
166  //Loop over loaders and variables to fill spectra.
167  for(unsigned int iload = 0; iload < loaders.size(); ++iload)
168  {
170  // SpillTruthVar weightST = kXSecCVWgt2017ST * kPPFXFluxCVWgtST;
171  // one vector per loader for all the vars
172 
173  //Total Preselection.
174  std::vector<ana::Spectrum*> vspec_totmc;
175 
176  //Denominator: True signal fiducial.
177  std::vector<ana::Spectrum*> vspec_Nsig_fid_only;
178 
179  //Signal numerator: Signal cuts.
180  std::vector<ana::Spectrum*> vspec_NSel_Cuts;
181 
182  //Preselection + !Signal
183  std::vector<ana::Spectrum*> vspec_Nbkg_Cuts;
184 
185  //Loop over axes.
186  for(int ivar = 0; ivar < nvar; ivar++)
187  {
188  //Total Preselection including RemID.
189  vspec_totmc.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection,kNoShift,weight));
190 
191  //Total signal.
192  vspec_NSel_Cuts.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Signal,kNoShift,weight));
193 
194  //True Signal with only true fiducial (Fiducial, pi0, NCurrent)
195  vspec_Nsig_fid_only.push_back(new Spectrum(*loaders[iload],taxes[ivar], kTrueFiducial && Signal,kNoShift,weight));
196 
197  //Total Background
198  vspec_Nbkg_Cuts.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Bkg,kNoShift,weight));
199 
200 
201  //Only build universes for the nominal loader.
202  if(iload == 0)
203  {
204  /*Define flux spectra.*/
205 
206  //Total flux.
207  std::vector<ana::Spectrum*> vspec_totmc_flux;
208 
209  //Total signal flux.
210  std::vector<ana::Spectrum*> vspec_Nsig_fid_only_flux;
211 
212  //Signal that passes preselection.
213  std::vector<ana::Spectrum*> vspec_NSel_Cuts_flux;
214 
215  //Total background.
216  std::vector<ana::Spectrum*> vspec_Nbkg_Cuts_flux;
217 
218  /*****************************************************************/
219  /*Define cross section spectra.*/
220 
221  //Total cross section.
222  std::vector<ana::Spectrum*> vspec_totmc_xsec;
223 
224  //Total signal cross section.
225  std::vector<ana::Spectrum*> vspec_Nsig_fid_only_xsec;
226 
227  //Signal that passes preselection.
228  std::vector<ana::Spectrum*> vspec_NSel_Cuts_xsec;
229 
230  //Total background.
231  std::vector<ana::Spectrum*> vspec_Nbkg_Cuts_xsec;
232 
233  //Loop over all 100 universes.
234  for(int iuniv = 0; iuniv < 100; iuniv++)
235  {
236  /**********Fill spectra from ppfx universes.***********/
237  vspec_totmc_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection,kNoShift,ppfx_univ[iuniv]));
238 
239  vspec_NSel_Cuts_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Signal,kNoShift,ppfx_univ[iuniv]));
240 
241  vspec_Nsig_fid_only_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],kTrueFiducial && Signal,kNoShift,ppfx_univ[iuniv]));
242 
243  vspec_Nbkg_Cuts_flux.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Bkg,kNoShift,ppfx_univ[iuniv]));
244 
245  /**********Fill spectra from genie universes.**********/
246  vspec_totmc_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection,genie_shifts[iuniv],weight));
247 
248  vspec_NSel_Cuts_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Signal,genie_shifts[iuniv],weight) );
249 
250  vspec_Nsig_fid_only_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],kTrueFiducial && Signal,genie_shifts[iuniv],weight));
251 
252  vspec_Nbkg_Cuts_xsec.push_back(new Spectrum(*loaders[iload],taxes[ivar],Preselection && Bkg,genie_shifts[iuniv],weight));
253  }
254  spec_totmc_xsec.push_back(vspec_totmc_xsec);
255  spec_Nsig_fid_only_xsec.push_back(vspec_Nsig_fid_only_xsec);
256  spec_NSel_Cuts_xsec.push_back(vspec_NSel_Cuts_xsec);
257  spec_Nbkg_Cuts_xsec.push_back(vspec_Nbkg_Cuts_xsec);
258 
259  spec_totmc_flux.push_back(vspec_totmc_flux);
260  spec_Nsig_fid_only_flux.push_back(vspec_Nsig_fid_only_flux);
261  spec_NSel_Cuts_flux.push_back(vspec_NSel_Cuts_flux);
262  spec_Nbkg_Cuts_flux.push_back(vspec_Nbkg_Cuts_flux);
263  }
264  }
265  spec_totmc.push_back(vspec_totmc);
266  spec_Nsig_fid_only.push_back(vspec_Nsig_fid_only);
267  spec_NSel_Cuts.push_back(vspec_NSel_Cuts);
268  spec_Nbkg_Cuts.push_back(vspec_Nbkg_Cuts);
269  }
270  /************************************************************************/
271 
272  for(unsigned int iload = 0; iload < loaders.size(); ++iload)
273  {
274  loaders[iload]->Go();
275  }
276 
277  //Loop over loaders
278  for(unsigned int iload = 0; iload < loaders.size(); ++iload)
279  {
280  TDirectory* d = fout.mkdir(dataset_labels[iload].c_str());
281 
282  //Loop over variables.
283  for(int ivar = 0; ivar < nvar; ++ivar)
284  {
285  std::cout<<"loader "<<dataset_labels[iload]<<" var "<<taxes_labels[ivar]<<"\n";
286  spec_totmc[iload][ivar]->SaveTo(d->mkdir( ("totmc_"+taxes_labels[ivar]).c_str()));
287  spec_Nsig_fid_only[iload][ivar]->SaveTo(d->mkdir( ("Nsig_fid_only_"+taxes_labels[ivar]).c_str()));
288  spec_NSel_Cuts[iload][ivar]->SaveTo(d->mkdir( ("NSel_Cuts_"+taxes_labels[ivar]).c_str()));
289  spec_Nbkg_Cuts[iload][ivar]->SaveTo(d->mkdir( ("Nbkg_Cuts_"+taxes_labels[ivar]).c_str()));
290 
291  if(iload == 0)
292  {
293  TDirectory* dftmc = fout.mkdir(("flux_totmc_"+taxes_labels[ivar]).c_str());
294  TDirectory* dft = fout.mkdir(("flux_Nsig_fid_only_"+taxes_labels[ivar]).c_str());
295  TDirectory* dfp = fout.mkdir(("flux_NSel_Cuts_"+taxes_labels[ivar]).c_str());
296  TDirectory* dfb = fout.mkdir(("flux_Nbkg_Cuts_"+taxes_labels[ivar]).c_str());
297 
298  TDirectory* dxtmc = fout.mkdir(("xsec_totmc_"+taxes_labels[ivar]).c_str());
299  TDirectory* dxt = fout.mkdir(("xsec_Nsig_fid_only_"+taxes_labels[ivar]).c_str());
300  TDirectory* dxp = fout.mkdir(("xsec_NSel_Cuts_"+taxes_labels[ivar]).c_str());
301  TDirectory* dxb = fout.mkdir(("xsec_Nbkg_Cuts_"+taxes_labels[ivar]).c_str());
302 
303  for(int iuniv = 0; iuniv < 100; iuniv++)
304  {
305  spec_totmc_flux[ivar][iuniv]->SaveTo(dftmc->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
306  spec_Nsig_fid_only_flux[ivar][iuniv]->SaveTo(dft->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
307  spec_NSel_Cuts_flux[ivar][iuniv]->SaveTo(dfp->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
308  spec_Nbkg_Cuts_flux[ivar][iuniv]->SaveTo(dfb->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
309  spec_totmc_xsec[ivar][iuniv]->SaveTo(dxtmc->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
310  spec_Nsig_fid_only_xsec[ivar][iuniv]->SaveTo(dxt->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
311  spec_NSel_Cuts_xsec[ivar][iuniv]->SaveTo(dxp->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
312  spec_Nbkg_Cuts_xsec[ivar][iuniv]->SaveTo(dxb->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
313  }
314  }
315  }
316  }
317  std::cout << "\nOutput saved to BDTGOutput_optimization.root" <<std::endl;
318 }
319 
const Var kNCPi0BDTID
Definition: NDNCPi0Xsec.h:23
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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 Cut Signal
Definition: ncpi0Cuts.h:811
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
TFile fout(fname.c_str(),"RECREATE")
void SetSpillCut(const SpillCut &cut)
const Cut Preselection
std::vector< Var > GetkPPFXFluxUnivWgt()
Definition: PPFXWeights.cxx:53
const SpillCut kTightBeamQualityCuts([](const caf::SRSpillProxy *s){if(s->ismc) return true; if(s->trigger==2) return true;if(s->spilltimesec==0 &&s->deltaspilltimensec==0 &&s->widthx==0) return bool(s->isgoodspill);if(std::abs(s->deltaspilltimensec) > 0.5e9) return false;if(s->spillpot< 2e12) return false;if(s->hornI< -202|| s->hornI >-198) return false;if(s->posx< -2.00|| s->posx >+2.00) return false;if(s->posy< -2.00|| s->posy >+2.00) return false;return kBeamWidthCut(s);})
Definition: SpillCuts.h:10
const XML_Char * s
Definition: expat.h:262
std::vector< std::map< const ISyst *, double > > ShiftTables()
Return the tables of knob shift values.
void uncertainty1png()
const Cut kTrueFiducial
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
const Cut Bkg
Float_t d
Definition: plot.C:236
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:570
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
Base class for the various types of spectrum loader.
std::string fname
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
std::vector< Loaders * > loaders
Definition: syst_header.h:386
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
std::vector< Dataset > datasets
Definition: Periods.h:3
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
std::vector< std::string > parse_as_vector(std::string s)