43 TLegend *
ContourLegend(
int hie,
bool fillContour,
bool fccorr,
44 Int_t kFillColor1, Int_t kFillColor2,
45 Int_t kFillColor3, Int_t kDarkColor,
47 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
48 const Int_t kFillColor1,
const Int_t kFillColor2,
const Int_t kFillColor3,
49 const Int_t kDarkColor,
const TString surftype);
52 double maxx = 2,
double miny = 0,
double maxy = 1){
55 c1->SetLeftMargin(0.14);
56 c1->SetBottomMargin(0.15);
57 TH2*
axes =
new TH2F();
59 if(surfName.Contains(
"delta_th23")) title=
";#delta_{CP};sin^{2}#theta_{23}";
60 if(surfName.Contains(
"th13_delta")) title=
";sin^{2}2#theta_{13};#delta_{CP}/#pi";
61 if(surfName.Contains(
"th23_dm32")) title=
";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
62 axes =
new TH2F(
"",title,100,minx,maxx,100,miny,
maxy);
70 axes->GetXaxis()->SetTitleOffset(0.8);
71 axes->GetYaxis()->SetTitleOffset(0.8);
72 TGaxis::SetMaxDigits(3);
73 if(surfName.Contains(
"th23_dm32")) {
74 axes->GetYaxis()->SetDecimals() ;
75 axes->GetYaxis()->SetTitleOffset(0.85);
76 axes->GetYaxis()->SetLabelOffset(0.002);
77 axes->GetYaxis()->SetLabelSize(0.058);
78 axes->GetYaxis()->SetTitleSize(0.078);
85 bool createFile=
false,
86 bool corrSysts =
false,
87 TString
options=
"joint_fakeNHLO",
88 bool dmsqSurf =
false,
89 bool th13Surf =
false,
92 bool fillContour =
false,
96 bool nueOnly =
options.Contains(
"nueOnly");
97 bool numuOnly =
options.Contains(
"numuOnly");
98 bool joint =
options.Contains(
"joint");
99 assert (nueOnly || numuOnly || joint);
101 bool fakeNHLO =
options.Contains(
"fakeNHLO");
102 bool fakeNHUO =
options.Contains(
"fakeNHUO");
103 bool fake2017 =
options.Contains(
"fake2017");
104 bool realData =
options.Contains(
"realData");
106 bool onlyNH =
options.Contains(
"onlyNH");
107 bool onlyIH =
options.Contains(
"onlyIH");
111 if(dmsqSurf) suffix+=
"_dmsq";
112 if(th13Surf) suffix+=
"_th13";
113 if(corrSysts) suffix+=
"_systs";
114 else suffix+=
"_stats";
116 TString
outdir=
"/nova/ana/nu_e_ana/Ana2017/Results/contours/";
117 if(th13Surf) outdir +=
"th13_delta/";
118 else if (dmsqSurf) outdir +=
"th23_dmsq/";
119 else outdir +=
"delta_th23/";
123 TString
outfilename (outdir +
"hist_contours_2017_" + suffix);
124 TString outFCfilename (outdir +
"contours_FCInput_2017_" +
decomp +
".root");
132 std::vector <const IPrediction * > preds;
133 std::vector <std::pair <TH1D*, double > >
cosmics;
134 std::vector <Spectrum * >
data;
135 std::vector <const IExperiment * > expts;
138 if(fakeNHLO)
SetFakeCalc(calc_fake, 0.404, 2.7
e-3, 1.48);
139 else if(fakeNHUO)
SetFakeCalc(calc_fake, 0.623, 2.7
e-3, 0.74);
140 else if(fake2017)
SetFakeCalc(calc_fake, 0.56, 2.43
e-3, 1.21);
141 else if(!realData) {
std::cerr <<
"need setting for data\n";
exit(1);}
150 preds.insert(preds.end(),numu_preds.begin(), numu_preds.end());
152 cosmics.insert(cosmics.end(),numu_cosmics.begin(), numu_cosmics.end());
160 data.insert(data.end(),numu_data.begin(), numu_data.end());
163 for(
int i = 0;
i <
int(preds.size()); ++
i){
166 if(realData) thisdata = data[
i];
168 cosmics[i].first,cosmics[i].
second));
173 cosmics[
i].first->SetMarkerStyle(34);
175 cosmics[
i].first->Draw(
"hist p same");
176 gPad->Print(outdir +
"debug_predictions_" + suffix +
"_" +
std::to_string(i) +
".pdf");
185 std::cout <<
"Adding WorldReactorConstraint2017\n";
189 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
192 std::cout <<
"Creating Multiexperiment with a total of " 193 << expts.size() <<
" experiments\n\n" <<
std::endl;
199 std::cout <<
"Systematics for the fit:\n";
205 std::cout <<
"\n\nSystematic correlations...\n";
206 if(!nueOnly && ! numuOnly && corrSysts){
209 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
228 std::vector <SystShifts> seedShifts = {};
230 for (
double systshift:{+0.5,-0.5}){
232 seedShifts.push_back(tempShifts);
236 double minchi23 = 1E20;
237 for(
double thseed:{0.45,0.55}){
240 << (thseed<0.5?
"LO ":
"UO ")
241 << (hie>0?
"NH " :
"IH ")
243 auto thisminchi = fit23.
Fit(calc,auxShifts ,{
249 minchi23=
min(minchi23, thisminchi);
257 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
260 TFile * outFCfile =
new TFile(outFCfilename+
".root",
"recreate");
265 v.Write(
"overall_minchi");
270 std::cerr <<
"\n WARNING Using 20 bins for dmsq surface\n\n";
274 for(
int hie: {-1, +1}){
276 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
278 std::cout <<
"Starting surface " << (hie>0?
"NH ":
"IH") <<
"\n\n";
281 if(!th13Surf && ! dmsqSurf){
291 auto surf1=surf23.
ToTH2(minchi23);
293 surf1->Write((TString)
"delta_th23_"+hieStr);
294 surf23.SaveTo(outfile, (TString)
"surf_delta_th23_"+hieStr);
297 std::vector<TH2*> profhists = surf23.GetProfiledHists();
299 profhists[0]->Write((hieStr+
"_DmSq32").c_str());
300 profhists[1]->Write((hieStr+
"_SinSq2Theta13").c_str());
302 for (
int i = 0;
i <
int(systs.size());
i++){
303 profhists[
i+2]->Write((hieStr+
"_"+systs[
i]->ShortName()).c_str());
306 if(!th13Surf && dmsqSurf){
312 &kFitDmSq32Scaled, steps,
313 (hie>0?2.1:-2.8),(hie>0?2.8:-2.1),
316 auto surf6=surf23m.
ToTH2(minchi23);
318 surf6->Write((TString)
"th23_dm32_"+(hie>0?
"NH":
"IH"));
319 surf23m.SaveTo(outfile, (TString)
"surf_th23_dm32_"+(hie>0?
"NH":
"IH"));
322 std::vector<TH2*> profhists = surf23m.GetProfiledHists();
324 profhists[0]->Write((hieStr+
"_SinSq2Theta13").c_str());
325 profhists[1]->Write((hieStr+
"_DeltaCPpi").c_str());
327 for (
int i = 0;
i <
int(systs.size());
i++){
328 profhists[
i+2]->Write((hieStr+
"_"+systs[
i]->ShortName()).c_str());
337 &kFitSinSq2Theta13, steps, 0, 0.35,
341 auto surf4 = surf13.
ToTH2(minchi23);
343 surf4->Write((TString)
"th13_delta_"+(hie>0?
"NH":
"IH"));
344 surf13.SaveTo(outfile, (TString)
"surf_th13_delta_"+(hie>0?
"NH":
"IH"));
347 std::vector<TH2*> profhists = surf13.GetProfiledHists();
349 profhists[0]->Write((hieStr+
"_SinSqTheta23").c_str());
350 profhists[1]->Write((hieStr+
"_DmSq32").c_str());
352 for (
int i = 0;
i <
int(systs.size());
i++){
353 profhists[
i+2]->Write((hieStr+
"_"+systs[
i]->ShortName()).c_str());
368 TFile *
infile =
new TFile (outfilename+
".root",
"read");
370 if(corrSysts &&
options.Contains(
"joint") && fccorr && 0) {
371 official =
new TFile(
"NOvA_nue_official_results.root",
"recreate");
372 TH2*
axes =
new TH2F(
"axes",
";#delta_{CP}/#pi;sin^{2}#theta_{23}",100,0,2,100,0.225,0.725);
377 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
378 double minchi23 = mins[0];
380 std::vector <TString> surfNames;
381 if(!th13Surf && ! dmsqSurf)surfNames = {
"delta_th23_"};
382 if (dmsqSurf) surfNames.push_back(
"th23_dm32_");
383 if (th13Surf) surfNames.push_back(
"th13_delta_");
384 for(TString surfName:surfNames){
386 for (
int hie:{-1,+1}){
388 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
391 if(surfName.Contains(
"delta_th23")) {
395 if(surfName.Contains(
"th23_dm32")) {
399 if(surfName.Contains(
"th13_delta")) {
411 if(!joint || !corrSysts ||
decomp !=
"combo" || dmsqSurf || th13Surf)
414 TString fcFolder=
"/nova/ana/nu_e_ana/Ana2017/FC/";
415 TString fcFileName ((TString)
"FCCol_"+(hie>0?
"NH":
"IH")+
"_surf_ssth23delta.root");
424 surf_1Sigma->Smooth();
425 surf_2Sigma->Smooth();
426 surf_3Sigma->Smooth();
431 Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.70);
432 Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.35);
433 Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.20);
434 Int_t kDarkColor = kGray+3;
437 auto g1 =
surf.GetGraphs(surf_1Sigma,minchi23);
438 auto g2 =
surf.GetGraphs(surf_2Sigma,minchi23);
439 auto g3 =
surf.GetGraphs(surf_3Sigma,minchi23);
441 kFillColor1, kFillColor2, kFillColor3,
442 kDarkColor, suffix + surfName + (hie>0?
"NH":
"IH") + (fccorr?
"_fccorr":
""));
445 surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
447 surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
454 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
457 leg->SetFillStyle(1001);
459 leg->SetX1(0.76);
leg->SetX2(0.89);
460 leg->SetY1(0.52);
leg->SetY2(0.85);
461 leg->SetHeader(
"No FC");
466 for(TString
ext: {
".pdf",
".root"}){
467 gPad->Print(outdir +
"contour_" + suffix +
"_" +
469 (hie>0?
"NH":
"IH") + (fccorr?
"_fccorr":
"") +
470 ((fccorr && smooth)?
"_smooth":
"") +
471 (zoomIn?
"_zoom":
"") +
475 if(debug && !fccorr){
477 int nSysts = corrSysts ? 30 : 0;
478 for (
int i = 0;
i < nSysts; ++
i){
479 hists.push_back((TH2*) infile->Get((TString)
"surf_" + surfName +
480 (hie > 0 ?
"NH" :
"IH") +
487 if(TString(
hists[i]->GetTitle()).Contains(
"delta")) {
489 hists[
i]->GetZaxis()->SetRangeUser(0,2);
493 double minz =
hists[
i]->GetBinContent(
hists[i]->GetMinimumBin());
494 double maxz =
hists[
i]->GetBinContent(
hists[i]->GetMaximumBin());
496 hists[
i]->GetZaxis()->SetRangeUser(-lim,+lim);
498 if(i==(
hists.size()-1)) {
500 hists[
i]->GetZaxis()->SetRangeUser(0.1,20);
505 hists[
i]->GetYaxis()->SetTitleOffset(0.8);
506 gPad->SetRightMargin(0.14);
507 surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
508 surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
509 surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
510 gPad->Print(outdir +
"debug_contour_" + suffix +
"_" +
512 (hie > 0 ?
"NH":
"IH") + (fccorr?
"_fccorr":
"") +
513 (zoomIn ?
"_zoom" :
"") +
525 TLatex * ltxNH =
new TLatex(0.16,0.89,
"#bf{NH}");
528 ltxNH->SetTextAlign(22);
530 if(th13) ltxNH->Draw();
531 else if(!vert) ltxNH->DrawLatex(.85,.20,
"#bf{NH}");
532 else ltxNH->DrawLatex(.85,.85,
"#bf{NH}");
535 TLatex * ltxIH =
new TLatex (0.16,0.46,
"#bf{IH}");
538 ltxIH->SetTextAlign(22);
540 if(th13) ltxIH->Draw();
541 else if(!vert) ltxIH->DrawLatex(.85,.20,
"#bf{IH}");
542 else ltxIH->DrawLatex(.85,.85,
"#bf{IH}");
547 Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
549 TLegend *
leg =
new TLegend(0.16,fccorr?0.19:0.17,0.65,fccorr?0.26:0.28);
551 leg->SetFillStyle(0);
552 leg->SetMargin(1.3*leg->GetMargin());
555 TGraph *
dummy =
new TGraph;
556 dummy->SetLineWidth(3);
559 dummy->SetLineColor(kDarkColor);
560 dummy->SetFillColor(kFillColor1);
561 leg->AddEntry(dummy->Clone(),
"1 #sigma",
"f");
562 dummy->SetFillColor(kFillColor2);
563 leg->AddEntry(dummy->Clone(),
"2 #sigma",
"f");
564 dummy->SetFillColor(kFillColor3);
565 leg->AddEntry(dummy->Clone(),
"3 #sigma",
"f");
568 leg->SetMargin(1.4*leg->GetMargin());
569 dummy->SetLineColor(kFillColor1);
570 dummy->SetLineStyle(kSolid);
571 leg->AddEntry(dummy->Clone(),
"1 #sigma",
"l");
573 dummy->SetLineColor(kFillColor2);
574 leg->AddEntry(dummy->Clone(),
"2 #sigma",
"l");
575 dummy->SetLineStyle(kSolid);
576 dummy->SetLineColor(kFillColor3);
577 leg->AddEntry(dummy->Clone(),
"3 #sigma",
"l");
583 dummy->SetMarkerStyle(kFullCircle);
584 dummy->SetMarkerColor(kDarkColor);
585 dummy->SetLineColor(kDarkColor);
586 dummy->SetLineStyle(7);
587 if(hie>0) leg->AddEntry(dummy->Clone(),
"Best Fit",
"p");
588 if(hie<0) leg->AddEntry((TObject*)0,
"#color[0]{Best Fit}",
"");
589 dummy->SetLineStyle(1);
591 if(!fccorr) leg->SetHeader(
"No Feldman-Cousins");
595 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
596 const Int_t kFillColor1,
const Int_t kFillColor2,
const Int_t kFillColor3,
597 const Int_t kDarkColor,
const TString surftype)
599 bool isJoint = surftype.Contains(
"joint");
600 bool isNH = surftype.Contains(
"NH");
601 bool isSysts = surftype.Contains(
"systs");
602 bool isFccorr = surftype.Contains(
"fccorr");
604 if (surftype.Contains(
"delta_th23")){
605 if(isJoint && isNH ){
606 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
608 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f");
609 if(isSysts && !isFccorr)
JoinGraphs(g1[0], g1[1], kWhite)->Draw(
"f");
610 if(isSysts && !isFccorr)
JoinGraphs(g1[0], g1[1], kFillColor1)->Draw(
"f");
611 if(!isSysts || isFccorr)
for (
auto &
g:g1) {
g->SetFillColor(kWhite);
g->DrawClone(
"f");}
612 if(!isSysts || isFccorr)
for (
auto &
g:g1) {
g->SetFillColor(kFillColor1);
g->Draw(
"f");}
615 if(isJoint && !isNH){
617 g->SetFillColor(kFillColor3);
g->Draw(
"f");}
619 g->SetFillColor(kWhite);
g->DrawClone(
"f");
620 g->SetFillColor(kFillColor2);
g->Draw(
"f");
625 if (surftype.Contains(
"th23_dm32")){
627 g->SetFillColor(kFillColor3);
g->Draw(
"f");}
629 g->SetFillColor(kWhite);
g->DrawClone(
"f");
630 g->SetFillColor(kFillColor2);
g->Draw(
"f");
633 g->SetFillColor(kWhite);
g->DrawClone(
"f");
634 g->SetFillColor(kFillColor1);
g->Draw(
"f");
638 if (surftype.Contains(
"th13_delta")) {
639 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
641 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f");
643 JoinGraphs(g1[0], g1[1], kFillColor1)->Draw(
"f");
647 for (
auto &
gs:{g3,g2,g1}){
649 g->SetLineColor(kDarkColor);
const Color_t kPrimColorIH
T max(const caf::Proxy< T > &a, T b)
const Color_t kPrimColorNH
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
fvar< T > fabs(const fvar< T > &x)
const Color_t kCosmicBackgroundColor
const FitDmSq32 kFitDmSq32
Simple record of shifts applied to systematic parameters.
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
const Color_t k2SigmaConfidenceColorNH
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
void CenterTitles(TH1 *histo)
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
virtual void SetDmsq32(const T &dmsq32)=0
const XML_Char const XML_Char * data
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
void SetPaletteBlueRedWhite()
void DrawHieTag(int hie, bool th13=true, bool vert=false)
void SetPaletteWhiteBlueDark()
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
void SetPaletteBlueRedCyclic()
virtual T GetDmsq32() const
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.
static float min(const float a, const float b, const float c)
TH2 * ToTH2(double minchi=-1) const
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
Combine multiple component experiments.
const Style_t k2SigmaConfidenceStyle
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
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)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const Float_t kBlessedLabelFontSize
const FitSinSq2Theta13 kFitSinSq2Theta13
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
std::string to_string(ModuleType mt)
void joint_fit_2017_contours(std::string decomp="prop", bool createFile=false, bool corrSysts=false, TString options="joint_fakeNHLO", bool dmsqSurf=false, bool th13Surf=false, bool fccorr=false, bool zoomIn=true, bool fillContour=false, bool smooth=true)
TH2 * SurfaceForSignificance(double sig)
const Float_t kBlessedTitleFontSize
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)
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
std::string UniqueName()
Return a different string each time, for creating histograms.
const Color_t k2SigmaConfidenceColorIH
Pseudo-experiments, binned to match a log-likelihood surface.
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.