syst_table_fit.C
Go to the documentation of this file.
1 #include "joint_fit_2018_tools.h"
6 #include "CAFAna/Fit/Fit.h"
8 #include "CAFAna/Systs/Systs.h"
11 #include "CAFAna/Analysis/Plots.h"
12 
13 using namespace ana;
14 
15 float GetLimits(TH1F *plot, float &hi, float &low) {
16 
17  low = 0.0;
18  hi = 0.0;
19 
20  for(int bin = 1; bin <= plot->GetNbinsX(); bin++) {
21  if(plot->GetBinContent(bin+1) < plot->GetBinContent(bin)) {
22  if(plot->GetBinContent(bin+1) < 1.0 && plot->GetBinContent(bin) > 1.0 && bin < plot->GetNbinsX()) {
23  float step = plot->GetBinContent(bin) - 1.0;
24  float inbin = plot->GetBinContent(bin) - plot->GetBinContent(bin+1);
25  float frac = step/inbin;
26  low = plot->GetXaxis()->GetXmin() + plot->GetXaxis()->GetBinWidth(1)/2.0 + (float(bin-1)+frac)*plot->GetXaxis()->GetBinWidth(1);
27  std::cout << "Crossed chisq of 1, descending, with x value of: " << low << std::endl;
28  }
29  }
30  else {
31  if(plot->GetBinContent(bin+1) > 1.0 && plot->GetBinContent(bin) < 1.0 && bin < plot->GetNbinsX()) {
32  float step = 1.0 - plot->GetBinContent(bin);
33  float inbin = plot->GetBinContent(bin+1) - plot->GetBinContent(bin);
34  float frac = step/inbin;
35  hi = plot->GetXaxis()->GetXmin() + plot->GetXaxis()->GetBinWidth(1)/2.0 + (float(bin-1)+frac)*plot->GetXaxis()->GetBinWidth(1);
36  std::cout << "Crossed chisq of 1, ascending, with x value of: " << hi << std::endl;
37  }
38  }
39  }
40 
41  return fabs(hi-low);
42 }
43 
45 {
46  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
47  prelim->SetTextColor(kAzure);
48  prelim->SetNDC();
49  prelim->SetTextSize(1.5/30.);
50  prelim->SetTextAlign(32);
51  prelim->Draw();
52 }
53 
54 void syst_table_fit(bool full = false)
55 {
56 
57  // 2017 Best Fit
58  // const double dCP = 1.21;
59  // const double ssth23 = 0.558;
60  // const double dmsq32 = 2.44e-3;
61  // 2018 Best Fit
62  const double dCP = 0.186;
63  const double ssth23 = 0.585;
64  const double dmsq32 = 2.501e-3;
65 
66  auto calc = DefaultOscCalc();
67  SetFakeCalc(calc, ssth23, dmsq32, dCP);
68 
69  // All the numus
70  const int nqs = 4;
71 
72  std::vector <const IPrediction*> preds;
73  std::vector <std::pair<TH1D*, double>> cosmics;
74  std::vector <Spectrum*> data;
75  std::vector <const IExperiment*> expts;
76 
77  // Get Predictions
78  preds.push_back(GetNuePrediction2018("combo", DefaultOscCalc(), true, "fhc", false));
79  preds.push_back(GetNuePrediction2018("prop", DefaultOscCalc(), true, "rhc", false));
80  auto numu_preds_fhc = GetNumuPredictions2018(nqs, true, "fhc");
81  preds.insert(preds.end(),numu_preds_fhc.begin(), numu_preds_fhc.end());
82  auto numu_preds_rhc = GetNumuPredictions2018(nqs, true, "rhc");
83  preds.insert(preds.end(),numu_preds_rhc.begin(), numu_preds_rhc.end());
84 
85  // Load the cosmics
86  cosmics.push_back(GetNueCosmics2018("fhc"));
87  cosmics.push_back(GetNueCosmics2018("rhc"));
88  auto numu_cosmics_fhc = GetNumuCosmics2018(nqs, "fhc");
89  cosmics.insert(cosmics.end(),numu_cosmics_fhc.begin(), numu_cosmics_fhc.end());
90  auto numu_cosmics_rhc = GetNumuCosmics2018(nqs, "rhc");
91  cosmics.insert(cosmics.end(),numu_cosmics_rhc.begin(), numu_cosmics_rhc.end());
92 
93  // make fake data:
94  for(unsigned int i = 0; i < preds.size(); ++i){
95  double POT;
96  if(i == 0 || (i >= 2 && i < 6)) POT = kAna2018FHCPOT;
97  if(i == 1 || (i >= 6 && i < 10)) POT = kAna2018RHCPOT;
98 
99  Spectrum* thisdata = GetFakeData (preds[i], calc, POT, cosmics[i].first);
100  expts.push_back(new SingleSampleExperiment(preds[i],*thisdata, cosmics[i].first,cosmics[i].second));
101  }
102 
103  expts.push_back(WorldReactorConstraint2017());
104 
105  std::cout << "\nCreating multiexperiment\n" << std::endl;
106  std::cout << "Creating Multiexperiment with a total of "
107  << expts.size() << " experiments\n\n" << std::endl;
108  auto expt = new MultiExperiment(expts);
109 
110  // Setting Systematics
111  // Full list
112  std::vector<const ISyst*> totSysts = GetJointFitSystematicList(true);
113 
114  // Categories
115  std::vector<const ISyst*> xsecSysts;
116  AddJointAna2018XSecSysts(xsecSysts, kJointAna2018, true);
117  std::vector<const ISyst*> fluxSysts;
119  std::vector<const ISyst*> calibSysts;
120  calibSysts.push_back(&kAnaCalibrationSyst);
121  calibSysts.push_back(&kAnaCalibShapeSyst);
122  calibSysts.push_back(&kAnaRelativeCalibSyst);
123  std::vector<const ISyst*> scintSysts;
124  scintSysts.push_back(&kAnaLightlevelSyst);
125  scintSysts.push_back(&kAnaCherenkovSyst);
126  std::vector<const ISyst*> normSysts;
127  normSysts.push_back(&kAna2018NormFHC);
128  normSysts.push_back(&kAna2018NormRHC);
129  std::vector<const ISyst*> absMuESysts;
130  absMuESysts.push_back(&kMuEScaleSyst2017);
131  absMuESysts.push_back(&kRelMuEScaleSyst2017);
132  std::vector<const ISyst*> neutronSysts;
133  neutronSysts.push_back(&kNeutronVisEScalePrimariesSyst2018);
134  std::vector<const ISyst*> nearFarSysts;
135  nearFarSysts.push_back(&kNueAcceptSystSignalKin2018FHC);
136  nearFarSysts.push_back(&kNueAcceptSystBkg2018FHC);
137  nearFarSysts.push_back(&kNueAcceptSystSignalKin2018RHC);
138  nearFarSysts.push_back(&kNueAcceptSystBkg2018RHC);
139  nearFarSysts.push_back(&kMichelTaggingSyst2018);
140 
141  std::map< std::string, std::vector<const ISyst*> > mapSystVecs;
142  mapSystVecs["Muon Energy Scale"] = absMuESysts;
143  mapSystVecs["Detector Calibration"] = calibSysts;
144  mapSystVecs["Normalization"] = normSysts;
145  mapSystVecs["Neutron Uncertainty"] = neutronSysts;
146  mapSystVecs["Neutrino Cross Sections"] = xsecSysts;
147  mapSystVecs["Beam Flux"] = fluxSysts;
148  mapSystVecs["Near-Far Differences"] = nearFarSysts;
149  mapSystVecs["Detector Response"] = scintSysts;
150 
151  // Set Correlations
152  expt->SetSystCorrelations(0, GetCorrelations(true, true));
153  expt->SetSystCorrelations(1, GetCorrelations(true, false));
154  auto notnumufhc = GetCorrelations(false, true);
155  auto notnumurhc = GetCorrelations(false, false);
156  for(int i = 2; i < 6; ++i) expt->SetSystCorrelations(i, notnumufhc);
157  for(int i = 6; i < 10; ++i) expt->SetSystCorrelations(i, notnumurhc);
158 
159  // Define what to profile over
160  const int steps = 50;
161  struct ProfDef{
162  const IFitVar * fitvar;
163  double minx;
164  double maxx;
165  std::vector<const IFitVar *> profvars;
166  double bf;
167  int oom;
170  std::string suffix;
171  };
172 
173  std::vector<ProfDef> profiles = {
175  "$\\mathrm{sin}^2\\theta_{23}$ ($\\times 10^{-3}$)", "Uncertainty in sin^{2}#theta_{23} (#times10^{-3})", "th23"},
177  "$\\delta_{CP}/\\pi$", "Uncertainty in #delta_{CP}/#pi", "dcp"},
178  {&kFitDmSq32, 2e-3, 3e-3, {&kFitDeltaInPiUnits, &kFitSinSqTheta23}, dmsq32, 3,
179  "$\\Delta m^2_{32}$ ($\\times 10^{-3}~\\mathrm{eV}^2$)", "Uncertainty in #Deltam^{2}_{32} (#times10^{-3} eV^{2})", "dmsq"}
180  };
181 
182  // outputs
183  TFile* fOut = new TFile("syst_table_profiles.root","RECREATE");
184 
185  ofstream txt;
186  txt.open("syst_table.txt");
187 
188  txt<<std::setprecision(2);
189 
190  // Profiles n things
191  SetFakeCalc(calc, ssth23, dmsq32, dCP);
192 
193  double range;
194  float hi, low;
195  std::map<std::string, double> statup;
196  std::map<std::string, double> statdn;
197  std::map<std::string,std::map<std::string,std::pair<double,double>>> bars;
198 
199  txt<< "Source of Uncertainty";
200  for(auto prof: profiles) txt<< " & " << prof.latex;
201  txt<< " \\\\" << std::endl;
202  txt<< "\\hline" << std::endl;
203 
204  // Stat only
205  for(auto prof: profiles){
206  SetFakeCalc(calc, ssth23, dmsq32, dCP);
207 
208  TH1F* plot = (TH1F*)Profile(expt, calc,
209  prof.fitvar,
210  steps, prof.minx, prof.maxx, 0,
211  prof.profvars);
212  plot -> Write((prof.suffix+"_stat").c_str());
213 
214  range = GetLimits(plot, hi, low);
215  if(prof.suffix=="dcp") low-=2;
216 
217  double up = hi-prof.bf;
218  double dn = prof.bf-low;
219 
220  statup[prof.suffix] = up;
221  statdn[prof.suffix] = dn;
222  }
223 
224  // Systs
225  if(!full){
226  for(auto const &systVec: mapSystVecs) {
227  txt<< systVec.first ;
228  for(auto prof: profiles){
229  SetFakeCalc(calc, ssth23, dmsq32, dCP);
230 
231  TH1F* systplot = (TH1F*)Profile(expt, calc,
232  prof.fitvar,
233  steps, prof.minx, prof.maxx, 0,
234  prof.profvars,
235  systVec.second);
236  systplot -> Write((prof.suffix+"_" + systVec.first).c_str());
237 
238  range = GetLimits(systplot, hi, low);
239  if(prof.suffix=="dcp") low-=2;
240 
241  double plus = sqrt( fabs( (hi-prof.bf)*(hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) * pow(10,prof.oom);
242  double minus = sqrt( fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) * pow(10,prof.oom);
243 
244  txt<< " & +"
245  << plus
246  << " / -"
247  << minus ;
248 
249  bars[prof.suffix][systVec.first] = {minus,plus};
250  }
251  txt<< " \\\\" <<std::endl;
252  }
253  }
254 
255 
256  if(full){
257  for(const ISyst* syst: totSysts) {
258  txt<< syst->LatexName() ;
259  for(auto prof: profiles){
260  SetFakeCalc(calc, ssth23, dmsq32, dCP);
261 
262  TH1F* systplot = (TH1F*)Profile(expt, calc,
263  prof.fitvar,
264  steps, prof.minx, prof.maxx, 0,
265  prof.profvars,
266  {syst});
267  systplot -> Write((prof.suffix+"_" + syst->ShortName()).c_str());
268 
269  range = GetLimits(systplot, hi, low);
270  if(prof.suffix=="dcp") low-=2;
271 
272  txt<< " & +"
273  << sqrt( fabs( (hi-prof.bf)*(hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) * pow(10,prof.oom)
274  << " / -"
275  << sqrt( fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) * pow(10,prof.oom) ;
276  }
277  txt<< " \\\\" <<std::endl;
278  }
279  }
280 
281  txt<< "\\hline" <<std::endl;
282 
283  // Total uncertainty
284  txt<< "Systematic Uncertainty" ;
285  for(auto prof: profiles){
287 
288  TH1F* totplot = (TH1F*)Profile(expt, calc,
289  prof.fitvar,
290  steps, prof.minx, prof.maxx, 0,
291  prof.profvars,
292  totSysts);
293  totplot -> Write((prof.suffix+"_totError").c_str());
294 
295  range = GetLimits(totplot, hi, low);
296  if(prof.suffix=="dcp") low-=2;
297 
298  txt<< " & +"
299  << sqrt( fabs( (hi-prof.bf)*(hi-prof.bf) - (statup[prof.suffix]*statup[prof.suffix]) ) ) * pow(10,prof.oom)
300  << " / -"
301  << sqrt( fabs( (prof.bf-low)*(prof.bf-low) - (statdn[prof.suffix]*statdn[prof.suffix]) ) ) * pow(10,prof.oom) ;
302  }
303  txt<< " \\\\" <<std::endl;
304 
305  txt<< "Statistical Uncertainty" ;
306  for(auto prof: profiles){
307  txt<< " & +"
308  << statup[prof.suffix] * pow(10,prof.oom)
309  << " / -"
310  << statdn[prof.suffix] * pow(10,prof.oom) ;
311  }
313 
314  fOut -> Close();
315  delete fOut;
316 
317  std::string sumStr = full ? "full" : "short";
318 
319  for(auto prof: profiles){
320  TCanvas *c;
321  ErrorBarChart(bars[prof.suffix],{statdn[prof.suffix]*pow(10,prof.oom),statup[prof.suffix]*pow(10,prof.oom)},prof.root);
322  gPad->SetLeftMargin(0.33);
323  gPad->SetBottomMargin(0.15);
324  preliminary();
325  gPad->Print(("plots/syst_table_"+prof.suffix+"_"+sumStr+".pdf").c_str());
326  gPad->Print(("plots/syst_table_"+prof.suffix+"_"+sumStr+".C").c_str());
327  }
328 
329 }
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NueAcceptSystSignalKin2018RHC kNueAcceptSystSignalKin2018RHC
double ssth23
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
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)
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 AddJointAna2018BeamSysts(std::vector< const ISyst * > &systs, const EAnaType2018 ana)
std::vector< const IPrediction * > GetNumuPredictions2018(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const DummyAnaSyst kAnaRelativeCalibSyst("RelativeCalib","RelCalib")
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 double kAna2018RHCPOT
Definition: Exposures.h:208
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
const XML_Char const XML_Char * data
Definition: expat.h:268
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
TSpline3 hi("hi", xhi, yhi, 18,"0")
void preliminary()
expt
Definition: demo5.py:34
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
Definition: lz4.cxx:387
double dCP
root
Link up the nodes tree #####.
const NueAcceptSystBkg2018RHC kNueAcceptSystBkg2018RHC
RHC.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const NueAcceptSystSignalKin2018FHC kNueAcceptSystSignalKin2018FHC
double frac(double x)
Fractional part.
float bin[41]
Definition: plottest35.C:14
const NueAcceptSystBkg2018FHC kNueAcceptSystBkg2018FHC
FHC.
TLatex * prelim
Definition: Xsec_final.C:133
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
c1 Close()
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
Combine multiple component experiments.
const DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
double dmsq32
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
const MichelTaggingSyst2018 kMichelTaggingSyst2018
void syst_table_fit(bool full=false)
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)
Interface definition for fittable variables.
Definition: IFitVar.h:16
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const double kAna2018FHCPOT
Definition: Exposures.h:207
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
const DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
const DummyAnaSyst kAnaCalibShapeSyst("CalibShape","CalibShape")
T minus(const T &x)
Definition: minus.hpp:15
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()
enum BeamMode string