barchart_extrapcomparison.C
Go to the documentation of this file.
2
3 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Systs/Systs.h"
7
9
12
13 #include "Utilities/func/MathUtil.h"
14
15 #include "TCanvas.h"
16 #include "TH2.h"
17 #include "THStack.h"
18 #include "TLegend.h"
19 #include "TH1.h"
20 #include "TLine.h"
21 #include "TFile.h"
22
23 #include <iostream>
24 #include <cmath>
25
26 using namespace ana;
27
28
29 void DrawCompBarChart(bool isSig,
30  const std::vector<std::pair<double,std::string> >& ne,
31  const std::vector<std::pair<double,std::string> >& pb,
32  const std::vector<std::pair<double,std::string> >& cb={})
33 {
34  TCanvas *canBar = new TCanvas();
35
36  int nbins = pb.size()+1;
37  // Systematics are on top, total systematic and (optional) statistical
38  // error bar are at the bottom.
39  int firstSystBin = 1;
40
41  // set up bin boundaries that depend on whether we're plotting sig / bkgd
42  double ne_x1 = isSig ? 0.5 : 0.65;
43  double ne_x2 = 0.95;
44
45  double pb_x1 = isSig ? 0.05 : 0.35;
46  double pb_x2 = isSig ? 0.5 : 0.65;
47
48  double cb_x1 = 0.05;
49  double cb_x2 = 0.35;
50
51  // Sum all the systematics in qaudrature
52  double ntot = 0;
53  for (std::pair<double, std::string> pair : ne) ntot += pair.first*pair.first;
54  ntot = sqrt(ntot);
55  double ctot = 0;
56  for (std::pair<double, std::string> pair : cb) ctot += pair.first*pair.first;
57  ctot = sqrt(ctot);
58  double ptot = 0;
59  for (std::pair<double, std::string> pair : pb) ptot += pair.first*pair.first;
60  ptot = sqrt(ptot);
61
62  // x-axis range is a multiple of 5% with at least 2% padding
63  double xrange = 20;
64  TH2* axes = new TH2F(UniqueName().c_str(), ";Background uncertainty (%)",
65  10, -xrange, +xrange, nbins, 0, nbins);
66  if (isSig)
67  axes->GetXaxis()->SetTitle("Signal uncertainty (%)");
68  axes->GetXaxis()->CenterTitle();
69
70  TAxis *yax = axes->GetYaxis();
71  for(unsigned int i = 0; i < pb.size(); ++i)
72  yax->SetBinLabel(i+1+1, pb[i].second.c_str());
73  yax->SetBinLabel(1, "Total syst. error");
74  // Make the labels approximately match the x axis title
75  yax->SetLabelSize(.06);
76
77  axes->Draw();
78
79  // Blue boxes for each systematic
80  for(unsigned int i = 0; i < ne.size(); ++i){
81  TBox* nbox = new TBox(-ne[i].first, i+1+ne_x1,
82  +ne[i].first, i+1+ne_x2);
83  nbox->SetFillColor(kGreen+2);
84  nbox->Draw();
85  TBox* pbox = new TBox(-pb[i].first, i+1+pb_x1,
86  +pb[i].first, i+1+pb_x2);
87  pbox->SetFillColor(kAzure-8);
88  pbox->Draw();
89  if (!isSig){ // Only plot combo if we're doing background
90  TBox* cbox = new TBox(-cb[i].first, i+1+cb_x1,
91  +cb[i].first, i+1+cb_x2);
92  cbox->SetFillColor(kRed-4);
93  cbox->Draw();
94  }
95  }
96
97  // Total systematic error
98  TBox* nsyst = new TBox(-ntot, ne_x1,
99  +ntot, ne_x2);
100  nsyst->SetFillColor(kGreen+2);
101  nsyst->Draw();
102  TBox* psyst = new TBox(-ptot, pb_x1,
103  +ptot, pb_x2);
104  psyst->SetFillColor(kAzure-8);
105  psyst->Draw();
106  if (!isSig){
107  TBox* csyst = new TBox(-ctot, cb_x1,
108  +ctot, cb_x2);
109  csyst->SetFillColor(kRed-4);
110  csyst->Draw();
111  }
112  if (isSig)
113  std::cout << "Totals: " << ntot << " " << ptot << std::endl;
114  else
115  std::cout << "Totals: " << ntot << " " << ptot << " " << ctot << std::endl;
116
117
118  // Separate individual systematics from the total and statistical errors
119  TLine* div = new TLine(-xrange, 1, +xrange, 1);
120  div->SetLineWidth(2);
121  div->Draw();
122
123  // Vertical line marking the symmetry point
124  TLine* zero = new TLine(0, 0, 0, nbins);
125  zero->SetLineStyle(7);
126  zero->SetLineWidth(2);
127  zero->Draw();
128
129  // Leave enough space for the systematic labels
131
132  TH1* ncolor = new TH1F(UniqueName().c_str(), "", 10, -xrange, +xrange);
133  TH1* pcolor = new TH1F(UniqueName().c_str(), "", 10, -xrange, +xrange);
134  TH1* ccolor = new TH1F(UniqueName().c_str(), "", 10, -xrange, +xrange);
135  ncolor->SetFillColor(kGreen+2);
136  pcolor->SetFillColor(kAzure-8);
137  ccolor->SetFillColor(kRed-4);
138
139  double yleghi = isSig ? 0.50 : 0.60;
140  TLegend *leg = new TLegend(0.32,0.32,0.48,yleghi);
142  if (isSig)
144  else{
147  }
148  leg->Draw();
149
150  std::string name = "Fig_NueExtrapCompare";
151  if (isSig) name += "Sig";
152  else name += "Bkgd";
153  canBar->SaveAs((name+".png").c_str());
154  canBar->SaveAs((name+".pdf").c_str());
155 }
156
158 {
162  Spectrum nc = pred->PredictComponent(&calc, Flavors::kAll,
164  Spectrum ne = pred->PredictComponent(&calc, Flavors::kNuEToNuE,
166  double pot = 9.48e20; // Absolute norm not important
167  TH1 *h = (cc+nc+ne).ToTH1(pot);
168  return h->Integral();
169 }
170
172 {
176  double pot = 9.48e20; // Absolute norm not important
177  TH1 *h = sig.ToTH1(pot);
178  return h->Integral();
179 }
180
181 double GetDif(TH1 *pred, TH1 *fdmc)
182 {
183  double diff = 0;
184  int count = 0;
185  for (int i = 1; i <= pred->GetNbinsX(); i++){
186  if (fdmc->GetBinContent(i) == 0) continue;
187  count ++;
188  double cur = pred->GetBinContent(i)-fdmc->GetBinContent(i);
189  diff += std::abs(cur);
190  //diff += util::sqr(cur/fdmc->GetBinContent(i));
191  }
192  return diff/count;
193 }
194
196  std::vector< std::pair<double, std::string> > &extrap_sig,
197  std::vector< std::pair<double, std::string> > &mc_sig,
198  std::vector< std::pair<double, std::string> > &combo_bg,
199  std::vector< std::pair<double, std::string> > &prop_bg,
200  std::vector< std::pair<double, std::string> > &mc_bg)
201 {
202  std::string n_nxpUp = "pred_nxp_" +name+"/plusOne";
203  std::string n_nxpDn = "pred_nxp_" +name+"/minusOne";
204  std::string n_propUp = "pred_xp_prop_" +name+"/plusOne";
205  std::string n_propDn = "pred_xp_prop_" +name+"/minusOne";
206  std::string n_comboUp = "pred_xp_combo_"+name+"/plusOne";
207  std::string n_comboDn = "pred_xp_combo_"+name+"/minusOne";
208
209  IPrediction* mcUp = LoadFromFile<IPrediction>(fname,n_nxpUp) .release();
210  IPrediction* mcDn = LoadFromFile<IPrediction>(fname,n_nxpDn) .release();
211  IPrediction* propUp = LoadFromFile<IPrediction>(fname,n_propUp) .release();
212  IPrediction* propDn = LoadFromFile<IPrediction>(fname,n_propDn) .release();
215
216  double valextrap_sig =
217  (GetSig(propUp)-GetSig(propDn)) / (GetSig(propUp)+GetSig(propDn));
218  double valmc_sig =
219  (GetSig(mcUp)-GetSig(mcDn)) / (GetSig(mcUp)+GetSig(mcDn));
220  double valcombo_bg =
221  (GetBkgd(comboUp)-GetBkgd(comboDn)) / (GetBkgd(comboUp)+GetBkgd(comboDn));
222  double valprop_bg =
223  (GetBkgd(propUp)-GetBkgd(propDn)) / (GetBkgd(propUp)+GetBkgd(propDn));
224  double valmc_bg =
225  (GetBkgd(mcUp)-GetBkgd(mcDn)) / (GetBkgd(mcUp)+GetBkgd(mcDn));
226
227
228  std::pair<double,std::string> cextrap_sig(100*valextrap_sig,name);
229  std::pair<double,std::string> cmc_sig(100*valmc_sig,name);
230  std::pair<double,std::string> ccombo_bg(100*valcombo_bg,name);
231  std::pair<double,std::string> cprop_bg(100*valprop_bg,name);
232  std::pair<double,std::string> cmc_bg(100*valmc_bg,name);
233
234  extrap_sig.push_back(cextrap_sig);
235  mc_sig.push_back(cmc_sig);
236  combo_bg.push_back(ccombo_bg);
237  prop_bg.push_back(cprop_bg);
238  mc_bg.push_back(cmc_bg);
239
240  std::cout << name << " " << valcombo_bg << " " << valprop_bg << " " << valmc_bg << std::endl;
241  std::cout << name << " " << valextrap_sig << " " << valmc_sig << std::endl;
242 }
243
244
245 void barchart_extrapcomparison(std::string fname = "/nova/ana/nu_e_ana/Ana2017/Predictions/provisional/pred_nue_reduced_v4.root")
246 {
247  double res = 0;
248
249  std::vector<std::string> names =
250  {"MaCCRES","MaNCRES","MvCCRES","DISvnCC1pi","MECq0Shape",
251  "RPAShapeRES","MaCCQE_reduced","ppfx_pc00","ppfx_pc01"};
252
253
254
255  std::vector< std::pair<double, std::string> > extrap_sig, mc_sig;
256  std::vector< std::pair<double, std::string> > combo_bg, prop_bg, mc_bg;
257
258  std::cout << "Numbers for total res: " << std::endl;
259  for (std::string name : names){
261  extrap_sig, mc_sig,
262  combo_bg, prop_bg, mc_bg);
263  }
264
265  DrawCompBarChart(true,mc_sig,extrap_sig);
266  DrawCompBarChart(false,mc_bg,prop_bg,combo_bg);
267
268  /*
269  TCanvas *can = new TCanvas();
270  DrawCompBarChart(ntot,ptot,ctot);
271  can->SaveAs("FDBkgdRes_PropCombo_comparison.png");
272  */
273 }
const XML_Char * name
Definition: expat.h:151
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
void DrawCompBarChart(bool isSig, const std::vector< std::pair< double, std::string > > &ne, const std::vector< std::pair< double, std::string > > &pb, const std::vector< std::pair< double, std::string > > &cb={})
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
double GetDif(TH1 *pred, TH1 *fdmc)
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
double GetSig(IPrediction *pred)
static const double pb
Definition: Units.h:90
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
#define pot
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void GetBiases(std::string fname, std::string name, std::vector< std::pair< double, std::string > > &extrap_sig, std::vector< std::pair< double, std::string > > &mc_sig, std::vector< std::pair< double, std::string > > &combo_bg, std::vector< std::pair< double, std::string > > &prop_bg, std::vector< std::pair< double, std::string > > &mc_bg)
void cc()
Definition: test_ana.C:28
Neutral-current interactions.
Definition: IPrediction.h:40
double GetBkgd(IPrediction *pred)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
All neutrinos, any flavor.
Definition: IPrediction.h:26
float count
Definition: make_cosmics.C:27
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
h
Definition: demo3.py:41
void barchart_extrapcomparison(std::string fname="/nova/ana/nu_e_ana/Ana2017/Predictions/provisional/pred_nue_reduced_v4.root")