Classes | Enumerations | Functions | Variables
systematics_summary_from_pred_interp.C File Reference
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "OscLib/IOscCalc.h"
#include "OscLib/OscCalcDumb.h"
#include "3FlavorAna/Ana2018/nue/joint_fit_2018_tools.h"
#include "3FlavorAna/Systs/JointAna2018Systs.h"
#include "TCanvas.h"
#include <iomanip>

Go to the source code of this file.

Classes

struct  ColumnDef
 

Enumerations

enum  NueSystType {
  kNFDiffSys, kGenieSys, kCalibSys, kBeamSys,
  kBirksSys, kNeutronSys, kMuScaleSys, kNormSys,
  kNumSystTypes, kNFDiffSys, kGenieSys, kCalibSys,
  kBeamSys, kBirksSys, kNeutronSys, kMuScaleSys,
  kNormSys, kNumSystTypes, kNFUncorSys, kGenieSys,
  kCalibSys, kBeamSys, kBirksSys, kNeutronSys,
  kLeptonRecoSys, kNumSystTypes
}
 

Functions

void drawLabel (std::string s, double x=0.3, double y=0.93)
 
double Integral (const Spectrum &s, double pot)
 
void preliminary ()
 
std::map< std::string, std::pair< double, double > > SumSysts (std::map< std::string, std::pair< double, double >> bars)
 
void systematics_summary_from_pred_interp (std::string beam="rhc", std::string extrapStr="prop", bool sum=true)
 

Variables

std::string systTypes [kNumSystTypes]
 

Enumeration Type Documentation

Enumerator
kNFDiffSys 
kGenieSys 
kCalibSys 
kBeamSys 
kBirksSys 
kNeutronSys 
kMuScaleSys 
kNormSys 
kNumSystTypes 
kNFDiffSys 
kGenieSys 
kCalibSys 
kBeamSys 
kBirksSys 
kNeutronSys 
kMuScaleSys 
kNormSys 
kNumSystTypes 
kNFUncorSys 
kGenieSys 
kCalibSys 
kBeamSys 
kBirksSys 
kNeutronSys 
kLeptonRecoSys 
kNumSystTypes 

Definition at line 31 of file systematics_summary_from_pred_interp.C.

Function Documentation

void drawLabel ( std::string  s,
double  x = 0.3,
double  y = 0.93 
)

Definition at line 53 of file systematics_summary_from_pred_interp.C.

References tex, submit_syst::x, and submit_syst::y.

Referenced by systematics_summary_from_pred_interp().

53  {
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 }
const XML_Char * s
Definition: expat.h:262
TLatex * tex
Definition: f2_nu.C:499
double Integral ( const Spectrum s,
double  pot 
)

Definition at line 61 of file systematics_summary_from_pred_interp.C.

References ana::Spectrum::Integral().

Referenced by systematics_summary_from_pred_interp().

62 {
63  return s.Integral(pot);
64 }
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:277
#define pot
void preliminary ( )

Definition at line 66 of file systematics_summary_from_pred_interp.C.

References prelim, string, and SumSysts().

Referenced by systematics_summary_from_pred_interp().

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 }
TLatex * prelim
Definition: Xsec_final.C:133
std::map<std::string, std::pair<double,double> > SumSysts ( std::map< std::string, std::pair< double, double >>  bars)

Definition at line 270 of file systematics_summary_from_pred_interp.C.

References om::cout, allTimeWatchdog::endl, MECModelEnuComparisons::i, it, kBeamSys, kBirksSys, kCalibSys, kGenieSys, kMuScaleSys, kNeutronSys, kNFDiffSys, kNormSys, kNumSystTypes, moon_position_table_new3::second, std::sqrt(), and systTypes.

Referenced by preliminary(), and systematics_summary_from_pred_interp().

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 }
set< int >::iterator it
std::string systTypes[kNumSystTypes]
T sqrt(T number)
Definition: d0nt_math.hpp:156
OStream cout
Definition: OStream.cxx:6
void systematics_summary_from_pred_interp ( std::string  beam = "rhc",
std::string  extrapStr = "prop",
bool  sum = true 
)

Definition at line 78 of file systematics_summary_from_pred_interp.C.

References std::asin(), POTSpillRate::beam, col, cols, om::cout, dCP, ana::DefaultOscCalc(), dmsq32, drawLabel(), allTimeWatchdog::endl, ana::ErrorBarChart(), ana::getAllAna2018Systs(), GetNuePrediction2018(), MECModelEnuComparisons::i, Integral(), ana::Flavors::kAll, ana::Flavors::kAllNuMu, ana::Flavors::kAllNuTau, ana::kAna2018FHCPOT, ana::kAna2018RHCPOT, ana::Sign::kAntiNu, ana::Current::kBoth, ana::Sign::kBoth, ana::Current::kCC, ana::kFHC, ana::Current::kNC, ana::Sign::kNu, ana::kNueAna2018, ana::Flavors::kNuEToNuE, kNumSystTypes, ana::Flavors::kNuMuToNuE, ana::kRHC, cet::sqlite::max(), std::max(), min(), std::min(), pot, plot_validation_datamc::pred, ana::IPrediction::PredictComponent(), ana::IPrediction::PredictComponentSyst(), preliminary(), moon_position_table_new3::second, osc::_IOscCalcAdjustable< T >::SetdCP(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh23(), std::sqrt(), ssth23, string, sum, SumSysts(), systs, and systTypes.

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 }
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
double ssth23
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
(&#39; appearance&#39;)
Definition: IPrediction.h:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
std::string systTypes[kNumSystTypes]
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
Definition: IPrediction.cxx:79
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const double kAna2018RHCPOT
Definition: Exposures.h:208
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
Definition: IPrediction.h:39
const int cols[3]
double dCP
Int_t col[ntarg]
Definition: Style.C:29
#define pot
void drawLabel(std::string s, double x=0.3, double y=0.93)
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
double dmsq32
double Integral(const Spectrum &s, double pot)
virtual void SetTh23(const T &th23)=0
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:800
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
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 max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string

Variable Documentation

Initial value:
=
{"Near-Far Differences",
"Neutrino Cross Sections",
"Detector Calibration",
"Beam Flux",
"Detector Response",
"Neutron Uncertainty",
"Muon Energy Scale",
"Normalization"
}

Definition at line 42 of file systematics_summary_from_pred_interp.C.

Referenced by SumSysts(), and systematics_summary_from_pred_interp().