multiverse_efficiency_plot.C
Go to the documentation of this file.
1 /*
2 #include "CAFAna/Core/Binning.h"
3 #include "CAFAna/Cuts/Cuts.h"
4 
5 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/SpectrumLoader.h"
7 
8 #include "CAFAna/Core/Var.h"
9 #include "CAFAna/Core/MultiVar.h"
10 #include "CAFAna/Core/HistAxis.h"
11 #include "CAFAna/Core/SystShifts.h"
12 
13 #include "CAFAna/Vars/Vars.h"
14 #include "CAFAna/Vars/NueVarsExtra.h"
15 
16 #include "NDAna/numucc_1Pi/NumuCC1PiVars.h"
17 #include "NDAna/numucc_1Pi/NumuCC1PiCuts.h"
18 
19 #include "G4RwgtPlots/Geant4WeightVars.h"
20 
21 #include "CAFAna/XSec/Multiverse.h"
22 #include "CAFAna/XSec/SystematicDef.h"
23 //#include "CAFAna/XSec/CutOptimization.h"
24 
25 // for POT
26 #include "CAFAna/Analysis/Exposures.h"
27 
28 // for plotting
29 #include "CAFAna/Analysis/Plots.h"
30 
31 // to get datetime in plot names
32 #include "/nova/app/users/csweeney/xsec/code/Utilities/misc.h"
33 
34 #include "StandardRecord/Proxy/SRProxy.h"
35 
36 #include <vector>
37 #include <string>
38 #include <math.h>
39 
40 #include "TFile.h"
41 #include "TH1.h"
42 #include "TCanvas.h"
43 #include "TGraphAsymmErrors.h"
44 #include "TString.h"
45 #include "Math/ProbFunc.h"
46 */
47 
48 
49 #include "CAFAna/Core/Binning.h"
50 #include "CAFAna/Cuts/Cuts.h"
51 
52 #include "CAFAna/Core/Spectrum.h"
54 
55 #include "CAFAna/Core/Var.h"
56 
57 #include "CAFAna/Vars/Vars.h"
58 
59 #include "CAFAna/XSec/Multiverse.h"
61 
63 
64 #include <vector>
65 #include <string>
66 
67 #include "TFile.h"
68 #include "TH1.h"
69 #include "TCanvas.h"
70 #include "TGraphAsymmErrors.h"
71 #include "TString.h"
72 #include "Math/ProbFunc.h"
73 
76 
77 #include "G4RwgtPlots/Geant4WeightVars.h"
78 
79 // for POT
81 
82 // for plotting
83 #include "CAFAna/Analysis/Plots.h"
84 
85 // to get datetime in plot names
86 #include "/nova/app/users/csweeney/xsec-TR/Utilities/misc.h"
87 
88 // for kVtxX
90 
91 
92 
93 namespace ana
94 {
95 
96  void PlotDebug(std::vector<TH1 *> universes,
97  TH1 * h_nominal,
100  //Scavenged from CAFAna/XSec/CutOptimization.h
101  {
102  TDirectory * tmp = gDirectory;
103  TCanvas * c = new TCanvas();
104  c->SetLeftMargin(0.15);
105  c->SetRightMargin(0.01);
106  c->cd();
107  double max_val = -1e9;
108  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
109  max_val = std::max(max_val, universes[ihist]->GetMaximum());
110  }
111  h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
112  h_nominal->SetTitle(title.c_str());
113  h_nominal->Draw("hist");
114  for(auto ihist = 0u; ihist < universes.size(); ihist++) {
115  universes[ihist]->Draw("hist same");
116  }
117  h_nominal->Draw("hist same");
118 
119  c->Print(name.c_str());
120  tmp->cd();
121  }
122 
123 
124  Var kUnity([](const caf::SRProxy *sr){
125  // To be used as "common_weight" for Multiverse constructor
126  return 1.f;
127  });
128 
129 
130  //Scavenged fromn CAFAna/XSec/Multiverse.cxx
131  double BinSigma(std::vector<double> events,
132  double nsigma,
133  double pivot)
134  {
135  int pivotbin = 0;
136  double pivotbincenter = 0;
137  std::sort(events.begin(), events.end());
138  for(int i = 0; i < (int) events.size() - 1; i++) {
139  if(pivot >= events[i] && pivot < events[i+1]) {
140  pivotbin = i;
141  break;
142  }
143  }
144  pivotbincenter = pivotbin+0.5;
145  double count_fraction = 2.0 * (ROOT::Math::normal_cdf(nsigma) - ROOT::Math::normal_cdf(0));
146  int nsideevents = 0;
147  int lastbinindex = (int) events.size() - 1;
148  if(nsigma >= 0) nsideevents = lastbinindex - pivotbin;
149  else nsideevents = pivotbin;
150  int boundIdx = pivotbincenter + count_fraction*(double)nsideevents;
151  int index = 0;
152  if(nsigma >= 0) index = std::min(boundIdx, (int)events.size() - 1);
153  else index = std::max(boundIdx, 0);
154  return events.at(index);
155  }
156 
157 
158  //Scavenged fromn CAFAna/XSec/Multiverse.cxx
159  TH1 * GetNSigmaShift(TH1 * nom_hist, std::vector<TH1*> mv_hists, double nsigma)
160  {
161 
162  // transform to 1D hists
163  TH1 * hret = (TH1*) nom_hist->Clone();
164 
165  /*
166  for(int i = 0; i < (int) fSpectra.size(); i++) {
167  // verify histograms are compatable
168  assert(nom_hist->GetNbinsX() == mv_hists[i]->GetNbinsX());
169  assert(nom_hist->GetNbinsY() == mv_hists[i]->GetNbinsY());
170  assert(nom_hist->GetNbinsZ() == mv_hists[i]->GetNbinsZ());
171  }
172  */
173 
174  // for each bin, calculate the shift
175  for(int i = 1; i <= nom_hist->GetNbinsX(); i++) {
176  for(int j = 1; j <= nom_hist->GetNbinsY(); j++) {
177  for(int k = 1; k <= nom_hist->GetNbinsZ(); k++) {
178  std::vector<double> vals;
179  for(int ihist = 0; ihist < (int) mv_hists.size(); ihist++) {
180  vals.push_back(mv_hists[ihist]->GetBinContent(i, j, k));
181  }
182  double cv = hret->GetBinContent(i, j, k);
183  hret->SetBinContent(i, j, k, BinSigma(vals, nsigma, cv));
184  }
185  }
186  }
187  return hret;
188  }
189 
190 
191 
192 }
193 
194 
195 
196 using namespace ana;
197 
199 {
200 
201  std::string inName = "./10_noMuonID_outFile_vtxX.root";
202  std::string outName = "10_noMuonID_vtxX_";
203  bool hack = true;
204 
205 
206  TFile* inFile = TFile::Open(inName.c_str(), "read");
207 
208 
209 
210 
211  std::vector<std::string> myLabel_vec{
212  "Vertex in X Direction (cm)"//, "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
213  //"Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)",
214  //"Vertex in X Direction (cm)", "Vertex in Y Direction (cm)", "Vertex in Z Direction (cm)"
215  };
216 
217  std::vector<std::string> myName_vec = {
218  "piplus_vtxX"//, "piplus_vtxY", "piplus_vtxZ",
219  // "piminus_vtxX", "piminus_vtxY", "piminus_vtxZ",
220  //"total_vtxX", "total_vtxY", "total_vtxZ"
221  };
222 
223 
224 /*
225  std::vector<std::string> myLabel_vec{ "True neutrino energy (GeV)" };
226  std::vector<std::string> myName_vec = { "piplus_nuE" };
227 */
228 
229 
230  for(uint i=0; i< myName_vec.size(); ++i){
231 
232  std::string base_dir = myName_vec[i];
233 
234  Multiverse * mySig_MV = Multiverse::LoadFrom( inFile->GetDirectory(base_dir.c_str()) , "mv_Sig" ).release();
235  // Spectrum * sSig = Spectrum::LoadFrom( , base_dir + "/spec_Sig").c_str() ) ).release();
236  Spectrum * sSig = Spectrum::LoadFrom( inFile->GetDirectory(base_dir.c_str()) , "spec_Sig" ).release();
237 
238 
239  Multiverse * mySelsig_MV = Multiverse::LoadFrom( inFile->GetDirectory(base_dir.c_str()), "mv_Selsig" ).release();
240  Spectrum * sSelsig = Spectrum::LoadFrom( inFile->GetDirectory(base_dir.c_str()), "spec_Selsig" ).release();
241 
242 
243  std::string myName = myName_vec[i];
244  std::string myLabel = myLabel_vec[i];
245 
246 
247 
248  //double pot = kAna2020FHCPOT;
249  double pot = sSig->POT();
250 
251 
252  TH1 * hSig = sSig->ToTH1(pot);
253  const std::vector<TH1*> mv_Sig_hists = mySig_MV->GetUniversesHist(kCyan);
254 
255  //PlotDebug(mv_signal_hists, hSig, "this is a plot", "plot.png");
256 
257  TH1 * hSelsig = sSelsig->ToTH1(pot);
258  const std::vector<TH1*> mv_Selsig_hists = mySelsig_MV->GetUniversesHist(kCyan);
259 
260  TH1* hNom_eff = (TH1*) hSelsig->Clone();
261  hNom_eff->Divide(hSig);
262  hNom_eff->GetYaxis()->SetTitle("Efficiency");
263 
264 
265 
266  // EXTREMELY hacky way of getting plots to look right. Values are nan at edges
267  if(hack == true){
268  hNom_eff->SetBinContent(1, 0.001);
269  hNom_eff->SetBinContent(hNom_eff->GetNbinsX(), 0.001);
270  }
271 
272  std::vector<TH1*> mv_eff_hists;
273  for(uint i=0; i< mv_Sig_hists.size(); ++i){
274 
275  TH1* eff = (TH1*) mv_Selsig_hists[i]->Clone();
276  eff->Divide(mv_Sig_hists[i]);
277 
278  mv_eff_hists.push_back(eff);
279  }
280 
281 
282  PlotDebug(mv_eff_hists, hNom_eff, "this is a plot", (outName + "big_eff_plot.png").c_str());
283 
284 
285 
286 
287  //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
288  TH1* hUp = GetNSigmaShift(hNom_eff, mv_eff_hists, 1);
289  std::vector<TH1*> hUp_vec{hUp};
290 
291  TH1* hDown = GetNSigmaShift(hNom_eff, mv_eff_hists, -1);
292  std::vector<TH1*> hDown_vec{hDown};
293 
294  // Get fractional uncertainty (take larger of +/-)
295  int nBins = hNom_eff->GetNbinsX();
296  double xLow = hNom_eff->GetXaxis()->GetXmin();
297  double xHigh = hNom_eff->GetXaxis()->GetXmax();
298 
299  TH1D * fracUncert = new TH1D("frac", "", nBins, xLow, xHigh);
300  for(int i=1; i <= nBins; ++i){
301  double up = hUp->GetBinContent(i);
302  double down = hDown->GetBinContent(i);
303  double nom = hNom_eff->GetBinContent(i);
304 
305  double abs_shift = std::max( std::abs(up - nom), std::abs(down - nom) );
306  double frac = abs_shift / nom;
307  //std::cout << "{ " << up << ", " << down << " }" << std::endl;
308  //std::cout << nom << " : " << abs_shift << " : " << frac << std::endl;
309 
310  if( isnan(frac) ) fracUncert->SetBinContent(i, 0.f);
311  else fracUncert->SetBinContent(i, frac);
312 
313  if(hack == true){
314  fracUncert->SetBinContent(1, 0.001);
315  fracUncert->SetBinContent(fracUncert->GetNbinsX(), 0.001);
316  }
317 
318  }
319 
320  fracUncert->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
321  fracUncert->GetXaxis()->SetTitle(myLabel.c_str());
322 
323 
324 
325  //TCanvas * canv2 = new TCanvas("canv2", "canv2", 800, 600);
326  //fracUncert->Draw();
327 
328  hNom_eff->SetTitle("");
329 
330  hNom_eff->SetLineColor(4);
331  hNom_eff->SetLineWidth(1);
332 
333  TGraphAsymmErrors * graph = PlotWithSystErrorBand(hNom_eff, hUp_vec, hDown_vec, 4, 2);
334 
335  TString rCname = "rC";
336  TCanvas* rC = new TCanvas(rCname, rCname);
337 
338  rC -> SetBottomMargin(0.);
339  double Spl = 0.3;
340  TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
341  TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
342  P2 -> SetRightMargin (.03);
343  P2 -> SetTopMargin (.0);
344  P2 -> SetBottomMargin(.3);
345  P2 -> SetLeftMargin (.13);
346  P2 -> Draw();
347  P1 -> SetRightMargin (.03);
348  P1 -> SetLeftMargin (.13);
349  P1 -> SetTopMargin (.02);
350  P1 -> SetBottomMargin(.00);
351  P1 -> Draw();
352  // Set some label sizes.
353  double Lb1 = 0.07;
354  double Lb2 = 0.15;
355  // --- First, draw the fracUncert so cd onto Pad2
356  P2 -> cd();
357 
358  // Set axis ranges etc.
359  fracUncert->GetYaxis()->SetTitleSize( Lb2 );
360  fracUncert->GetYaxis()->SetTitleOffset(0.4);
361  fracUncert->GetYaxis()->SetLabelSize( Lb2 );
362  fracUncert->GetXaxis()->SetTitleSize( Lb2 );
363  fracUncert->GetXaxis()->SetLabelSize( Lb2 );
364  fracUncert->Draw("");
365 
366  //fracUncert->Draw("axis same");
367  P1->cd();
368  // graph goes on P1
369  hNom_eff->GetYaxis()->SetTitleSize( Lb1 );
370  hNom_eff->GetYaxis()->SetLabelSize( Lb1 );
371  hNom_eff->GetYaxis()->SetTitleOffset( 0.7 );
372  // Remove the x axis labels
373  hNom_eff->GetXaxis()->SetLabelSize (0 );
374  hNom_eff->GetXaxis()->SetLabelOffset(99);
375  hNom_eff->SetLineColor(4);
376  hNom_eff->DrawCopy("hist");
377 
378  hNom_eff->Draw("hist ][");
379  graph->Draw("e2 same");
380  hNom_eff->Draw("hist ][ same");
381 
382  rC->SaveAs( (outName + "big_eff_ratio.pdf").c_str() );
383  rC->SaveAs( (outName + "big_eff_ratio.png").c_str() );
384 
385 
386  inFile->Close();
387  }
388 
389 
390 }
391 
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
int nBins
Definition: plotROC.py:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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:148
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype, double alpha)
Plot prediction with +/-1sigma error band.
Definition: Plots.cxx:304
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
void multiverse_efficiency_plot()
Float_t tmp
Definition: plot.C:36
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
ifstream inFile
Definition: AnaPlotMaker.h:34
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const std::string cv[Ncv]
#define pot
caf::StandardRecord * sr
void PlotDebug(std::vector< TH1 * > universes, TH1 *h_nominal, std::string title, std::string name)
const double j
Definition: BetheBloch.cxx:29
double frac(double x)
Fractional part.
double POT() const
Definition: Spectrum.h:227
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Multiverse.cxx:573
double BinSigma(std::vector< double > events, double nsigma, double pivot)
std::vector< TH1D * > nom_hist
const std::vector< TH1 * > GetUniversesHist()
Definition: Multiverse.cxx:120
TH1 * GetNSigmaShift(TH1 *nom_hist, std::vector< TH1 * > mv_hists, double nsigma)
T min(const caf::Proxy< T > &a, T b)
Var kUnity([](const caf::SRProxy *sr){return 1.f;})
return_type< T_y, T_loc, T_scale >::type normal_cdf(const T_y &y, const T_loc &mu, const T_scale &sigma)
Definition: normal_cdf.hpp:39
void events(int which)
Definition: Cana.C:52
c cd(1)
enum BeamMode string
unsigned int uint