59 float minx = 0.35, maxx = 0.65;
60 float minx_zoom = 0.58, maxx_zoom = 0.60;
61 float miny = 2.30,
maxy = 2.70;
62 float miny_zoom = 2.43, maxy_zoom = 2.50;
67 minx = 0.55; maxx = 0.70;
68 minx_zoom = 0.62; maxx_zoom = 0.68;
69 miny = 2.20;
maxy = 3.00;
70 miny_zoom = 2.555;maxy_zoom = 2.65;
79 minx = 0.35; maxx = 0.65;
80 minx_zoom = 0.48; maxx_zoom = 0.55;
81 miny = 2.20;
maxy = 2.70;
82 miny_zoom = 2.40; maxy_zoom = 2.55;
91 std::vector< const ISyst* >
systs;
94 std::map< TString, std::vector<const ISyst*> > mapSystLists;
105 std::string inFile =
"../results/surfprof_realdata_cosmics_"+anaName+
"__extrap_stats_2018calc_nh__numu2018.root";
106 TFile* fSurfSens =
new TFile(inFile.c_str());
108 fSurfSens ->
Close();
110 double bestFitDm = surfaceSens->GetBestFitY();
111 double bestFitSsq = surfaceSens->GetBestFitX();
112 std::cout<<
"Best fit from result contour: ssqth23: " << bestFitSsq
113 <<
", dm^2_32: " << bestFitDm <<
std::endl;
118 TString dout_name =
"../plots_blessing/";
119 std::string ddata_name =
"/nova/ana/nu_mu_ana/Ana2018/Data/";
120 std::string dcosm_name =
"/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
121 std::string dpred_name =
"/nova/ana/nu_mu_ana/Ana2018/Predictions/";
127 const std::string dpred_name_fhc= dpred_name+
"pred_numuconcat_fhc__numu2018.root";
128 const std::string dpred_name_rhc= dpred_name+
"pred_numuconcat_rhc__numu2018.root";
131 if(anaName==
"fhc") dpred_name_one = dpred_name_fhc;
136 for(
int quant = 1; quant <
totquant+1; quant++){
145 std::vector<Spectrum> s_cosmics;
152 TFile fcosmics_one(fcosmics_name_one.c_str());
153 auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get(
"cosmics_q1");
154 auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get(
"cosmics_q2");
155 auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get(
"cosmics_q3");
156 auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get(
"cosmics_q4");
165 TFile fcosmics_two(fcosmics_name_two.c_str());
166 auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get(
"cosmics_q1");
167 auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get(
"cosmics_q2");
168 auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get(
"cosmics_q3");
169 auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get(
"cosmics_q4");
180 for(
auto const& thisSystList: mapSystLists){
183 TCanvas* canvas_zoom =
new TCanvas(
"canvas_zoom", (NameZoom).c_str());
184 TString out_zoom = dout_name+NameZoom.c_str();
187 TCanvas*
canvas =
new TCanvas(
"canvas", (Name).c_str());
188 TString
out = dout_name+Name.c_str();
190 TH2F *axes_zoom =
new TH2F(
"setaxes_zoom",
"", 30, minx_zoom, maxx_zoom, 50, miny_zoom, maxy_zoom);
191 axes_zoom->GetXaxis()->SetTitle(
"sin^{2}#theta_{23}");
192 axes_zoom->GetYaxis()->SetTitle(
"#Deltam^{2}_{32} (10^{-3} eV^{2})");
193 axes_zoom->GetXaxis()->SetTitleOffset(0.95);
194 axes_zoom->GetYaxis()->SetTitleOffset(0.85);
195 axes_zoom->GetXaxis()->SetTitleSize(0.055);
196 axes_zoom->GetYaxis()->SetTitleSize(0.055);
197 axes_zoom->GetXaxis()->SetLabelSize(0.04);
198 axes_zoom->GetYaxis()->SetLabelSize(0.04);
199 axes_zoom->GetXaxis()->CenterTitle();
200 axes_zoom->GetYaxis()->CenterTitle();
201 axes_zoom->SetTitle(
"");
202 axes_zoom->GetYaxis()->SetDecimals();
204 TH2F *
axes =
new TH2F(
"setaxes",
"" ,38 ,minx, maxx, 52, miny,
maxy);
205 axes->GetXaxis()->SetTitle(
"sin^{2}#theta_{23}");
206 axes->GetYaxis()->SetTitle(
"#Deltam^{2}_{32} (10^{-3} eV^{2})");
207 axes->GetXaxis()->SetTitleOffset(0.95);
208 axes->GetYaxis()->SetTitleOffset(0.85);
209 axes->GetXaxis()->SetTitleSize(0.055);
210 axes->GetYaxis()->SetTitleSize(0.055);
211 axes->GetXaxis()->SetLabelSize(0.04);
212 axes->GetYaxis()->SetLabelSize(0.04);
213 axes->GetXaxis()->CenterTitle();
214 axes->GetYaxis()->CenterTitle();
216 axes->GetYaxis()->SetDecimals();
222 gPad->SetFillStyle(0);
223 gPad->SetTopMargin(0.08);
224 gPad->SetRightMargin(0.08);
225 gPad->SetBottomMargin(0.12);
226 gPad->SetLeftMargin(0.12);
227 gStyle->SetTitleOffset(0.95,
"x");
228 gStyle->SetTitleOffset(0.85,
"y");
231 surfaceSens->DrawBestFit(kBlack, kFullStar);
246 gPad->SetFillStyle(0);
247 gPad->SetTopMargin(0.08);
248 gPad->SetRightMargin(0.08);
249 gPad->SetBottomMargin(0.12);
250 gPad->SetLeftMargin(0.12);
251 gStyle->SetTitleOffset(0.95,
"x");
252 gStyle->SetTitleOffset(0.85,
"y");
255 surfaceSens->DrawBestFit(kBlack, kFullStar);
256 DrawLegendBFNull(surfaceSens->GetBestFitX(), surfaceSens->GetBestFitY(), kBlack, kFullStar);
258 TLegend *
leg =
new TLegend(0.40, 0.65, 0.70, 0.85);
259 TLegendEntry *
entry = leg->AddEntry(
"NULL",
"NOvA NH #nu_{#mu}+#bar{#nu}_{#mu} Systematic",
"h");
269 std::vector<const ISyst*> systList = thisSystList.second;
271 for(
const ISyst*
s: systList){
273 TGraph*
g =
new TGraph;
274 TGraph* gLO =
new TGraph;
283 std::vector<Spectrum> s_fakedata;
284 for(
int quant = 0; quant <
totquant; quant++){
288 s_fakedata.push_back( predictions[quant]->PredictSyst(calc, systListShift).
FakeData(
pot) );
292 std::vector <const IExperiment*> experiments;
294 for(
int quant = 0; quant <
totquant; quant++){
296 s_fakedata[quant] += s_cosmics[quant];
297 experiments.push_back(
new SingleSampleExperiment(predictions[quant], s_fakedata[quant], s_cosmics[quant], 0.1));
304 s_fakedata[quantN-1] += s_cosmics[quantN-1];
305 experiments.push_back(
new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1], s_cosmics[quantN-1], 0.1));
319 fitter.Fit(calc, IFitter::kQuiet);
321 g->SetPoint(g->GetN(),
323 kFitDmSq32Scaled.GetValue(calc));
330 fitterLO.Fit(calc, IFitter::kQuiet);
332 gLO->SetPoint(gLO->GetN(),
334 kFitDmSq32Scaled.GetValue(calc));
341 g->SetLineColor(
color[graphId]);
343 gLO->SetLineColor(
color[graphId]);
347 g->SetTitle((
s->LatexName()).c_str());
349 g->SetMarkerStyle(0);
351 gLO->SetTitle((
s->LatexName()).c_str());
352 gLO->SetFillColor(10);
353 gLO->SetMarkerStyle(0);
354 gLO->SetLineWidth(3);
356 leg->
AddEntry(g, (
s->LatexName()).c_str(),
"l");
363 surfaceSens->DrawBestFit(kBlack,kFullStar);
368 surfaceSens->DrawBestFit(kBlack,kFullStar);
373 canvas->SaveAs((out+
".root"));
374 canvas->SaveAs((out+
".pdf"));
375 canvas->SaveAs((out+
".png"));
377 file.open (out +
".txt");
378 file <<
"90% confidence level sensitivity at NOvA's combined neutrino + antineutrino 2018 stats only best fit. Each point in the coloured lines represent the new best fit when a single systematic is applied to the fit: from start to end, a shift to the systematic is applied to cover a +-1$\sigma$ range.";
384 canvas_zoom->SaveAs((out_zoom+
".root"));
385 canvas_zoom->SaveAs((out_zoom+
".pdf"));
386 canvas_zoom->SaveAs((out_zoom+
".png"));
387 std::ofstream file_zoom;
388 file_zoom.open (out_zoom +
".txt");
389 file_zoom <<
"Zoom in the 90% confidence level sensitivity at NOvA's combined neutrino + antineutrino 2018 stats only best fit. Each point in the coloured lines represent the new best fit when a single systematic is applied to the fit: from start to end, a shift to the systematic is applied to cover a +-1$\sigma$ range.";
391 canvas_zoom->Close();
std::vector< const ISyst * > MyMECSysts
std::vector< const ISyst * > MyFluxSysts
std::vector< const ISyst * > MyGeniePCASysts
bin1_2sigma SetFillColor(3)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Loads shifted spectra from files.
Simple record of shifts applied to systematic parameters.
const double livetime_fhc
std::vector< const ISyst * > MyCalibSysts
TCanvas * canvas(const char *nm, const char *ti, const char *a)
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Encapsulate code to systematically shift a caf::SRProxy.
ntuple SetFillStyle(1001)
Representation of a spectrum in any variable, with associated POT.
std::vector< const ISyst * > MyOtherSysts
Log-likelihood scan across two parameters.
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
std::vector< const ISyst * > MyGenieRwgtSysts
leg AddEntry(GRdata,"data","p")
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
std::vector< const ISyst * > MyGenieRwgtFrSysts
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
Combine multiple component experiments.
const FitSinSqTheta23UO kFitSinSqTheta23UO
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
std::vector< Style_t > styleline
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
A simple ascii-art progress bar.
std::vector< const ISyst * > MyXSecSysts
const FitSinSqTheta23LO kFitSinSqTheta23LO
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
void SetShift(const ISyst *syst, double shift, bool force=false)
std::vector< Color_t > color
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.