makeSystTable.C
Go to the documentation of this file.
1 // This version marginalizes over a systematic (or group of them), gets the resultant 1D profile and gets the 1-sigma limits
2 // These limits are compared to the stats-only limits to see how much effect the systematic(s) have (asymmetric, look each direction)
3 // usage tips:
4 // change the evaluated at point by changing theta and msq
5 // I think it's right to still evaluate denom for fractional uncertainty at 0.514
6 // run this with a ana18 pred, send it to a log file. Then grep ++ <logfile> > <textfile> to get relevant lines
7 
8 #include "includes.h"
9 #include "Utilities/rootlogon.C"
10 
11 using namespace ana;
12 
13 
14 class SystEntry{
17  std::vector< const ISyst* > systVec;
18  const int xBins, yBins;
19  const double xLow, xHi, yLow, yHi;
20  const double thetaSeed, msqSeed;
21  TH1F *hSinUncertainty, *hDmsqUncertainty, *hSinDcpUncertainty, *hDmsqDcpUncertainty;
23  ofstream* fOut;
24 public:
25  // constructor
26  SystEntry(osc::IOscCalcAdjustable* kCal, MultiExperiment* kExpt, std::vector< const ISyst*> kSysts,
27  int kXBins, int kYBins, double kXLow, double kXHi, double kYLow, double kYHi, double kThetaSeed, double kMsqSeed,
28  std::string kSystName, ofstream* kFOut) :
29  calc(kCal), expt(kExpt), systVec(kSysts), xBins(kXBins), yBins(kYBins),
30  xLow(kXLow), xHi(kXHi), yLow(kYLow), yHi(kYHi), thetaSeed(kThetaSeed), msqSeed(kMsqSeed),
31  systName(kSystName), fOut(kFOut) {}
32 
33  // destructor
34  ~SystEntry () {delete calc; delete expt;
35  for(auto & sys:systVec) delete sys;
36  }
37 
38  void setUncertaintyHists();
39  void setLinearUncertaintyHists();
40  TH1F* getSin() { return hSinUncertainty;}
41  TH1F* getDm() { return hDmsqUncertainty;}
42  TH1F* getSinDcp() { return hSinDcpUncertainty;}
43  TH1F* getDmDcp() { return hDmsqDcpUncertainty;}
44  TH1F* hSin();
45  TH1F* hDm();
46  TH1F* hSinDcp();
47  TH1F* hDmDcp();
48 
49  std::map< std::string, double > GetLimits(TH1F* hist);
50 
51  // void printUncertainty();
52  // void printTotalUncertaintyLinearDcp(std::map< std::string, TH1F*> statMap);
53  // void printUncertaintyLinear(std::map< std::string, TH1F*> statMap);
54  // void printUncertaintySyst(std::map< std::string, TH1F*> statMap);
55  // void printUncertaintySystLinearDcp(std::map< std::string, TH1F*> statMap);
56  std::map< std::string, double > printUncertainty();
57  std::map< std::string, double > printTotalUncertaintyLinearDcp(std::map< std::string, TH1F*> statMap);
58  std::map< std::string, double > printUncertaintyLinear(std::map< std::string, TH1F*> statMap);
59  std::map< std::string, double > printUncertaintySyst(std::map< std::string, TH1F*> statMap);
60  std::map< std::string, double > printUncertaintySystLinearDcp(std::map< std::string, TH1F*> statMap);
61 
62 
63  void enterUncertaintyLine(double sin, double dmup, double dmdn, TString name);
64  void printTexFilePreamble();
65  void printTexFileHLine();
66 
67 };
68 
70  hSinUncertainty = (TH1F*)this -> hSin();
71  hDmsqUncertainty = (TH1F*)this -> hDm();
72 };
73 
75  hSinDcpUncertainty = (TH1F*)this -> hSinDcp();
76  hDmsqDcpUncertainty = (TH1F*)this -> hDmDcp();
77 };
78 
80  kFitSinSqTheta23.SetValue(calc,thetaSeed);
83  std::vector<const IFitVar*> oscparDM = {&kFitDmSq32Scaled};
84  TH1F *xplot = new TH1F("xplot","xplot",xBins,xLow,xHi);
85  xplot = (TH1F*)Profile(expt, calc, &kFitSinSqTheta23,
86  xBins, xLow, xHi, 0,
87  oscparDM, systVec);
88  return xplot;
89 }
90 
92  kFitSinSqTheta23.SetValue(calc,thetaSeed);
95  std::vector<const IFitVar*> oscparDM = {&kFitDmSq32Scaled, &kFitDeltaInPiUnits};
96  TH1F *xplot = new TH1F("xplot","xplot",xBins,xLow,xHi);
97  xplot = (TH1F*)Profile(expt, calc, &kFitSinSqTheta23,
98  xBins, xLow, xHi, 0,
99  oscparDM, {},
100  { {&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}, });
101  return xplot;
102 }
103 
105  kFitSinSqTheta23.SetValue(calc,thetaSeed);
106  kFitDmSq32Scaled.SetValue(calc,msqSeed);
108  std::vector<const IFitVar*> oscparSin = {&kFitSinSqTheta23};
109  TH1F *yplot = new TH1F("yplot","yplot",yBins,yLow,yHi);
110  yplot = (TH1F*)Profile(expt, calc, &kFitDmSq32Scaled,
111  yBins, yLow, yHi, 0,
112  oscparSin, systVec);
113  return yplot;
114 }
115 
117  kFitSinSqTheta23.SetValue(calc,thetaSeed);
118  kFitDmSq32Scaled.SetValue(calc,msqSeed);
120  std::vector<const IFitVar*> oscparSin = {&kFitSinSqTheta23, &kFitDeltaInPiUnits};
121  TH1F *yplot = new TH1F("yplot","yplot",yBins,yLow,yHi);
122  yplot = (TH1F*)Profile(expt, calc, &kFitDmSq32Scaled,
123  yBins, yLow, yHi, 0,
124  oscparSin, {},
125  { {&kFitSinSqTheta23, {0.45,0.55}},
126  {&kFitDeltaInPiUnits,{0.,0.5,1.,1.5}}, });
127  return yplot;
128 }
129 
130 
131 std::map< std::string, double > SystEntry::GetLimits(TH1F* hist){
132  double low = 0.0;
133  double hi = 0.0;
134  for(int bin = 1; bin<hist->GetNbinsX()+1; bin++) {
135  if(hist->GetBinContent(bin+1) < hist->GetBinContent(bin)) {
136  if(hist->GetBinContent(bin+1) < 1.0 && hist->GetBinContent(bin) > 1.0 && bin < hist->GetNbinsX()) {
137  float step = hist->GetBinContent(bin) - 1.0;
138  float inbin = hist->GetBinContent(bin) - hist->GetBinContent(bin+1);
139  float frac = step/inbin;
140  float temp = hist->GetXaxis()->GetXmin() + hist->GetXaxis()->GetBinWidth(1)/2.0 + (float(bin-1)+frac)*hist->GetXaxis()->GetBinWidth(1);
141  std::cout << "Crossed chisq of 1, descending, with x value of: " << temp << std::endl;
142  low = temp; // this is tricky with two dips. This should keep the min and max for the second dip which is fine
143  }
144  }
145  else {
146  if(hist->GetBinContent(bin+1) > 1.0 && hist->GetBinContent(bin) < 1.0 && bin < hist->GetNbinsX()) {
147  float step = 1.0 - hist->GetBinContent(bin);
148  float inbin = hist->GetBinContent(bin+1) - hist->GetBinContent(bin);
149  float frac = step/inbin;
150  hi = hist->GetXaxis()->GetXmin() + hist->GetXaxis()->GetBinWidth(1)/2.0 + (float(bin-1)+frac)*hist->GetXaxis()->GetBinWidth(1);
151  std::cout << "Crossed chisq of 1, ascending, with x value of: " << hi << std::endl;
152  }
153  }
154  }
155  std::map< std::string, double > map;
156  map["low"] = low;
157  map["hi"] = hi;
158  map["range"] = fabs(hi-low);
159  return map;
160 }
161 
162 
163 // print abs uncertainty
164 std::map< std::string, double > SystEntry::printUncertainty(){
165 //void SystEntry::printUncertainty(){
166 
167  TH1F* xplot = (TH1F*)hSinUncertainty -> Clone();
168  TH1F* yplot = (TH1F*)hDmsqUncertainty -> Clone();
169 
170  std::map< std::string, double > mapSin = this -> GetLimits(xplot);
171  std::map< std::string, double > mapDmsq = this -> GetLimits(yplot);
172  double stattheta = mapSin["range"];
173  double statmsqup = mapDmsq["hi"] - msqSeed;
174  double statmsqdn = msqSeed - mapDmsq["low"];
175 
176  double sinUnc = stattheta * pow(10,3) / 2;
177  double dmUncUp = statmsqup * pow(10,3);
178  double dmUncDn = statmsqdn * pow(10,3);
179 
180  this->enterUncertaintyLine(sinUnc, dmUncUp, dmUncDn, "Statistical");
181 
182  std::map< std::string, double > uncMap;
183  uncMap["sin"] = sinUnc;
184  uncMap["dmdn"] = dmUncDn;
185  uncMap["dmup"] = dmUncUp;
186  return uncMap;
187 }
188 
189 std::map< std::string, double > SystEntry::printTotalUncertaintyLinearDcp(std::map< std::string, TH1F*> statMap){
190  // void SystEntry::printTotalUncertaintyLinearDcp(std::map< std::string, TH1F*> statMap){
191 
192  TH1F* xplot = (TH1F*)hSinUncertainty -> Clone();
193  TH1F* yplot = (TH1F*)hDmsqUncertainty -> Clone();
194  std::map< std::string, double > systMapSin = this -> GetLimits(xplot);
195  std::map< std::string, double > systMapDmsq = this -> GetLimits(yplot);
196  double systtheta = systMapSin["range"];
197  double systmsqup = systMapDmsq["hi"] - msqSeed;
198  double systmsqdn = msqSeed - systMapDmsq["low"];
199 
200  TH1F* xplotDcp = (TH1F*)hSinDcpUncertainty -> Clone();
201  TH1F* yplotDcp = (TH1F*)hDmsqDcpUncertainty -> Clone();
202  std::map< std::string, double > dcpMapSin = this -> GetLimits(xplotDcp);
203  std::map< std::string, double > dcpMapDmsq = this -> GetLimits(yplotDcp);
204  double dcpTheta = dcpMapSin["range"];
205  double dcpDmsqUp = dcpMapDmsq["hi"] - msqSeed;
206  double dcpDmsqDn = msqSeed - dcpMapDmsq["low"];
207 
208  std::map< std::string, double > mapStatSin = GetLimits(statMap["sinStat"]);
209  std::map< std::string, double > mapStatDm = GetLimits(statMap["dmStat"]);
210  double stattheta = mapStatSin["range"];
211  double statmsqup = mapStatDm["hi"] - msqSeed;
212  double statmsqdn = msqSeed - mapStatDm["low"];
213 
214  double sinUncDcp = (dcpTheta - stattheta) * pow(10,3) / 2;
215  double dmUncUpDcp = (dcpDmsqUp - statmsqup) * pow(10,3);
216  double dmUncDnDcp = (dcpDmsqDn - statmsqdn) * pow(10,3);
217 
218  double sinUnc = (systtheta * pow(10,3) / 2) + sinUncDcp;
219  double dmUncUp = (systmsqup * pow(10,3)) + dmUncUpDcp;
220  double dmUncDn = (systmsqdn * pow(10,3)) + dmUncDnDcp;
221  this->enterUncertaintyLine(sinUnc, dmUncUp, dmUncDn, "Total uncer + linear $\\delta_{CP}$");
222 
223  std::map< std::string, double > uncMap;
224  uncMap["sin"] = sinUnc;
225  uncMap["dmdn"] = dmUncDn;
226  uncMap["dmup"] = dmUncUp;
227  return uncMap;
228 }
229 
230 std::map< std::string, double > SystEntry::printUncertaintySystLinearDcp(std::map< std::string, TH1F*> statMap){
231 
232  TH1F* xplot = (TH1F*)hSinUncertainty -> Clone();
233  TH1F* yplot = (TH1F*)hDmsqUncertainty -> Clone();
234 
235  std::map< std::string, double > mapSystSin = this -> GetLimits(xplot);
236  std::map< std::string, double > mapSystDmsq = this -> GetLimits(yplot);
237  double systTheta = mapSystSin["range"];
238  double systMsqUp = mapSystDmsq["hi"] - msqSeed;
239  double systMsqDn = msqSeed - mapSystDmsq["low"];
240 
241  std::map< std::string, double > mapStatSin = GetLimits(statMap["sinStat"]);
242  std::map< std::string, double > mapStatDm = GetLimits(statMap["dmStat"]);
243  double stattheta = mapStatSin["range"];
244  double statmsqup = mapStatDm["hi"] - msqSeed;
245  double statmsqdn = msqSeed - mapStatDm["low"];
246 
247  TH1F* xplotDcp = (TH1F*)hSinDcpUncertainty -> Clone();
248  TH1F* yplotDcp = (TH1F*)hDmsqDcpUncertainty -> Clone();
249  std::map< std::string, double > dcpMapSin = this -> GetLimits(xplotDcp);
250  std::map< std::string, double > dcpMapDmsq = this -> GetLimits(yplotDcp);
251  double dcpTheta = dcpMapSin["range"];
252  double dcpDmsqUp = dcpMapDmsq["hi"] - msqSeed;
253  double dcpDmsqDn = msqSeed - dcpMapDmsq["low"];
254 
255  double sinUncDcp = (dcpTheta - stattheta) * pow(10,3) / 2;
256  double dmUncUpDcp = (dcpDmsqUp - statmsqup) * pow(10,3);
257  double dmUncDnDcp = (dcpDmsqDn - statmsqdn) * pow(10,3);
258 
259  double sinUnc = (sqrt(systTheta*systTheta - stattheta*stattheta) * pow(10,3) / 2) + sinUncDcp;
260  double dmUncUp = (sqrt(systMsqUp*systMsqUp - statmsqup*statmsqup) * pow(10,3)) + dmUncUpDcp;
261  double dmUncDn = (sqrt(systMsqDn*systMsqDn - statmsqdn*statmsqdn) * pow(10,3)) + dmUncDnDcp;
262  this->enterUncertaintyLine(sinUnc, dmUncUp, dmUncDn, "Total syst. + linear $\\delta_{CP}$");
263 
264  sinUnc = (sqrt(systTheta*systTheta - stattheta*stattheta) * pow(10,3) / 2);
265  dmUncUp = (sqrt(systMsqUp*systMsqUp - statmsqup*statmsqup) * pow(10,3));
266  dmUncDn = (sqrt(systMsqDn*systMsqDn - statmsqdn*statmsqdn) * pow(10,3));
267  this->enterUncertaintyLine(sinUnc, dmUncUp, dmUncDn, "Testing... Total syst$");
268 
269  std::map< std::string, double > uncMap;
270  uncMap["sin"] = sinUnc;
271  uncMap["dmdn"] = dmUncDn;
272  uncMap["dmup"] = dmUncUp;
273  return uncMap;
274 }
275 
276 // print uncertainty relative to stats only
277 std::map< std::string, double > SystEntry::printUncertaintySyst(std::map< std::string, TH1F*> statMap){
278  //void SystEntry::printUncertaintySyst(std::map< std::string, TH1F*> statMap){
279 
280  TH1F* xplot = (TH1F*)hSinUncertainty -> Clone();
281  TH1F* yplot = (TH1F*)hDmsqUncertainty -> Clone();
282 
283  std::map< std::string, double > mapSystSin = this -> GetLimits(xplot);
284  std::map< std::string, double > mapSystDmsq = this -> GetLimits(yplot);
285  double systTheta = mapSystSin["range"];
286  double systMsqUp = mapSystDmsq["hi"] - msqSeed;
287  double systMsqDn = msqSeed - mapSystDmsq["low"];
288 
289  std::map< std::string, double > mapStatSin = GetLimits(statMap["sinStat"]);
290  std::map< std::string, double > mapStatDm = GetLimits(statMap["dmStat"]);
291  double stattheta = mapStatSin["range"];
292  double statmsqup = mapStatDm["hi"] - msqSeed;
293  double statmsqdn = msqSeed - mapStatDm["low"];
294 
295  double sinUnc = sqrt(systTheta*systTheta - stattheta*stattheta) * pow(10,3) / 2;
296  double dmUncUp = sqrt(systMsqUp*systMsqUp - statmsqup*statmsqup) * pow(10,3);
297  double dmUncDn = sqrt(systMsqDn*systMsqDn - statmsqdn*statmsqdn) * pow(10,3);
298  this->enterUncertaintyLine(sinUnc, dmUncUp, dmUncDn, systName);
299 
300  std::map< std::string, double > uncMap;
301  uncMap["sin"] = sinUnc;
302  uncMap["dmdn"] = dmUncDn;
303  uncMap["dmup"] = dmUncUp;
304  return uncMap;
305 }
306 
307 std::map< std::string, double > SystEntry::printUncertaintyLinear(std::map< std::string, TH1F*> statMap){
308 //void SystEntry::printUncertaintyLinear(std::map< std::string, TH1F*> statMap){
309 
310  TH1F* xplot = (TH1F*)hSinDcpUncertainty -> Clone();
311  TH1F* yplot = (TH1F*)hDmsqDcpUncertainty -> Clone();
312 
313  std::map< std::string, double > mapDcpSin = this -> GetLimits(xplot);
314  std::map< std::string, double > mapDcpDmsq = this -> GetLimits(yplot);
315  double dcpTheta = mapDcpSin["range"];
316  double dcpDmsqUp = mapDcpDmsq["hi"] - msqSeed;
317  double dcpDmsqDn = msqSeed - mapDcpDmsq["low"];
318 
319  std::map< std::string, double > mapStatSin = GetLimits(statMap["sinStat"]);
320  std::map< std::string, double > mapStatDm = GetLimits(statMap["dmStat"]);
321  double stattheta = mapStatSin["range"];
322  double statmsqup = mapStatDm["hi"] - msqSeed;
323  double statmsqdn = msqSeed - mapStatDm["low"];
324 
325  double sinUnc = (dcpTheta - stattheta) * pow(10,3) / 2;
326  double dmUncUp = (dcpDmsqUp - statmsqup) * pow(10,3);
327  double dmUncDn = (dcpDmsqDn - statmsqdn) * pow(10,3);
328  this->enterUncertaintyLine(sinUnc, dmUncUp, dmUncDn, "$\\Delta_{CP} (0-2\\pi$)");
329 
330  std::map< std::string, double > uncMap;
331  uncMap["sin"] = sinUnc;
332  uncMap["dmdn"] = dmUncDn;
333  uncMap["dmup"] = dmUncUp;
334  return uncMap;
335 }
336 
337 
338 void SystEntry::enterUncertaintyLine(double sin, double dmup, double dmdn, TString name){
339  TString stOut = name + Form(" & $\\pm$%.2g & +%.2g / -%.2g \\\\ \n", sin, dmup, dmdn);
340  std::string stringOut = (std::string) stOut;
341  std::cout<< "stringOut: " << stringOut <<std::endl;
342  *(fOut) << stringOut;
343 }
344 
346  std::string preamble = "\\documentclass{article} \n\\usepackage[english]{babel} \n\\usepackage{amsmath} \n\\usepackage{multirow} \n";
347  preamble += "\\usepackage[dvipsnames]{xcolor} \n\n\\begin{document} \n\n\\begin{table*}[t] \n\\centering \n";
348  preamble += "\\begin{tabular}{c c c} \n\\hline \n\\multirow{2}{*}{Source of uncertainty} & Uncertainty in & Uncertainty in \\\\";
349  preamble += "& $\\sin^2\\!\\theta_{23} (\\times 10^{-3})$ & $\\Delta m^2_{32} \\left(10^{-6}\\text{eV}^{2}\\right)$\\\\ \n\\hline \n";
350  *(fOut) << preamble;
351 }
352 
354  std::string hline = "\\hline \n";
355  *(fOut) << hline;
356 }
357 
358 void DrawHCLabel(TString hcName = "fhc")
359 {
360  TLatex* label = new TLatex(.35, .95, hcName);
361  label->SetTextColor(kBlue);
362  label->SetNDC();
363  label->SetTextSize(2/30.);
364  label->SetTextAlign(32);
365  label->Draw();
366 }
367 
368 void makeSystTable(TString current = "both", bool makeFullSystTable = false)
369 {
370 
371  std::string OUTDIR = "systTables/";
372 
373  // set current to "fhc" or "rhc" to make table for only one horn current
374  osc::IOscCalcAdjustable* calc = DefaultOscCalc(); // full systs normal h
375  osc::IOscCalcAdjustable* calcfake = DefaultOscCalc(); // full systs normal h
376 
377  //Absolute muon energy scale (`2%)
378  std::vector<const ISyst*> absMuESyst;
379  absMuESyst.push_back(&kMuEScaleSyst2017);
380  // Relative muon energy scale (`2%)
381  std::vector<const ISyst*> relMuESyst;
382  relMuESyst.push_back(&kRelMuEScaleSyst2017);
383  // Absolute hadronic energy scale (`5%)
384  std::vector<const ISyst*> absHadESyst;
385  absHadESyst.push_back(&kAnaCalibrationSyst);
386  absHadESyst.push_back(&kAnaCalibShapeSyst);
387  // Relative hadronic energy scale (`5%)
388  std::vector<const ISyst*> relHadESyst;
389  relHadESyst.push_back(&kAnaRelativeCalibSyst);
390  // Normalization (`5%)
391  std::vector< const ISyst* > normSysts;
392  normSysts.push_back(&kAna2018NormFHC);
393  normSysts.push_back(&kAna2018NormRHC);
394  std::vector<const ISyst*> xsecSysts;
395  AddJointAna2018XSecSysts(xsecSysts, kNumuAna2018, 1);
396  // Neutrino flux
397  std::vector<const ISyst*> fluxSysts;
399  // Neutron syst
400  std::vector< const ISyst* > neutronSyst;
401  neutronSyst.push_back(&kNeutronVisEScalePrimariesSyst2018);
402  // Scintillation model
403  std::vector<const ISyst*> scintSysts;
404  scintSysts.push_back(&kAnaLightlevelSyst);
405  scintSysts.push_back(&kAnaCherenkovSyst);
406  // Total systematic uncertainty
407 
409  // systs = absMuESyst; // for quick testing
410  std::vector<const ISyst*> totSysts = getAllAna2018Systs(kNumuAna2018,false);
411  std::vector<const ISyst*> noSysts = {};
412 
413  std::map< std::string, std::vector<const ISyst*> > mapSystVecs;
414  mapSystVecs["Absolute muon energy scale"] = absMuESyst;
415  mapSystVecs["Relative muon energy scale"] = relMuESyst;
416  mapSystVecs["Absolute hadronic energy scale"] = absHadESyst;
417  mapSystVecs["Relative hadronic energy scale"] = relHadESyst;
418  mapSystVecs["Normalisation"] = normSysts;
419  mapSystVecs["Neutron uncertainty"] = neutronSyst;
420  mapSystVecs["Cross sections and final-state interaction"] = xsecSysts;
421  mapSystVecs["Neutrino flux"] = fluxSysts;
422  mapSystVecs["Detector response"] = scintSysts;
423 
424  float theta_real = 0.58;
425  float msq = 2.51;
426  std::cout<< "\n++theta_real: " << theta_real
427  << "\n++msq: " << msq
428  <<std::endl;
429  kFitSinSqTheta23.SetValue(calc,theta_real);
430  kFitDmSq32Scaled.SetValue(calc,msq);
431  kFitDeltaInPiUnits.SetValue(calc,1.5);
432 
433  const std::string fakeNDPredFileNameFHC="/nova/ana/nu_mu_ana/Ana2018/Predictions/pred_fakeNDData_numuconcat_fhc__numu2018.root";
434  const std::string fakeNDPredFileNameRHC="/nova/ana/nu_mu_ana/Ana2018/Predictions/pred_fakeNDData_numuconcat_rhc__numu2018.root";
435  std::cout<< "FHC pred file: " << fakeNDPredFileNameFHC
436  << "\nRHC pred file: " << fakeNDPredFileNameRHC
437  <<std::endl;
438 
439 
440  ENu2018ExtrapType extrap = kNuMu;
441  PredictionSystJoint2018 prediction_Q1_FHC(extrap, calc, "fhc", 1, systs, fakeNDPredFileNameFHC);
442  PredictionSystJoint2018 prediction_Q2_FHC(extrap, calc, "fhc", 2, systs, fakeNDPredFileNameFHC);
443  PredictionSystJoint2018 prediction_Q3_FHC(extrap, calc, "fhc", 3, systs, fakeNDPredFileNameFHC);
444  PredictionSystJoint2018 prediction_Q4_FHC(extrap, calc, "fhc", 4, systs, fakeNDPredFileNameFHC);
445  PredictionSystJoint2018 prediction_Q1_RHC(extrap, calc, "rhc", 1, systs, fakeNDPredFileNameRHC);
446  PredictionSystJoint2018 prediction_Q2_RHC(extrap, calc, "rhc", 2, systs, fakeNDPredFileNameRHC);
447  PredictionSystJoint2018 prediction_Q3_RHC(extrap, calc, "rhc", 3, systs, fakeNDPredFileNameRHC);
448  PredictionSystJoint2018 prediction_Q4_RHC(extrap, calc, "rhc", 4, systs, fakeNDPredFileNameRHC);
449 
450  // make fake data:
451  Spectrum fakeQ1_FHC = prediction_Q1_FHC.Predict(calc).FakeData(kAna2018FHCPOT);
452  Spectrum fakeQ2_FHC = prediction_Q2_FHC.Predict(calc).FakeData(kAna2018FHCPOT);
453  Spectrum fakeQ3_FHC = prediction_Q3_FHC.Predict(calc).FakeData(kAna2018FHCPOT);
454  Spectrum fakeQ4_FHC = prediction_Q4_FHC.Predict(calc).FakeData(kAna2018FHCPOT);
455  Spectrum fakeQ1_RHC = prediction_Q1_RHC.Predict(calc).FakeData(kAna2018RHCPOT);
456  Spectrum fakeQ2_RHC = prediction_Q2_RHC.Predict(calc).FakeData(kAna2018RHCPOT);
457  Spectrum fakeQ3_RHC = prediction_Q3_RHC.Predict(calc).FakeData(kAna2018RHCPOT);
458  Spectrum fakeQ4_RHC = prediction_Q4_RHC.Predict(calc).FakeData(kAna2018RHCPOT);
459 
460  // get cosmics:
461  // Get the cosmic distributions:
462  DontAddDirectory guard;
463  TFile* fCosFHC = new TFile("Files/cosmics_fhc__numu2018.root");
464  TFile* fCosRHC = new TFile("Files/cosmics_rhc__numu2018.root");
465  TH1D* hCosFHC_Q1 = (TH1D*) fCosFHC->Get("cosmics_q1");
466  TH1D* hCosFHC_Q2 = (TH1D*) fCosFHC->Get("cosmics_q2");
467  TH1D* hCosFHC_Q3 = (TH1D*) fCosFHC->Get("cosmics_q3");
468  TH1D* hCosFHC_Q4 = (TH1D*) fCosFHC->Get("cosmics_q4");
469  TH1D* hCosRHC_Q1 = (TH1D*) fCosRHC->Get("cosmics_q1");
470  TH1D* hCosRHC_Q2 = (TH1D*) fCosRHC->Get("cosmics_q2");
471  TH1D* hCosRHC_Q3 = (TH1D*) fCosRHC->Get("cosmics_q3");
472  TH1D* hCosRHC_Q4 = (TH1D*) fCosRHC->Get("cosmics_q4");
473  fCosFHC->Close();
474  fCosRHC->Close();
475  delete fCosFHC;
476  delete fCosRHC;
477  Spectrum specCosFHC_Q1(hCosFHC_Q1, kAna2018FHCPOT, kAna2018FHCLivetime);
478  Spectrum specCosFHC_Q2(hCosFHC_Q2, kAna2018FHCPOT, kAna2018FHCLivetime);
479  Spectrum specCosFHC_Q3(hCosFHC_Q3, kAna2018FHCPOT, kAna2018FHCLivetime);
480  Spectrum specCosFHC_Q4(hCosFHC_Q4, kAna2018FHCPOT, kAna2018FHCLivetime);
481  Spectrum specCosRHC_Q1(hCosRHC_Q1, kAna2018RHCPOT, kAna2018RHCLivetime);
482  Spectrum specCosRHC_Q2(hCosRHC_Q2, kAna2018RHCPOT, kAna2018RHCLivetime);
483  Spectrum specCosRHC_Q3(hCosRHC_Q3, kAna2018RHCPOT, kAna2018RHCLivetime);
484  Spectrum specCosRHC_Q4(hCosRHC_Q4, kAna2018RHCPOT, kAna2018RHCLivetime);
485 
486  // Add cosmic bkg to fake FD data
487  fakeQ1_FHC += specCosFHC_Q1;
488  fakeQ2_FHC += specCosFHC_Q2;
489  fakeQ3_FHC += specCosFHC_Q3;
490  fakeQ4_FHC += specCosFHC_Q4;
491  fakeQ1_RHC += specCosRHC_Q1;
492  fakeQ2_RHC += specCosRHC_Q2;
493  fakeQ3_RHC += specCosRHC_Q3;
494  fakeQ4_RHC += specCosRHC_Q4;
495 
496  // // make expt for each quantile:
497  SingleSampleExperiment* exptQ1_FHC = new SingleSampleExperiment(&prediction_Q1_FHC, fakeQ1_FHC, specCosFHC_Q1);
498  SingleSampleExperiment* exptQ2_FHC = new SingleSampleExperiment(&prediction_Q2_FHC, fakeQ2_FHC, specCosFHC_Q2);
499  SingleSampleExperiment* exptQ3_FHC = new SingleSampleExperiment(&prediction_Q3_FHC, fakeQ3_FHC, specCosFHC_Q3);
500  SingleSampleExperiment* exptQ4_FHC = new SingleSampleExperiment(&prediction_Q4_FHC, fakeQ4_FHC, specCosFHC_Q4);
501  SingleSampleExperiment* exptQ1_RHC = new SingleSampleExperiment(&prediction_Q1_RHC, fakeQ1_RHC, specCosRHC_Q1);
502  SingleSampleExperiment* exptQ2_RHC = new SingleSampleExperiment(&prediction_Q2_RHC, fakeQ2_RHC, specCosRHC_Q2);
503  SingleSampleExperiment* exptQ3_RHC = new SingleSampleExperiment(&prediction_Q3_RHC, fakeQ3_RHC, specCosRHC_Q3);
504  SingleSampleExperiment* exptQ4_RHC = new SingleSampleExperiment(&prediction_Q4_RHC, fakeQ4_RHC, specCosRHC_Q4);
505 
506  std::vector<const IExperiment*> exptsAllFake;
507  if(current=="fhc" || current=="both"){
508  exptsAllFake.push_back(exptQ1_FHC);
509  exptsAllFake.push_back(exptQ2_FHC);
510  exptsAllFake.push_back(exptQ3_FHC);
511  exptsAllFake.push_back(exptQ4_FHC);
512  }
513  if(current=="rhc" || current=="both"){
514  exptsAllFake.push_back(exptQ1_RHC);
515  exptsAllFake.push_back(exptQ2_RHC);
516  exptsAllFake.push_back(exptQ3_RHC);
517  exptsAllFake.push_back(exptQ4_RHC);
518  }
519  exptsAllFake.push_back(&kSolarConstraintsPDG2017);
520  exptsAllFake.push_back(WorldReactorConstraint2017());
521  MultiExperiment exptFakeCC(exptsAllFake);
522 
523  std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam();
524  std::vector<const IFitVar*> oscparDM = {&kFitDmSq32Scaled};
525  std::vector<const IFitVar*> oscparSin = {&kFitSinSqTheta23};
526 
527  int kXBins = 100;
528  int kYBins = 100;
529  // int kXBins = 4;
530  // int kYBins = 4;
531  double xLowEdge = 0.3;
532  double xHiEdge = 0.7;
533  double yLowEdge = 2.1;
534  double yHiEdge = 2.8;
535 
536  double range = 0;
537  float stattheta = 0, statmsq = 0;
538  float hi, low;
539  float hix, lowx;
540 
541  TString texFileName = "systTable_" + current;
542  if(makeFullSystTable) texFileName += (TString)"_fullTable.tex";
543  if(!makeFullSystTable) texFileName += (TString)"_summaryTable.tex";
544  ofstream *fTexOut = new ofstream(OUTDIR + (std::string)texFileName);
545 
546  // std::cout<< "Get the stats only uncertainty" <<std::endl;
547  SystEntry* entStat = new SystEntry(calcfake, &exptFakeCC, {}, kXBins, kYBins, xLowEdge, xHiEdge, yLowEdge, yHiEdge, theta_real, msq, "Stats. only", fTexOut);
548  entStat->setUncertaintyHists();
549  entStat->printTexFilePreamble();
550 
551  std::map<std::string, std::pair<double, double>> ss23SystErrors;
552  std::map<std::string, std::pair<double, double>> dm32SystErrors;
553 
554  std::map< std::string, TH1F* > mapStatHists;
555  mapStatHists["sinStat"] = (TH1F*) entStat->getSin();
556  mapStatHists["dmStat"] = (TH1F*) entStat->getDm();
557 
558  if(!makeFullSystTable){
559  for(auto const& thisSystVec: mapSystVecs){
560  SystEntry* entSyst = new SystEntry(calcfake, &exptFakeCC, thisSystVec.second, kXBins, kYBins, xLowEdge, xHiEdge, yLowEdge, yHiEdge,
561  theta_real, msq, (std::string)thisSystVec.first, fTexOut);
562  entSyst->setUncertaintyHists();
563 
564  std::map< std::string, double > mapUnc = entSyst->printUncertaintySyst(mapStatHists);
565  ss23SystErrors[(std::string)thisSystVec.first] = {mapUnc["sin"], mapUnc["sin"]};
566  dm32SystErrors[(std::string)thisSystVec.first] = {mapUnc["dmdn"], mapUnc["dmup"]};
567 
568  //delete entSyst;
569  }
570  }
571 
572 
573  entStat->setLinearUncertaintyHists();
574  std::map< std::string, double > mapDcpUnc = entStat->printUncertaintyLinear(mapStatHists);
575  ss23SystErrors["Value of #delta_{CP}"] = {mapDcpUnc["sin"], mapDcpUnc["sin"]};
576  dm32SystErrors["Value of #delta_{CP}"] = {mapDcpUnc["dmdn"], mapDcpUnc["dmup"]};
577 
578  entStat->printTexFileHLine();
579  std::map< std::string, double > mapStatUnc = entStat->printUncertainty();
580  std::pair<double, double> ss23StatError {mapStatUnc["sin"], mapStatUnc["sin"]};
581  std::pair<double, double> dm32StatError {mapStatUnc["dmdn"], mapStatUnc["dmup"]};
582 
583  //delete entStat;
584 
585  TString plotFileName = "systTable_" + current;
586  if(makeFullSystTable) plotFileName += (TString)"_full_plot";
587  if(!makeFullSystTable) plotFileName += (TString)"_summary_plot";
588 
589  TString plotCurrentLabel = "";
590  if(current == "fhc") plotCurrentLabel = "#nu beam";
591  if(current == "rhc") plotCurrentLabel = "#bar{#nu} beam";
592  if(current == "both") plotCurrentLabel = "#nu + #bar{#nu} beam";
593 
594  TCanvas c;
595  c.SetBottomMargin(0.15); // otherwise axis titles get chopped off
596  ana::ErrorBarChart(ss23SystErrors, ss23StatError, "Uncertainty on sin^{2}#theta_{23} (10^{-3})");
597  c.SetLeftMargin(0.5); // otherwise axis titles get chopped off
598  Preliminary();
599  DrawHCLabel(plotCurrentLabel);
600  c.SaveAs( (OUTDIR).c_str() + (TString)"sin_" + plotFileName + ".pdf");
601  c.SaveAs( (OUTDIR).c_str() + (TString)"sin_" + plotFileName + ".root");
602  c.Clear();
603  ana::ErrorBarChart(dm32SystErrors, dm32StatError, "Uncertainty on #Deltam_{32}^{2} (10^{-6} eV^{2})");
604  c.SetLeftMargin(0.5); // otherwise axis titles get chopped off
605  Preliminary();
606  DrawHCLabel(plotCurrentLabel);
607  c.SaveAs( (OUTDIR).c_str() + (TString)"dm_" + plotFileName + ".pdf");
608  c.SaveAs( (OUTDIR).c_str() + (TString)"dm_" + plotFileName + ".root");
609 
610 
611 
612  SystEntry* entTotSyst = new SystEntry(calcfake, &exptFakeCC, systs, kXBins, kYBins, xLowEdge, xHiEdge, yLowEdge, yHiEdge,
613  theta_real, msq, "Total Systematic", fTexOut);
614  entTotSyst->setUncertaintyHists();
615  entTotSyst->setLinearUncertaintyHists();
616  std::map< std::string, double > mapSystWithDcpUnc = entTotSyst->printUncertaintySystLinearDcp(mapStatHists);
617  std::map< std::string, double > mapTotalWithDcpUnc = entTotSyst->printTotalUncertaintyLinearDcp(mapStatHists);
618 
619  //delete entTotSyst;
620 
621  std::string postamble = "\\hline \n\\end{tabular} \n\\end{table*} \n\n\\end{document} \n";
622  *(fTexOut) << postamble;
623  fTexOut -> close();
624  delete fTexOut;
625 
626 
627 
628 }
const XML_Char * name
Definition: expat.h:151
std::map< std::string, double > printUncertaintySyst(std::map< std::string, TH1F * > statMap)
TH1F * hDmDcp()
std::map< std::string, double > printUncertaintyLinear(std::map< std::string, TH1F * > statMap)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
TH1F * hDm()
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void setLinearUncertaintyHists()
Definition: makeSystTable.C:74
void enterUncertaintyLine(double sin, double dmup, double dmdn, TString name)
const int yBins
Definition: makeSystTable.C:18
std::map< std::string, double > printUncertaintySystLinearDcp(std::map< std::string, TH1F * > statMap)
Loads shifted spectra from files.
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
TH1F * hSin()
Definition: makeSystTable.C:79
float GetLimits(TH1F *plot, float &hi, float &low)
MultiExperiment * expt
Definition: makeSystTable.C:16
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void AddJointAna2018BeamSysts(std::vector< const ISyst * > &systs, const EAnaType2018 ana)
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
std::map< std::string, double > GetLimits(TH1F *hist)
const DummyAnaSyst kAnaRelativeCalibSyst("RelativeCalib","RelCalib")
osc::IOscCalcAdjustable * calc
Definition: makeSystTable.C:15
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const double kAna2018RHCPOT
Definition: Exposures.h:208
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
const char * label
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
const double thetaSeed
Definition: makeSystTable.C:20
TH1F * getDmDcp()
Definition: makeSystTable.C:43
void makeSystTable(TString current="both", bool makeFullSystTable=false)
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
TSpline3 hi("hi", xhi, yhi, 18,"0")
expt
Definition: demo5.py:34
TH1F * getDm()
Definition: makeSystTable.C:41
TH1F * getSin()
Definition: makeSystTable.C:40
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
std::string systName
Definition: makeSystTable.C:22
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
TH1F * hSinDcp()
Definition: makeSystTable.C:91
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void printTexFilePreamble()
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
std::vector< const ISyst * > systVec
Definition: makeSystTable.C:17
double frac(double x)
Fractional part.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
float bin[41]
Definition: plottest35.C:14
void Preliminary()
OStream cout
Definition: OStream.cxx:6
TH1F * hSinUncertainty
Definition: makeSystTable.C:21
void DrawHCLabel(TString hcName="fhc")
std::map< std::string, double > printUncertainty()
Spectrum Predict(osc::IOscCalc *calc) const override
TH1F * getSinDcp()
Definition: makeSystTable.C:42
void setUncertaintyHists()
Definition: makeSystTable.C:69
Combine multiple component experiments.
SystEntry(osc::IOscCalcAdjustable *kCal, MultiExperiment *kExpt, std::vector< const ISyst * > kSysts, int kXBins, int kYBins, double kXLow, double kXHi, double kYLow, double kYHi, double kThetaSeed, double kMsqSeed, std::string kSystName, ofstream *kFOut)
Definition: makeSystTable.C:26
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
T sin(T number)
Definition: d0nt_math.hpp:132
const DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
General interface to any calculator that lets you set the parameters.
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:798
void AddJointAna2018XSecSysts(std::vector< const ISyst * > &systs, const EAnaType2018 ana, bool smallgenie)
const double yLow
Definition: makeSystTable.C:19
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: Utilities.h:62
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::map< std::string, double > printTotalUncertaintyLinearDcp(std::map< std::string, TH1F * > statMap)
ofstream * fOut
Definition: makeSystTable.C:23
const DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
const double kAna2018FHCLivetime
Definition: Exposures.h:213
procfile close()
const DummyAnaSyst kAnaCalibShapeSyst("CalibShape","CalibShape")
const double kAna2018RHCLivetime
Definition: Exposures.h:214
Compare a single data spectrum to the MC + cosmics expectation.
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
void printTexFileHLine()