84 anaName2=
"Antineutrino beam";
92 anaName2=
"Neutrino beam";
110 auto calcNoOsc = osc.
Copy();
120 TString dout_name =
"";
121 std::string ddata_name =
"/nova/ana/nu_mu_ana/Ana2018/Data/";
122 std::string dcosm_name =
"/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
123 std::string dpred_name =
"/nova/ana/nu_mu_ana/Ana2018/Predictions/";
127 TFile* fdata_one =
new TFile((ddata_name +
"fd_dataspectrum_" +
horn_one +
"_full__numu2018.root").c_str());
128 TFile* fdata_two =
new TFile((ddata_name +
"fd_dataspectrum_" +
horn_two +
"_full__numu2018.root").c_str());
129 std::vector<Spectrum*> s_realdata;
131 for(
int quant=1; quant<
totquant+1; quant++){
132 if(quant<5) s_realdata.push_back(LoadFrom<Spectrum>(fdata_one, Form(
"fd_data_q%d", quant) ).
release());
133 else s_realdata.push_back(LoadFrom<Spectrum>(fdata_two, Form(
"fd_data_q%d", quant-4) ).
release());
137 double pot1here = s_realdata[0]->POT();
142 std::vector<Spectrum> s_cosmics;
143 std::vector<TH1*> hCosmics;
147 TFile fcosmics_one(fcosmics_name_one.c_str());
148 hCosmics.push_back((TH1*)fcosmics_one.Get(
"cosmics_q0"));
149 hCosmics.push_back((TH1*)fcosmics_one.Get(
"cosmics_q1"));
150 hCosmics.push_back((TH1*)fcosmics_one.Get(
"cosmics_q2"));
151 hCosmics.push_back((TH1*)fcosmics_one.Get(
"cosmics_q3"));
152 hCosmics.push_back((TH1*)fcosmics_one.Get(
"cosmics_q4"));
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");
163 TFile fcosmics_two(fcosmics_name_two.c_str());
164 hCosmics.push_back((TH1*)fcosmics_two.Get(
"cosmics_q1"));
165 hCosmics.push_back((TH1*)fcosmics_two.Get(
"cosmics_q2"));
166 hCosmics.push_back((TH1*)fcosmics_two.Get(
"cosmics_q3"));
167 hCosmics.push_back((TH1*)fcosmics_two.Get(
"cosmics_q4"));
168 auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get(
"cosmics_q1");
169 auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get(
"cosmics_q2");
170 auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get(
"cosmics_q3");
171 auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get(
"cosmics_q4");
178 for(
int quantId = 5; quantId < totquant+1; quantId++){
179 hCosmics[0]->Add(hCosmics[quantId]);
185 const std::string dpred_name_fhc= dpred_name+
"pred_numuconcat_fhc__numu2018.root";
186 const std::string dpred_name_rhc= dpred_name+
"pred_numuconcat_rhc__numu2018.root";
189 if(anaName==
"fhc") dpred_name_one = dpred_name_fhc;
191 std::vector< const ISyst* >
systs;
193 int systSize = systs.size();
198 for(
int quant = 1; quant < totquant+1; quant++){
200 if(quant<5) predictions.push_back(
new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
201 else predictions.push_back(
new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
205 std::vector<TH1*> hDataUnscaled;
206 std::vector<TH1*> hData;
207 std::vector<TH1*> hNC;
208 std::vector<TH1*> hCC;
209 std::vector<TH1*> hNuMu;
210 std::vector<TH1*> hNoOsc;
211 std::vector<TH1*> hPred;
212 std::vector<TH1*> hBkg;
213 std::vector<TH1*> hAllBkg;
214 std::vector<TH1*> hWS;
218 hData.push_back(s_realdata[0]->ToTH1(s_realdata[0]->
POT()));
219 hDataUnscaled.push_back(s_realdata[0]->ToTH1(s_realdata[0]->
POT()));
223 hBkg.push_back(predictions[0]->Predict(calc).ToTH1(
pot));
224 hAllBkg.push_back(predictions[0]->Predict(calc).ToTH1(
pot));
225 hPred.push_back(predictions[0]->Predict(calc).ToTH1(
pot));
226 hNoOsc.push_back(predictions[0]->Predict(calcNoOsc).ToTH1(
pot));
229 for(
int quantId = 0; quantId <
totquant; quantId++){
231 hData.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->
POT()));
232 hDataUnscaled.push_back(s_realdata[quantId]->ToTH1(s_realdata[0]->
POT()));
236 hBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(
pot));
237 hAllBkg.push_back(predictions[quantId]->Predict(calc).ToTH1(
pot));
238 hPred.push_back(predictions[quantId]->Predict(calc).ToTH1(
pot));
239 hNoOsc.push_back(predictions[quantId]->Predict(calcNoOsc).ToTH1(
pot));
242 hDataUnscaled[0] ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->
POT()));
243 hData[0] ->Add(s_realdata[quantId]->ToTH1(s_realdata[0]->
POT()));
247 hBkg[0] ->Add(predictions[quantId]->Predict(calc).ToTH1(
pot));
248 hAllBkg[0]->Add(predictions[quantId]->Predict(calc).ToTH1(
pot));
249 hPred[0] ->Add(predictions[quantId]->Predict(calc).ToTH1(
pot));
250 hNoOsc[0] ->Add(predictions[quantId]->Predict(calcNoOsc).ToTH1(
pot));
257 for(
int quantId = 0; quantId < totquant+1; quantId++){
258 hBkg[quantId]->Add(hNuMu[quantId], -1);
259 hAllBkg[quantId]->Add(hNuMu[quantId], -1);
260 hAllBkg[quantId]->Add(hCosmics[quantId]);
261 hWS[quantId]->Add(hAllBkg[quantId]);
262 hPred[quantId]->Add(hCosmics[quantId]);
266 std::cout <<
"\n\nIntegrated events = " << hData[0]->Integral() <<
"\n\n" <<
std::endl;
270 for(
int quantId = 0; quantId < totquant+1; quantId++){
271 hData[quantId] ->Scale(0.1,
"width");
272 hCosmics[quantId]->Scale(0.1,
"width");
273 hNC[quantId] ->Scale(0.1,
"width");
274 hCC[quantId] ->Scale(0.1,
"width");
275 hBkg[quantId] ->Scale(0.1,
"width");
276 hAllBkg[quantId] ->Scale(0.1,
"width");
277 hNuMu[quantId] ->Scale(0.1,
"width");
278 hPred[quantId] ->Scale(0.1,
"width");
279 hNoOsc[quantId] ->Scale(0.1,
"width");
280 hWS[quantId] ->Scale(0.1,
"width");
285 std::vector<TGraph*> grData;
286 for(
int quantId = 0; quantId < totquant+1; quantId++){
293 std::vector<Spectrum> sPrediction;
294 for(
int quantId = 0; quantId < totquant+1; quantId++){
298 sPrediction.push_back(prediction);
303 sPrediction.push_back(prediction);
310 std::vector<Spectrum> ups[totquant+1], dns[totquant+1];
314 for(
int quantId = 0; quantId < totquant+1; quantId++){
318 for(
const ISyst* syst: systs){
324 TH1 * htempup = predictions[0]->PredictSyst(calc, shifts).ToTH1(
pot);
325 for(
int i = 1;
i <
totquant;
i++){htempup->Add(predictions[
i]->PredictSyst(calc, shifts).ToTH1(
pot));}
326 htempup->Scale(0.1,
"width");
327 htempup->Add(hCosmics[0]);
329 ups[quantId].push_back(specup);
332 TH1 * htempdn = predictions[0]->PredictSyst(calc, shifts).ToTH1(
pot);
333 for(
int i = 1;
i <
totquant;
i++){htempdn->Add(predictions[
i]->PredictSyst(calc, shifts).ToTH1(
pot));}
334 htempdn->Scale(0.1,
"width");
335 htempdn->Add(hCosmics[0]);
337 dns[quantId].push_back(specdn);
343 TH1* htempup = predictions[quantId-1]->PredictSyst(calc, shifts).ToTH1(
pot);
344 htempup->Scale(0.1,
"width");
345 htempup->Add(hCosmics[quantId]);
347 ups[quantId].push_back(specup);
349 TH1 * htempdn = predictions[quantId-1]->PredictSyst(calc, shifts).ToTH1(
pot);
350 htempdn->Scale(0.1,
"width");
351 htempdn->Add(hCosmics[quantId]);
353 dns[quantId].push_back(specdn);
362 for(
int quantId = 0; quantId< 1; quantId++){
373 grData[quantId]->SetMarkerSize(1.);
377 for(
int quantId = 1; quantId< totquant+1; quantId++){
388 grData[quantId]->SetMarkerSize(1.);
393 TLatex* xTitle =
new TLatex(0.5, 0.03,
"Reconstructed Energy (GeV)");
395 xTitle->SetTextAlign(22);
397 TLatex* yTitle =
new TLatex(0.02, 0.5,
"Events / 0.1 GeV");
399 yTitle->SetTextAlign(22);
400 yTitle->SetTextAngle(90.);
403 const int cx = 960, cy = 800;
408 TString
outName = dout_name +
"fd_dataprediction_spectra__" + anaName;
411 TCanvas *
canvas =
new TCanvas(
"canvas",
"canvas", cx, cy);
414 hPred[0]->GetYaxis()->SetTitle(
"Events / 0.1 GeV");
415 hPred[0]->GetYaxis()->SetRangeUser(0.0,
maxy_pad0);
416 hPred[0]->GetXaxis()->SetLabelSize(hPred[0]->
GetXaxis()->GetLabelSize());
417 hPred[0]->GetYaxis()->SetLabelSize(hPred[0]->
GetYaxis()->GetLabelSize());
418 hPred[0]->Draw(
"hist");
421 else hPred[quantId]->Draw(
"hist same");
422 hWS[quantId]->Draw(
"hist same");
423 hAllBkg[quantId]->Draw(
"hist same");
424 hCosmics[quantId]->Draw(
"hist same");
425 hPred[quantId]->Draw(
"hist same");
426 grData[quantId]->Draw(
"p0 same");
427 grData[quantId]->Draw(
"p0 same");
428 TLegend *
legend =
new TLegend(0.58,0.54,0.89,0.89);
429 legend->AddEntry(hData[quantId],
"FD Data",
"lep");
430 legend->AddEntry(hPred[quantId],
"Prediction",
"l");
431 TLegendEntry *
entry=legend->AddEntry(
"error",
"1-#sigma syst. range",
"bf");
432 entry->SetFillStyle(1001);
435 if(anaName==
"fhc")legend->AddEntry(hWS[quantId],
"Wrong Sign:#bar{#nu}_{#mu}CC",
"bf");
436 else legend->AddEntry(hWS[quantId],
"Wrong Sign:#nu_{#mu}CC",
"bf");
437 legend->AddEntry(hAllBkg[quantId],
"Total bkg.",
"bf");
438 legend->AddEntry(hCosmics[quantId],
"Cosmic bkg.",
"bf");
439 legend->SetTextSize(0.04);
440 legend->SetBorderSize(0);
441 legend->SetFillStyle(0);
450 canvas->SaveAs(outName +
"__"+calcName+
".png");
451 canvas->SaveAs(outName +
"__"+calcName+
".pdf");
452 canvas->SaveAs(outName +
"__"+calcName+
".root");
453 canvas->SaveAs(outName +
"__"+calcName+
".C");
454 file.open (outName +
"__"+calcName+
".txt");
455 file<<
"FD spectra showing data, prediction at best fit with systematics band about the best fit with systematics, and beam and cosmic backgrounds for the NuMI beam in " << anaName <<
" mode.";
469 TString outName = dout_name +
"fd_dataprediction_spectra__" + anaName;
471 for(
int quantId:{3,4,1,2}){
472 if(quantId==1) pad1->cd();
473 if(quantId==2) pad2->cd();
474 if(quantId==3) pad3->cd();
475 if(quantId==4) pad4->cd();
479 hPred[quantId]->GetXaxis()->SetLabelSize(hPred[quantId]->
GetXaxis()->GetLabelSize()/
paperscale);
480 hPred[quantId]->GetYaxis()->SetLabelSize(hPred[quantId]->
GetYaxis()->GetLabelSize()/
paperscale);
481 hPred[quantId]->Draw(
"hist");
483 TGaxis* gax1 =
new TGaxis(0., 0., 4., 0., 0.,4.,4,
"");
484 gax1->SetLabelSize(hPred[quantId]->
GetXaxis()->GetLabelSize());
485 gax1->SetLabelOffset(hPred[quantId]->
GetXaxis()->GetLabelOffset());
486 gax1->SetLabelFont(hPred[quantId]->
GetXaxis()->GetLabelFont());
487 hPred[quantId]->GetXaxis()->SetLabelSize(0.);
488 hPred[quantId]->Draw(
"hist");
493 else hPred[quantId]->Draw(
"hist same");
494 hWS[quantId]->Draw(
"hist same");
495 hAllBkg[quantId]->Draw(
"hist same");
496 hCosmics[quantId]->Draw(
"hist same");
497 hPred[quantId]->Draw(
"hist same");
498 grData[quantId]->Draw(
"p0 same");
502 TLegend *legend =
new TLegend(0.09,0.73,0.27,0.89,
"",
"NDC");
503 legend->AddEntry(hData[quantId],
"FD Data",
"lep");
504 legend->AddEntry(hPred[quantId],
"Prediction",
"l");
505 TLegendEntry *entry=legend->AddEntry(
"error",
"1-#sigma syst. range",
"bf");
506 entry->SetFillStyle(1001);
509 if(anaName==
"fhc")legend->AddEntry(hWS[quantId],
"Wrong Sign:#bar{#nu}_{#mu}CC",
"bf");
510 else legend->AddEntry(hWS[quantId],
"Wrong Sign:#nu_{#mu}CC",
"bf");
511 legend->AddEntry(hAllBkg[quantId],
"Total bkg.",
"bf");
512 legend->AddEntry(hCosmics[quantId],
"Cosmic bkg.",
"bf");
513 legend->SetTextSize(0.025);
514 legend->SetBorderSize(0);
515 legend->SetFillStyle(0);
532 canvas->SaveAs(outName +
"_quant14__"+calcName+
".png");
533 canvas->SaveAs(outName +
"_quant14__"+calcName+
".pdf");
534 canvas->SaveAs(outName +
"_quant14__"+calcName+
".root");
535 canvas->SaveAs(outName +
"_quant14__"+calcName+
".C");
536 file.open (outName +
"_quant14__"+calcName+
".txt");
537 file<<
"FD spectra showing data, prediction at best fit with systematics band about the best fit with systematics, and beam and cosmic backgrounds for the NuMI beam in " << anaName <<
" mode, separated per quartiles.";
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
Pass neutrinos through unchanged.
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
correl_xv GetYaxis() -> SetDecimals()
Loads shifted spectra from files.
Simple record of shifts applied to systematic parameters.
const double livetime_fhc
void CenterTitles(TH1 *histo)
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.
Representation of a spectrum in any variable, with associated POT.
void PimpHist(TH1 *hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize)
Pimp histogram once and for all.
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
Charged-current interactions.
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
correl_xv GetXaxis() -> SetDecimals()
virtual _IOscCalc< T > * Copy() const override
void HornSettings(int quantile, double &pot, double &livetime, Sign::Sign_t &wrongs, double &maxy)
Oscillation probability calculators.
std::vector< double > POT
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float maxy, bool newaxis)
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
All neutrinos, any flavor.
void SetShift(const ISyst *syst, double shift, bool force=false)
void FillWithDimColor(TH1 *h, bool usealpha, float dim)