uncertainty2png.C
Go to the documentation of this file.
1 
2 
3 
4 #ifdef __CINT__
5 
6 void uncertainty2png()
7 {
8  std::cout << "Sorry, you must run in compiled mode" << std::endl;
9 }
10 #else
11 
12 
16 #include "CAFAna/Core/HistAxis.h"
17 #include "CAFAna/Vars/Vars.h"
18 #include "CAFAna/Core/Binning.h"
19 #include "CAFAna/Core/Cut.h"
20 #include "CAFAna/Core/Ratio.h"
21 #include "CAFAna/Cuts/Cuts.h"
22 #include "CAFAna/Cuts/SpillCuts.h"
24 #include "CAFAna/Cuts/TruthCuts.h"
25 
28 
30 
31 #include "TFile.h"
32 #include "CAFAna/Core/Spectrum.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "/nova/app/users/acedeno/tag_releasesS18-06-14/CAFAna/Vars/NDNCPi0Xsec.h"
40 
41 
42 #include "CAFAna/Analysis/Plots.h"
43 #include "TLegend.h"
44 #include "TCanvas.h"
45 #include "TColor.h"
46 
47 #include <stdio.h>
48 #include <string.h>
49 
50 using namespace ana;
51 
52 
53 // declare vectors of all Spectra that we want to generate
54 
55 std::string fname = "PNG1OPTI_newpngcvn_nominalWFlux.root";
56 
57 TFile fout(fname.c_str(),"RECREATE");
58 
59 
60 std::vector<std::string> parse_as_vector(std::string s)
61 {
62  std::vector<std::string> ret;
63 
64  char * pch;
65  char *cstr = &s[0u];
66  // strtok splits at delimiter, but input has to be char *
67  pch = strtok (cstr,",");
68 
69  while (pch != NULL)
70  {
71  // cout << pch << endl;
72  ret.push_back(pch);
73  pch = strtok (NULL, ",");
74  }
75 
76  return ret;
77 }
78 
79 //////////////////////////////////////////////////////////////
80 
82 {
83  const int dataIdx = 0;
84  const int nominalIdx = 1;
85 
86  // const std::string& datasetString = "prod_caf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1 with limit 10,prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1 with limit 10,prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1 with limit 10,prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1 with limit 10,prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-neg-offset_v1 with limit 10,prod_caf_R17-03-01-prod3reco.h_nd_genie_nonswap_fhc_nova_v08_periods1235_calib-shift-nd-xyview-pos-offset_v1 with limit 10,prod_caf_R17-03-01-prod3reco.j_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1 with limit 10";
87 
88 
89 
90  const std::string& datasetString = "CVN_Prong_Prod4_ND_NuMI_MC_2view_Genielike_FlatFlux_CAFs,dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightdown-calibup_v1,dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_lightmodel-lightup-calibdown_v1,dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.l_nd_genie_nonswap_fhc_nova_v08_full_ckv-proton-shift-down_v1,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,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,dataset_def_name_newest_snapshot prod_caf_R17-03-01-prod3reco.j_nd_genie_nonswap_fhc_nova_v08_full_calib-shift-nd-func_v1";
91 
92  std::vector<std::string> dataset_labels = {"nominal", "lightdown", "lightup", "ckv", "calibneg", "calibpos", "calibshape"};
93 
94  std::vector<ana::SpectrumLoaderBase*> loaders;
95 
96  std::vector<std::string> datasets =
97  parse_as_vector(datasetString);
98 
99 
100 
101 
102  // Print all the inputs for easy checking
103  std::cout<<"\nDataset names:\n";
104  for(auto &id: datasets)
105  std::cout<<id<<"\n";
106 
107  // vector of spec of different vars for different loaders
108  std::vector<std::vector<ana::Spectrum*> > spec_totmc;
109  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only; // true signal (No Cuts)
110  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts; // true signal that passes (Signal Cuts)
111  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts; // true signal that passes all cuts (Total Background)
112 
113  // vector of spec of different vars for nominal loader, but different flux univs
114  std::vector<std::vector<ana::Spectrum*> > spec_totmc_flux;
115  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only_flux; // true signal (No Cuts)
116  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts_flux; // true signal that passes (Signal Cuts)
117  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts_flux; // true signal that passes all cuts (Total Background)
118 
119  // vector of spec of different vars for nominal loader, but different xsec univs
120  std::vector<std::vector<ana::Spectrum*> > spec_totmc_xsec;
121  std::vector<std::vector<ana::Spectrum*> > spec_Nsig_fid_only_xsec; // true signal (No Cuts)
122  std::vector<std::vector<ana::Spectrum*> > spec_NSel_Cuts_xsec; // true signal that passes (Signal Cuts)
123  std::vector<std::vector<ana::Spectrum*> > spec_Nbkg_Cuts_xsec; // true signal that passes all cuts (Total Background)
124 
125 
126  // set spill cuts on all spectra
127  for(unsigned int idataset = 0; idataset < datasets.size(); ++idataset){
128 
130  loader = new SpectrumLoader(datasets[idataset]);
131 
133  loaders.push_back(loader);
134  }// end loop over loaders
135 
136  const int nvar = 1;
137 
138  HistAxis taxes[nvar] = {
139  HistAxis ("BDTG",BDTGBins,kNCPi0BDTID)
140  };
141 
142  std::vector<std::string> taxes_labels = {"BDTG"};
143 
144  std::vector<ana::Var> ppfx_univ = ana::GetkPPFXFluxUnivWgt();
146  //GenieMultiverseParameters verse(100, "/nova/app/users/acedeno/tag_releasesS18-06-14/CAFAna/XSec/Utilities/knob_config.txt");
147  std::vector<std::map<const ISyst*, double> > genie_shifts = verse.ShiftTables();
148 
149  // loop over loaders and vars to declare spectra
150  for(unsigned int iload = 0; iload < loaders.size(); ++iload){
151 
153  // SpillTruthVar weightST = kXSecCVWgt2017ST * kPPFXFluxCVWgtST;
154 
155  // one vector per loader for all the vars
156 
157  std::vector<ana::Spectrum*> vspec_totmc; // Total PRESel
158  std::vector<ana::Spectrum*> vspec_Nsig_fid_only; // denominator: (true signal fiducial)
159  std::vector<ana::Spectrum*> vspec_NSel_Cuts; //signal(NSel) numerator: (Signal Cuts)
160  std::vector<ana::Spectrum*> vspec_Nbkg_Cuts; //kIsPreSel + !kPi0 + !kNCurrent (Nbkg)
161 
162  for(int ivar = 0; ivar < nvar; ivar++){
163 
164  vspec_totmc.push_back(new Spectrum(*loaders[iload], ////Total PreSel : kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNpCVN1gammaID;kNpCVN2gammaID for optimizing prong 2;
165  taxes[ivar],
166  kIsPreSel ,
167  kNoShift,
168  weight) );
169 
170  vspec_NSel_Cuts.push_back(new Spectrum(*loaders[iload], ///Signal Presel : kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && kNCurrent && kPi0 (0.3geV) ; (Signal PreSel Cuts)
171 
172  taxes[ivar],
174  kNoShift,
175  weight) );
176 
177  vspec_Nsig_fid_only.push_back(new Spectrum(*loaders[iload], ///true signal fiducial (Nsig,fid only): kVtxInFiducial && kNCurrent && kPi0 (0.3geV);no 2 prong...delete, delete.
178  taxes[ivar],
180  kNoShift,
181  weight) );
182 
183 
184 
185  vspec_Nbkg_Cuts.push_back(new Spectrum(*loaders[iload], ////Selected Bakcground (Nbkg):kVtxInFiducial && kContainment && k2Prongs && kPi0RemIDND && !kPi0 && !kNCurrent; (Not Signal Background)
186  taxes[ivar],
188  kNoShift,
189  weight));
190 
191 
192  // let's loop over all the flux and genie universes
193  // there are 100 of them so we can do it together
194  if(iload == 0){
195  std::vector<ana::Spectrum*> vspec_totmc_flux;
196  std::vector<ana::Spectrum*> vspec_Nsig_fid_only_flux; // true signal (No Cuts)
197  std::vector<ana::Spectrum*> vspec_NSel_Cuts_flux; // true signal that passes (Signal Cuts)
198  std::vector<ana::Spectrum*> vspec_Nbkg_Cuts_flux; // true signal that passes all cuts (Total Background)
199 
200  std::vector<ana::Spectrum*> vspec_totmc_xsec;
201  std::vector<ana::Spectrum*> vspec_Nsig_fid_only_xsec; // true signal (No Cuts)
202  std::vector<ana::Spectrum*> vspec_NSel_Cuts_xsec; // true signal that passes (Signal Cuts)
203  std::vector<ana::Spectrum*> vspec_Nbkg_Cuts_xsec; // true signal that passes all cuts (Total Background)
204 
205  for(int iuniv = 0; iuniv < 100; iuniv++){
206 
207  // ppfx universes
208 
209 
210  vspec_totmc_flux.push_back(new Spectrum(*loaders[iload],
211  taxes[ivar],
212  kIsPreSel,
213  kNoShift,
214  ppfx_univ[iuniv]) );
215 
216  vspec_NSel_Cuts_flux.push_back(new Spectrum(*loaders[iload],
217  taxes[ivar],
219  kNoShift,
220  ppfx_univ[iuniv]) );
221 
222  vspec_Nsig_fid_only_flux.push_back(new Spectrum(*loaders[iload],
223  taxes[ivar],
225  kNoShift,
226  ppfx_univ[iuniv]) );
227 
228 
229  vspec_Nbkg_Cuts_flux.push_back(new Spectrum(*loaders[iload],
230  taxes[ivar],
231  kIsPreSel && !kIsSignal01L ,
232  kNoShift,
233  ppfx_univ[iuniv]) );
234 
235  // genie universes
236 
237  vspec_totmc_xsec.push_back(new Spectrum(*loaders[iload],
238  taxes[ivar],
239  kIsPreSel,
240  genie_shifts[iuniv],
241  weight) );
242  vspec_NSel_Cuts_xsec.push_back(new Spectrum(*loaders[iload],
243  taxes[ivar],
245  genie_shifts[iuniv],
246  weight) );
247  vspec_Nsig_fid_only_xsec.push_back(new Spectrum(*loaders[iload],
248  taxes[ivar],
250  genie_shifts[iuniv],
251  weight) );
252 
253  vspec_Nbkg_Cuts_xsec.push_back(new Spectrum(*loaders[iload],
254  taxes[ivar],
256  genie_shifts[iuniv],
257  weight) );
258 
259  }// end loop over univeres
260 
261  spec_totmc_xsec.push_back(vspec_totmc_xsec);
262  spec_Nsig_fid_only_xsec.push_back(vspec_Nsig_fid_only_xsec);
263  spec_NSel_Cuts_xsec.push_back(vspec_NSel_Cuts_xsec);
264  spec_Nbkg_Cuts_xsec.push_back(vspec_Nbkg_Cuts_xsec);
265 
266  spec_totmc_flux.push_back(vspec_totmc_flux);
267  spec_Nsig_fid_only_flux.push_back(vspec_Nsig_fid_only_flux);
268  spec_NSel_Cuts_flux.push_back(vspec_NSel_Cuts_flux);
269  spec_Nbkg_Cuts_flux.push_back(vspec_Nbkg_Cuts_flux);
270 
271  }// only build universes for nominal loader
272 
273  }// end loop over axes
274 
275  spec_totmc.push_back(vspec_totmc);
276  spec_Nsig_fid_only.push_back(vspec_Nsig_fid_only);
277  spec_NSel_Cuts.push_back(vspec_NSel_Cuts);
278  spec_Nbkg_Cuts.push_back(vspec_Nbkg_Cuts);
279 
280  }// end loop over loaders
281 
282  //----------------------------------------------
283 
284 
285  for(unsigned int iload = 0; iload < loaders.size(); ++iload){
286  loaders[iload]->Go();
287  }
288 
289 
290  for(unsigned int iload = 0; iload < loaders.size(); ++iload){
291 
292  TDirectory* d = fout.mkdir(dataset_labels[iload].c_str());
293 
294  for(int ivar = 0; ivar < nvar; ++ivar){
295  std::cout<<"loader "<<dataset_labels[iload]<<" var "<<taxes_labels[ivar]<<"\n";
296  spec_totmc[iload][ivar]->SaveTo(d->mkdir( ("totmc_"+taxes_labels[ivar]).c_str()));
297  spec_Nsig_fid_only[iload][ivar]->SaveTo(d->mkdir( ("Nsig_fid_only_"+taxes_labels[ivar]).c_str()));
298  spec_NSel_Cuts[iload][ivar]->SaveTo(d->mkdir( ("NSel_Cuts_"+taxes_labels[ivar]).c_str()));
299  spec_Nbkg_Cuts[iload][ivar]->SaveTo(d->mkdir( ("Nbkg_Cuts_"+taxes_labels[ivar]).c_str()));
300 
301  if(iload == 0){
302 
303  TDirectory* dftmc = fout.mkdir(("flux_totmc_"+taxes_labels[ivar]).c_str());
304  TDirectory* dft = fout.mkdir(("flux_Nsig_fid_only_"+taxes_labels[ivar]).c_str());
305  TDirectory* dfp = fout.mkdir(("flux_NSel_Cuts_"+taxes_labels[ivar]).c_str());
306  TDirectory* dfb = fout.mkdir(("flux_Nbkg_Cuts_"+taxes_labels[ivar]).c_str());
307 
308  TDirectory* dxtmc = fout.mkdir(("xsec_totmc_"+taxes_labels[ivar]).c_str());
309  TDirectory* dxt = fout.mkdir(("xsec_Nsig_fid_only_"+taxes_labels[ivar]).c_str());
310  TDirectory* dxp = fout.mkdir(("xsec_NSel_Cuts_"+taxes_labels[ivar]).c_str());
311  TDirectory* dxb = fout.mkdir(("xsec_Nbkg_Cuts_"+taxes_labels[ivar]).c_str());
312 
313  for(int iuniv = 0; iuniv < 100; iuniv++){
314  // std::cout<<"- universe "<<iuniv<<"\n";
315  //std::cout<<"---------------------------------------\n";
316  spec_totmc_flux[ivar][iuniv]->SaveTo(dftmc->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
317  //std::cout<<"---------------------------------------\n";
318  spec_Nsig_fid_only_flux[ivar][iuniv]->SaveTo(dft->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
319  //std::cout<<"---------------------------------------\n";
320 
321  spec_NSel_Cuts_flux[ivar][iuniv]->SaveTo(dfp->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
322  //std::cout<<"---------------------------------------\n";
323 
324  spec_Nbkg_Cuts_flux[ivar][iuniv]->SaveTo(dfb->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
325  //std::cout<<"---------------------------------------\n";
326 
327  spec_totmc_xsec[ivar][iuniv]->SaveTo(dxtmc->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
328  //std::cout<<"---------------------------------------\n";
329 
330  spec_Nsig_fid_only_xsec[ivar][iuniv]->SaveTo(dxt->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
331  //std::cout<<"---------------------------------------\n";
332 
333  spec_NSel_Cuts_xsec[ivar][iuniv]->SaveTo(dxp->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
334  //std::cout<<"---------------------------------------\n";
335 
336  spec_Nbkg_Cuts_xsec[ivar][iuniv]->SaveTo(dxb->mkdir( ("spec"+std::to_string(iuniv)).c_str()));
337  //std::cout<<"---------------------------------------\n";
338 
339  }
340  }
341 
342  }//end loop over vars
343  }// end loop over loaders
344 std::cout << "Succesfull****** " <<std::endl;
345 }//end of void
346 
347 #endif
const Var kNCPi0BDTID
Definition: NDNCPi0Xsec.h:23
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
const Cut TrueSignalVtxFid
FIDUCIAL OPTIMIZATION ANALYSYS 2 PRONG NO CUTS 0.3 GeV THRESHOLD///.
Definition: ncpi0Cuts.h:1127
void SetSpillCut(const SpillCut &cut)
const Cut kIsSignal01L
FIDUCIAL OPTIMIZATION ANALYSIS 2 PRONG TOTAL PRESELECTION CUTS///.
Definition: ncpi0Cuts.h:1092
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.
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
TFile fout(fname.c_str(),"RECREATE")
Float_t d
Definition: plot.C:236
std::vector< std::string > parse_as_vector(std::string s)
loader
Definition: demo0.py:10
std::vector< float > Spectrum
Definition: Constants.h:570
const Cut kIsPreSel
Definition: ncpi0Cuts.h:1099
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
Base class for the various types of spectrum loader.
::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
void uncertainty2png()
std::string fname
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const Binning BDTGBins
Definition: ncpi0Bins.h:55
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
std::vector< Dataset > datasets
Definition: Periods.h:3
std::vector< const ISyst * > getAllXsecSysts_2017()
Get master XSec syst list for 2017 analyses.