20 #include "../joint_fit_2020_loader_tools.h" 21 #include "../joint_fit_2020_datarelease_tools.h" 46 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
47 const Int_t kFillColor1,
const Int_t kFillColor2,
const Int_t kFillColor3,
48 const Int_t kDarkColor,
const TString surftype);
55 Double_t
white[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
56 Double_t Red[9] = { 61./255., 99./255., 136./255., 181./255., 213./255., 225./255., 198./255., 99./255., 61./255.};
57 Double_t Green[9] = { 149./255., 140./255., 96./255., 83./255., 132./255., 178./255., 190./255., 140./255., 149./255};
58 Double_t Blue[9] = { 214./255., 203./255., 168./255., 135./255., 110./255., 100./255., 111./255., 203./255., 214./255.};
59 Double_t Length[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
61 TColor::CreateGradientColorTable(9, Length, Red, Green, Blue, 255, 1.0);
62 gStyle->SetNumberContours(255);
67 bool runOnGrid =
false,
68 TString
options=
"joint_fake2019_both_onlyIH",
69 bool dmsqSurf =
false,
70 bool th13Surf =
false,
73 bool fillContour =
false,
75 bool fc_overlay =
false)
78 bool nueOnly =
options.Contains(
"nueOnly");
79 bool numuOnly =
options.Contains(
"numuOnly");
80 bool joint =
options.Contains(
"joint");
81 assert (nueOnly || numuOnly || joint);
83 bool FHCOnly =
options.Contains(
"FHCOnly");
84 bool RHCOnly =
options.Contains(
"RHCOnly");
85 bool both =
options.Contains(
"both");
86 assert (FHCOnly || RHCOnly || both);
88 bool fake2019 =
options.Contains(
"fake2019");
89 bool realData =
options.Contains(
"realData");
91 bool onlyNH =
options.Contains(
"onlyNH");
92 bool onlyIH =
options.Contains(
"onlyIH");
95 if(dmsqSurf) suffix+=
"_dmsq";
96 if(th13Surf) suffix+=
"_th13";
97 if(corrSysts) suffix+=
"_systs";
98 else suffix+=
"_stats";
100 TString
outdir =
"/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/";
101 if (runOnGrid) outdir =
"./";
102 TString outFCfilename (outdir +
"contours_FCInput_2020_" + suffix);
103 TString
outfilename (outdir +
"hist_contours_2020_" + suffix);
115 TString infilename =
outdir;
116 infilename +=
"hist_contours_2020_";
118 if (dmsqSurf) infilename +=
"_dmsq";
119 infilename += corrSysts?
"_systs.root":
"_stats.root";
121 TFile *
infile =
new TFile (infilename,
"read");
123 TString officialpreffix =
"./";
124 bool makedatarelease =
true;
126 if(!makedatarelease && corrSysts &&
options.Contains(
"joint") && fccorr && 0) {
127 official =
new TFile(
"NOvA_nue_official_results.root",
"recreate");
128 TH2*
axes =
new TH2F(
"axes",
";#delta_{CP}/#pi;sin^{2}#theta_{23}",100,0,2,100,0.225,0.725);
133 if (makedatarelease){
138 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
139 double minchi23 = mins[0];
141 std::vector <TString> surfNames;
142 if(!th13Surf && ! dmsqSurf)surfNames = {
"delta_th23_"};
143 if (dmsqSurf) surfNames.push_back(
"th23_dm32_");
144 if (th13Surf) surfNames.push_back(
"th13_delta_");
146 for(TString surfName:surfNames){
147 for (
int hie:{-1,+1}){
149 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
152 if(surfName.Contains(
"delta_th23")) {
153 if (zoomIn)
DrawContourCanvas(surfName, 0., 2., (RHCOnly?0:0.275), (RHCOnly?1.0:0.725));
156 if(surfName.Contains(
"th23_dm32")) {
157 if (hie>0)
DrawContourCanvas(surfName, ((RHCOnly)?0.2:0.3), ((RHCOnly)?0.8:0.7), (RHCOnly?2.0:2.075), (RHCOnly?3.3:2.725));
158 else DrawContourCanvas(surfName, ((RHCOnly)?0.2:0.3), ((RHCOnly)?0.8:0.7), (RHCOnly?-3.0:-2.8), (RHCOnly?-2.0:-2.2));
160 if(surfName.Contains(
"th13_delta")) {
166 TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?
"NH":
"IH"));
170 TH2 * fc_surf_1Sigma =
new TH2F();
171 TH2 * fc_surf_2Sigma =
new TH2F();
172 TH2 * fc_surf_3Sigma =
new TH2F();
173 TH2 * fc_surf_90CL =
new TH2F();
178 TString fcFolder=
"/nova/ana/nu_e_ana/Ana2020/FC/";
181 if(dmsqSurf) fcFileName =
"contours/ssth23dmsq/ssth23dmsq_";
182 else if(!dmsqSurf && !th13Surf)
183 fcFileName =
"contours/deltassth23/deltassth23_";
184 fcFileName += hie > 0?
"nh":
"ih";
185 fcFileName +=
".root";
189 int NX = hsurf->GetXaxis()->GetNbins();
190 int NY = hsurf->GetYaxis()->GetNbins();
191 double minX = hsurf->GetXaxis()->GetBinLowEdge(1);
192 double maxX = hsurf->GetXaxis()->GetBinUpEdge(NX);
193 double minY = hsurf->GetYaxis()->GetBinLowEdge(1);
194 double maxY = hsurf->GetYaxis()->GetBinUpEdge(NY);
197 std::cout <<
"-------------------------------\n";
204 std::cout <<
"\n------------------------------\n";
206 fcXH =
new FCSurface(NX, minX, maxX, NY, minY, maxY);
208 if(dmsqSurf) plot_name =
"ssth23dmsq32";
209 if(th13Surf) plot_name =
"deltath13";
211 fcXH->
Add(*fccol, plot_name);
217 fc_surf_1Sigma->Smooth();
218 fc_surf_2Sigma->Smooth();
219 fc_surf_90CL->Smooth();
220 if(dmsqSurf || th13Surf)
221 fc_surf_3Sigma->Smooth();
226 Color_t k2SigmaConfidenceColor = ((hie>0)?
cnh:
cih);
227 Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.99);
228 Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.5);
229 Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.2);
230 Int_t kFCFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.65: 0.75);
231 Int_t kFCFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.4: 0.38);
232 Int_t kFCFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, fc_overlay? 0.18: 0.18);
233 Int_t kDarkColor = kGray+3;
237 auto g1 =
surf.GetGraphs(surf_1Sigma,minchi23);
238 auto g2 =
surf.GetGraphs(surf_2Sigma,minchi23);
239 auto g3 =
surf.GetGraphs(surf_3Sigma,minchi23);
241 kFillColor1, kFillColor2, kFillColor3,
242 kDarkColor, suffix + surfName + (hie>0?
"NH":
"IH"));
245 auto fc_g1 =
surf.GetGraphs(fc_surf_1Sigma,minchi23);
246 auto fc_g2 =
surf.GetGraphs(fc_surf_2Sigma,minchi23);
247 auto fc_g3 =
surf.GetGraphs(fc_surf_3Sigma,minchi23);
249 kFillColor1, kFillColor2, kFillColor3,
250 kDarkColor, suffix + surfName + (hie>0?
"NH":
"IH")+
"fccorr");
254 cout<<
"minchi from file is "<<minchi23<<
endl;
255 if(!fccorr || fc_overlay){
258 auto g1 =
surf.GetGraphs(surf_1Sigma,minchi23);
259 auto g2 =
surf.GetGraphs(surf_2Sigma,minchi23);
260 auto g3 =
surf.GetGraphs(surf_3Sigma,minchi23);
262 kFillColor1, kFillColor2, kFillColor3,
263 kDarkColor, suffix + surfName + (hie>0?
"NH":
"IH"));
265 surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
267 surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
270 surf.DrawContour(fc_surf_1Sigma, kSolid, kFCFillColor1, minchi23);
272 surf.DrawContour(fc_surf_3Sigma, kSolid, kFCFillColor3, minchi23);
281 if(bestfit && both && corrSysts && hie>0){
282 TFile * bffile =
new TFile(
"/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/bestfits_joint_realData_both_systs.root",
"READ");
283 auto calc_tmp = LoadFrom<osc::IOscCalc>(bffile,
"NH_UO_calc").
release();
291 else{delta = 0.678516; th13 = 0.0957894;
cout<<
"NB! Old 2019 values"<<
endl;}
292 bf.SetMarkerColor(kDarkColor);
294 bf.SetMarkerStyle(43);
296 if(!dmsqSurf && ! th13Surf){
313 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
318 leg->SetFillStyle(1001);
320 leg->SetX1(0.7);
leg->SetX2(0.8);
321 leg->SetY1(0.4);
leg->SetY2(0.85);
322 if(!fccorr)
leg->SetHeader(
"No FC");
324 TGraph *
dummy =
new TGraph;
325 dummy->SetLineWidth(2);
326 dummy->SetLineColor(kGray+2);
327 dummy->SetFillColor(kGray);
328 leg->AddEntry(dummy,
"#splitline{Reactor}{68%C.L.}",
"f");
335 for(TString
ext: {
".png",
".pdf",
".root"}){
336 gPad->Print(
"./contour_" + suffix +
"_" +
338 (hie>0?
"NH":
"IH") + (fccorr?
"_fccorr":
"") +
339 ((fccorr && smooth)?
"_smooth":
"") +
340 (zoomIn?
"_zoom":
"") +
341 (fc_overlay?
"_overlay":
"") +
346 auto fc_g1 =
surf.GetGraphs(fc_surf_1Sigma,minchi23);
347 auto fc_g2 =
surf.GetGraphs(fc_surf_2Sigma,minchi23);
348 auto fc_g3 =
surf.GetGraphs(fc_surf_3Sigma,minchi23);
349 auto fc_g9 =
surf.GetGraphs(fc_surf_90CL,minchi23);
359 if(debug && !fccorr){
361 int nSysts = corrSysts ? 70 : 0;
362 for (
int i = 0;
i < nSysts; ++
i){
363 hists.push_back((TH2*) infile->Get((TString)
"surf_" + surfName +
364 (hie > 0 ?
"NH" :
"IH") +
371 if(TString(
hists[i]->GetTitle()).Contains(
"delta")) {
373 hists[
i]->GetZaxis()->SetRangeUser(0,2);
375 else gStyle->SetPalette(87);
377 double minz =
hists[
i]->GetBinContent(
hists[i]->GetMinimumBin());
378 double maxz =
hists[
i]->GetBinContent(
hists[i]->GetMaximumBin());
380 hists[
i]->GetZaxis()->SetRangeUser(-lim,+lim);
382 if(i==(
hists.size()-1)) {
383 gStyle->SetPalette(56);
384 hists[
i]->GetZaxis()->SetRangeUser(0.,11.83);
388 hists[
i]->GetYaxis()->SetTitleOffset(0.8);
389 gPad->SetRightMargin(0.14);
390 surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
391 surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
392 surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
393 gPad->Print(
"./debug_contour_" + suffix +
"_" +
395 (hie > 0 ?
"NH":
"IH") + (fccorr?
"_fccorr":
"") +
396 (zoomIn ?
"_zoom" :
"") +
408 TLatex * ltxNH =
new TLatex(0.16,0.89,
"#bf{NH}");
409 ltxNH->SetTextColor(
cnh);
411 ltxNH->SetTextAlign(22);
413 if(th13) ltxNH->Draw();
414 else if(!vert) ltxNH->DrawLatex(.85,.20,
"#bf{NH}");
415 else ltxNH->DrawLatex(.85,.85,
"#bf{NH}");
418 TLatex * ltxIH =
new TLatex (0.16,0.46,
"#bf{IH}");
421 ltxIH->SetTextAlign(22);
423 if(th13) ltxIH->Draw();
424 else if(!vert) ltxIH->DrawLatex(.85,.20,
"#bf{IH}");
425 else ltxIH->DrawLatex(.85,.85,
"#bf{IH}");
430 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
431 const Int_t kFillColor1,
const Int_t kFillColor2,
const Int_t kFillColor3,
432 const Int_t kDarkColor,
const TString surftype)
434 bool isJoint = surftype.Contains(
"joint");
435 bool isBoth = surftype.Contains(
"both");
436 bool isRHC = surftype.Contains(
"RHCOnly");
437 bool isNH = surftype.Contains(
"NH");
438 bool isSysts = surftype.Contains(
"systs");
439 bool isFccorr = surftype.Contains(
"fccorr");
441 if (surftype.Contains(
"delta_th23")){
443 if(isJoint && isNH ){
446 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
449 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f same");
450 JoinGraphs(g2[1], g2[2], kFillColor2)->Draw(
"f same");
453 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
456 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f");
461 for (
auto &
g:g1) {
g->SetFillColor(kWhite);
g->DrawClone(
"c f");}
462 for (
auto &
g:g1) {
g->SetFillColor(kFillColor1);
g->Draw(
"c f");}
463 for (
auto &
g:g1) {
g->Draw(
"c");}
466 JoinGraphs(g1[0], g1[1], kFillColor1)->Draw(
"c f");
472 if(isJoint && !isNH){
473 if(isRHC)
JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
476 g3[1]->SetFillColor(kFillColor3); g3[1]->Draw(
"f c");
477 g3[0]->SetFillColor(kFillColor3); g3[0]->Draw(
"f c");
484 g->SetFillColor(kWhite);
g->DrawClone(
"f c");
485 g->SetFillColor(kFillColor2);
g->Draw(
"f c");
486 g->SetLineColor(kDarkColor);
491 g->SetFillColor(kWhite);
g->DrawClone(
"f c");
492 g->SetFillColor(kFillColor1);
g->Draw(
"f c");
493 g->SetLineColor(kDarkColor);
497 for(
auto &
gs:{g3, g2, g1}){
499 g->SetLineColor(kDarkColor);
507 if (surftype.Contains(
"th23_dm32")){
509 g->SetFillColor(kFillColor3);
g->Draw(
"cf");
512 g->SetFillColor(kWhite);
g->DrawClone(
"cf");
513 g->SetFillColor(kFillColor2);
g->Draw(
"cf");
516 g->SetFillColor(kWhite);
g->DrawClone(
"cf");
517 g->SetFillColor(kFillColor1);
g->Draw(
"cf");
519 for (
auto &
gs:{g3,g2,g1}){
521 g->SetLineColor(kDarkColor);
529 if (surftype.Contains(
"th13_delta")) {
531 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
534 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f");
537 for (
auto &
g:g2) {
g->SetFillColor(kWhite);
g->DrawClone(
"f");}
538 for (
auto &
g:g2) {
g->SetFillColor(kFillColor2);
g->Draw(
"f");}
540 for (
auto &
g:g1) {
g->SetFillColor(kWhite);
g->DrawClone(
"f");}
541 for (
auto &
g:g1) {
g->SetFillColor(kFillColor1);
g->Draw(
"f");}
547 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f c");
548 JoinGraphs(g2[0], g2[1], kWhite)->Draw(
"f c");
549 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f c");
554 JoinGraphs(g1[0], g1[1], kWhite)->Draw(
"f c");
555 JoinGraphs(g1[0], g1[1], kFillColor1)->Draw(
"f c");
558 for (
auto &
gs:{g3,g2,g1}){
560 g->SetLineColor(kDarkColor);
571 TBox* box =
new TBox(mean-err, 0, mean+err, 2);
572 box->SetFillColorAlpha(kGray,0.7);
574 TLine* linLo =
new TLine(mean-err, 0, mean-err, 2);
575 linLo->SetLineColor(kGray+2);
576 linLo->SetLineWidth(2);
578 TLine* linHi =
new TLine(mean+err, 0, mean+err, 2);
579 linHi->SetLineColor(kGray+2);
580 linHi->SetLineWidth(2);
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
fvar< T > fabs(const fvar< T > &x)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var's GetValue()
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
void DrawHieTag(int hie, bool th13=true, bool vert=false)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void SetPaletteDeltaCyclicNew()
string outfilename
knobs that need extra care
void FillGraphs(std::vector< TGraph * > g1, std::vector< TGraph * > g2, std::vector< TGraph * > g3, const Int_t kFillColor1, const Int_t kFillColor2, const Int_t kFillColor3, const Int_t kDarkColor, const TString surftype)
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
void plot_joint_fit_2020_contours(bool corrSysts=false, bool runOnGrid=false, TString options="joint_fake2019_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false, bool fccorr=false, bool zoomIn=true, bool fillContour=false, bool smooth=true, bool fc_overlay=false)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void Add(const FCPoint &pt, std::string plot)
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
const Style_t k2SigmaConfidenceStyle
void DrawReactorConstraint(double mean, double err)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
const Float_t kBlessedLabelFontSize
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
TH2 * SurfaceForSignificance(double sig)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
virtual IOscCalc * Copy() const override
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
Pseudo-experiments, binned to match a log-likelihood surface.