plot_3flavor_withsysts.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TColor.h"
3 #include "TGraph.h"
4 #include "TH1D.h"
5 #include "TH2D.h"
6 #include "THashList.h"
7 #include "TLegend.h"
8 #include "TLine.h"
9 #include "TMarker.h"
10 #include "TMath.h"
11 #include "TProfile.h"
12 
13 #include "CAFAna/Analysis/Plots.h"
16 #include "CAFAna/Fit/MCMCSamples.h"
17 
18 #include "OscLib/OscCalcAnalytic.h"
19 
20 #include "Utilities/func/MathUtil.h"
21 
24 
25 #include "CAFAna/3flavor/Ana2020/joint_fit_2020_loader_tools.h"
26 
27 namespace mcmc
28 {
29  const std::map<const ana::IFitVar*, std::pair<double, double>> fitVarDrawRanges
30  {
32  {&ana::fitDmSq32Scaled_UniformPrior, {2.2, 2.8}},
34  };
35 
37  {
38  const int NRGBs = 9;
39  const int n_color_contours = 999;
40  static bool initialized=false;
41  static int* colors=new int[n_color_contours];
42  static int colmin = 0;
43 
44  if(!initialized){
45  // White -> red
46  Double_t stops[NRGBs] = { 0.00, 0.125, 0.250, 0.375, 0.500, 0.625, 0.750, 0.875, 1.000};
47  Double_t red[NRGBs] = { 1.00, 1.00, 0.99, 0.99, 0.98, 0.94, 0.80, 0.65, 0.40 };
48  Double_t green[NRGBs] = { 0.96, 0.88, 0.73, 0.57, 0.42, 0.23, 0.09, 0.06, 0.00 };
49  Double_t blue[NRGBs] = { 0.94, 0.82, 0.63, 0.45, 0.29, 0.17, 0.11, 0.08, 0.05 };
50  colmin=TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, n_color_contours);
51  for(uint i=0; i<n_color_contours; ++i) colors[i]=colmin+i;
52 
53  initialized=true;
54  }
55 
56  return colmin;
57  }
58 
59 }
60 
61 void plot_3flavor_withsysts(const std::string & samplesFilename,
62  const std::string & fakeDataFilename,
63  const std::string& outDirPrefix,
64  bool drawFakeDataPredComparison=true,
65  bool drawSurfaces=true,
66  bool draw1Dmarginals=true,
67  bool draw2Dmarginals=true,
68  bool drawPulls=true)
69 {
70  const auto TFText = [](bool yes) { return yes ? "YES" : "NO"; };
71  std::cout << "Drawing best fit - fake data comparison: " << TFText(drawFakeDataPredComparison) << std::endl;
72  std::cout << "Drawing osc param surfaces: " << TFText(drawSurfaces) << std::endl;
73  std::cout << "Drawing 1D marginals: " << TFText(draw1Dmarginals) << std::endl;
74  std::cout << "Drawing 2D marginals: " << TFText(draw2Dmarginals) << std::endl;
75  std::cout << "Drawing pulls: " << TFText(drawPulls) << std::endl;
76 
77  // load the systs so that they're in memory when we try to reload the systTruePulls below
79 
80  mcmc::FileContents fileContents = mcmc::LoadFromFile(samplesFilename);
81  std::unique_ptr<ana::MCMCSamples> warmup = std::move(fileContents.warmup);
82  std::unique_ptr<ana::MCMCSamples> samples = std::move(fileContents.samples);
83 
84  std::map<std::string, ana::Spectrum> fakeDataSpectra = mcmc::LoadFakeData(fakeDataFilename, "spectra");
85 
86  fileContents = mcmc::LoadFromFile(fakeDataFilename);
87  std::unique_ptr<ana::SystShifts> systTruePulls = std::move(fileContents.trueSystPulls);
88  std::unique_ptr<osc::IOscCalcAdjustable> calcTruth = std::move(fileContents.trueOscParams);
89 
90  std::cout << "samples: " << samples.get() << std::endl;
91  std::cout << "syst truth: " << systTruePulls.get() << std::endl;
92  std::cout << "calc truth: " << calcTruth.get() << std::endl;
93  assert (samples && systTruePulls && calcTruth);
94 
95  // todo: what to do here? really want the maximum of the joint marginal distribution over all the variables,
96  // but the single point with the best LL isn't a foolproof estimator of that.
97  auto bestFitIdx = samples->BestFitSampleIdx();
98  auto calc = std::make_unique<osc::OscCalcAnalytic>();
100  for (const auto & var : samples->Vars())
101  var->SetValue(calc.get(), samples->SampleValue(var, bestFitIdx));
102  auto shifts = std::make_unique<ana::SystShifts>();
103  for (const auto & syst : samples->Systs())
104  shifts->SetShift(syst, samples->SampleValue(syst, bestFitIdx));
105 
106  //auto bfDir = dynamic_cast<TDirectory*>(inf.Get("bestfit"));
107  //calc.reset(dynamic_cast<osc::OscCalcPMNSOpt*>(LoadFrom<osc::IOscCalc>(dynamic_cast<TDirectory*>(bfDir->Get("osc"))).release()));
108  //shifts = LoadFrom<SystShifts>(dynamic_cast<TDirectory*>(bfDir->Get("systs")));
109  //assert(shifts);
110 
111  if (!fakeDataSpectra.empty() && drawFakeDataPredComparison)
112  {
113  // draw a comparison of the fake data and the best fit prediction
114  auto preds = mcmc::LoadPreds();
115  TCanvas c;
116  for (const auto & predBundle : preds)
117  {
118  c.Clear();
119 
120  // calc we passed to the surface now contains best fit params
121  ana::DataMCComparison(fakeDataSpectra.at(predBundle.name),
122  predBundle.pred, // missing the cosmics for the moment...
123  calc.get(),
124  *shifts,
125  kBinDensity);
126  c.SaveAs(mcmc::FullFilename(outDirPrefix, "bestfitpred." + predBundle.name + ".systs.png").c_str());
127  }
128  }
129 
131  TCanvas c;
132 
133  if (drawSurfaces)
134  {
135  // draw the contours in oscillation parameter space overlaid on the LL surface.
136  // first, atmospheric
137  ana::BayesianSurface surf = samples->MarginalizeTo(&ana::fitSsqTh23_UniformPriorSsqTh23, 100, .3, .7,
138  &ana::fitDmSq32Scaled_UniformPrior, 100, 2.2, 2.8);
139  surf.Draw();
140  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian1Sigma), 7, kGreen + 2); // dashed
141  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian2Sigma), 7, kBlack); // dashed
142  surf.DrawBestFit(kGray);
143  TMarker marker(util::sqr(sin(calcTruth->GetTh23())), calcTruth->GetDmsq32() * 1e3, kFullStar);
144  marker.SetMarkerColor(kBlue);
145  marker.SetMarkerSize(3);
146  marker.Draw();
147  c.SaveAs(mcmc::FullFilename(outDirPrefix, "surface_contour_ssth23-dm32.png").c_str());
148 
149  // now the 'reactor' contour
150  c.Clear();
151  surf = samples->MarginalizeTo(&ana::fitDeltaInPiUnits_UniformPriordCP, 75, 0, 2,
152  &ana::fitSsqTh23_UniformPriorSsqTh23, 75, 0.3, 0.7);
153  surf.Draw();
154  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian1Sigma), 7, kGreen + 2); // dashed
155  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian2Sigma), 7, kBlack); // dashed
156  surf.DrawBestFit(kGray);
157  marker = TMarker(calcTruth->GetdCP() / TMath::Pi(), util::sqr(sin(calcTruth->GetTh23())), kFullStar);
158  marker.SetMarkerColor(kBlue);
159  marker.SetMarkerSize(3);
160  marker.Draw();
161  c.SaveAs(mcmc::FullFilename(outDirPrefix, "surface_contour_ssth23-dcp.png").c_str());
162  }
163 
164  if (draw1Dmarginals)
165  {
166  // draw the marginals and the pulls
167  for (const auto &fitVarPtr : samples->Vars())
168  {
169  c.Clear();
170  ana::Bayesian1DMarginal marg = samples->MarginalizeTo(fitVarPtr);
171  auto h = marg.ToTH1(ana::Binning::Simple(50, samples->MinValue(fitVarPtr), samples->MaxValue(fitVarPtr)));
172  h.Draw("hist");
173 
174  double truthVal = fitVarPtr->GetValue(calcTruth.get());
175  TLine l(truthVal, h.GetMinimum(), truthVal, h.GetMaximum());
176  l.SetLineColor(kRed);
177  l.SetLineWidth(2);
178  l.SetLineStyle(kDashed);
179  l.Draw();
180 
181  c.SaveAs(Form(mcmc::FullFilename(outDirPrefix, "marg_%s.systs.png").c_str(), fitVarPtr->ShortName().c_str()));
182  }
183  for (const auto &systPtr : shifts->ActiveSysts())
184  {
185  c.Clear();
186  ana::Bayesian1DMarginal marg = samples->MarginalizeTo(systPtr);
187  auto h = marg.ToTH1(ana::Binning::Simple(50, samples->MinValue(systPtr), samples->MaxValue(systPtr)));
188  h.Draw("hist");
189 
190  double truthVal = systTruePulls->GetShift(systPtr);
191  TLine l(truthVal, h.GetMinimum(), truthVal, h.GetMaximum());
192  l.SetLineColor(kRed);
193  l.SetLineWidth(2);
194  l.SetLineStyle(kDashed);
195  l.Draw();
196 
197  c.SaveAs(Form(mcmc::FullFilename(outDirPrefix, "marg_%s.systs.png").c_str(), systPtr->ShortName().c_str()));
198  }
199  }
200 
201  if (draw2Dmarginals)
202  {
203  // pairs of fit vars
204  for (const auto fitVar1 : samples->Vars())
205  {
206  for (auto itVar2 = samples->Vars().rbegin(); itVar2 != samples->Vars().rend(); itVar2++)
207  {
208  const ana::IFitVar *fitVar2 = *itVar2;
209  if (fitVar2 == fitVar1)
210  break;
211 
212  c.Clear();
213  std::pair<double, double> drawRange1 = mcmc::fitVarDrawRanges.at(fitVar1);
214  std::pair<double, double> drawRange2 = mcmc::fitVarDrawRanges.at(fitVar2);
215  ana::BayesianSurface surf = samples->MarginalizeTo(fitVar1, 100, drawRange1.first, drawRange1.second,
216  fitVar2, 100, drawRange2.first, drawRange2.second);
217  surf.Draw();
218  surf.DrawContour(surf.QuantileSurface(ana::Quantile::kGaussian1Sigma), 7, kGreen + 2); // dashed
219  surf.DrawBestFit(kGray);
220  TMarker marker(fitVar1->GetValue(calcTruth.get()), fitVar2->GetValue(calcTruth.get()), kFullStar);
221  marker.SetMarkerColor(kMagenta);
222  marker.SetMarkerSize(3);
223  marker.Draw();
224  c.SaveAs(Form(mcmc::FullFilename(outDirPrefix, "surface_contour_%s-%s.systs.png").c_str(),
225  fitVar1->ShortName().c_str(),
226  fitVar2->ShortName().c_str()));
227  }
228  }
229 
230  // 2D syst space
231  auto systs = shifts->ActiveSysts();
232  std::set<std::string> systNames;
233  std::transform(systs.begin(), systs.end(),
234  std::inserter(systNames, systNames.end()), [](const ana::ISyst *s)
235  { return s->ShortName(); });
236  assert(systNames.size() == systs.size());
237  for (const auto &systName1 : systNames)
238  {
239  auto systPtr1 = *std::find_if(systs.begin(), systs.end(),
240  [&systName1](const ana::ISyst *s)
241  { return s->ShortName() == systName1; });
242  for (auto itSystName2 = systNames.rbegin(); itSystName2 != systNames.rend(); itSystName2++)
243  {
244  if (*itSystName2 == systName1)
245  break;
246 
247  const ana::ISyst *systPtr2 = *std::find_if(systs.begin(), systs.end(),
248  [&itSystName2](const ana::ISyst *s)
249  { return s->ShortName() == *itSystName2; });
250 
251  c.Clear();
252  ana::BayesianSurface marg = samples->MarginalizeTo(systPtr1, 30, samples->MinValue(systPtr1),
253  samples->MaxValue(systPtr1),
254  systPtr2, 30, samples->MinValue(systPtr2),
255  samples->MaxValue(systPtr2));
256  auto h2 = marg.ToTH2();
257  h2->Draw("colz");
258  marg.DrawContour(marg.QuantileSurface(ana::Quantile::kGaussian1Sigma), 7, kGreen + 2); // dashed
259 
260  TMarker marker(systTruePulls->GetShift(systPtr1), systTruePulls->GetShift(systPtr2), kFullStar);
261  marker.SetMarkerColor(kBlack);
262  marker.SetMarkerSize(3);
263  marker.Draw();
264 
265  std::string outf = Form(mcmc::FullFilename(outDirPrefix, "marg_%s-%s.allSysts.png").c_str(),
266  systPtr1->ShortName().c_str(),
267  systPtr2->ShortName().c_str());
268  c.SaveAs(outf.c_str());
269  }
270  }
271  }
272 
273  if (drawPulls && !samples->Systs().empty())
274  { // show that the pulls were what we expect
275  const ana::Binning margBins = ana::Binning::Simple(100, -2.95,
276  2.95); // then labels don't collide in the ratio plot
277  TH2D h_fittedPulls("fitted_pulls", ";Systematic;Pull (#sigma)",
278  shifts->ActiveSysts().size(), 0, shifts->ActiveSysts().size(),
279  margBins.NBins(), margBins.Min(), margBins.Max());
280  TH1D h_truePulls("true_pulls", ";Systematic;Pull (#sigma)", shifts->ActiveSysts().size(), 0,
281  shifts->ActiveSysts().size());
282  std::size_t systIdx = 0;
283  double maxShift = 0;
284  auto systsSorted = shifts->ActiveSysts();
285  std::sort(systsSorted.begin(),
286  systsSorted.end(),
287  [](const auto syst1, const auto syst2)
288  {
289  return syst1->ShortName() < syst2->ShortName();
290  });
291  for (const auto &syst : systsSorted)
292  {
293  ++systIdx;
294  auto systMarginal = samples->MarginalizeTo(syst);
295  TH1D h_systMarginal = systMarginal.ToTH1(margBins);
296  for (int binIdx = 0; binIdx <= h_systMarginal.GetNbinsX(); binIdx++)
297  h_fittedPulls.SetBinContent(systIdx, binIdx, h_systMarginal.GetBinContent(binIdx));
298  maxShift = std::max(maxShift, std::abs(shifts->GetShift(syst)));
299  h_truePulls.SetBinContent(systIdx, systTruePulls->GetShift(syst));
300  maxShift = std::max(maxShift, std::abs(systTruePulls->GetShift(syst)));
301  for (auto &h : {dynamic_cast<TH1 *>(&h_fittedPulls), dynamic_cast<TH1 *>(&h_truePulls)})
302  h->GetXaxis()->SetBinLabel(systIdx, syst->ShortName().c_str());
303  }
304  maxShift = std::max(3., 1.2 * maxShift);
305  h_truePulls.GetYaxis()->SetRangeUser(-maxShift, maxShift);
306  h_fittedPulls.GetYaxis()->SetRangeUser(-maxShift, maxShift);
307 
308  std::unique_ptr<TProfile> h_fittedPulls_prof(h_fittedPulls.ProfileX(ana::UniqueName().c_str(), 1, -1, "s"));
309  h_fittedPulls_prof->SetLineColor(kBlue);
310 
311  c.Clear();
312  c.SetCanvasSize(900, 600);
313  TAxis ax(*h_fittedPulls.GetXaxis());
314  h_fittedPulls.GetXaxis()->GetLabels()->Clear();
315 
316  TPad upper("", "", 0, 0, 1.0, 1);
317  upper.cd();
318  h_fittedPulls.SetMaximum(0.25); // so that more of the color winds up red
320  h_fittedPulls.SetFillColor(kRed + 1); // for the legend
321  h_fittedPulls.GetZaxis()->SetTitle("Fraction of MCMC samples");
322  h_fittedPulls.GetZaxis()->RotateTitle();
323  h_fittedPulls.Draw("colz");
324  h_truePulls.Draw("hist x same");
325  h_fittedPulls_prof->Draw("pe same");
326  upper.SetLeftMargin(0.15);
327  upper.SetRightMargin(0.15);
328  upper.SetBottomMargin(0.4);
329 
330  TPad lower("", "", 0, 0, 1.0, 1);
331  lower.SetFillStyle(0);
332  lower.cd();
333 
334  // a "pull" plot
335  TGraph fitPulls;
336  for (int bin = 1; bin <= h_fittedPulls_prof->GetNbinsX(); bin++)
337  {
338  double pull = 0;
339  if (h_fittedPulls_prof->GetBinError(bin) != 0)
340  pull = (h_fittedPulls_prof->GetBinContent(bin) - h_truePulls.GetBinContent(bin))
341  / h_fittedPulls_prof->GetBinError(bin);
342  fitPulls.SetPoint(fitPulls.GetN(), double(bin) - 0.5, pull);
343  }
344  fitPulls.Draw("al");
345  *fitPulls.GetXaxis() = ax;
346  fitPulls.SetTitle(";;#frac{True - Fit mean}{Fit RMS}");
347  fitPulls.GetYaxis()->SetRangeUser(-0.95, 0.95);
348  lower.SetLeftMargin(0.15);
349  lower.SetRightMargin(0.15);
350  lower.SetTopMargin(0.6);
351  lower.SetBottomMargin(0.3);
352 
353  c.cd();
354  upper.Draw();
355  lower.Draw();
356  TLegend leg(0.2, 0.9, 0.8, 1.0);
357  leg.SetFillStyle(0);
358  leg.SetBorderSize(0);
359  leg.SetNColumns(3);
360  leg.AddEntry(&h_truePulls, "True", "l");
361  leg.AddEntry(&h_fittedPulls, "MCMC samples", "f");
362  leg.AddEntry(h_fittedPulls_prof.get(), "MCMC (mean+RMS)", "lpe");
363  leg.SetColumnSeparation(0.15);
364  leg.Draw();
365 
366 // c.SetBottomMargin(0.5);
367  c.SaveAs(mcmc::FullFilename(outDirPrefix, "pulls.systs.png").c_str());
368  }
369 }
T max(const caf::Proxy< T > &a, T b)
std::unique_ptr< osc::IOscCalcAdjustable > trueOscParams
Definition: MCMC3FShared.h:40
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
void plot_3flavor_withsysts(const std::string &samplesFilename, const std::string &fakeDataFilename, const std::string &outDirPrefix, bool drawFakeDataPredComparison=true, bool drawSurfaces=true, bool draw1Dmarginals=true, bool draw2Dmarginals=true, bool drawPulls=true)
std::string FullFilename(const std::string &dir, std::string file)
Definition: MCMC3FShared.h:53
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP
TH2 * QuantileSurface(Quantile quantile) const
int GetRedHeatPalette()
Divide bin contents by bin widths.
Definition: UtilsExt.h:31
std::unique_ptr< ana::MCMCSamples > samples
Definition: MCMC3FShared.h:38
std::unique_ptr< ana::MCMCSamples > warmup
Definition: MCMC3FShared.h:37
osc::OscCalcDumb calc
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
float abs(float number)
Definition: d0nt_math.hpp:39
const XML_Char * s
Definition: expat.h:262
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
int colors[6]
Definition: tools.h:1
double Min() const
Definition: Binning.h:30
tuple blue
Definition: rootlogon.py:65
int NRGBs
Definition: rootlogon.py:77
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
const std::map< const ana::IFitVar *, std::pair< double, double > > fitVarDrawRanges
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior
TFile * outf
Definition: testXsec.C:51
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23
double Max() const
Definition: Binning.h:31
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T sin(T number)
Definition: d0nt_math.hpp:132
std::unique_ptr< ana::SystShifts > trueSystPulls
Definition: MCMC3FShared.h:39
const std::string & ShortName() const
Definition: IFitVar.h:36
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
int NBins() const
Definition: Binning.h:29
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
TH1D ToTH1(const Binning &bins) const
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
Definition: Plots.cxx:35
void ResetCalculator(osc::IOscCalcAdjustable &calc)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
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)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
std::vector< ana::predictions > LoadPreds(bool extrap)
Just wrap up all the args somewhere centralized for the moment.
FileContents LoadFromFile(const std::string &filename)
surf
Definition: demo5.py:35
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0
unsigned int uint