makeSystTable_reduced.C
Go to the documentation of this file.
1 // This version marginalizes over a systematic (or group of them), gets the resultant contour, then projects to 1D 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 SA pred, send it to a log file. Then grep ++ <logfile> > <textfile> to get relevant lines
7 
8 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
11 
12 #include "CAFAna/Core/Loaders.h"
13 #include "CAFAna/Core/Binning.h"
14 //JPD for dontadddirectory guard
15 #include "CAFAna/Core/Utilities.h"
16 
19 #include "CAFAna/Prediction/PredictionNumuFAHadE.h"
23 
24 #include "CAFAna/Systs/Systs.h"
26 #include "CAFAna/Systs/BeamSysts.h"
30 
32 
34 #include "CAFAna/Analysis/Calcs.h"
35 #include "CAFAna/Fit/Fit.h"
37 #include "CAFAna/Analysis/Plots.h"
38 
43 
44 #include "CAFAna/Vars/FitVars.h"
46 
47 #include "CAFAna/FC/FCCollection.h"
48 #include "CAFAna/FC/FCPoint.h"
49 #include "CAFAna/FC/FCSurface.h"
50 
51 #include "OscLib/OscCalcGeneral.h"
52 #include "OscLib/OscCalcPMNS.h"
53 #include "OscLib/OscCalc.h"
54 #include "OscLib/IOscCalc.h"
55 
56 
57 #include "TCanvas.h"
58 #include "TFile.h"
59 #include "TGraph.h"
60 #include "TH1.h"
61 #include "TH2.h"
62 #include "TMath.h"
63 #include "TGaxis.h"
64 #include "TMultiGraph.h"
65 #include "TLegend.h"
66 #include "TLatex.h"
67 #include "TStyle.h"
68 #include "THStack.h"
69 #include "TPaveText.h"
70 #include "TGraphAsymmErrors.h"
71 #include "TLegendEntry.h"
72 #include "TMarker.h"
73 
74 #include <cmath>
75 #include <iostream>
76 #include <vector>
77 #include <sstream>
78 #include <string>
79 #include <fstream>
80 #include <iomanip>
81 
82 #define ERROR std::cerr << __FUNCTION__ << " : " << __LINE__ << " " <<std::endl
83 
84 using namespace ana;
85 
86 
87 float GetLimits(TH1F *plot, float &hi, float &low) {
88 
89  low = 0.0;
90  hi = 0.0;
91 
92  for(int bin = 1; bin<plot->GetNbinsX()+1; bin++) {
93  if(plot->GetBinContent(bin+1) < plot->GetBinContent(bin)) {
94  if(plot->GetBinContent(bin+1) < 1.0 && plot->GetBinContent(bin) > 1.0 && bin < plot->GetNbinsX()) {
95  float step = plot->GetBinContent(bin) - 1.0;
96  float inbin = plot->GetBinContent(bin) - plot->GetBinContent(bin+1);
97  float frac = step/inbin;
98  float temp = plot->GetXaxis()->GetXmin() + plot->GetXaxis()->GetBinWidth(1)/2.0 + (float(bin-1)+frac)*plot->GetXaxis()->GetBinWidth(1);
99  std::cout << "Crossed chisq of 1, descending, with x value of: " << temp << std::endl;
100  low = temp; // this is tricky with two dips. This should keep the min and max for the second dip which is fine
101  }
102  }
103  else {
104  if(plot->GetBinContent(bin+1) > 1.0 && plot->GetBinContent(bin) < 1.0 && bin < plot->GetNbinsX()) {
105  float step = 1.0 - plot->GetBinContent(bin);
106  float inbin = plot->GetBinContent(bin+1) - plot->GetBinContent(bin);
107  float frac = step/inbin;
108  hi = plot->GetXaxis()->GetXmin() + plot->GetXaxis()->GetBinWidth(1)/2.0 + (float(bin-1)+frac)*plot->GetXaxis()->GetBinWidth(1);
109  std::cout << "Crossed chisq of 1, ascending, with x value of: " << hi << std::endl;
110  }
111  }
112  }
113 
114  return fabs(hi-low);
115 }
116 
117 
118 
119 void makeSystTable_reduced(bool makeFullSystTable = false)
120 {
121 
122  osc::IOscCalcAdjustable* calc = DefaultOscCalc(); // full systs normal h
123  osc::IOscCalcAdjustable* calcfake = DefaultOscCalc(); // full systs normal h
124 
125  //Absolute muon energy scale (`2%)
126  std::vector<const ISyst*> absMuESyst;
127  absMuESyst.push_back(&kMuEScaleSyst2017);
128  // Relative muon energy scale (`2%)
129  std::vector<const ISyst*> relMuESyst;
130  relMuESyst.push_back(&kRelMuEScaleSyst2017);
131  // Absolute hadronic energy scale (`5%)
132  std::vector<const ISyst*> absHadESyst;
133  absHadESyst.push_back(&kAna2017CalibrationSyst);
134  absHadESyst.push_back(&kAna2017CalibShapeSyst);
135  // Relative hadronic energy scale (`5%)
136  std::vector<const ISyst*> relHadESyst;
137  relHadESyst.push_back(&kAna2017RelativeCalibSyst);
138  // Normalization (`5%)
139  std::vector<const ISyst*> normSyst;
140  normSyst.push_back(&kAna2017NormBoth);
141  // Cross sections and final-state interaction
142  std::vector<const ISyst*> xsecSysts;
143  AddJointAna2017XSecSysts(xsecSysts, kNumuAna2017, 1);
144  // Neutrino flux
145  std::vector<const ISyst*> fluxSysts;
147  // Beam background normalization (`100%)
148  // Scintillation model
149  std::vector<const ISyst*> scintSysts;
150  scintSysts.push_back(&kAna2017LightlevelSyst);
151  scintSysts.push_back(&kAna2017CherenkovSyst);
152  // Total systematic uncertainty
153  std::vector<const ISyst*> totSysts = getAllAna2017Systs(kNumuAna2017);
154 
155 
156  std::map< std::string, std::vector<const ISyst*> > mapSystVecs;
157  mapSystVecs["Absolute muon energy scale"] = absMuESyst;
158  mapSystVecs["Relative muon energy scale"] = relMuESyst;
159  mapSystVecs["Absolute hadronic energy scale"] = absHadESyst;
160  mapSystVecs["Relative hadronic energy scale"] = relHadESyst;
161  mapSystVecs["Normalisation"] = normSyst;
162  mapSystVecs["Cross sections and final-state interaction"] = xsecSysts;
163  mapSystVecs["Neutrino flux"] = fluxSysts;
164  mapSystVecs["Detector response"] = scintSysts;
165 
166 
167  //mapSystVecs["Total systematic uncertainty"] = totSysts;
168 
169 
170  float theta_real = 0.468;
171  float msq = 2.44;
172  std::cout<< "\n++theta_real: " << theta_real
173  << "\n++msq: " << msq
174  <<std::endl;
175  kFitSinSqTheta23.SetValue(calc,theta_real);
176  kFitDmSq32Scaled.SetValue(calc,msq);
177  kFitDeltaInPiUnits.SetValue(calc,1.5);
178 
179  //PredictionSystNumu2017 (ENumu2017ExtrapType extrap, osc::IOscCalc *osc, int WhichQuant=1, const std::vector< const DummyNumu2017Syst * > &systs=kNumu2017Systs, const std::string fname="/nova/ana/nu_mu_ana/Ana2017/Predictions/provisional/pred_numu_reduced_v5.root")
180  const std::string fakeNDPredFileName = "/pnfs/nova/persistent/users/lvinton1/numu/pd/numu/3A/newSystPreds/fakeNDData/pred_full_fakeNDData.root";
181  ENumu2017ExtrapType extrap = kNuMu;
182  //PredictionSystNumu2017 prediction_Q1(extrap, calc, 1);
183  // PredictionSystNumu2017 prediction_Q2(extrap, calc, 2);
184  // PredictionSystNumu2017 prediction_Q3(extrap, calc, 3);
185  // PredictionSystNumu2017 prediction_Q4(extrap, calc, 4);
186  PredictionSystNumu2017 prediction_Q1(extrap, calc, 1, kNumu2017Systs, fakeNDPredFileName);
187  PredictionSystNumu2017 prediction_Q2(extrap, calc, 2, kNumu2017Systs, fakeNDPredFileName);
188  PredictionSystNumu2017 prediction_Q3(extrap, calc, 3, kNumu2017Systs, fakeNDPredFileName);
189  PredictionSystNumu2017 prediction_Q4(extrap, calc, 4, kNumu2017Systs, fakeNDPredFileName);
190 
191 
192 
193  // make fake data:
194  Spectrum fakeQ1 = prediction_Q1.Predict(calc).FakeData(kAna2017POT);
195  Spectrum fakeQ2 = prediction_Q2.Predict(calc).FakeData(kAna2017POT);
196  Spectrum fakeQ3 = prediction_Q3.Predict(calc).FakeData(kAna2017POT);
197  Spectrum fakeQ4 = prediction_Q4.Predict(calc).FakeData(kAna2017POT);
198 
199  // get cosmics:
200  // Get the cosmic distributions:
201  DontAddDirectory guard;
202  TFile* fCosmic = new TFile((TString)"/nova/app/users/bays/dev_xchecks/work/done5/cos_bkg_hists_double.root");
203  TH1D* hCos_Q1 = (TH1D*) fCosmic -> Get("q1all");
204  TH1D* hCos_Q2 = (TH1D*) fCosmic -> Get("q2all");
205  TH1D* hCos_Q3 = (TH1D*) fCosmic -> Get("q3all");
206  TH1D* hCos_Q4 = (TH1D*) fCosmic -> Get("q4all");
207  fCosmic -> Close();
208  delete fCosmic;
209 
210  std::cout<< "hCos_Q1 -> Integral(): " << hCos_Q1 -> Integral() <<std::endl;
211 
212  // add cosmics to fake data:
213  Spectrum specCos_Q1(hCos_Q1, kAna2017POT, kAna2017Livetime);
214  fakeQ1 += specCos_Q1;
215  Spectrum specCos_Q2(hCos_Q2, kAna2017POT, kAna2017Livetime);
216  fakeQ2 += specCos_Q2;
217  Spectrum specCos_Q3(hCos_Q3, kAna2017POT, kAna2017Livetime);
218  fakeQ3 += specCos_Q3;
219  Spectrum specCos_Q4(hCos_Q4, kAna2017POT, kAna2017Livetime);
220  fakeQ4 += specCos_Q4;
221 
222 
223  // make expt for each quantile:
224  SingleSampleExperiment* exptQ1 = new SingleSampleExperiment(&prediction_Q1, fakeQ1, specCos_Q1);
225  SingleSampleExperiment* exptQ2 = new SingleSampleExperiment(&prediction_Q2, fakeQ2, specCos_Q2);
226  SingleSampleExperiment* exptQ3 = new SingleSampleExperiment(&prediction_Q3, fakeQ3, specCos_Q3);
227  SingleSampleExperiment* exptQ4 = new SingleSampleExperiment(&prediction_Q4, fakeQ4, specCos_Q4);
228 
229  std::vector<const IExperiment*> exptsAllFake;
230  exptsAllFake.push_back(exptQ1);
231  exptsAllFake.push_back(exptQ2);
232  exptsAllFake.push_back(exptQ3);
233  exptsAllFake.push_back(exptQ4);
234  exptsAllFake.push_back(&kSolarConstraintsPDG2017);
235  exptsAllFake.push_back(WorldReactorConstraint2017());
236  MultiExperiment exptFakeCC(exptsAllFake);
237 
238  std::vector<const IFitVar*> oscpar = getNumuMarginalisedOscParam();
239  std::vector<const IFitVar*> oscparDM = {&kFitDmSq32Scaled};
240  std::vector<const IFitVar*> oscparSin = {&kFitSinSqTheta23};
241 
242 
243  int kXBins = 200;
244  int kYBins = 200;
245  double xLowEdge = 0.3;
246  double xHiEdge = 0.7;
247  double yLowEdge = 2.1;
248  double yHiEdge = 2.8;
249 
250  TH1F *xplot = new TH1F("xplot","xplot",kXBins,xLowEdge,xHiEdge);
251  TH1F *yplot = new TH1F("yplot","yplot",kYBins,yLowEdge,yHiEdge);
252 
253  double range = 0;
254  float stattheta = 0, statmsq = 0;
255  float hi, low;
256  float hix, lowx;
257 
258 
259  TFile* fOut = new TFile("systMagsHists_reducedTable_newPredsFakeNDData.root","RECREATE");
260 
261  kFitSinSqTheta23.SetValue(calc,theta_real);
262  kFitDmSq32Scaled.SetValue(calc,msq);
263  kFitDeltaInPiUnits.SetValue(calc,1.5);
264  kFitSinSqTheta23.SetValue(calcfake,theta_real);
265  kFitDmSq32Scaled.SetValue(calcfake,msq);
266  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
267 
268  xplot = (TH1F*)Profile(&exptFakeCC, calcfake,
270  kXBins, xLowEdge, xHiEdge, 0,
271  oscparDM);
272  //{&kFitDmSq32Scaled});
273  xplot -> Write("ssq_stat");
274  TH1F* hXStat = (TH1F*) xplot -> Clone();
275 
276  kFitSinSqTheta23.SetValue(calc,theta_real);
277  kFitDmSq32Scaled.SetValue(calc,msq);
278  kFitDeltaInPiUnits.SetValue(calc,1.5);
279  kFitSinSqTheta23.SetValue(calcfake,theta_real);
280  kFitDmSq32Scaled.SetValue(calcfake,msq);
281  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
282 
283  yplot = (TH1F*)Profile(&exptFakeCC, calcfake,
285  kYBins, yLowEdge, yHiEdge, 0,
286  oscparSin);
287  //{&kFitSinSqTheta23});
288  yplot -> Write("msq_stat");
289  TH1F* hYStat = (TH1F*) yplot -> Clone();
290 
291 
292  range = GetLimits(xplot, hi, low);
293  stattheta = range;
294  float statthetaup = hi-theta_real;
295  float statthetadn = theta_real-low;
296  std::cout<< std::setprecision(2)
297  << "Raw. sin^2 low: " << low << ", high: " << hi
298  << ", updated bf: " << xplot->GetBinCenter(xplot->GetMinimumBin())
299  <<std::endl;
300 
301 
302  range = GetLimits(yplot, hi, low);
303  statmsq = range;
304  float statmsqup = hi-msq;
305  float statmsqdn = msq-low;
306  std::cout<< std::setprecision(2)
307  << "Raw. dm^2 low: " << low << ", high: " << hi
308  << ", updated bf: " << yplot->GetBinCenter(yplot->GetMinimumBin())
309  <<std::endl;
310 
311  std::cout << "++stats only delta m^2 range is: " << range << " " << range/2.0 << " " << range/(2*msq) << " " << (hi-msq)/msq << " " << (msq-low)/msq << " " << (hi-msq) << " " << (msq-low) << std::endl;
312 
313  std::cout << "LatexFormatLegend: in sin (10^-3) , dM (10^-6)" << std::endl;
314  std::cout<< std::setprecision(2)
315  << "LatexFormat stats only & +" << statthetaup * pow(10,3) << " / -"
316  << statthetadn * pow(10,3) << " & +" << statmsqup * pow(10,3) << " / - "
317  << statmsqdn * pow(10,3) << "\\\\"
318  <<std::endl;
319 
320  std::cout<< "LatexFormatFracSin Legend: in sin (10^-3) , dM (10^-6)" << std::endl;
321  std::cout<< std::setprecision(2)
322  << "LatexFormatFracSin stats only & $\\pm$ " << stattheta * pow(10,3) / 2
323  << " & +" << statmsqup * pow(10,3) << " / - "
324  << statmsqdn * pow(10,3) << "\\\\"
325  <<std::endl;
326 
327 
328 
329  if(!makeFullSystTable){
330  for(auto const& thisSystVec: mapSystVecs){
331 
332  kFitSinSqTheta23.SetValue(calc,theta_real);
333  kFitDmSq32Scaled.SetValue(calc,msq);
334  kFitDeltaInPiUnits.SetValue(calc,1.5);
335 
336  kFitSinSqTheta23.SetValue(calcfake,theta_real);
337  kFitDmSq32Scaled.SetValue(calcfake,msq);
338  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
339  xplot = (TH1F*)Profile(&exptFakeCC, calcfake,
341  kXBins, xLowEdge, xHiEdge, 0,
342  oscparDM,
343  //{&kFitDmSq32Scaled},
344  thisSystVec.second);
345  xplot -> Write("ssq_" + TString(thisSystVec.first));
346  TH1F* hXDiff = (TH1F*) hXStat -> Clone();
347  hXDiff -> Add(xplot, -1);
348  hXDiff -> Write("ssqDiff_" + TString(thisSystVec.first));
349 
350  kFitSinSqTheta23.SetValue(calc,theta_real);
351  kFitDmSq32Scaled.SetValue(calc,msq);
352  kFitDeltaInPiUnits.SetValue(calc,1.5);
353 
354  kFitSinSqTheta23.SetValue(calcfake,theta_real);
355  kFitDmSq32Scaled.SetValue(calcfake,msq);
356  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
357  yplot = (TH1F*)Profile(&exptFakeCC, calcfake,
359  kYBins, yLowEdge, yHiEdge, 0,
360  oscparSin,
361  //{&kFitSinSqTheta23},
362  thisSystVec.second);
363  yplot -> Write("msq_" + TString(thisSystVec.first));
364  TH1F* hYDiff = (TH1F*) hYStat -> Clone();
365  hYDiff -> Add(yplot, -1);
366  hYDiff -> Write("msqDiff_" + TString(thisSystVec.first));
367 
368  range = GetLimits(yplot, hi, low);
369  range = GetLimits(xplot, hix, lowx);
370 
371  // print out for table
372  std::cout << std::setprecision(2)
373  << "LatexFormatFracSin Frac. uncer. " << thisSystVec.first << " & $\\pm$ "
374  << sqrt( (hix-lowx)*(hix-lowx) - (stattheta*stattheta) ) * pow(10,3) / 2
375  << " & +"
376  << sqrt((hi-msq)*(hi-msq)-statmsqup*statmsqup) * pow(10,3) << " / -"
377  << sqrt((msq-low)*(msq-low)-statmsqdn*statmsqdn) * pow(10,3) << "\\\\" <<std::endl;
378 
379  // print out for checks
380  std::cout<< std::setprecision(9)
381  << "Raw. sin^2 low: " << lowx
382  << ", high: " << hix
383  << ", updated bf: " << xplot->GetBinCenter(xplot->GetMinimumBin())
384  <<std::endl;
385  std::cout<< std::setprecision(9)
386  << "Raw. dm^2 low: " << low
387  << ", high: " << hi
388  << ", updated bf: " << yplot->GetBinCenter(yplot->GetMinimumBin())
389  <<std::endl;
390 
391  }
392  }
393 
394 
395  if(makeFullSystTable){
396  for(const ISyst* syst: totSysts) {
397 
398  kFitSinSqTheta23.SetValue(calc,theta_real);
399  kFitDmSq32Scaled.SetValue(calc,msq);
400  kFitDeltaInPiUnits.SetValue(calc,1.5);
401  kFitSinSqTheta23.SetValue(calcfake,theta_real);
402  kFitDmSq32Scaled.SetValue(calcfake,msq);
403  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
404 
405  xplot = (TH1F*)Profile(&exptFakeCC, calcfake,
407  kXBins, xLowEdge, xHiEdge, 0,
408  oscparDM,
409  //{&kFitDmSq32Scaled},
410  {syst});
411  xplot -> Write("ssq_" + TString(syst->ShortName()));
412  TH1F* hXDiff = (TH1F*) hXStat -> Clone();
413  hXDiff -> Add(xplot, -1);
414  hXDiff -> Write("ssqDiff_" + TString(syst->ShortName()));
415 
416  kFitSinSqTheta23.SetValue(calc,theta_real);
417  kFitDmSq32Scaled.SetValue(calc,msq);
418  kFitDeltaInPiUnits.SetValue(calc,1.5);
419  kFitSinSqTheta23.SetValue(calcfake,theta_real);
420  kFitDmSq32Scaled.SetValue(calcfake,msq);
421  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
422 
423  yplot = (TH1F*)Profile(&exptFakeCC, calcfake,
425  kYBins, yLowEdge, yHiEdge, 0,
426  oscparSin,
427  //{&kFitSinSqTheta23},
428  {syst});
429  yplot -> Write("msq_" + TString(syst->ShortName()));
430 
431  TH1F* hYDiff = (TH1F*) hYStat -> Clone();
432  hYDiff -> Add(yplot, -1);
433  hYDiff -> Write("msqDiff_" + TString(syst->ShortName()));
434 
435  range = GetLimits(yplot, hi, low);
436  range = GetLimits(xplot, hix, lowx);
437 
438 
439  std::cout << std::setprecision(2)
440  << "LatexFormatFracSin Frac. uncer. " << syst->ShortName() << " & $\\pm$ "
441  << sqrt( (hix-lowx)*(hix-lowx) - (stattheta*stattheta) ) * pow(10,3) / 2
442  << " & +"
443  << sqrt((hi-msq)*(hi-msq)-statmsqup*statmsqup) * pow(10,3) << " / -"
444  << sqrt((msq-low)*(msq-low)-statmsqdn*statmsqdn) * pow(10,3) << "\\\\" <<std::endl;
445 
446 
447  std::cout<< std::setprecision(9)
448  << "Raw. sin^2 low: " << lowx
449  << ", high: " << hix
450  << ", updated bf: " << xplot->GetBinCenter(xplot->GetMinimumBin())
451  <<std::endl;
452  std::cout<< std::setprecision(9)
453  << "Raw. dm^2 low: " << low
454  << ", high: " << hi
455  << ", updated bf: " << yplot->GetBinCenter(yplot->GetMinimumBin())
456  <<std::endl;
457 
458  }
459  }
460 
461 
462 
463  // get effect of deltaCP
464  xplot = (TH1F*)Profile(&exptFakeCC, calcfake,
466  kXBins, xLowEdge, xHiEdge, 0,
468  xplot -> Write("ssq_dcp");
469 
470  yplot = (TH1F*)Profile(&exptFakeCC, calcfake,
472  kYBins, yLowEdge, yHiEdge, 0,
474  yplot -> Write("msq_dcp");
475 
476  range = GetLimits(yplot, hi, low);
477  range = GetLimits(xplot, hix, lowx);
478 
479  std::cout<<setprecision(4)
480  << "\n\n\nrange: " << range
481  << ", statrange: " << stattheta
482  << "\nhi: " << hi
483  << ", low: " << low
484  << ", msq: " << msq
485  << "\n\n\n" <<std::endl;
486 
487  //float statmsqup = hi-msq;
488  double deltaCP_sinSq = (range-stattheta)*pow(10,3)/2;
489  double deltaCP_dmSqUp = ( (hi-msq) - (statmsqup) )*pow(10,3);
490  double deltaCP_dmSqDn = ( (msq-low) - (statmsqdn) )*pow(10,3);
491 
492  std::cout<< "Effect of deltaCP: "
493  << "ssqth23: (range - statrange)*pow(10,3)/2: " << deltaCP_sinSq
494  << "\n" << "dmsq32 up: " << deltaCP_dmSqUp
495  << "\n" << "dmsq32 dn: " << deltaCP_dmSqDn
496  <<std::endl;
497 
498  std::cout << std::setprecision(2)
499  << "LatexFormatFracSin Frac. uncer. $\\delta_{CP}$ (0-$2\\pi$) & $\\pm$ "
500  << deltaCP_sinSq
501  << " & +"
502  << deltaCP_dmSqUp
503  << " / -"
504  << deltaCP_dmSqDn
505  << "\\\\" <<std::endl;
506 
507 
508 
509  std::cout<< std::setprecision(9)
510  << "Raw. sin^2 low: " << lowx << ", high: " << hix
511  << ", updated bf: " << xplot->GetBinCenter(xplot->GetMinimumBin())
512  <<std::endl;
513  std::cout<< "Raw. dm^2 low: " << low << ", high: " << hi
514  << ", updated bf: " << yplot->GetBinCenter(yplot->GetMinimumBin())
515  <<std::endl;
516 
517 
518 
519 
520  // get total uncertainty
521  kFitSinSqTheta23.SetValue(calc,theta_real);
522  kFitDmSq32Scaled.SetValue(calc,msq);
523  kFitDeltaInPiUnits.SetValue(calc,1.5);
524 
525  kFitSinSqTheta23.SetValue(calcfake,theta_real);
526  kFitDmSq32Scaled.SetValue(calcfake,msq);
527  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
528  xplot = (TH1F*)Profile(&exptFakeCC, calcfake,
530  kXBins, xLowEdge, xHiEdge, 0,
531  oscparDM,
532  //{&kFitDmSq32Scaled},
533  totSysts);
534  xplot -> Write("ssq_totSysts");
535  TH1F* hXDiff = (TH1F*) hXStat -> Clone();
536  hXDiff -> Add(xplot, -1);
537  hXDiff -> Write("ssqDiff_totSysts");
538 
539  kFitSinSqTheta23.SetValue(calc,theta_real);
540  kFitDmSq32Scaled.SetValue(calc,msq);
541  kFitDeltaInPiUnits.SetValue(calc,1.5);
542 
543  kFitSinSqTheta23.SetValue(calcfake,theta_real);
544  kFitDmSq32Scaled.SetValue(calcfake,msq);
545  kFitDeltaInPiUnits.SetValue(calcfake,1.5);
546  yplot = (TH1F*)Profile(&exptFakeCC, calcfake,
548  kYBins, yLowEdge, yHiEdge, 0,
549  oscparSin,
550  //{&kFitSinSqTheta23},
551  totSysts);
552  yplot -> Write("msq_totSysts");
553  TH1F* hYDiff = (TH1F*) hYStat -> Clone();
554  hYDiff -> Add(yplot, -1);
555  hYDiff -> Write("msqDiff_totSysts");
556 
557  range = GetLimits(yplot, hi, low);
558  range = GetLimits(xplot, hix, lowx);
559 
560  // print out for table
561  std::cout << std::setprecision(2)
562  << "LatexFormatFracSin Frac. uncer. total syst uncertainty & $\\pm$ "
563  << sqrt( (hix-lowx)*(hix-lowx) - (stattheta*stattheta) ) * pow(10,3) / 2
564  << " & +"
565  << sqrt((hi-msq)*(hi-msq)-statmsqup*statmsqup) * pow(10,3) << " / -"
566  << sqrt((msq-low)*(msq-low)-statmsqdn*statmsqdn) * pow(10,3) << "\\\\" <<std::endl;
567 
568 
569  // print out for checks
570  std::cout<< std::setprecision(9)
571  << "Raw. sin^2 low: " << lowx << ", high: " << hix
572  << ", updated bf: " << xplot->GetBinCenter(xplot->GetMinimumBin())
573  <<std::endl;
574  std::cout<< std::setprecision(9)
575  << "Raw. dm^2 low: " << low << ", high: " << hi
576  << ", updated bf: " << yplot->GetBinCenter(yplot->GetMinimumBin())
577  <<std::endl;
578 
579  // print out for table
580  std::cout << std::setprecision(2)
581  << "LatexFormatFracSin Frac. uncer. total syst uncertainty with dCP linear & $\\pm$ "
582  << (sqrt( (hix-lowx)*(hix-lowx) - (stattheta*stattheta) ) * pow(10,3) / 2) + deltaCP_sinSq
583  << " & +"
584  << (sqrt((hi-msq)*(hi-msq)-statmsqup*statmsqup) * pow(10,3)) + deltaCP_dmSqUp << " / -"
585  << (sqrt((msq-low)*(msq-low)-statmsqdn*statmsqdn) * pow(10,3)) + deltaCP_dmSqDn << "\\\\"
586  <<std::endl;
587 
588 
589  std::cout<< std::setprecision(2)
590  << "LatexFormatFracSin Total uncertainty & $\\pm$ " << ( (hix-lowx) * pow(10,3) / 2)
591  << " & +" << ((hi-msq) * pow(10,3)) << " / - "
592  << ((msq-low) * pow(10,3)) << "\\\\"
593  <<std::endl;
594 
595  std::cout<< std::setprecision(2)
596  << "LatexFormatFracSin Total uncertainty including dCP & $\\pm$ " << ( (hix-lowx) * pow(10,3) / 2) + deltaCP_sinSq
597  << " & +" << ((hi-msq) * pow(10,3)) + deltaCP_dmSqUp << " / - "
598  << ((msq-low) * pow(10,3)) + deltaCP_dmSqDn << "\\\\"
599  <<std::endl;
600 
601 
602 
603  fOut -> Close();
604  delete fOut;
605 
606 
607 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double Integral(const Spectrum &s, const double pot, int cvnregion)
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
float GetLimits(TH1F *plot, float &hi, float &low)
const double kAna2017POT
Definition: Exposures.h:174
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
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
osc::OscCalcDumb calc
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
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
const std::vector< const DummyNumu2017Syst * > kNumu2017Systs
TSpline3 hi("hi", xhi, yhi, 18,"0")
const double kAna2017Livetime
Definition: Exposures.h:200
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:341
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Loads shifted spectra from files.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
double frac(double x)
Fractional part.
void AddJointAna2017XSecSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana, bool smallgenie)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
float bin[41]
Definition: plottest35.C:14
c1 Close()
OStream cout
Definition: OStream.cxx:6
Spectrum Predict(osc::IOscCalc *calc) const override
Combine multiple component experiments.
::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)
void AddJointAna2017BeamSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
void makeSystTable_reduced(bool makeFullSystTable=false)
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
const DummyAnaSyst kAna2017NormBoth("NormBoth","Sig+Bkg Norm.")
Compare a single data spectrum to the MC + cosmics expectation.
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
gm Write()
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9