systematics_summary_from_pred_interp.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Macro for making the nue systematics bar plots
4 //
5 ////////////////////////////////////////////////////////////////////////////////
6 
10 
13 
16 
17 #include "TCanvas.h"
18 
19 #include <iomanip>
20 
21 using namespace ana;
22 
23 struct ColumnDef
24 {
26  Flavors::Flavors_t flav;
29 };
30 
40 };
41 
43  {"Near-Far Differences",
44  "Neutrino Cross Sections",
45  "Detector Calibration",
46  "Beam Flux",
47  "Detector Response",
48  "Neutron Uncertainty",
49  "Muon Energy Scale",
50  "Normalization"
51 };
52 
53 void drawLabel(std::string s,double x=0.3,double y=0.93){
54  TLatex *tex = new TLatex(x,y,s.c_str());
55  tex->SetTextSize(1./30);
56  tex->SetNDC();
57  tex->SetTextAlign(12);
58  tex->Draw();
59 }
60 
61 double Integral(const Spectrum& s, double pot)
62 {
63  return s.Integral(pot);
64 }
65 
67 {
68  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
69  prelim->SetTextColor(kAzure);
70  prelim->SetNDC();
71  prelim->SetTextSize(2/30.);
72  prelim->SetTextAlign(32);
73  prelim->Draw();
74 }
75 
76 std::map<std::string, std::pair<double,double>> SumSysts(std::map<std::string, std::pair<double,double>> bars);
77 
78 void systematics_summary_from_pred_interp(std::string beam = "rhc", std::string extrapStr = "prop", bool sum = true)
79 {
80  std::vector<ColumnDef> cols;
81  if( beam == "fhc" )
82  cols = {{"$\\nu_e$ Signal", Flavors::kNuMuToNuE, Current::kCC, Sign::kNu},
84  {"$\\bar\\nu_e$ WS", Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu},
85  {"beam $\\nu_e$ CC", Flavors::kNuEToNuE, Current::kCC, Sign::kBoth},
86  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
87  {"$\\nu_\\mu$ CC", Flavors::kAllNuMu, Current::kCC, Sign::kBoth},
88  {"$\\nu_\\tau$ CC", Flavors::kAllNuTau, Current::kCC, Sign::kBoth}
89  };
90  else
91  cols = {{"$\\bar\\nu_e$ Signal",Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu},
92  {"Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
93  {"$\\nu_e$ WS", Flavors::kNuMuToNuE, Current::kCC, Sign::kNu},
94  {"beam $\\nu_e$ CC", Flavors::kNuEToNuE, Current::kCC, Sign::kBoth},
95  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
96  {"$\\nu_\\mu$ CC", Flavors::kAllNuMu, Current::kCC, Sign::kBoth},
97  {"$\\nu_\\tau$ CC", Flavors::kAllNuTau, Current::kCC, Sign::kBoth}
98  };
99 
100  std::cout << std::fixed;
101  std::cout << std::setprecision(2);
102 
103  // 2017 Best Fit
104  const double dCP = 1.21;
105  const double ssth23 = 0.556;
106  const double dmsq32 = 2.44e-3;
107 
109  calcosc->SetdCP ( dCP*(TMath::Pi()) );
110  calcosc->SetTh23 ( asin(sqrt(ssth23)) );
111  calcosc->SetDmsq32( dmsq32 );
112 
113  auto systs = getAllAna2018Systs(kNueAna2018,true,beam=="fhc" ? kFHC : kRHC);
114  //std::vector<const ISyst*> systs;
115  //AddJointAna2018XSecSysts(systs,kNueAna2018,false);
116 
117  double pot = beam == "fhc" ? kAna2018FHCPOT : kAna2018RHCPOT;
118 
119  const IPrediction *pred = GetNuePrediction2018(extrapStr,calcosc,true,beam,true);
120 
121  std::vector<Spectrum> noms;
122  for(ColumnDef col: cols) noms.push_back(pred->PredictComponent(calcosc, col.flav, col.curr, col.sign));
123 
124  noms[1] -= noms[0];
125 
126  double sigErr = 100/sqrt(Integral(noms[0], pot));
127  double bkgErr = 100/sqrt(Integral(noms[1], pot));
128 
129  std::cout<<"\\makebox[\\textwidth]{"<<std::endl;
130  std::cout<<"\\scriptsize"<<std::endl;
131  std::cout<<"\\begin{tabular}{l c c c c c c c c c c c c c c}"<<std::endl;
132  std::cout<<"\\multicolumn{15}{c}{\\bf{Systematic Summary (\\%)}} \\\\"<<std::endl;
133  std::cout << "{Systematic}";
134  for(ColumnDef col: cols) std::cout << " & \\multicolumn{2}{c}{" << col.name << "}";
135  std::cout << " \\\\" << std::endl;
136  std::cout << "\\hline" << std::endl;
137 
138  std::pair<double,double> quads[cols.size()];
139  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) quads[colIdx] = {0,0};
140 
141  std::map<std::string, std::pair<double,double>> barsSig, barsBkg;
142 
143  for(const ISyst* syst: systs){
144  std::cout << syst->LatexName();
145 
146  std::vector<Spectrum> shiftsUp, shiftsDn;
147  for(ColumnDef col: cols){
148  shiftsUp.push_back(pred->PredictComponentSyst(calcosc,
149  SystShifts(syst, +1),
150  col.flav, col.curr, col.sign));
151  shiftsDn.push_back(pred->PredictComponentSyst(calcosc,
152  SystShifts(syst, -1),
153  col.flav, col.curr, col.sign));
154  }
155 
156  shiftsUp[1] -= shiftsUp[0];
157  shiftsDn[1] -= shiftsDn[0];
158 
159  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
160  const Spectrum& nom = noms[colIdx];
161 
162  const double ratioUp = Integral(shiftsUp[colIdx], pot) / Integral(nom, pot);
163  const double errUp = 100*(ratioUp-1);
164 
165  const double ratioDn = Integral(shiftsDn[colIdx], pot) / Integral(nom, pot);
166  const double errDn = 100*(ratioDn-1);
167 
168  double min = std::min(errUp,errDn);
169  double max = std::max(errUp,errDn);
170  if(max < 0.01) max = 0;
171  if(min > -0.01) min = 0;
172 
173  quads[colIdx].first += min*min;
174  quads[colIdx].second += max*max;
175 
176  std::string minstr = "white";
177  std::string maxstr = "white";
178  if( fabs(min) >= 1 && fabs(min) < 2) minstr="Apricot";
179  if( fabs(min) >= 2 && fabs(min) < 3) minstr="YellowOrange";
180  if( fabs(min) >= 3 ) minstr="RedOrange";
181 
182  if( fabs(max) >= 1 && fabs(max) < 2) maxstr="Apricot";
183  if( fabs(max) >= 2 && fabs(max) < 3) maxstr="YellowOrange";
184  if( fabs(max) >= 3 ) maxstr="RedOrange";
185 
186  std::cout << " & \\cellcolor{" << minstr << "}" << (min == 0 ? "-" : "") << min << " & \\cellcolor{"<< maxstr << "}+" << max;
187 
188  min = fabs(min);
189  max = fabs(max);
190  if(cols[colIdx].name == "$\\nu_e$ Signal" || cols[colIdx].name == "$\\bar\\nu_e$ Signal")
191  barsSig[syst->LatexName()] = {min,max};
192  if(cols[colIdx].name == "Beam Bkg") barsBkg[syst->LatexName()] = {min,max};
193  } // end for colIdx
194 
195  std::cout << " \\\\" << std::endl;
196  } // end for syst
197 
198  std::cout << "\\hline" << std::endl;
199 
200  std::cout << "Quadrature";
201  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) std::cout << " & -" << sqrt(quads[colIdx].first) << " & " << sqrt(quads[colIdx].second);
202  std::cout << std::endl;
203  std::cout<<"\\end{tabular}"<<std::endl;
204  std::cout<<"}"<<std::endl;
205 
206  if(sum){
207  barsSig = SumSysts(barsSig);
208  barsBkg = SumSysts(barsBkg);
209 
211  std::cout<<"Summing..."<<std::endl;
213 
214  std::cout<<"\\begin{tabular}{l S@{/}S S@{/}S}"<<std::endl;
215  std::cout<<"\\multicolumn{5}{c}{\\bf{Systematic Summary (\\%)}} \\\\"<<std::endl;
216  std::cout<<"{Systematic}";
217  for(int colIdx = 0; colIdx < 2; ++colIdx) std::cout << " & \\multicolumn{2}{c}{" << cols[colIdx].name << "}";
218  std::cout<<" \\\\"<<std::endl;
219  std::cout << "\\hline" << std::endl;
220 
221  for(int i = 0;i < kNumSystTypes; ++i){
222  std::cout << systTypes[i] << " & "
223  << "-" << barsSig[systTypes[i]].first << " & " << barsSig[systTypes[i]].second << " & "
224  << "-" << barsBkg[systTypes[i]].first << " & " << barsBkg[systTypes[i]].second <<" \\\\"
225  << std::endl;
226  }
227 
228  std::cout << "\\hline" << std::endl;
229  std::cout << "Quadrature";
230  for(int colIdx = 0; colIdx < 2; ++colIdx) std::cout << " & -" << sqrt(quads[colIdx].first) << " & " << sqrt(quads[colIdx].second);
232  std::cout<<"\\end{tabular}"<<std::endl;
234  }
235 
236  std::string beamLab = beam == "fhc" ? "#nu Beam" : "#bar{#nu} Beam";
237  std::string sumStr = sum == true ? "short" : "full";
238 
239  TCanvas *cSig = new TCanvas("cSig","cSig");
240  ErrorBarChart(barsSig,{sigErr,sigErr},"Signal Uncertainty (%)");
241  gPad->SetLeftMargin(.3);
242  drawLabel(beamLab);
243  preliminary();
244  gPad->Print(("plots/syst_summary_"+sumStr+"_"+beam+"_sig.pdf").c_str());
245 
246  TCanvas *cSig2 = new TCanvas("cSig2","cSig2");
247  ErrorBarChart(barsSig,{0,0},"Signal Uncertainty (%)");
248  gPad->SetLeftMargin(.3);
249  drawLabel(beamLab);
250  preliminary();
251  gPad->Print(("plots/syst_summary_"+sumStr+"_"+beam+"_sig_nostat.pdf").c_str());
252 
253  TCanvas *cBkg = new TCanvas("cBkg","cBkg");
254  ErrorBarChart(barsBkg,{bkgErr,bkgErr},"Background Uncertainty (%)");
255  gPad->SetLeftMargin(.3);
256  drawLabel(beamLab);
257  preliminary();
258  gPad->Print(("plots/syst_summary_"+sumStr+"_"+beam+"_bkg.pdf").c_str());
259 
260  TCanvas *cBkg2 = new TCanvas("cBkg2","cBkg2");
261  ErrorBarChart(barsBkg,{0,0},"Background Uncertainty (%)");
262  gPad->SetLeftMargin(.3);
263  drawLabel(beamLab);
264  preliminary();
265  gPad->Print(("plots/syst_summary_"+sumStr+"_"+beam+"_bkg_nostat.pdf").c_str());
266 
267 
268 }
269 
270 std::map<std::string, std::pair<double,double>> SumSysts(std::map<std::string, std::pair<double,double>> bars)
271 {
272  std::map<std::string, std::pair<double,double>> summedbars;
273  std::pair<double,double> summedsyst[kNumSystTypes];
274  for(int i = 0; i < kNumSystTypes; ++i) summedsyst[i] = {0,0};
275 
276  for(auto it: bars){
277  if( it.first == "Acceptance Bkg RHC" ||
278  it.first == "Acceptance ND to FD Kinematics Signal RHC" ||
279  it.first == "Acceptance Bkg FHC" ||
280  it.first == "Acceptance ND to FD Kinematics Signal FHC" ||
281  it.first == "Michel Electrons Tagging Uncertainty" )
282  {
283  summedsyst[kNFDiffSys].first += it.second.first*it.second.first;
284  summedsyst[kNFDiffSys].second+= it.second.second*it.second.second;
285  continue;
286  }
287 
288  if( it.first == "Absolute Muon Energy Scale 2017" )
289  {
290  summedsyst[kMuScaleSys].first += it.second.first*it.second.first;
291  summedsyst[kMuScaleSys].second+= it.second.second*it.second.second;
292  continue;
293  }
294 
295  if( it.first == "FHC. Norm." ||
296  it.first == "RHC. Norm." )
297  {
298  summedsyst[kNormSys].first += it.second.first*it.second.first;
299  summedsyst[kNormSys].second+= it.second.second*it.second.second;
300  continue;
301  }
302 
303  if( it.first == "AbsCalib" ||
304  it.first == "RelCalib" ||
305  it.first == "CalibShape" )
306  {
307  summedsyst[kCalibSys].first += it.second.first*it.second.first;
308  summedsyst[kCalibSys].second+= it.second.second*it.second.second;
309  continue;
310  }
311 
312  if( it.first == "Neutron visible energy systematic 2018" )
313  {
314  summedsyst[kNeutronSys].first += it.second.first*it.second.first;
315  summedsyst[kNeutronSys].second += it.second.second*it.second.second;
316  continue;
317  }
318 
319  if( it.first == "#nu_{#tau} Scale" ||
320  it.first == "CCQEPauliSupViaKF" ||
321  it.first == "Coherent CC Scale" ||
322  it.first == "Coherent NC Scale" ||
323  it.first == "FormZone" ||
324  it.first == "FrAbs_N" ||
325  it.first == "FrCEx_N" ||
326  it.first == "FrElas_N" ||
327  it.first == "FrInel_pi" ||
328  it.first == "Genie PC 0" ||
329  it.first == "Genie PC 1" ||
330  it.first == "Genie PC 2" ||
331  it.first == "Genie PC 3" ||
332  it.first == "Genie PC 4" ||
333  it.first == "MEC 2018 shape, anti-neutrinos" ||
334  it.first == "MEC 2018 shape, neutrinos" ||
335  it.first == "MEC E_{#nu} shape, anti-neutrinos" ||
336  it.first == "MEC E_{#nu} shape, neutrinos" ||
337  it.first == "MEC initial state np fraction, anti-neutrinos" ||
338  it.first == "MEC initial state np fraction, neutrinos" ||
339  it.first == "MaCCQE" ||
340  it.first == "MaCCRES" ||
341  it.first == "MaNCRES" ||
342  it.first == "MvCCRES" ||
343  it.first == "RPA shape: enh2018" ||
344  it.first == "RPA shape: resonance events" ||
345  it.first == "Radiative corrections for #bar{#nu}_{e}" ||
346  it.first == "Radiative corrections for #nu_{e}" ||
347  it.first == "Second class currents" ||
348  it.first == "DIS vnCC1pi" )
349  {
350  summedsyst[kGenieSys].first += it.second.first*it.second.first;
351  summedsyst[kGenieSys].second+= it.second.second*it.second.second;
352  continue;
353  }
354 
355  if( it.first == "Cherenkov" || it.first == "Lightlevel" )
356  {
357  summedsyst[kBirksSys].first += it.second.first*it.second.first;
358  summedsyst[kBirksSys].second+= it.second.second*it.second.second;
359  continue;
360  }
361 
362  if( it.first == "Flux Component 00" ||
363  it.first == "Flux Component 01" ||
364  it.first == "Flux Component 02" ||
365  it.first == "Flux Component 03" ||
366  it.first == "Flux Component 04" )
367  {
368  summedsyst[kBeamSys].first += it.second.first*it.second.first;
369  summedsyst[kBeamSys].second+= it.second.second*it.second.second;
370  continue;
371  }
372 
373  std::cout<<"it.first == \""<<it.first<<"\" || "<<std::endl;
374  }
375 
376  for (int i = 0; i < kNumSystTypes; ++i){
377  summedbars[systTypes[i]] = {sqrt(summedsyst[i].first),sqrt(summedsyst[i].second)};
378  }
379 
380  return summedbars;
381 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
double ssth23
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
const TString systTypes[kNumSystTypes]
virtual void SetdCP(const T &dCP)=0
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const double kAna2018RHCPOT
Definition: Exposures.h:208
virtual Spectrum PredictComponentSyst(osc::IOscCalculator *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
const XML_Char * s
Definition: expat.h:262
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
Definition: IPrediction.h:39
double Integral(const Spectrum &s, double pot)
General interface to any calculator that lets you set the parameters.
Interactions of both types.
Definition: IPrediction.h:42
double dCP
Int_t col[ntarg]
Definition: Style.C:29
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalculator *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
void drawLabel(std::string s, double x=0.3, double y=0.93)
#define pot
const int cols[3]
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
TLatex * prelim
Definition: Xsec_final.C:133
Neutrinos-only.
Definition: IPrediction.h:49
def sign(x)
Definition: canMan.py:204
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void systematics_summary_from_pred_interp(std::string beam="rhc", std::string extrapStr="prop", bool sum=true)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double dmsq32
string syst
Definition: plotSysts.py:176
TLatex * tex
Definition: f2_nu.C:499
Neutral-current interactions.
Definition: IPrediction.h:40
void ErrorBarChart(const std::map< std::string, std::pair< double, double >> &systErrors, const std::pair< double, double > &statErr, const std::string &label)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:730
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const double kAna2018FHCPOT
Definition: Exposures.h:207
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
Double_t sum
Definition: plot.C:31
T asin(T number)
Definition: d0nt_math.hpp:60