Macros | Functions
makeSystTable_reduced.C File Reference
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionNumuFAHadE.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Prediction/PredictionSystNumu2017.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "3FlavorAna/Systs/JointAna2017Systs.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/SolarConstraints.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Vars/FitVarsNumu.h"
#include "CAFAna/FC/FCCollection.h"
#include "CAFAna/FC/FCPoint.h"
#include "CAFAna/FC/FCSurface.h"
#include "OscLib/OscCalcGeneral.h"
#include "OscLib/OscCalcPMNS.h"
#include "OscLib/OscCalc.h"
#include "OscLib/IOscCalc.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TGaxis.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TStyle.h"
#include "THStack.h"
#include "TPaveText.h"
#include "TGraphAsymmErrors.h"
#include "TLegendEntry.h"
#include "TMarker.h"
#include <cmath>
#include <iostream>
#include <vector>
#include <sstream>
#include <string>
#include <fstream>
#include <iomanip>

Go to the source code of this file.

Macros

#define ERROR   std::cerr << __FUNCTION__ << " : " << __LINE__ << " " <<std::endl
 

Functions

float GetLimits (TH1F *plot, float &hi, float &low)
 
void makeSystTable_reduced (bool makeFullSystTable=false)
 

Macro Definition Documentation

#define ERROR   std::cerr << __FUNCTION__ << " : " << __LINE__ << " " <<std::endl

Definition at line 82 of file makeSystTable_reduced.C.

Referenced by rwgt::MakeGENIEReweightTable::beginJob(), and genie::Messenger::PriorityFromString().

Function Documentation

float GetLimits ( TH1F *  plot,
float &  hi,
float &  low 
)

Definition at line 87 of file makeSystTable_reduced.C.

References bin, om::cout, allTimeWatchdog::endl, stan::math::fabs(), check_time_usage::float, util::frac(), fillBadChanDBTables::step, and APDHVSetting::temp.

Referenced by SystEntry::getDmDcp(), makeSystTable_reduced(), SystEntry::printTotalUncertaintyLinearDcp(), SystEntry::printUncertainty(), SystEntry::printUncertaintyLinear(), SystEntry::printUncertaintySyst(), and SystEntry::printUncertaintySystLinearDcp().

87  {
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 }
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
hi
TSpline3 hi("hi", xhi, yhi, 18,"0")
check_time_usage.float
float
Definition: check_time_usage.py:7
util::frac
double frac(double x)
Fractional part.
Definition: RemoveBeamSpillsPrecise_module.cc:83
bin
float bin[41]
Definition: plottest35.C:14
om::cout
OStream cout
Definition: OStream.cxx:6
plot
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
Definition: plot_resolution.C:184
void makeSystTable_reduced ( bool  makeFullSystTable = false)

Definition at line 119 of file makeSystTable_reduced.C.

