2 #include "TDirectory.h" 23 #include "CAFAna/Extrap/ExtrapSterile.h" 61 const IFitVar*
var, std::vector<const IFitVar*> profVars,
62 TDirectory* rootOut,
int nbins,
double min,
double max,
136 TDirectory* rootOut, FILE* textOFS,
154 TDirectory* rootOut, FILE* textOFS,
162 TH1::AddDirectory(0);
169 std::string loadLocation = folder + filenm +
".root";
170 std::string saveLocation = folder + filenm +
"Ana.root";
171 std::string textLocation = folder + filenm +
"Ana.txt";
174 TFile* rootL =
new TFile(loadLocation.c_str(),
"UPDATE");
175 TDirectory* tmpL = gDirectory;
176 TDirectory* loadDir = gDirectory;
178 loadDir->cd((loadLocation +
":/decompNC").c_str());
180 loadDir->cd((loadLocation +
":/prediction").c_str());
182 loadDir->cd((loadLocation +
":/sData1Yr").c_str());
184 loadDir->cd((loadLocation +
":/sData3Yr").c_str());
186 loadDir->cd((loadLocation +
":/sCos1Yr").c_str());
188 loadDir->cd((loadLocation +
":/sCos3Yr").c_str());
190 loadDir->cd((loadLocation +
":/pFDPENum").c_str());
192 loadDir->cd((loadLocation +
":/pFDPDen").c_str());
194 loadDir->cd((loadLocation +
":/pFDEDen").c_str());
196 loadDir->cd((loadLocation +
":/sNDPENum").c_str());
198 loadDir->cd((loadLocation +
":/sNDPDen").c_str());
200 loadDir->cd((loadLocation +
":/sNDEDen").c_str());
211 calc4f->
SetDm(4, 0.5);
215 std::string labelEReco =
"Calorimetric Energy (GeV)";
218 TFile* rootF =
new TFile(saveLocation.c_str(),
"RECREATE");
220 textF = fopen(textLocation.c_str(),
"w");
228 PlotStack(&pred, sCosmic1y, calc3f, calc4f, rootF, textF,
"FDSpectra1yr", titleHelper, 6);
229 PlotStack(&decompNC, rootF, textF,
"NDSpectra1yr", titleHelper, 6);
232 titleHelper =
" Spectra for 18e20 POT";
233 PlotStack(&pred, sCosmic3y, calc3f, calc4f, rootF, textF,
"FDSpectra3yr", titleHelper, 18);
234 PlotStack(&decompNC, rootF, textF,
"NDSpectra3yr", titleHelper, 18);
237 PlotPurEff(&pFDPurEffNum, &pFDPurityDen, &pFDEfcncyDen,
238 sNDPurEffNum, sNDPurityDen, sNDEfcncyDen,
239 calc3f, sCosmic1y, rootF, textF,
"PurEff", labelEReco);
254 std::vector<const IFitVar*> fitvars;
281 const Color_t kFitColor =
kRed;
289 rootF, avals,
"1", kFitColor
298 rootF, avals,
"3", kFitColor
329 const IFitVar*
var, std::vector<const IFitVar*> profVars,
330 TDirectory* rootOut,
int nbins,
double min,
double max,
335 TH1*
h = (profVars.size() > 0 ?
336 Profile(expt, calc, var, nbins, min, max, -1, profVars) :
337 Slice(expt, calc, var, nbins, min, max));
342 h->SetLineColor(color);
345 h->SetName(name.c_str());
346 rootOut->WriteTObject(h);
349 TLatex *
tex =
new TLatex();
351 tex->SetTextFont(42);
352 tex->SetTextSize(0.037);
353 tex->SetTextAlign(11);
356 double xTex = 0.92, y68 = 0.25, y90 = 0.52;
359 TLine* line68 =
new TLine(min, 1, max, 1);
360 TLine* line90 =
new TLine(min, 2.71, max, 2.71);
364 TLine* line5 =
new TLine(min, 5, max, 5);
367 line68->SetLineStyle(2);
368 line90->SetLineStyle(2);
371 line68->SetLineWidth(2);
372 line90->SetLineWidth(2);
373 line5 ->SetLineWidth(2);
377 std::string cname =
"c1DSlice" + name.substr(1, name.size()-1);
383 std::string cpotstr = ( name.find(
"1") != std::string::npos ?
384 "6e20 POT" :
"18e20 POT" );
387 std::string ctitle =
"Delta chi^2 for Theta" + indices;
389 ctitle += ( name.find(
"Prof") != std::string::npos ?
390 " (Profiling Over Other Angles) for " + cpotstr :
394 TCanvas*
c =
new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
396 line68->DrawClone(
"same");
397 line90->DrawClone(
"same");
398 line5 ->DrawClone(
"same");
399 tex->DrawLatex(xTex, y68,
"68%");
400 tex->DrawLatex(xTex, y90,
"90%");
402 rootOut->WriteTObject(c);
418 fullname = name +
"23";
421 fullname = name +
"34";
424 fullname = name +
"24";
432 fullname = name +
"23Prof";
435 fullname = name +
"34Prof";
438 fullname = name +
"24Prof";
492 Plot2DSlices(&s3424Prof, &s2324Prof, &s3423Prof, rootOut, years,
"Prof");
506 fullname = name +
"34Th24" + prof;
507 TH1* h3424 = s3424->
ToTH2();
508 h3424->SetName(fullname.c_str());
509 rootOut->WriteTObject(h3424);
511 fullname = name +
"23Th24" + prof;
512 TH1* h2324 = s2324->
ToTH2();
513 h2324->SetName(fullname.c_str());
514 rootOut->WriteTObject(h2324);
516 fullname = name +
"34Th23" + prof;
517 TH1* h3423 = s3423->
ToTH2();
518 h3423->SetName(fullname.c_str());
519 rootOut->WriteTObject(h3423);
522 std::string cname =
"c2DSlice" + years +
"yr" + prof;
524 std::string potstr = ( name.find(
"1") != std::string::npos ?
525 "6e20 POT" :
"18e20 POT" );
529 "Surfaces (Profiling Over Other Angles) for " + potstr :
530 "Surfaces for " + potstr
534 TCanvas*
c =
new TCanvas(cname.c_str(), ctitle.c_str(), 800, 800);
552 rootOut->WriteTObject(c);
563 double realPOT = (double)POT*1e20;
566 sprintf(potBuff,
"%d", POT);
570 TH1* hAll = spectra[0].
ToTH1(realPOT);
571 TH1* hNC = spectra[1].
ToTH1(realPOT);
572 TH1* hNumu = spectra[2].
ToTH1(realPOT);
573 TH1* hNue = spectra[3].
ToTH1(realPOT);
574 TH1* hCos = spectra[4].
ToTH1(realPOT);
576 TH1* hStrl = (sterile ? sterile ->
ToTH1(realPOT) :
nullptr);
577 TH1* hStrlNC = (sterileNC ? sterileNC->
ToTH1(realPOT) :
nullptr);
582 std::string yLabel = (det.compare(
"ND") == 0 ?
"10^{3} Events" :
"Events");
583 yLabel +=
"/" + potStr +
" / 0.25 GeV";
584 std::string histTitle =
";" + xLabel +
";" + yLabel;
587 fprintf(textOFS,
"Event numbers at the %2.2s for %2de20 POT:\n", det.c_str(),
POT);
588 if(hStrlNC) { fprintf(textOFS,
"%12.12s, ",
"NC (Sterile)"); }
589 else { fprintf(textOFS,
"%14.14s",
""); }
590 fprintf(textOFS,
"%12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
591 "NC (3 Flav)",
"All",
"Numu",
"Nue",
"Cosmic",
"FOM");
594 double nNC3F = hNC ->Integral();
595 double nAll = hAll ->Integral();
596 double nNumu = hNumu->Integral();
597 double nNue = hNue ->Integral();
598 double nCos = hCos ->Integral();
599 double fom = nNC3F/
sqrt(nAll);
601 double nNCSt = hStrlNC->Integral();
602 fprintf(textOFS,
"%10E, ", nNCSt);
605 fprintf(textOFS,
"%14.14s",
"");
607 fprintf(textOFS,
"%10E, %10E, %10E, %10E, %10E, %10E\n\n",
608 nNC3F, nAll, nNumu, nNue, nCos, fom);
617 hAll ->SetMinimum(0);
619 hNumu->SetMinimum(0);
620 hNue ->SetMinimum(0);
621 hCos ->SetMinimum(0);
622 if(hStrlNC) { hStrl->SetMinimum(0); }
629 hCos ->SetLineColor(kCosmicBackgroundColor);
632 hStrlNC->SetLineStyle(2);
635 hAll ->SetTitle(histTitle.c_str());
636 hNC ->SetTitle(histTitle.c_str());
637 hNumu->SetTitle(histTitle.c_str());
638 hNue ->SetTitle(histTitle.c_str());
639 hCos ->SetTitle(histTitle.c_str());
640 if(hStrlNC) { hStrlNC->SetTitle(histTitle.c_str()); }
644 TH1* hCosStack = (TH1*)hCos->Clone();
647 TH1* hNumuStack = (TH1*)hCosStack->Clone();
648 hNumuStack->Add(hNumu);
651 TH1* hNueStack = (TH1*)hNumuStack->Clone();
652 hNueStack->Add(hNue);
655 TH1* hNCStack = (TH1*)hNC->Clone();
670 hCosStack ->SetLineColor(kCosmicBackgroundColor);
671 hCosStack ->SetFillColor(kCosmicBackgroundColor);
673 hNCStack ->SetMinimum(0);
674 hNumuStack->SetMinimum(0);
675 hNueStack ->SetMinimum(0);
676 hCosStack ->SetMinimum(0);
680 double maxvalstack = 0.;
681 for(
int i = 1,
n = hNC->GetNbinsX();
i <=
n; ++
i) {
683 hNC ->SetBinError(
i, error);
684 hNCStack->SetBinError(
i, error);
686 maxvalstack =
std::max(maxvalstack, hNCStack->GetBinContent(
i) +
error);
690 if(det.compare(
"ND") == 0) {
691 double ndScale = 1./1000.;
692 hAll ->Scale(ndScale);
693 hNC ->Scale(ndScale);
694 hNue ->Scale(ndScale);
695 hNumu->Scale(ndScale);
696 hNCStack ->Scale(ndScale);
697 hNueStack ->Scale(ndScale);
698 hNumuStack->Scale(ndScale);
700 maxvalstack *= ndScale;
703 TLatex *
tex =
new TLatex();
705 tex->SetTextFont(42);
706 tex->SetTextSize(0.037);
707 tex->SetTextAlign(11);
710 TCanvas*
c =
new TCanvas(name.c_str(), title.c_str(), 800, 500);
711 hNue->SetMaximum(maxval*1.1);
712 hNue->GetXaxis()->SetNdivisions(-406);
715 double xL = 0.5, xR = 0.7, yB = 0.59, yT = 0.84;
716 if(det.compare(
"ND") == 0) { yB += 0.05; }
717 TLegend*
leg =
new TLegend(xL, yB, xR, yT);
718 leg->SetBorderSize(0);
719 leg->SetFillColor(kWhite);
720 leg->SetFillStyle(0);
721 leg->SetFillStyle(0);
722 leg->SetTextFont(42);
723 leg->SetTextSize(0.037);
724 leg->AddEntry(hAll,
"Total Prediction",
"l");
725 leg->AddEntry(hNC,
"NC 3 Flavor Prediction",
"le");
726 leg->AddEntry(hNue,
"#nu_{e} CC Background",
"l");
727 leg->AddEntry(hNumu,
"#nu_{#mu} CC Background",
"l");
728 if(det.compare(
"FD") == 0) {
729 leg->AddEntry(hCos,
"Cosmic Background",
"l");
733 hNumu->Draw(
"hist same");
734 if(det.compare(
"FD") == 0) { hCos->Draw(
"hist same"); }
735 hAll ->Draw(
"hist same");
738 if(det.compare(
"FD") == 0) {
739 tex->DrawLatex(0.505, yB - 0.045,
"#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
740 tex->DrawLatex(0.505, yB - 0.100,
"#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
741 tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
748 rootOut->WriteTObject(c);
751 TCanvas* cStack =
new TCanvas((name +
"Stack").c_str(), title.c_str(), 800, 500);
752 hNueStack->SetMaximum(maxvalstack*1.1);
753 hNueStack->GetXaxis()->SetNdivisions(-406);
755 if(det.compare(
"ND") == 0) { yB += 0.05; }
756 TLegend* legStack =
new TLegend(xL, yB, xR, yT);
757 legStack->SetBorderSize(0);
758 legStack->SetFillColor(kWhite);
759 legStack->SetFillStyle(0);
760 legStack->SetFillStyle(0);
761 legStack->SetTextFont(42);
762 legStack->SetTextSize(0.037);
763 legStack->AddEntry(hNCStack,
"NC 3 Flavor Prediction",
"le");
764 if(hStrlNC) { legStack->AddEntry(hStrlNC,
"NC 4 Flavor Prediction",
"l"); }
765 legStack->AddEntry(hNueStack,
"#nu_{e} CC Background",
"f");
766 legStack->AddEntry(hNumuStack,
"#nu_{#mu} CC Background",
"f");
767 if(det.compare(
"FD") == 0) {
768 legStack->AddEntry(hCosStack,
"Cosmic Background",
"f");
771 hNueStack ->Draw(
"hist");
772 hNumuStack->Draw(
"hist same");
773 if(det.compare(
"FD") == 0) { hCosStack->Draw(
"hist same"); }
774 hNCStack ->Draw(
"same");
775 if(hStrlNC) { hStrlNC->Draw(
"hist same"); }
777 if(det.compare(
"FD") == 0) {
778 tex->DrawLatex(0.505, yB - 0.045,
"#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
779 tex->DrawLatex(0.505, yB - 0.100,
"#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
780 tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
786 rootOut->WriteTObject(cStack);
810 PlotStack(spectraND, rootOut, textOFS, name, fulltitle, det, POT);
818 TDirectory* rootOut, FILE* textOFS,
834 PlotStack(spectraFD, rootOut, textOFS, name, fulltitle, det, POT, &sterile, &sterileNC);
843 TDirectory* rootOut, FILE* textOFS,
858 TH1* hFDEff = (TH1*)hFDNum->Clone();
859 hFDEff->Divide(hFDEffDen);
861 TH1* hNDEff = (TH1*)hNDNum->Clone();
862 hNDEff->Divide(hNDEffDen);
864 TH1* hFDPur = (TH1*)hFDNum->Clone();
865 hFDPur->Divide(hFDPurDen);
867 TH1* hNDPur = (TH1*)hNDNum->Clone();
868 hNDPur->Divide(hNDPurDen);
870 hFDEff->SetTitle((
";" + title +
";Efficiency, Purity / 0.25 GeV").c_str());
871 hFDPur->SetTitle((
";" + title +
";Efficiency, Purity / 0.25 GeV").c_str());
872 hNDEff->SetTitle((
";" + title +
";Efficiency, Purity / 0.25 GeV").c_str());
873 hNDPur->SetTitle((
";" + title +
";Efficiency, Purity / 0.25 GeV").c_str());
875 hFDEff->SetLineColor(
kBlue);
876 hNDEff->SetLineColor(
kBlue);
878 hFDPur->SetLineColor(
kRed);
879 hNDPur->SetLineColor(
kRed);
881 hNDEff->SetLineStyle(2);
882 hNDPur->SetLineStyle(2);
884 hFDEff->SetMaximum(1.);
885 hFDEff->SetMinimum(0.);
886 hFDEff->GetXaxis()->SetNdivisions(-406);
888 hNDEff->SetMaximum(1.);
889 hNDEff->SetMinimum(0.);
891 hFDPur->SetMaximum(1.);
892 hFDPur->SetMinimum(0.);
894 hNDPur->SetMaximum(1.);
895 hNDPur->SetMinimum(0.);
902 TLegend*
leg =
new TLegend(0.6, 0.35, 0.85, 0.55);
903 leg->AddEntry(hFDEff,
"FD MC Efficiency",
"L");
904 leg->AddEntry(hNDEff,
"ND MC Efficiency",
"L");
905 leg->AddEntry(hFDPur,
"FD MC Purity",
"L");
906 leg->AddEntry(hNDPur,
"ND MC Purity",
"L");
908 TCanvas*
c =
new TCanvas(name.c_str(),
"Purity and Efficiency", 800, 500);
909 hFDEff->Draw(
"hist");
910 hNDEff->Draw(
"hist same");
911 hFDPur->Draw(
"hist same");
912 hNDPur->Draw(
"hist same");
918 rootOut->WriteTObject(c);
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string det, double POT, double potEquiv, PlotOptions opt, Spectrum *dataspec)
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
void SetNFlavors(int nflavors)
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual Spectrum AntiNueComponent() const =0
void PlotPurEff(Spectrum pureff[], Spectrum effcyD[], TDirectory *out, std::string name, std::string title)
const Color_t kCosmicBackgroundColor
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
virtual Spectrum NumuComponent() const =0
A simple Gaussian constraint on an arbitrary IFitVar.
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
General interface to oscillation calculators.
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
void Plot1DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, std::string years, Color_t color)
static std::unique_ptr< PredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
void Plot1DSlice(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var, std::vector< const IFitVar * > profVars, TDirectory *rootOut, int nbins, double min, double max, std::string name, Color_t color)
void CenterTitles(TH1 *histo)
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
const Color_t kNumuBackgroundColor
static std::unique_ptr< ProportionalDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Representation of a spectrum in any variable, with associated POT.
const Color_t kNueSignalColor
Log-likelihood scan across two parameters.
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Charged-current interactions.
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
void Draw() const
Draw the surface itself.
void Plot2DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, std::string years)
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
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.
TGraph * Slice(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, MinuitFitter::FitOpts opts)
scan in one variable, holding all others constant
virtual Spectrum AntiNumuComponent() const =0
std::vector< double > POT
static float min(const float a, const float b, const float c)
TH2 * ToTH2(double minchi=-1) const
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
Splits Data proportionally according to MC.
void SetDm(int i, double dm)
void ResetAngles(osc::OscCalcSterile *calc)
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
virtual Spectrum NCTotalComponent() const
Standard interface to all decomposition techniques.
Base class defining interface for experiments.
Neutral-current interactions.
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Both neutrinos and antineutrinos.
A prediction object compatible with sterile oscillations.
Standard interface to all prediction techniques.
const Color_t kTotalMCColor
All neutrinos, any flavor.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const Color_t kNCBackgroundColor
Compare a data spectrum to MC expectation using only the event count.
virtual Spectrum NueComponent() const =0
Perform MINUIT fits in one or two dimensions.