3 #include "CAFAna/Analysis/Fit.h" 5 #include "CAFAna/Analysis/NuePlotStyle.h" 6 #include "CAFAna/Analysis/Surface.h" 13 #include "CAFAna/Prediction/PredictionSystJoint2018.h" 15 #include "CAFAna/Systs/NumuSysts.h" 20 #include "CAFAna/nue/Ana2019/joint_fit_2019_loader_tools.h" 43 std::vector <Spectrum * >
GetMockData(TString ana_options);
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);
56 Double_t
white[9] = { 1, 1, 1, 1, 1, 1, 1, 1, 1 };
57 Double_t Red[9] = { 61./255., 99./255., 136./255., 181./255., 213./255., 225./255., 198./255., 99./255., 61./255.};
58 Double_t Green[9] = { 149./255., 140./255., 96./255., 83./255., 132./255., 178./255., 190./255., 140./255., 149./255};
59 Double_t Blue[9] = { 214./255., 203./255., 168./255., 135./255., 110./255., 100./255., 111./255., 203./255., 214./255.};
60 Double_t Length[9] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
62 TColor::CreateGradientColorTable(9, Length, Red, Green, Blue, 255, 1.0);
63 gStyle->SetNumberContours(255);
69 bool corrSysts =
true,
70 TString
options=
"joint_realData_both_onlyIH",
71 bool dmsqSurf =
false,
72 bool th13Surf =
false,
76 bool fillContour =
false)
80 bool nueOnly =
options.Contains(
"nueOnly");
81 bool numuOnly =
options.Contains(
"numuOnly");
82 bool joint =
options.Contains(
"joint");
83 assert (nueOnly || numuOnly || joint);
85 bool FHCOnly =
options.Contains(
"FHCOnly");
86 bool RHCOnly =
options.Contains(
"RHCOnly");
87 bool both =
options.Contains(
"both");
88 assert (FHCOnly || RHCOnly || both);
90 bool fake2019 =
options.Contains(
"fake2019");
91 bool realData =
options.Contains(
"realData");
93 bool onlyNH =
options.Contains(
"onlyNH");
94 bool onlyIH =
options.Contains(
"onlyIH");
97 if(dmsqSurf) suffix+=
"_dmsq";
98 if(th13Surf) suffix+=
"_th13";
99 if(corrSysts) suffix+=
"_systs";
100 else suffix+=
"_stats";
114 TString
outfilename (outdir +
"hist_contours_2019_" + suffix);
131 std::pair <TH1D*, double>
cos;
135 std::vector <predictions> preds;
136 std::vector <Spectrum*>
data;
137 std::vector <const IExperiment * > expts;
140 if(fake2019)
SetFakeCalc(calc_fake, 0.565, 2.48
e-3, 1.99);
143 if(FHCOnly || both ) {
146 if(RHCOnly || both ) {
153 if(FHCOnly || both ) {
156 for (
int i = 0;
i< nnumu;
i++ ){
160 if(RHCOnly || both ) {
163 for (
int i = 0;
i< nnumu;
i++ ){
173 cout<<
"\n\n Predictions used for the fit:\n";
177 if(FHCOnly || both ) {
179 data.insert(data.end(),
temp.begin(),
temp.end());
181 if(RHCOnly || both ) {
183 data.insert(data.end(),
temp.begin(),
temp.end());
187 if(FHCOnly || both ){
189 data.insert(data.end(),
temp.begin(),
temp.end());}
190 if(RHCOnly || both ){
192 data.insert(data.end(),
temp.begin(),
temp.end());}
196 for(
int i = 0;
i <
int(preds.size()); ++
i){
197 double POT = preds[
i].pot;
199 if(realData) thisdata = data[
i];
200 cout<<i<<
" "<<preds[
i].name<<
" POT "<<POT<<
" tot MC "<<preds[
i].pred->Predict(calc_fake).Integral(POT)<<
" cos "<<preds[
i].cos.first->Integral()<<
" cos er "<<preds[
i].cos.second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
216 std::cout <<
"Adding WorldReactorConstraint2017\n";
220 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
224 std::cout <<
"Creating Multiexperiment with a total of " 225 << expts.size() <<
" experiments\n\n" <<
std::endl;
236 std::vector< const ISyst * >
systs = {};
251 std::vector <SystShifts> seedShifts = {};
259 std::vector<double> seeds;
264 std::vector <th23helper> th23seeds;
271 std::vector<double> th23_all_seeds_NH = {0.45, 0.5, 0.55};
273 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
274 double seed_dmsq = 2.5e-3;
280 double minchi23 = 1E20;
283 for (
auto & thseed: th23seeds){
286 << (hie>0?
" NH " :
" IH ")
288 for (
auto const &
s:thseed.seeds)
std::cout <<
s <<
", ";
295 Fitter fit23(exptThis, fitvars, systs);
296 auto thisminchi = fit23.Fit(calc,auxShifts ,{
298 {thseed.var, thseed.seeds},
304 minchi23=
min(minchi23, thisminchi);
312 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
317 cout<<
"minchi wrote into the file "<<minchi23<<
endl;
318 v.Write(
"overall_minchi");
323 std::cerr <<
"\n WARNING Using 20 bins for dmsq surface\n\n";
328 for(
int hie: {-1, +1}){
330 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
332 std::cout <<
"Starting surface " << (hie>0?
"NH ":
"IH") <<
"\n\n";
333 cout<<
"the delta seeds are ";
334 for (
auto const &
s:delta_seeds)
std::cout <<
s <<
", ";
338 if(!th13Surf && ! dmsqSurf){
342 double low = (nueOnly||RHCOnly)?0:0.3;
343 double high = (nueOnly||RHCOnly)?1:0.7;
344 cout<<
"low "<<low<<
" high "<<high<<
endl;
346 Surface surf23(exptThis, calc,
353 auto surf1=surf23.ToTH2(minchi23);
355 surf1->Write((TString)
"delta_th23_"+hieStr);
356 surf23.SaveTo(outfile->mkdir((TString)
"surf_delta_th23_"+hieStr));
360 if(!th13Surf && dmsqSurf){
364 double low = (nueOnly||RHCOnly)?0.2:0.35;
365 double high = (nueOnly||RHCOnly)?1:0.7;
366 cout<<
"low "<<low<<
" high "<<high<<
endl;
370 Surface surf23m(exptThis, calc,
373 (hie>0?(RHCOnly?2.0:2.2):(RHCOnly?-3.5:-2.8)),(hie>0?(RHCOnly?3.2:2.85):(RHCOnly?-2.0:-2.3)),
376 auto surf6=surf23m.ToTH2(minchi23);
378 surf6->Write((TString)
"th23_dm32_"+(hie>0?
"NH":
"IH"));
379 surf23m.SaveTo(outfile->mkdir((TString)
"surf_th23_dm32_"+(hie>0?
"NH":
"IH")));
389 Surface surf13(exptThis, calc,
394 auto surf4 = surf13.ToTH2(minchi23);
396 surf4->Write((TString)
"th13_delta_"+(hie>0?
"NH":
"IH"));
397 surf13.SaveTo(outfile->mkdir((TString)
"surf_th13_delta_"+(hie>0?
"NH":
"IH")));
416 TFile *
infile =
new TFile (outfilename+
".root",
"read");
418 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
419 double minchi23 = mins[0];
421 std::vector <TString> surfNames;
422 if(!th13Surf && ! dmsqSurf)surfNames = {
"delta_th23_"};
423 if (dmsqSurf) surfNames.push_back(
"th23_dm32_");
424 if (th13Surf) surfNames.push_back(
"th13_delta_");
425 for(TString surfName:surfNames){
427 for (
int hie:{-1,+1}){
429 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
432 if(surfName.Contains(
"delta_th23")) {
433 if (zoomIn)
DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
436 if(surfName.Contains(
"th23_dm32")) {
437 if (hie>0)
DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?1.8:2.05), (RHCOnly?3.3:2.85));
438 else DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?-3.5:-2.95), (RHCOnly?-2.0:-2.25));
440 if(surfName.Contains(
"th13_delta")) {
446 TH2 * hsurf = (TH2*)infile->Get(surfName+(hie>0?
"NH":
"IH"));
452 Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
453 Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
454 Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
455 Int_t kDarkColor = kGray+3;
458 auto g1 =
surf.GetGraphs(surf_1Sigma,minchi23);
459 auto g2 =
surf.GetGraphs(surf_2Sigma,minchi23);
460 auto g3 =
surf.GetGraphs(surf_3Sigma,minchi23);
462 kFillColor1, kFillColor2, kFillColor3,
463 kDarkColor, suffix + surfName + (hie>0?
"NH":
"IH"));
466 cout<<
"minchi from file is "<<minchi23<<
endl;
468 surf.DrawContour(surf_1Sigma, kSolid, kFillColor1,minchi23);
470 surf.DrawContour(surf_3Sigma, kSolid, kFillColor3,minchi23);
477 kFillColor1,kFillColor2,kFillColor3,kDarkColor,
482 leg->SetFillStyle(1001);
484 leg->SetX1(0.7);
leg->SetX2(0.8);
485 leg->SetY1(0.4);
leg->SetY2(0.85);
486 leg->SetHeader(
"No FC");
488 TGraph *
dummy =
new TGraph;
489 dummy->SetLineWidth(2);
490 dummy->SetLineColor(kGray+2);
491 dummy->SetFillColor(kGray);
492 leg->AddEntry(dummy,
"#splitline{Reactor}{68%C.L.}",
"f");
499 for(TString
ext: {
".pdf",
".root"}){
500 gPad->Print(outdir +
"contour_" + suffix +
"_" +
503 (zoomIn?
"_zoom":
"") +
511 int nSysts = corrSysts ? 51 : 0;
512 for (
int i = 0;
i < nSysts; ++
i){
513 hists.push_back((TH2*) infile->Get((TString)
"surf_" + surfName +
514 (hie > 0 ?
"NH" :
"IH") +
521 if(TString(
hists[i]->GetTitle()).Contains(
"delta")) {
523 hists[
i]->GetZaxis()->SetRangeUser(0,2);
525 else gStyle->SetPalette(87);
527 double minz =
hists[
i]->GetBinContent(
hists[i]->GetMinimumBin());
528 double maxz =
hists[
i]->GetBinContent(
hists[i]->GetMaximumBin());
530 hists[
i]->GetZaxis()->SetRangeUser(-lim,+lim);
532 if(i==(
hists.size()-1)) {
533 gStyle->SetPalette(56);
534 hists[
i]->GetZaxis()->SetRangeUser(0.,11.83);
538 hists[
i]->GetYaxis()->SetTitleOffset(0.8);
539 gPad->SetRightMargin(0.14);
540 surf.DrawContour(surf_1Sigma, kSolid, kDarkColor,minchi23);
541 surf.DrawContour(surf_2Sigma, kSolid, kDarkColor,minchi23);
542 surf.DrawContour(surf_3Sigma, kSolid, kDarkColor,minchi23);
543 gPad->Print(outdir +
"debug_contour_" + suffix +
"_" +
545 (hie > 0 ?
"NH":
"IH") +
546 (zoomIn ?
"_zoom" :
"") +
559 TLatex * ltxNH =
new TLatex(0.16,0.89,
"#bf{NH}");
562 ltxNH->SetTextAlign(22);
564 if(th13) ltxNH->Draw();
565 else if(!vert) ltxNH->DrawLatex(.85,.20,
"#bf{NH}");
566 else ltxNH->DrawLatex(.85,.85,
"#bf{NH}");
569 TLatex * ltxIH =
new TLatex (0.16,0.46,
"#bf{IH}");
572 ltxIH->SetTextAlign(22);
574 if(th13) ltxIH->Draw();
575 else if(!vert) ltxIH->DrawLatex(.85,.20,
"#bf{IH}");
576 else ltxIH->DrawLatex(.85,.85,
"#bf{IH}");
582 void FillGraphs(std::vector<TGraph*> g1,std::vector<TGraph*> g2,std::vector<TGraph*> g3,
583 const Int_t kFillColor1,
const Int_t kFillColor2,
const Int_t kFillColor3,
584 const Int_t kDarkColor,
const TString surftype)
586 bool isJoint = surftype.Contains(
"joint");
587 bool isBoth = surftype.Contains(
"both");
588 bool isRHC = surftype.Contains(
"RHCOnly");
589 bool isNH = surftype.Contains(
"NH");
590 bool isSysts = surftype.Contains(
"systs");
591 bool isFccorr = surftype.Contains(
"fccorr");
593 if (surftype.Contains(
"delta_th23")){
595 if(isJoint && isNH ){
598 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
601 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f same");
602 JoinGraphs(g2[1], g2[2], kFillColor2)->Draw(
"f same");
605 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
607 JoinGraphs(g2[0], g2[1], kWhite)->Draw(
"c f");
608 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"c f");
612 for (
auto &
g:g1) {
g->SetFillColor(kWhite);
g->DrawClone(
"c f");}
613 for (
auto &
g:g1) {
g->SetFillColor(kFillColor1);
g->Draw(
"c f");}
614 for (
auto &
g:g1) {
g->Draw(
"c");}
618 if(isJoint && !isNH){
619 if(isRHC)
JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
621 JoinGraphs(g3[0], g3[2], kFillColor3)->Draw(
"f c");
622 g3[1]->SetFillColor(kFillColor3); g3[1]->Draw(
"f c");
623 g3[3]->SetFillColor(kFillColor3); g3[3]->Draw(
"f c");
630 g->SetFillColor(kWhite);
g->DrawClone(
"f c");
631 g->SetFillColor(kFillColor2);
g->Draw(
"f c");
632 g->SetLineColor(kDarkColor);
637 g->SetFillColor(kWhite);
g->DrawClone(
"f c");
638 g->SetFillColor(kFillColor1);
g->Draw(
"f c");
639 g->SetLineColor(kDarkColor);
643 for(
auto &
gs:{g3, g2, g1}){
645 g->SetLineColor(kDarkColor);
653 if (surftype.Contains(
"th23_dm32")){
655 g->SetFillColor(kFillColor3);
g->Draw(
"cf");
658 g->SetFillColor(kWhite);
g->DrawClone(
"cf");
659 g->SetFillColor(kFillColor2);
g->Draw(
"cf");
662 g->SetFillColor(kWhite);
g->DrawClone(
"cf");
663 g->SetFillColor(kFillColor1);
g->Draw(
"cf");
665 for (
auto &
gs:{g3,g2,g1}){
667 g->SetLineColor(kDarkColor);
675 if (surftype.Contains(
"th13_delta")) {
677 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f");
680 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f");
683 for (
auto &
g:g2) {
g->SetFillColor(kWhite);
g->DrawClone(
"f");}
684 for (
auto &
g:g2) {
g->SetFillColor(kFillColor2);
g->Draw(
"f");}
686 for (
auto &
g:g1) {
g->SetFillColor(kWhite);
g->DrawClone(
"f");}
687 for (
auto &
g:g1) {
g->SetFillColor(kFillColor1);
g->Draw(
"f");}
693 JoinGraphs(g3[0], g3[1], kFillColor3)->Draw(
"f c");
694 JoinGraphs(g2[0], g2[1], kWhite)->Draw(
"f c");
695 JoinGraphs(g2[0], g2[1], kFillColor2)->Draw(
"f c");
700 JoinGraphs(g1[0], g1[1], kWhite)->Draw(
"f c");
701 JoinGraphs(g1[0], g1[1], kFillColor1)->Draw(
"f c");
704 for (
auto &
gs:{g3,g2,g1}){
706 g->SetLineColor(kDarkColor);
717 TBox* box =
new TBox(mean-err, 0, mean+err, 2);
718 box->SetFillColorAlpha(kGray,0.7);
720 TLine* linLo =
new TLine(mean-err, 0, mean-err, 2);
721 linLo->SetLineColor(kGray+2);
722 linLo->SetLineWidth(2);
724 TLine* linHi =
new TLine(mean+err, 0, mean+err, 2);
725 linHi->SetLineColor(kGray+2);
726 linHi->SetLineWidth(2);
737 cout<<
"\n\n ============= You're loading Mock data from the file ============="<<
endl;
738 std::vector <Spectrum *> spec;
740 bool RHCmode = ana_options.Contains(
"rhc");
741 bool FHCmode = ana_options.Contains(
"fhc");
742 bool NueData = ana_options.Contains(
"nue");
743 bool NumuData = ana_options.Contains(
"numu");
745 if (FHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+
"mock_data_for_yn_tutorial.root",
"data_exp_0").
release());
746 if (RHCmode) spec.push_back(LoadFromFile<Spectrum>(Dir+
"mock_data_for_yn_tutorial.root",
"data_exp_1").
release());
750 for (
int i = 2;
i < 6; ++
i) {
752 spec.push_back(
temp);
756 for (
int i = 6;
i < 10; ++
i) {
758 spec.push_back(
temp);
const Color_t kPrimColorIH
T max(const caf::Proxy< T > &a, T b)
const Color_t kPrimColorNH
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)
void Tutorial2019FitContours(bool createFile=false, bool corrSysts=true, TString options="joint_realData_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false, bool zoomIn=true, bool fillContour=false)
Start the Demo script (it's a bit simplified nue / Ana2018 / FitandFC / joint_fit_2018_contours.C):
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)
Simple record of shifts applied to systematic parameters.
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
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.
static SystShifts Nominal()
std::vector< Spectrum * > GetMockData(TString ana_options)
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
string outfilename
knobs that need extra care
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
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.
virtual T GetDmsq32() const
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
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.
std::vector< double > POT
static float min(const float a, const float b, const float c)
virtual T GetTh13() 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.
void DrawReactorConstraint(double mean, double err)
const Style_t k2SigmaConfidenceStyle
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
double GetPOT(bool isFHC=true)
Standard interface to all prediction techniques.
const Float_t kBlessedLabelFontSize
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
std::pair< Spectrum *, double > cos
void DrawHieTag(int hie, bool th13=true, bool vert=false)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
void SetPaletteDeltaCyclicNew()
const Color_t k2SigmaConfidenceColorIH
Compare a single data spectrum to the MC + cosmics expectation.