Classes | Enumerations | Functions | Variables
systematics_extrap_comp_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 myBarChart (std::string xlabel, std::vector< std::string > ylabels, std::vector< std::vector< std::pair< double, double >>> errs, std::vector< std::pair< double, double >> sum, std::vector< std::string > legend)
 
std::map< std::string, std::pair< double, double > > SumSysts (std::map< std::string, std::pair< double, double >> bars)
 
void drawLabel (std::string s, double x=0.3, double y=0.93)
 
double Integral (const Spectrum &s, double pot)
 
void preliminary ()
 
void systematics_extrap_comp_from_pred_interp (std::string beam="fhc", 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 29 of file systematics_extrap_comp_from_pred_interp.C.

Function Documentation

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

Definition at line 60 of file systematics_extrap_comp_from_pred_interp.C.

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

Referenced by systematics_extrap_comp_from_pred_interp().

60  {
61  TLatex *tex = new TLatex(x,y,s.c_str());
62  tex->SetTextSize(1./30);
63  tex->SetNDC();
64  tex->SetTextAlign(12);
65  tex->Draw();
66 }
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 68 of file systematics_extrap_comp_from_pred_interp.C.

References ana::Spectrum::Integral().

Referenced by systematics_extrap_comp_from_pred_interp().

69 {
70  return s.Integral(pot);
71 }
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 myBarChart ( std::string  xlabel,
std::vector< std::string ylabels,
std::vector< std::vector< std::pair< double, double >>>  errs,
std::vector< std::pair< double, double >>  sum,
std::vector< std::string legend 
)

Definition at line 238 of file systematics_extrap_comp_from_pred_interp.C.

References file_size_ana::axes, colors, MECModelEnuComparisons::i, calib::j, kGreen, kRed, submit_hadd::l, cet::sqlite::max(), std::max(), nbins, moon_position_table_new3::second, gen_hdf5record::size, sum, and zero().

Referenced by systematics_extrap_comp_from_pred_interp(), and systematics_summary_from_pred_interp().

243 {
244  // add more colors as needed
245  std::vector<int> colors = {kAzure-8,kRed-7,kGreen+2};
246 
247  int nbins = ylabels.size() + 1;
248 
249  double xrange = 0;
250  double max = 0;
251  for(unsigned int i = 0; i < sum.size(); ++i)
252  if(std::max(sum[i].first,sum[i].second) > max)
253  max = std::max(sum[i].first,sum[i].second);
254  while(xrange < max + 2 ) xrange += 5;
255 
256  TH2* axes = new TH2F("", (";" + xlabel).c_str(),
257  10, -xrange, +xrange, nbins, 0, nbins);
258  axes->GetXaxis()->CenterTitle();
259 
260  TAxis* yax = axes->GetYaxis();
261  for(unsigned int i = 0; i < ylabels.size(); ++i)
262  yax->SetBinLabel(i+2, ylabels[i].c_str());
263  yax->SetBinLabel(1, "Systematic Uncertainty");
264  yax->SetLabelSize(.06);
265 
266  axes->Draw();
267 
268  TLegend *l = new TLegend(0.325,0.7,0.575,0.825);
269  l->SetFillStyle(0);
270 
271  for(unsigned int i = 0; i < errs.size(); ++i){
272  double size = 0.8/errs.size();
273  for(unsigned int j = 0; j < errs[i].size(); ++j){
274  TBox* box = new TBox(-errs[i][j].first, j+1+0.1+(errs.size()-i-1)*size,
275  +errs[i][j].second, j+1+0.1+(errs.size()-i)*size);
276  box->SetFillColor(colors[i]);
277  box->Draw();
278  }
279 
280  TBox* syst = new TBox(-sum[i].first, 0.1+(errs.size()-i-1)*size,
281  +sum[i].second, 0.1+(errs.size()-i)*size);
282  syst->SetFillColor(colors[i]);
283  syst->Draw();
284 
285  l->AddEntry(syst,legend[i].c_str(),"f");
286 
287  }
288 
289  l->Draw();
290 
291  TLine* div = new TLine(-xrange, 1, +xrange, 1);
292  div->SetLineWidth(2);
293  div->Draw();
294 
295  TLine* zero = new TLine(0, 0, 0, nbins);
296  zero->SetLineStyle(7);
297  zero->SetLineWidth(2);
298  zero->Draw();
299 
300  gPad->SetLeftMargin(.3);
301 }
T max(const caf::Proxy< T > &a, T b)
enum BeamMode kRed
const int nbins
Definition: cellShifts.C:15
int colors[6]
Definition: tools.h:1
const double j
Definition: BetheBloch.cxx:29
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
enum BeamMode kGreen
auto zero()
Definition: PMNS.cxx:55
void preliminary ( )

Definition at line 73 of file systematics_extrap_comp_from_pred_interp.C.

References prelim.

Referenced by systematics_extrap_comp_from_pred_interp().

74  {
75  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
76  prelim->SetTextColor(kAzure);
77  prelim->SetNDC();
78  prelim->SetTextSize(2/30.);
79  prelim->SetTextAlign(32);
80  prelim->Draw();
81  }
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 303 of file systematics_extrap_comp_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 systematics_extrap_comp_from_pred_interp().

304  {
305  std::map<std::string, std::pair<double,double>> summedbars;
306  std::pair<double,double> summedsyst[kNumSystTypes];
307  for(int i = 0; i < kNumSystTypes; ++i) summedsyst[i] = {0,0};
308 
309  for(auto it: bars){
310  if( it.first == "Acceptance Bkg RHC" ||
311  it.first == "Acceptance ND to FD Kinematics Signal RHC" ||
312  it.first == "Acceptance Bkg FHC" ||
313  it.first == "Acceptance ND to FD Kinematics Signal FHC" ||
314  it.first == "Michel Electrons Tagging Uncertainty" )
315  {
316  summedsyst[kNFDiffSys].first += it.second.first*it.second.first;
317  summedsyst[kNFDiffSys].second+= it.second.second*it.second.second;
318  continue;
319  }
320 
321  if( it.first == "Absolute Muon Energy Scale 2017" )
322  {
323  summedsyst[kMuScaleSys].first += it.second.first*it.second.first;
324  summedsyst[kMuScaleSys].second+= it.second.second*it.second.second;
325  continue;
326  }
327 
328  if( it.first == "FHC. Norm." ||
329  it.first == "RHC. Norm." )
330  {
331  summedsyst[kNormSys].first += it.second.first*it.second.first;
332  summedsyst[kNormSys].second+= it.second.second*it.second.second;
333  continue;
334  }
335 
336  if( it.first == "AbsCalib" ||
337  it.first == "RelCalib" ||
338  it.first == "CalibShape" )
339  {
340  summedsyst[kCalibSys].first += it.second.first*it.second.first;
341  summedsyst[kCalibSys].second+= it.second.second*it.second.second;
342  continue;
343  }
344 
345  if( it.first == "Neutron visible energy systematic 2018" )
346  {
347  summedsyst[kNeutronSys].first += it.second.first*it.second.first;
348  summedsyst[kNeutronSys].second+= it.second.second*it.second.second;
349  continue;
350  }
351 
352  if( it.first == "#nu_{#tau} Scale" ||
353  it.first == "CCQEPauliSupViaKF" ||
354  it.first == "Coherent CC Scale" ||
355  it.first == "Coherent NC Scale" ||
356  it.first == "FormZone" ||
357  it.first == "FrAbs_N" ||
358  it.first == "FrCEx_N" ||
359  it.first == "FrElas_N" ||
360  it.first == "FrInel_pi" ||
361  it.first == "Genie PC 0" ||
362  it.first == "Genie PC 1" ||
363  it.first == "Genie PC 2" ||
364  it.first == "Genie PC 3" ||
365  it.first == "Genie PC 4" ||
366  it.first == "MEC 2018 shape, anti-neutrinos" ||
367  it.first == "MEC 2018 shape, neutrinos" ||
368  it.first == "MEC E_{#nu} shape, anti-neutrinos" ||
369  it.first == "MEC E_{#nu} shape, neutrinos" ||
370  it.first == "MEC initial state np fraction, anti-neutrinos" ||
371  it.first == "MEC initial state np fraction, neutrinos" ||
372  it.first == "MaCCQE" ||
373  it.first == "MaCCRES" ||
374  it.first == "MaNCRES" ||
375  it.first == "MvCCRES" ||
376  it.first == "RPA shape: enh2018" ||
377  it.first == "RPA shape: resonance events" ||
378  it.first == "Radiative corrections for #bar{#nu}_{e}" ||
379  it.first == "Radiative corrections for #nu_{e}" ||
380  it.first == "Second class currents" ||
381  it.first == "DIS vnCC1pi" )
382  {
383  summedsyst[kGenieSys].first += it.second.first*it.second.first;
384  summedsyst[kGenieSys].second+= it.second.second*it.second.second;
385  continue;
386  }
387 
388  if( it.first == "Cherenkov" || it.first == "Lightlevel" )
389  {
390  summedsyst[kBirksSys].first += it.second.first*it.second.first;
391  summedsyst[kBirksSys].second+= it.second.second*it.second.second;
392  continue;
393  }
394 
395  if( it.first == "Flux Component 00" ||
396  it.first == "Flux Component 01" ||
397  it.first == "Flux Component 02" ||
398  it.first == "Flux Component 03" ||
399  it.first == "Flux Component 04" )
400  {
401  summedsyst[kBeamSys].first += it.second.first*it.second.first;
402  summedsyst[kBeamSys].second+= it.second.second*it.second.second;
403  continue;
404  }
405 
406  std::cout<<"it.first == \""<<it.first<<"\" || "<<std::endl;
407  }
408 
409  for (int i = 0; i < kNumSystTypes; ++i){
410  summedbars[systTypes[i]] = {sqrt(summedsyst[i].first),sqrt(summedsyst[i].second)};
411  }
412 
413  return summedbars;
414  }
set< int >::iterator it
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::string systTypes[kNumSystTypes]
OStream cout
Definition: OStream.cxx:6
void systematics_extrap_comp_from_pred_interp ( std::string  beam = "fhc",
bool  sum = true 
)

Definition at line 83 of file systematics_extrap_comp_from_pred_interp.C.

References std::asin(), POTSpillRate::beam, col, cols, dCP, ana::DefaultOscCalc(), dmsq32, drawLabel(), fSig, ana::getAllAna2018Systs(), GetNuePrediction2018(), MECModelEnuComparisons::i, Integral(), calib::j, ana::Flavors::kAll, ana::Sign::kAntiNu, ana::Current::kBoth, ana::Sign::kBoth, ana::Current::kCC, ana::kFHC, ana::Sign::kNu, ana::kNueAna2018, kNumSystTypes, ana::Flavors::kNuMuToNuE, ana::kRHC, PandAna.Demos.pi0_spectra::labels, cet::sqlite::max(), std::max(), min(), std::min(), myBarChart(), 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.

84 {
85  std::vector<ColumnDef> cols;
86  if( beam == "fhc" )
87  cols = {{"$\\nu_e$ Signal", Flavors::kNuMuToNuE, Current::kCC, Sign::kNu},
88  {"Total Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth}
89  };
90  else
91  cols = {{"$\\bar\\nu_e$ Signal", Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu},
92  {"Total Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth}
93  };
94  // 2017 Best Fit
95  const double dCP = 1.21;
96  const double ssth23 = 0.556;
97  const double dmsq32 = 2.44e-3;
98 
100  calcosc->SetdCP ( dCP*(TMath::Pi()) );
101  calcosc->SetTh23 ( asin(sqrt(ssth23)) );
102  calcosc->SetDmsq32( dmsq32 );
103 
104  auto systs = getAllAna2018Systs(kNueAna2018,true,beam=="fhc" ? kFHC : kRHC);
105 
106  std::vector<std::string> extrapStr;
107  std::vector<std::string> extrapLab;
108  if(beam == "fhc"){
109  extrapStr = {"noextrap","combo"};
110  extrapLab = {"Not Extrapolated","Extrapolated"};
111  }
112  else{
113  extrapStr = {"noextrap","prop"};
114  extrapLab = {"Not Extrapolated","Extrapolated"};
115  }
116  std::string beamLab = beam == "fhc" ? "#nu Beam" : "#bar{#nu} Beam";
117 
118  std::vector<std::map<std::string,std::pair<double,double>>> sigs, bkgs;
119  std::vector<std::pair<double,double>> sigQuad, bkgQuad;
120 
121  for(unsigned int i = 0; i < extrapStr.size(); ++i){
122  const IPrediction *pred = GetNuePrediction2018(extrapStr[i],calcosc,true,beam,true);
123 
124  std::vector<Spectrum> noms;
125  for(ColumnDef col: cols) noms.push_back(pred->PredictComponent(calcosc, col.flav, col.curr, col.sign));
126 
127  noms[1] -= noms[0];
128 
129  std::map<std::string, std::pair<double,double>> barsSig, barsBkg;
130 
131  std::pair<double,double> quads[cols.size()];
132  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) quads[colIdx] = {0,0};
133 
134  for(const ISyst* syst: systs){
135  std::vector<Spectrum> shiftsUp, shiftsDn;
136  for(ColumnDef col: cols){
137  shiftsUp.push_back(pred->PredictComponentSyst(calcosc,
138  SystShifts(syst, +1),
139  col.flav, col.curr, col.sign));
140  shiftsDn.push_back(pred->PredictComponentSyst(calcosc,
141  SystShifts(syst, -1),
142  col.flav, col.curr, col.sign));
143  }
144 
145  shiftsUp[1] -= shiftsUp[0];
146  shiftsDn[1] -= shiftsDn[0];
147 
148  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
149  const Spectrum& nom = noms[colIdx];
150 
151  // POT cancels, is irrelevant
152  const double ratioUp = Integral(shiftsUp[colIdx], 1e20) / Integral(nom, 1e20);
153  const double errUp = 100*(ratioUp-1);
154 
155  const double ratioDn = Integral(shiftsDn[colIdx], 1e20) / Integral(nom, 1e20);
156  const double errDn = 100*(ratioDn-1);
157 
158  double min = std::min(errUp,errDn);
159  double max = std::max(errUp,errDn);
160  if(max < 0) max = 0;
161  if(min > 0) min = 0;
162 
163  quads[colIdx].first += min*min;
164  quads[colIdx].second += max*max;
165 
166  min = fabs(min);
167  max = fabs(max);
168  if(cols[colIdx].name == "$\\nu_e$ Signal" || cols[colIdx].name == "$\\bar\\nu_e$ Signal")
169  barsSig[syst->LatexName()] = {min,max};
170  if(cols[colIdx].name == "Total Beam Bkg")
171  barsBkg[syst->LatexName()] = {min,max};
172  } // end for colIdx
173  } // end for syst
174 
175  if(sum){
176  barsSig = SumSysts(barsSig);
177  barsBkg = SumSysts(barsBkg);
178  }
179 
180  sigs.push_back(barsSig);
181  bkgs.push_back(barsBkg);
182 
183  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
184  quads[colIdx].first = sqrt(quads[colIdx].first);
185  quads[colIdx].second = sqrt(quads[colIdx].second);
186  }
187 
188  sigQuad.push_back(quads[0]);
189  bkgQuad.push_back(quads[1]);
190  } // end for extrap
191 
192  std::vector<std::string> labels;
193  std::vector<std::vector<std::pair<double,double>>> fSig;
194  std::vector<std::vector<std::pair<double,double>>> fBkg;
195 
196  for(unsigned int i = 0; i < extrapStr.size(); ++i){
197  std::vector<std::pair<double,double>> tempSig;
198  std::vector<std::pair<double,double>> tempBkg;
199 
200  if(sum){
201  for(int j = 0; j < kNumSystTypes; ++j){
202  if(i == 0) labels.push_back(systTypes[j]);
203  tempSig.push_back(sigs[i][systTypes[j]]);
204  tempBkg.push_back(bkgs[i][systTypes[j]]);
205  }
206  }
207  else{
208  for(const ISyst* syst: systs){
209  if(i == 0) labels.push_back(syst->LatexName());
210  tempSig.push_back(sigs[i][syst->LatexName()]);
211  tempBkg.push_back(bkgs[i][syst->LatexName()]);
212  }
213  }
214 
215  fSig.push_back(tempSig);
216  fBkg.push_back(tempBkg);
217  }
218 
219  new TCanvas;
220  myBarChart("Signal Uncertainty (%)", labels,
221  fSig, sigQuad, extrapLab);
222  drawLabel(beamLab);
223  preliminary();
224 
225  gPad->Print(("plots/syst_extrap_comp_"+beam+"_sig.pdf").c_str());
226 
227  new TCanvas;
228  myBarChart("Background Uncertainty (%)", labels,
229  fBkg, bkgQuad, extrapLab);
230  drawLabel(beamLab);
231  preliminary();
232 
233  gPad->Print(("plots/syst_extrap_comp_"+beam+"_bkg.pdf").c_str());
234 
235 }
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
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
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
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
std::string systTypes[kNumSystTypes]
const double j
Definition: BetheBloch.cxx:29
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
Neutrinos-only.
Definition: IPrediction.h:49
void myBarChart(std::string xlabel, std::vector< std::string > ylabels, std::vector< std::vector< std::pair< double, double >>> errs, std::vector< std::pair< double, double >> sum, std::vector< std::string > legend)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
double dmsq32
virtual void SetTh23(const T &th23)=0
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)
T min(const caf::Proxy< T > &a, T b)
double Integral(const Spectrum &s, double pot)
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
float fSig[xbins]
Definition: MakePlots.C:86
void drawLabel(std::string s, double x=0.3, double y=0.93)
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 40 of file systematics_extrap_comp_from_pred_interp.C.

Referenced by SumSysts(), SumSysts2017(), and systematics_extrap_comp_from_pred_interp().