References ana::AddJointAna2017BeamSysts(), ana::AddJointAna2017XSecSysts(), calc, plot_validation_datamc::Clone(), Close(), om::cout, ana::DefaultOscCalc(), allTimeWatchdog::endl, ana::Spectrum::FakeData(), makeBrightnessMap::fOut, Get, ana::getAllAna2017Systs(), GetLimits(), ana::getNumuMarginalisedOscParam(), hi(), Integral(), ana::kAna2017Livetime, ana::kAna2017NormBoth, ana::kAna2017POT, ana::kFitDeltaInPiUnits, ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, ana::kMuEScaleSyst2017, ana::kNuMu, ana::kNumu2017Systs, ana::kNumuAna2017, ana::kRelMuEScaleSyst2017, ana::kSolarConstraintsPDG2017, cet::pow(), ana::PredictionInterp::Predict(), ana::Profile(), PandAna.Demos.demo1::range, ana::FitVarWithPrior::SetValue(), ana::ConstrainedFitVarWithPrior::SetValue(), std::sqrt(), string, ana::WorldReactorConstraint2017(), and Write().

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;
146  AddJointAna2017BeamSysts(fluxSysts, kNumuAna2017);
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,
269  &kFitSinSqTheta23,
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,
284  &kFitDmSq32Scaled,
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,
340  &kFitSinSqTheta23,
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,
358  &kFitDmSq32Scaled,
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,
406  &kFitSinSqTheta23,
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,
424  &kFitDmSq32Scaled,
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,
465  &kFitSinSqTheta23,
466  kXBins, xLowEdge, xHiEdge, 0,
467  {&kFitDmSq32Scaled, &kFitDeltaInPiUnits});
468  xplot -> Write("ssq_dcp");
469 
470  yplot = (TH1F*)Profile(&exptFakeCC, calcfake,
471  &kFitDmSq32Scaled,
472  kYBins, yLowEdge, yHiEdge, 0,
473  {&kFitSinSqTheta23, &kFitDeltaInPiUnits});
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,
529  &kFitSinSqTheta23,
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,
547  &kFitDmSq32Scaled,
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 }
Integral
double Integral(const Spectrum &s, const double pot, int cvnregion)
Definition: syst_table_from_preds.C:46
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
GetLimits
float GetLimits(TH1F *plot, float &hi, float &low)
Definition: makeSystTable_reduced.C:87
ana::kAna2017POT
const double kAna2017POT
Definition: Exposures.h:174
ana::Profile
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
Get
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
ana::FitVarWithPrior::SetValue
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
Definition: FitVarWithPrior.h:51
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
ana::DefaultOscCalc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
ana::ISyst
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
ana::kMuEScaleSyst2017
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
Definition: EnergySysts2017.h:32
ana::kNumu2017Systs
const std::vector< const DummyNumu2017Syst * > kNumu2017Systs
Definition: PredictionSystNumu2017.h:70
hi
TSpline3 hi("hi", xhi, yhi, 18,"0")
ana::kAna2017Livetime
const double kAna2017Livetime
Definition: Exposures.h:200
ana::kRelMuEScaleSyst2017
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
Definition: EnergySysts2017.h:57
ana::Spectrum::FakeData
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
ana::ConstrainedFitVarWithPrior::SetValue
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVarWithPrior.h:104
ana::PredictionSystNumu2017
Loads shifted spectra from files.
Definition: PredictionSystNumu2017.h:115
ana::WorldReactorConstraint2017
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
Definition: ReactorExperiment.cxx:115
ana::AddJointAna2017XSecSysts
void AddJointAna2017XSecSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana, bool smallgenie)
Definition: JointAna2017Systs.cxx:91
ana::kFitDmSq32Scaled
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
Definition: VarsWithPriors.cxx:13
Close
c1 Close()
om::cout
OStream cout
Definition: OStream.cxx:6
ana::MultiExperiment
Combine multiple component experiments.
Definition: MultiExperiment.h:11
ana::kSolarConstraintsPDG2017
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
Definition: SolarConstraints.h:56
ana::AddJointAna2017BeamSysts
void AddJointAna2017BeamSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana)
Definition: JointAna2017Systs.cxx:102
ana::kFitSinSqTheta23
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Definition: VarsWithPriors.cxx:9
osc::_IOscCalcAdjustable< double >
ana::getAllAna2017Systs
std::vector< const ISyst * > getAllAna2017Systs(const EAnaType2017 ana, const bool smallgenie)
Definition: JointAna2017Systs.cxx:146
ana::DontAddDirectory
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
ana::kFitDeltaInPiUnits
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
Definition: VarsWithPriors.cxx:17
ana::kAna2017NormBoth
const DummyAnaSyst kAna2017NormBoth("NormBoth","Sig+Bkg Norm.")
Definition: JointAna2017Systs.h:36
ana::SingleSampleExperiment
Compare a single data spectrum to the MC + cosmics expectation.
Definition: SingleSampleExperiment.h:21
Write
gm Write()
ana::getNumuMarginalisedOscParam
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
string
enum BeamMode string
Generated on Sat Jan 23 2021 14:14:42 for NOvA by   doxygen 1.8.11