20 #include "Utilities/func/MathUtil.h" 25 #include "CAFAna/3flavor/Ana2020/joint_fit_2020_loader_tools.h" 29 const std::map<const ana::IFitVar*, std::pair<double, double>>
fitVarDrawRanges 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;
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;
64 bool drawFakeDataPredComparison=
true,
65 bool drawSurfaces=
true,
66 bool draw1Dmarginals=
true,
67 bool draw2Dmarginals=
true,
70 const auto TFText = [](
bool yes) {
return yes ?
"YES" :
"NO"; };
71 std::cout <<
"Drawing best fit - fake data comparison: " << TFText(drawFakeDataPredComparison) <<
std::endl;
81 std::unique_ptr<ana::MCMCSamples> warmup = std::move(fileContents.
warmup);
82 std::unique_ptr<ana::MCMCSamples> samples = std::move(fileContents.
samples);
84 std::map<std::string, ana::Spectrum> fakeDataSpectra =
mcmc::LoadFakeData(fakeDataFilename,
"spectra");
87 std::unique_ptr<ana::SystShifts> systTruePulls = std::move(fileContents.
trueSystPulls);
88 std::unique_ptr<osc::IOscCalcAdjustable> calcTruth = std::move(fileContents.
trueOscParams);
93 assert (samples && systTruePulls && calcTruth);
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));
111 if (!fakeDataSpectra.empty() && drawFakeDataPredComparison)
114 auto preds = mcmc::LoadPreds();
116 for (
const auto & predBundle : preds)
126 c.SaveAs(
mcmc::FullFilename(outDirPrefix,
"bestfitpred." + predBundle.name +
".systs.png").c_str());
143 TMarker marker(
util::sqr(
sin(calcTruth->GetTh23())), calcTruth->GetDmsq32() * 1e3, kFullStar);
144 marker.SetMarkerColor(
kBlue);
145 marker.SetMarkerSize(3);
147 c.SaveAs(
mcmc::FullFilename(outDirPrefix,
"surface_contour_ssth23-dm32.png").c_str());
157 marker = TMarker(calcTruth->GetdCP() / TMath::Pi(),
util::sqr(
sin(calcTruth->GetTh23())), kFullStar);
158 marker.SetMarkerColor(
kBlue);
159 marker.SetMarkerSize(3);
167 for (
const auto &fitVarPtr : samples->Vars())
174 double truthVal = fitVarPtr->GetValue(calcTruth.get());
175 TLine
l(truthVal,
h.GetMinimum(), truthVal,
h.GetMaximum());
176 l.SetLineColor(
kRed);
178 l.SetLineStyle(kDashed);
181 c.SaveAs(Form(
mcmc::FullFilename(outDirPrefix,
"marg_%s.systs.png").c_str(), fitVarPtr->ShortName().c_str()));
183 for (
const auto &systPtr : shifts->ActiveSysts())
190 double truthVal = systTruePulls->GetShift(systPtr);
191 TLine
l(truthVal,
h.GetMinimum(), truthVal,
h.GetMaximum());
192 l.SetLineColor(
kRed);
194 l.SetLineStyle(kDashed);
197 c.SaveAs(Form(
mcmc::FullFilename(outDirPrefix,
"marg_%s.systs.png").c_str(), systPtr->ShortName().c_str()));
204 for (
const auto fitVar1 : samples->Vars())
206 for (
auto itVar2 = samples->Vars().rbegin(); itVar2 != samples->Vars().rend(); itVar2++)
209 if (fitVar2 == fitVar1)
216 fitVar2, 100, drawRange2.first, drawRange2.second);
220 TMarker marker(fitVar1->GetValue(calcTruth.get()), fitVar2->
GetValue(calcTruth.get()), kFullStar);
221 marker.SetMarkerColor(kMagenta);
222 marker.SetMarkerSize(3);
224 c.SaveAs(Form(
mcmc::FullFilename(outDirPrefix,
"surface_contour_%s-%s.systs.png").c_str(),
225 fitVar1->ShortName().c_str(),
231 auto systs = shifts->ActiveSysts();
234 std::inserter(systNames, systNames.end()), [](
const ana::ISyst *
s)
237 for (
const auto &systName1 : systNames)
241 {
return s->
ShortName() == systName1; });
242 for (
auto itSystName2 = systNames.rbegin(); itSystName2 != systNames.rend(); itSystName2++)
244 if (*itSystName2 == systName1)
249 {
return s->
ShortName() == *itSystName2; });
253 samples->MaxValue(systPtr1),
254 systPtr2, 30, samples->MinValue(systPtr2),
255 samples->MaxValue(systPtr2));
260 TMarker marker(systTruePulls->GetShift(systPtr1), systTruePulls->GetShift(systPtr2), kFullStar);
261 marker.SetMarkerColor(kBlack);
262 marker.SetMarkerSize(3);
266 systPtr1->ShortName().c_str(),
268 c.SaveAs(outf.c_str());
273 if (drawPulls && !samples->Systs().empty())
277 TH2D h_fittedPulls(
"fitted_pulls",
";Systematic;Pull (#sigma)",
278 shifts->ActiveSysts().size(), 0, shifts->ActiveSysts().size(),
280 TH1D h_truePulls(
"true_pulls",
";Systematic;Pull (#sigma)", shifts->ActiveSysts().size(), 0,
281 shifts->ActiveSysts().size());
282 std::size_t systIdx = 0;
284 auto systsSorted = shifts->ActiveSysts();
285 std::sort(systsSorted.begin(),
287 [](
const auto syst1,
const auto syst2)
289 return syst1->ShortName() < syst2->ShortName();
291 for (
const auto &syst : systsSorted)
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));
299 h_truePulls.SetBinContent(systIdx, 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());
304 maxShift =
std::max(3., 1.2 * maxShift);
305 h_truePulls.GetYaxis()->SetRangeUser(-maxShift, maxShift);
306 h_fittedPulls.GetYaxis()->SetRangeUser(-maxShift, maxShift);
308 std::unique_ptr<TProfile> h_fittedPulls_prof(h_fittedPulls.ProfileX(
ana::UniqueName().c_str(), 1, -1,
"s"));
309 h_fittedPulls_prof->SetLineColor(
kBlue);
312 c.SetCanvasSize(900, 600);
313 TAxis
ax(*h_fittedPulls.GetXaxis());
314 h_fittedPulls.GetXaxis()->GetLabels()->Clear();
316 TPad upper(
"",
"", 0, 0, 1.0, 1);
318 h_fittedPulls.SetMaximum(0.25);
320 h_fittedPulls.SetFillColor(
kRed + 1);
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);
330 TPad lower(
"",
"", 0, 0, 1.0, 1);
331 lower.SetFillStyle(0);
336 for (
int bin = 1;
bin <= h_fittedPulls_prof->GetNbinsX();
bin++)
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);
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);
356 TLegend
leg(0.2, 0.9, 0.8, 1.0);
358 leg.SetBorderSize(0);
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);
T max(const caf::Proxy< T > &a, T b)
std::unique_ptr< osc::IOscCalcAdjustable > trueOscParams
Represent the binning of a Spectrum's x-axis.
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)
virtual const std::string & ShortName() const final
The name printed out to the screen.
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP
TH2 * QuantileSurface(Quantile quantile) const
Divide bin contents by bin widths.
std::unique_ptr< ana::MCMCSamples > samples
std::unique_ptr< ana::MCMCSamples > warmup
T sqr(T x)
More efficient square function than pow(x,2)
Encapsulate code to systematically shift a caf::SRProxy.
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
void Draw() const
Draw the surface itself.
const std::map< const ana::IFitVar *, std::pair< double, double > > fitVarDrawRanges
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23
TH2 * ToTH2(double minchi=-1) const
std::unique_ptr< ana::SystShifts > trueSystPulls
const std::string & ShortName() const
std::map< std::string, std::string > systNames
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
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)
void ResetCalculator(osc::IOscCalcAdjustable &calc)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
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.
FileContents LoadFromFile(const std::string &filename)
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0