systematics_extrap_comp_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 void myBarChart(std::string xlabel,
24  std::vector<std::string> ylabels,
25  std::vector<std::vector<std::pair<double,double>>> errs,
26  std::vector<std::pair<double,double>> sum,
27  std::vector<std::string> legend);
28 
38  };
39 
41  {"Near-Far Differences",
42  "Neutrino Cross Sections",
43  "Detector Calibration",
44  "Beam Flux",
45  "Detector Response",
46  "Neutron Uncertainty",
47  "Muon Energy Scale",
48  "Normalization"};
49 
50  std::map<std::string, std::pair<double,double>> SumSysts(std::map<std::string, std::pair<double,double>> bars);
51 
52 struct ColumnDef
53 {
55  Flavors::Flavors_t flav;
58 };
59 
60 void drawLabel(std::string s,double x=0.3,double y=0.93){
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 }
67 
68 double Integral(const Spectrum& s, double pot)
69 {
70  return s.Integral(pot);
71 }
72 
73  void preliminary()
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  }
82 
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 }
236 
237 
239  std::vector<std::string> ylabels,
240  std::vector<std::vector<std::pair<double,double>>> errs,
241  std::vector<std::pair<double,double>> sum,
242  std::vector<std::string> legend)
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 }
302 
303  std::map<std::string, std::pair<double,double>> SumSysts(std::map<std::string, std::pair<double,double>> bars)
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  }
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
void systematics_extrap_comp_from_pred_interp(std::string beam="fhc", bool sum=true)
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
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
sigs
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
const int nbins
Definition: cellShifts.C:15
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
Definition: IPrediction.h:39
int colors[6]
Definition: tools.h:1
General interface to any calculator that lets you set the parameters.
Interactions of both types.
Definition: IPrediction.h:42
float fSig[xbins]
Definition: MakePlots.C:86
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)
#define pot
const int cols[3]
const double j
Definition: BetheBloch.cxx:29
virtual void SetTh23(const T &th23)=0
virtual void SetDmsq32(const T &dmsq32)=0
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
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
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
::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
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
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
void drawLabel(std::string s, double x=0.3, double y=0.93)
T asin(T number)
Definition: d0nt_math.hpp:60