systematics_table_from_pred_interp.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Macro for making the nue and numu systematics bar plots
4 // Use quant = 0 for nue and quant = 1,2,3,4 for the 4 numu quantiles
5 //
6 ////////////////////////////////////////////////////////////////////////////////
7 
10 #include "CAFAna/Analysis/Plots.h"
11 using namespace ana;
12 
15 
18 
19 #include "TCanvas.h"
20 
21 #include <iomanip>
22 
23 struct ColumnDef
24 {
26  Flavors::Flavors_t flav;
29 };
30 
31 double Integral(const Spectrum& s, double pot)
32 {
33  return s.Integral(pot);
34 }
35 
37  bool shortchart = true,
38  std::string decomp = "combo")
39 {
40  std::vector<ColumnDef> cols;
41  if(quant == 0) cols =
42  {{"$\\nu_e$ Signal", Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth},
43  {"Total Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
44  {"beam $\\nu_e$ CC", Flavors::kNuEToNuE, Current::kCC, Sign::kBoth},
45  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
46  {"$\\nu_\\mu$ CC", Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth},
47  {"$\\nu_\\tau$ CC", Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth}};
48  else cols =
49  {{"$\\nu_\\mu$ Signal",Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth},
50  {"Total Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
51  {"$\\nu_e$ CC", Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth},
52  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
53  {"$\\nu_\\mu$ App", Flavors::kNuEToNuMu, Current::kCC, Sign::kBoth},
54  {"$\\nu_\\tau$ CC", Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth}};
55 
56  std::cout << std::setprecision(2);
57 
59  if(decomp == "noextrap") extrap = kNoExtrap;
60  if(decomp == "combo") extrap = kCombo;
61 
63 
64  osc::OscCalculatorDumb dumbosc;
65 
66  const std::vector<const ISyst*> systs = quant==0 ? GetJointFitSystematicList(true, true, false, true) :
67  GetJointFitSystematicList(true, false, true, true) ;
68 
69  const IPrediction *pred = quant==0 ? GetNuePrediction2017(decomp,&dumbosc,true) :
70  GetNumuPredictions2017()[quant-1];
71 
72  std::vector<Spectrum> noms;
73  for(ColumnDef col: cols) noms.push_back(pred->PredictComponent(&dumbosc, col.flav, col.curr, col.sign));
74 
75  noms[1] -= noms[0];
76 
77  const double sigErr = 100/sqrt(Integral(noms[0], kAna2017POT));
78  const double bkgErr = 100/sqrt(Integral(noms[1], kAna2017POT));
79 
80  std::cout << "Systematic";
81  for(ColumnDef col: cols) std::cout << " & " << col.name;
82  std::cout << " \\\\" << std::endl;
83  std::cout << "\\hline" << std::endl;
84 
85  std::vector<double> quads(cols.size(), 0);
86 
87  std::map<std::string, double> barsSig, barsBkg;
88 
89  std::vector<double> errs;
90 
91  for(const ISyst* syst: systs){
92  std::vector<Spectrum> shiftsUp, shiftsDn;
93  for(ColumnDef col: cols){
94  shiftsUp.push_back(pred->PredictComponentSyst(&dumbosc,
95  SystShifts(syst, +1),
96  col.flav, col.curr, col.sign));
97  shiftsDn.push_back(pred->PredictComponentSyst(&dumbosc,
98  SystShifts(syst, -1),
99  col.flav, col.curr, col.sign));
100  }
101 
102  shiftsUp[1] -= shiftsUp[0];
103  shiftsDn[1] -= shiftsDn[0];
104 
105  std::cout << syst->LatexName();
106 
107  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
108  const Spectrum& nom = noms[colIdx];
109 
110  // POT cancels, is irrelevant
111  const double ratioUp = Integral(shiftsUp[colIdx], 1e20) / Integral(nom, 1e20);
112  const double errUp = 100*fabs(ratioUp-1);
113 
114  const double ratioDn = Integral(shiftsDn[colIdx], 1e20) / Integral(nom, 1e20);
115  const double errDn = 100*fabs(ratioDn-1);
116 
117  const double err = .5*(errUp+errDn);
118 
119  quads[colIdx] += err*err;
120 
121  std::cout << " & " << err;
122 
123  if(cols[colIdx].name == "$\\nu_e$ Signal" || cols[colIdx].name == "$\\nu_\\mu$ Signal") barsSig[syst->LatexName()] = err;
124  if(cols[colIdx].name == "Total Beam Bkg") barsBkg[syst->LatexName()] = err;
125  } // end for colIdx
126 
127  std::cout << " \\\\" << std::endl;
128  } // end for syst
129 
130  std::cout << "\\hline" << std::endl;
131 
132  std::cout << "Quadrature";
133  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) std::cout << " & " << sqrt(quads[colIdx]);
134  std::cout << std::endl;
135 
136  if(shortchart){
137  std::cout<<"Reducing Systematic List..."<<std::endl;
138  barsSig=SumSysts2017(barsSig);
139  barsBkg=SumSysts2017(barsBkg);
140  }
141 
142  std::string shortStr = shortchart ? "short" : "full";
143  std::string quantStr = quant==0 ? "_nue_" : "_numu_q"+std::to_string(quant)+"_";
144  new TCanvas;
145  CountingExperimentErrorBarChart(barsSig, sigErr, true);
146  gPad->Print(("plots/bars_sig_"+decomp+quantStr+shortStr+".root").c_str());
147  gPad->Print(("plots/bars_sig_"+decomp+quantStr+shortStr+".pdf").c_str());
148  new TCanvas;
149  CountingExperimentErrorBarChart(barsBkg, bkgErr, false);
150  gPad->Print(("plots/bars_bkg_"+decomp+quantStr+shortStr+".root").c_str());
151  gPad->Print(("plots/bars_bkg_"+decomp+quantStr+shortStr+".pdf").c_str());
152 }
const XML_Char * name
Definition: expat.h:151
Oscillation analysis framework, runs over CAF files outside of ART.
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
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool intersection=true)
std::map< std::string, double > SumSysts2017(std::map< std::string, double > systbars)
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 double kAna2017POT
Definition: Exposures.h:174
double Integral(const Spectrum &s, double pot)
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
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
Charged-current interactions.
Definition: IPrediction.h:39
Interactions of both types.
Definition: IPrediction.h:42
std::vector< const IPrediction * > GetNumuPredictions2017(const int nq=4, bool useSysts=true, bool GetFromUPS=false)
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Int_t col[ntarg]
Definition: Style.C:29
#define pot
NuESystPID
Definition: NueSysts.h:12
void systematics_table_from_pred_interp(int quant=0, bool shortchart=true, std::string decomp="combo")
def sign(x)
Definition: canMan.py:204
OStream cout
Definition: OStream.cxx:6
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalculator *calc, bool corrSysts, bool GetFromUPS=false)
(&#39; survival&#39;)
Definition: IPrediction.h:19
void CountingExperimentErrorBarChart(const std::map< std::string, double > &systbars, double statErr, bool bkgdOrSig, bool shortchart)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:828
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
string syst
Definition: plotSysts.py:176
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const std::vector< ColumnDef > cols
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16