8 #include "NuXAna/Analysis/NusSystsMaker.h" 60 { v.push_back(v.back()+
b); }
65 std::vector<double>
v = { 0. };
66 for (
size_t i = 0;
i < 6; ++
i)
68 for (
size_t i = 0;
i < 2; ++
i)
94 std::vector<double>
v = { 0. };
111 std::vector<double>
v = { 0. };
122 std::vector<double>
v = { 0. };
124 for (
int i = 0;
i < 67; ++
i)
140 if (detector ==
"nd") use_nd =
true;
141 else if (detector ==
"fd") use_fd =
true;
142 else if (detector ==
"both") {
146 else throw std::runtime_error(
"Detector option \"" + detector +
"\" not understood.");
150 std::vector<Binning>
ret;
215 gStyle->SetCanvasDefW(1200);
216 gStyle->SetCanvasDefH(800);
217 gStyle->SetPadTopMargin(0.14);
218 gStyle->SetPadBottomMargin(0.14);
219 gStyle->SetPadLeftMargin(0.14);
220 gStyle->SetPadRightMargin(0.14);
230 if (detector !=
"neardet" && detector !=
"fardet")
231 throw std::runtime_error(
"Detector string must be either \"neardet\" or \"fardet\"!");
233 std::string dname = detector ==
"neardet" ?
"nd" :
"fd";
237 if (not rhc and detector ==
"neardet" and syst_type ==
"none") {
240 return LoadFromFile<PredictionInterp>(filename.c_str(),
241 "pred_interp_nus_fhc_neardet_none").
release();
246 +
"_" + syst_type +
".root";
252 if (!path)
assert(
false and
253 "NUSDATA_NUS18_PRED environment variable not set! Did you forget to set up nusdata?");
259 return LoadFromFile<PredictionInterp>(
filepath,
"pred_interp").
release();
268 bool fhc = opt.Contains(
"fhc", TString::kIgnoreCase);
269 bool rhc = opt.Contains(
"rhc", TString::kIgnoreCase);
272 bool syst = opt.Contains(
"syst", TString::kIgnoreCase);
275 filename <<
"pred_nus18_";
276 rhc ? filename <<
"rhc" : filename <<
"fhc";
277 filename <<
"_extrap";
278 if (syst) filename <<
"_systs";
282 if (!path)
throw std::runtime_error(
283 "NUSDATA_NUS18_PRED environment variable not set! Did you forget to set up nusdata?");
284 else std::cout <<
"Loading prediction from " << path <<
"/" << filename.str() <<
std::endl;
289 return LoadFromFile<IPrediction>(
filepath,
"FDNus18ExtrapPred").
release();
300 if (!path)
throw std::runtime_error(
301 "NUSDATA_NUS18_COVMX environment variable not set! Did you forget to set up nusdata?");
305 std::cout <<
"Loading covariance matrix " << name <<
" from " << filepath <<
std::endl;
317 if (detector !=
"nd" && detector !=
"fd")
throw std::runtime_error(
318 "Detector string must be \"nd\" or \"fd\" - " + detector +
" not understood!");
319 if (!path)
throw std::runtime_error(
320 "NUSDATA_NUS18_WEIGHTS environment variable not set! Did you forget to set up nusdata?");
365 if (!path)
throw std::runtime_error(
366 "NUSDATA_NUS18_FAKEDATA environment variable not set! Did you forget to set up nusdata?");
367 else std::cout <<
"Loading fake data spectrum universe_" << universe <<
" from " 368 << path <<
"/" << filename <<
".root" <<
std::endl;
373 TFile*
f =
new TFile(filepath.c_str(),
"read");
374 TH1D*
h = (TH1D*)f->Get(histname.c_str());
386 if (rhc) path =
"/nova/app/users/wallbank/test_ndosc/NuXAna/macros/Nus18_neutrino/box_opening/";
388 if (rhc)
return *LoadFromFile<Spectrum>(path +
"/fout_nus18_box_opening_rhc.root",
"spec_nus_vise_numi").
release();
390 if (nd)
return *LoadFromFile<Spectrum>(path +
"/files/fhc_NDDataMCLoadOut.root",
"RecoEPreselPtpData").
release();
391 else return *LoadFromFile<Spectrum>(path +
"/files/fout_nus18_box_opening_fhc.root",
"spec_nus_vise_numi").
release();
401 if (!path)
throw std::runtime_error(
402 "NUSDATA_NUS18_COSMIC environment variable not set! Did you forget to set up nusdata?");
403 else std::cout <<
"Loading cosmic data from " << path <<
"/nus18_cosmic_background.root" <<
std::endl;
405 TFile
f(Form(
"%s/nus18_cosmic_background.root", path));
406 return (TH1D*)f.Get(rhc?
"rhcnumi" :
"fhcnumi");
415 bool fhc = opt.Contains(
"fhc", TString::kIgnoreCase);
416 bool rhc = opt.Contains(
"rhc", TString::kIgnoreCase);
420 if (!path)
throw std::runtime_error(
421 "NUSDATA_NUS18_COSMIC environment variable not set! Did you forget to set up nusdata?");
422 else std::cout <<
"Loading cosmic data from " << path <<
"/nus18_cosmic_background.root" <<
std::endl;
424 TFile
f(Form(
"%s/nus18_cosmic_background.root", path));
425 return (TH1D*)f.Get(rhc?
"rhcnumi" :
"fhcnumi");
442 bool fhc = opt.Contains(
"fhc", TString::kIgnoreCase);
443 bool rhc = opt.Contains(
"rhc", TString::kIgnoreCase);
472 throw std::runtime_error(
"NUSDATA_NUS18_SYSTS environment variable not set!");
476 std::vector<const ISyst*>
systs;
478 TFile*
f =
new TFile(Form(
"%s/Nus18ISysts.root", path),
"READ");
479 TIter
next(f->GetListOfKeys());
481 while ((key = (TKey*)
next()))
494 char* ppfx_env =
std::getenv(
"NUSDATA_NUS18_WEIGHTS");
495 if (ppfx_env ==
nullptr)
throw std::runtime_error(
496 "NUSDATA_NUS18_WEIGHTS environment variable not set! Did you forget to set up nusdata?");
503 PredictionConcat* GetDefaultSimulation(bool rhc, std::string detector, std::string syst_type)
508 if (detector == "nd") use_nd = true;
509 else if (detector == "fd") use_fd = true;
510 else if (detector == "both") {
514 else throw std::runtime_error("Detector string must be \"nd\", \"fd\" or \"both\".");
517 if (syst_type != "none" && syst_type != "xsec" && syst_type != "nonxsec")
518 throw std::runtime_error("Syst type string must be \"none\", \"xsec\" or \"nonxsec\".");
521 double nd_pot = rhc ? kAna2018SensitivityRHCNDPOT : kAna2018SensitivityFHCNDPOT;
522 double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
523 double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
524 std::string polarity = rhc? "RHC" : "FHC";
526 // Vectors for predictions, POT & livetimes, systematic prefixes
527 std::vector<IPrediction*> predictions;
528 std::vector<double> pots;
529 std::vector<double> livetimes;
530 std::vector<std::string> samples;
531 std::vector<std::vector<const ISyst*>> systs;
535 samples.push_back("Near detector");
536 predictions.push_back(LoadPrediction("neardet", rhc, syst_type));
537 pots.push_back(nd_pot);
538 livetimes.push_back(0);
539 systs.push_back(GetNus18Systs(rhc, "nd", syst_type));
540 // syst_types.push_back(polarity + "_ND");
543 samples.push_back("Far detector");
544 predictions.push_back(LoadPrediction("fardet", rhc, syst_type));
545 pots.push_back(fd_pot);
546 livetimes.push_back(fd_livetime);
547 systs.push_back(GetNus18Systs(rhc, "fd", syst_type));
548 // syst_types.push_back(polarity + "_FD");
551 PredictionConcat* sim = new PredictionConcat(predictions, pots, livetimes);
552 sim->SetSysts(new SystConcat(samples, {}, systs));
553 sim->SetBinning(GetNus18Binning(rhc, detector));
564 if (detector ==
"nd") use_nd =
true;
565 else if (detector ==
"fd") use_fd =
true;
566 else if (detector ==
"both") {
570 else throw std::runtime_error(
"Detector string \"" + detector +
"\" not recognised!");
582 std::map<const ISyst*, CovarianceMatrix*>
m, std::vector<IPrediction*> preds,
585 if (m.size() == 0)
throw std::runtime_error(
"Matrix map is empty!");
586 else if (m.size() == 1)
throw std::runtime_error(
"Only one matrix in map; cannot add!");
588 gen->PredictCovMx(preds, pot, calc);
638 std::vector<TH1D*>
GetCosmics(std::vector<covmx::Sample> samples,
639 PredictionConcat*
sim)
643 std::map<covmx::Polarity, std::string>
pname;
647 std::map<covmx::Quantile, std::string>
qname;
653 assert(samples.size() == sim->GetNSamples() and
654 "Sample vector and concatenated prediction have different size!");
659 std::vector<TFile*>
f(4,
nullptr);
667 f[0] = TFile::Open(filepath.c_str(),
"read");
675 f[1] = TFile::Open(filepath.c_str(),
"read");
684 f[2] = TFile::Open(filepath.c_str(),
"read");
693 f[3] = TFile::Open(filepath.c_str(),
"read");
700 std::vector<TH1D*>
ret(samples.size(),
nullptr);
702 for (
size_t i = 0;
i < samples.size(); ++
i) {
718 ret[
i] = (TH1D*)f[0]->
Get(hname.c_str());
724 std::string sname =
"cosmic_spect_" + pname[samples[
i].polarity];
733 std::string hname =
"cosmics_" + qname[samples[
i].quantile];
734 ret[
i] = (TH1D*)f[2]->
Get(hname.c_str());
741 std::string hname =
"cosmics_" + qname[samples[
i].quantile];
742 ret[
i] = (TH1D*)f[3]->
Get(hname.c_str());
747 assert(
ret[
i] &&
"Found uninitialised cosmic spectrum!");
754 for (TFile*
tmp : f) {
769 assert(
false and
"Only NC FHC is supported right now!");
774 std::string filepath =
"/pnfs/nova/persistent/users/jhewes15/nus19/isysts/isysts_" + detector +
".root";
781 TFile*
f = TFile::Open(filepath.c_str(),
"read");
782 TDirectory*
dir = f->GetDirectory(name.c_str());
784 std::vector<const ISyst*>
systs;
786 TIter
next(f->GetListOfKeys());
788 while ((key = (TKey*)
next())) {
792 for (
auto syst : systs)
804 TLegend*
leg =
new TLegend(x1, y1, x2, y2);
805 leg->SetBorderSize(0);
806 leg->SetFillColor(kWhite);
807 leg->SetFillStyle(0);
808 leg->SetFillStyle(0);
809 leg->SetTextFont(42);
810 leg->SetTextSize(textsize);
818 TLatex*
tex =
new TLatex();
820 tex->SetTextFont(42);
821 tex->SetTextSize(0.037);
822 tex->SetTextAlign(11);
823 tex->DrawLatex(xpos, start,
"Neutrino Beam:");
824 tex->DrawLatex(xpos, start - 1*step,
"8.0#times10^{20} ND POT,");
825 tex->DrawLatex(xpos, start - 2*step,
"8.85#times10^{20} FD POT-equiv.");
826 tex->DrawLatex(xpos, start - 3*step,
"#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
827 tex->DrawLatex(xpos, start - 4*step,
"#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
828 tex->DrawLatex(xpos, start - 5*step,
"#Deltam^{2}_{41} = 0.5 eV^{2}");
835 TLatex*
tex =
new TLatex();
837 tex->SetTextFont(42);
838 tex->SetTextSize(0.037);
839 tex->SetTextAlign(11);
840 tex->DrawLatex(xpos, start,
"Neutrino Beam:");
841 tex->DrawLatex(xpos, start - 1*step,
"8.0#times10^{20} ND POT,");
842 tex->DrawLatex(xpos, start - 2*step,
"8.85#times10^{20} FD POT-equiv.");
843 tex->DrawLatex(xpos, start - 3*step,
"#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
844 tex->DrawLatex(xpos, start - 4*step,
"#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
845 tex->DrawLatex(xpos, start - 5*step, Form(
"#Deltam^{2}_{41} = %s eV^{2}",dm.c_str()));
852 TLatex*
tex =
new TLatex();
854 tex->SetTextFont(42);
855 tex->SetTextSize(0.037);
856 tex->SetTextAlign(11);
857 tex->DrawLatex(xpos, start,
"Neutrino beam:");
858 tex->DrawLatex(xpos, start - 1*step,
"8.0#times10^{20} ND POT,");
859 tex->DrawLatex(xpos, start - 2*step,
"8.85#times10^{20} FD POT-equiv.");
860 tex->DrawLatex(xpos, start - 3*step,
"#Deltam^{2}_{32} = 2.44#times10^{-3} eV^{2}");
861 tex->DrawLatex(xpos, start - 4*step,
"#theta_{13} = 8.3#circ, ^{}#theta_{23} = 48.4#circ");
868 std::vector<const IExperiment*>
ret;
889 std::vector<double>
ret;
900 std::vector<const IFitVar*>
ret;
910 int nbins = h->GetNbinsX();
911 std::vector<double> edges(nbins+1);
913 edges[
i] = h->GetXaxis()->GetBinLowEdge(
i+1);
920 int nbins = m->GetNbinsX();
921 std::vector<double> edges =
GetEdges((TH1*)m);
922 TH1D*
h =
new TH1D(
UniqueName().c_str(),
"",nbins,&edges[0]);
924 h->SetBinContent(
i,1);
925 double var = m->GetBinContent(
i,
i);
926 if (var > 0) h->SetBinError(
i,
sqrt(var));
927 else h->SetBinError(
i,0);
933 std::vector<std::tuple<ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t>>
GetComponents()
935 std::vector<std::tuple<ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t> >
components;
940 for (
auto & flavor: flavors) {
941 for (
auto &
sign: signs) {
952 TLatex *
l =
new TLatex(x, y, text);
954 l->SetTextSize(size);
955 l->SetTextColor(kBlack);
963 if (i==0)
return std::string(
"#nu_{e} #rightarrow #nu_{e}");
964 else if (i==1)
return std::string(
"#bar{#nu}_{e} #rightarrow #bar{#nu}_{e}");
965 else if (i==2)
return std::string(
"#nu_{e} #rightarrow #nu_{#mu}");
966 else if (i==3)
return std::string(
"#bar{#nu}_{e} #rightarrow #bar{#nu}_{#mu}");
967 else if (i==4)
return std::string(
"#nu_{e} #rightarrow #nu_{#tau}");
968 else if (i==5)
return std::string(
"#bar{#nu}_{e} #rightarrow #bar{#nu}_{#tau}");
969 else if (i==6)
return std::string(
"#nu_{#mu} #rightarrow #nu_{e}");
970 else if (i==7)
return std::string(
"#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{e}");
971 else if (i==8)
return std::string(
"#nu_{#mu} #rightarrow #nu_{#mu}");
972 else if (i==9)
return std::string(
"#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#mu}");
973 else if (i==10)
return std::string(
"#nu_{#mu} #rightarrow #nu_{#tau}");
974 else if (i==11)
return std::string(
"#bar{#nu}_{#mu} #rightarrow #bar{#nu}_{#tau}");
983 double last_edge = sim->GetBinningConcatenated().Edges().back();
984 double max = 13. * last_edge;
988 for (
int i = 0;
i < 13; ++
i) {
989 double pos = ((double)
i + 0.5) * last_edge;
990 TLine* lx =
new TLine(0, pos, max, pos);
991 lx->SetLineColor(kGray+1);
993 TLine* ly =
new TLine(pos, 0, pos, max);
994 ly->SetLineColor(kGray+1);
999 for (
int i = 0;
i <= 13; ++
i) {
1000 double pos =
i * last_edge;
1001 TLine* lx =
new TLine(0, pos, max, pos);
1003 TLine* ly =
new TLine(pos, 0, pos, max);
1008 for (
int i = 0;
i < 13; ++
i) {
1009 double pos = 0.16 + (0.72*(
i/13.));
1011 l1->SetTextAngle(38);
1013 l2->SetTextAlign(31);
1016 h->SetNdivisions(0,
"X");
1017 h->SetNdivisions(0,
"Y");
1023 if (opts.Contains(opt, TString::ECaseCompare::kIgnoreCase))
1031 size_t n_opts = options.size();
1033 std::vector<bool>
vec(n_opts);
1034 for (
size_t i = 0;
i < n_opts; ++
i)
1038 if (std::none_of(vec.begin(), vec.end(), [](
bool v) {
return v; })) {
1040 for (
size_t i = 0; i < n_opts-2; ++
i)
1041 err +=
"\"" + options[i] +
"\", ";
1042 err += options[n_opts-2] +
" or " + options[n_opts-1] +
".";
1043 throw std::runtime_error(err);
1048 for (
size_t i = 0; i < n_opts; ++
i)
1064 size_t pos = optstring.find(key);
1065 if (pos != std::string::npos) {
1066 if (optstring.substr(pos+key.size(), 1) !=
"=") {
1067 std::ostringstream err;
1068 err <<
"The " << key <<
" option must take the form \"" << key <<
"=xxx\".";
1069 throw std::runtime_error(err.str());
1071 size_t start = pos + key.size() + 1;
1072 size_t end = optstring.find(
"_", pos);
1073 if (end != std::string::npos)
1074 val =
std::stod(optstring.substr(start, end - start));
1076 val =
std::stod(optstring.substr(start));
1081 std::ostringstream err;
1082 err <<
"No value for " << key <<
" found in option string.";
1083 throw std::runtime_error(err.str());
1093 std::vector<TH1D*> h_split;
1095 for (
size_t i = 0;
i < b.size(); ++
i) {
1097 n_bins += b[
i].NBins();
1102 if (h->GetNbinsX() !=
int(n_bins))
1103 throw std::runtime_error(
"There are " +
std::to_string(h->GetNbinsX())
1104 +
" bins in merged histogram but a total of " +
std::to_string(n_bins)
1105 +
" bins across the binning vector provided.");
1109 for (
size_t i = 0;
i < b.size(); ++
i) {
1110 for (
int j = 1;
j <= h_split[
i]->GetNbinsX(); ++
j) {
1111 double val = h->GetBinContent(bin_count+
j);
1112 h_split[
i]->SetBinContent(
j, val);
1113 h_split[
i]->SetBinError(
j,
sqrt(val));
1118 bin_count += h_split[
i]->GetNbinsX();
1127 std::vector<std::pair<std::string, osc::IOscCalcAdjustable*>> calcs,
1128 TH1D* cosmic_merged,
1130 std::vector<double>
pot,
1131 std::vector<bool> is_nd,
1138 int cosmic_colour = kGray+2;
1142 if (calcs.empty())
throw std::runtime_error(
"Spectrum vector is empty.");
1143 size_t n_universes = calcs.size();
1144 if (n_universes > colours.size())
1145 throw std::runtime_error(
"Only " +
std::to_string(colours.size())
1146 +
" universes per plot are supported!");
1147 size_t n_detectors = sim->GetNSamples();
1148 if (pot.size() != n_detectors)
1149 throw std::runtime_error(
std::to_string(n_detectors) +
" detectors but only " 1151 if (is_nd.size() != n_detectors)
1152 throw std::runtime_error(
std::to_string(n_detectors) +
" detectors but only " 1157 std::vector<std::vector<Spectrum>>
specs;
1158 for (
size_t i = 0;
i < n_universes; ++
i)
1159 specs.push_back(sim->GetSpectra(calcs[
i].second));
1168 TCanvas*
c =
new TCanvas(
UniqueName().c_str(),
"c", 500 * n_detectors, 800);
1169 c->Divide(n_detectors, 3);
1175 std::vector<float> margin = { 0.1, 0.05 };
1177 gStyle->SetPadLeftMargin(margin[0]);
1178 gStyle->SetPadRightMargin(margin[1]);
1180 gStyle->SetLabelFont(43,
"xyz");
1181 gStyle->SetLabelSize(14,
"xyz");
1183 for (
size_t d = 0;
d < n_detectors; ++
d)
1193 gPad->SetPad(xmin, 0.45, xmax, 1);
1197 c->cd(n_detectors+
d+1);
1198 gPad->SetPad(xmin, 0.25, xmax, 0.45);
1202 c->cd((2*n_detectors)+
d+1);
1203 gPad->SetPad(xmin, 0.05, xmax, 0.25);
1207 TLegend*
l =
new TLegend(0.5, 0.65, 0.85, 0.85);
1212 std::vector<std::vector<TH1D*>>
hists, ratio_hists;
1214 hists.resize(n_universes);
1215 for (std::vector<TH1D*>&
h : hists)
1216 h.resize(n_detectors);
1218 ratio_hists.resize(n_universes);
1219 for (std::vector<TH1D*>&
h : ratio_hists)
1220 h.resize(n_detectors);
1222 std::vector<TGraph*>
graphs;
1226 std::vector<double>
ymax, rmax;
1227 std::vector<bool> rescaled;
1231 std::vector<TH1D*> data_ratio;
1232 data_ratio.resize(n_detectors);
1233 for (
size_t d = 0;
d < n_detectors; ++
d)
1235 data_ratio[
d] = (TH1D*)data[
d]->
Clone();
1236 for (
int b = 1;
b <= data_ratio[
d]->GetNbinsX(); ++
b)
1238 double val = data[
d]->GetBinContent(
b);
1239 double err = data[
d]->GetBinError(
b);
1240 data_ratio[
d]->SetBinContent(
b, 1);
1243 data_ratio[
d]->SetBinError(
b, 0);
1245 data_ratio[
d]->SetBinError(
b, err / val);
1252 for (
size_t d = 0;
d < n_detectors; ++
d)
1256 std::vector<double>
m,
r;
1260 if (is_nd[
d]) data[
d]->Scale(1
e-4);
1262 for (
size_t u = 0;
u < n_universes; ++
u)
1264 hists[
u][
d] = (specs[
u][
d].ToTH1(pot[d]));
1268 if (is_nd[d]) hists[
u][
d]->Scale(1
e-4);
1273 ratio_hists[
u][
d] = (TH1D*)hists[
u][d]->
Clone();
1274 TH1D*
h = ratio_hists[
u][
d];
1275 for (
int b = 1;
b <= h->GetNbinsX(); ++
b)
1277 double val = (hists[
u][
d]->GetBinContent(
b) + cosmic[
d]->GetBinContent(
b)) / data[d]->GetBinContent(
b);
1279 h->SetBinContent(
b, 1);
1280 h->SetBinError(
b, 0);
1283 h->SetBinContent(
b, val);
1284 r.push_back(
fabs(1-val));
1293 int maxdata = data[
d]->GetMaximumBin();
1294 m.push_back(data[d]->GetBinContent(maxdata) + data[d]->GetBinError(maxdata));
1298 for (
int b = 1;
b <= data_ratio[
d]->GetNbinsX(); ++
b)
1299 r.push_back(
fabs(data_ratio[d]->GetBinError(
b)));
1301 ymax.push_back(*std::max_element(m.begin(), m.end()));
1302 rmax.push_back(*std::max_element(r.begin(), r.end()));
1312 for (
size_t d = 0;
d < n_detectors; ++
d) {
1317 TH1D* hc = cosmic[
d];
1318 hc->SetLineColor(cosmic_colour);
1319 hc->SetLineWidth(2);
1321 hc->SetMaximum(1.1 * ymax[
d]);
1322 hc->Draw(
"hist same");
1328 for (
size_t u = 0;
u < n_universes; ++
u)
1332 for (
size_t d = 0;
d < n_detectors; ++
d)
1337 TH1D*
h = hists[
u][
d];
1339 h->SetLineColor(colours[
u]);
1342 h->SetMaximum(1.1 * ymax[d]);
1345 h->SetTitle(
";Energy deposited in detector (GeV);10^{4} events");
1347 h->SetTitle(
";Energy deposited in detector (GeV);Events");
1348 h->GetXaxis()->SetTitleSize(ts);
1349 h->GetYaxis()->SetTitleSize(ts);
1351 h->Draw(
"hist same");
1355 h = ratio_hists[
u][
d];
1356 c->cd(n_detectors+d+1);
1357 h->SetLineColor(colours[u]);
1359 h->SetMinimum(1 - (1.1 * rmax[d]));
1360 h->SetMaximum(1 + (1.1 * rmax[d]));
1362 h->SetTitle(
";;Ratio");
1363 h->GetYaxis()->SetTitleSize(3*ts);
1365 h->Draw(
"hist same");
1369 c->cd((2*n_detectors)+d+1);
1371 std::vector<float>
x,
y;
1372 double l = calcs[
u].second->GetL();
1373 if (is_nd[d]) calcs[
u].second->SetL(1);
1374 else calcs[
u].second->SetL(810);
1375 for (
size_t k = 0; k < 200; ++k)
1377 x.push_back(0.1 * (
float(k) + 0.5));
1378 y.push_back(calcs[u].
second->P(14, 0, x.back()));
1380 calcs[
u].second->SetL(l);
1382 TGraph*
g =
new TGraph(x.size(), &x[0], &y[0]);
1385 g->SetLineColor(colours[u]);
1387 g->SetTitle(
";;P(#nu_{#mu} #rightarrow #nu_{s})");
1388 g->GetYaxis()->SetTitleSize(3*ts);
1391 graphs.push_back(g);
1395 l->AddEntry(hists[
u][0], calcs[u].first.c_str(),
"l");
1399 l->AddEntry(cosmic[0],
"Cosmic background",
"l");
1401 if (data.size() != n_detectors)
1402 throw std::runtime_error(
"Number of data hists does not match number of samples!");
1406 for (
size_t d = 0;
d < n_detectors; ++
d)
1409 data[
d]->SetMarkerSize(2);
1410 data[
d]->Draw(
"e0 same");
1412 c->cd(n_detectors+
d+1);
1413 data_ratio[
d]->SetMarkerSize(2);
1414 data_ratio[
d]->Draw(
"e0 same");
1418 l->AddEntry(data[0],
"Data",
"lep");
1422 dir->WriteTObject(c, name.c_str());
1431 for (
auto g : graphs)
1433 for (
auto h_vec : hists)
1434 for (
auto h : h_vec)
1436 for (
auto h_vec : ratio_hists)
1437 for (
auto h : h_vec)
1439 for (
auto h : data_ratio)
std::vector< double > GetNus18SeedValues(osc::IOscCalcAdjustable *calc)
Seed values for Nus18 marginalisation parameters.
int isinf(const stan::math::var &a)
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
double GetDelta(int i, int j) const
IPrediction * LoadExtrapPrediction(TString opt)
Spectrum LoadCosmicSpectrum(bool rhc=false)
Function to return cosmic spectrum.
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
std::map< std::string, double > xmax
fvar< T > fabs(const fvar< T > &x)
Float_t y1[n_points_granero]
double stod(const std::string &val)
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
void FormatFullCovMxPlot(TH2D *h, PredictionConcat *sim)
bool RunningOnGrid()
Is this a grid (condor) job?
const std::string kFDMCRHCProd4
A simple Gaussian constraint on an arbitrary IFitVar.
const double kNus18Delta13
std::vector< TH1D * > SplitFakeData(TH1D *h, std::vector< Binning > b)
Float_t x1[n_points_granero]
void PlotTextNus18(const double xpos, const double start, const double step)
Add text to plot with delta m^2 41 at 0.5 eV^2.
TH1D * MakeTH1D(const char *name, const char *title, const Binning &bins)
Utility function to avoid need to switch on bins.IsSimple()
void SetDelta(int i, int j, double delta)
std::vector< const ISyst * > LoadISysts()
Function to load ISysts.
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
std::string GetNus18ComponentName(int i)
const int kOfficialNus18FHCColor
TH1D * LoadCosmicHist(bool rhc=false)
Function to load cosmic histogram.
std::string pnfs2xrootd(std::string loc, bool unauth)
void PrintOscParams(osc::OscCalcSterile *calc)
const std::string kFDDataFHC
std::vector< const IExperiment * > GetNus18Constraints()
Gaussian constrains on atmospheric parameters.
void DrawSurfacePoint(PredictionConcat *sim, std::vector< std::pair< std::string, osc::IOscCalcAdjustable * >> calcs, TH1D *cosmic_merged, TH1D *data_merged, std::vector< double > pot, std::vector< bool > is_nd, TDirectory *dir, std::string name)
const double kNus18Th23Sigma
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
void FormatCanvas()
Canvas formatting utility.
std::string FindCAFAnaDir()
void SetAngles(osc::OscCalcSterile *calc)
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Representation of a spectrum in any variable, with associated POT.
int isnan(const stan::math::var &a)
const XML_Char const XML_Char * data
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
const int kOfficialNus18RHCColor
CovarianceMatrix * LoadCovMx(std::string file, std::string name)
Function to load covariance matrix.
Charged-current interactions.
const Binning kNus18EnergyBinning
Energy binnings for Nus18 for nus analysis.
Sum up livetimes from individual cosmic triggers.
TLegend * MakeLegendNus18(double x1, double y1, double x2, double y2, double textsize=0.03)
Default legend.
const std::string kFDDataRHC
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
std::string getenv(std::string const &name)
Spectrum LoadData(bool nd, bool rhc=false)
Function to load NuS data file.
std::vector< std::string > flavors
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
IPrediction * LoadPPFX(std::string detector, int universe)
Function to load PPFX universe simulations.
void AddBin(std::vector< double > &v, double b)
const std::string kFDCosmicRHCProd4
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
static Binning Custom(const std::vector< double > &edges)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
std::vector< float > Spectrum
void PlotTextNus18WithoutDm41(const double xpos, const double start, const double step)
Add text to plot without delta m^2 41 value.
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
std::string GetPPFXWeightsDir()
Function to get location of PPFX weights.
bool CheckOption(TString opts, TString opt)
std::vector< Binning > GetNus18Binning(bool rhc, std::string detector)
TH1D * DiagonalErrors(TH2D *m)
void SetAngle(int i, int j, double th)
const std::string kFDCosmicFHCProd4
void SetDm(int i, double dm)
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
const std::string kFDMCFHCProd4
void PlotTextNus18Dm41(const double xpos, const double start, const double step, std::string dm)
Add text to plot with delta m^2 41 at custom value.
TLatex * MiscText(double x, double y, double size, TString text)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Neutral-current interactions.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::vector< TH1D * > GetCosmics(bool rhc, std::string detector)
Function to get default cosmics.
assert(nhit_max >=nhit_nbins)
Both neutrinos and antineutrinos.
double GetDm(int i) const
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Standard interface to all prediction techniques.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
::xsd::cxx::tree::qname< char, simple_type, uri, ncname > qname
std::string to_string(ModuleType mt)
std::vector< const IFitVar * > GetNus18FitVars()
All neutrinos, any flavor.
Prevent histograms being added to the current directory.
Spectrum LoadFakeData(std::string filename, int universe, double pot, double livetime)
Function to load the unblinded data spectrum.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
std::vector< double > GetEdges(TH1 *h)
std::string UniqueName()
Return a different string each time, for creating histograms.
const double kNus18Dm32Sigma
const double kNus18Delta24
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
double GetAngle(int i, int j) const
const double kAna2018RHCLivetime
size_t GetOptionValue(TString opt, std::string key)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void CombineMatrices(CovarianceMatrix *gen, std::map< const ISyst *, CovarianceMatrix * > m, std::vector< IPrediction * > preds, std::vector< double > pot, osc::IOscCalcAdjustable *calc)
Function to combine covariance matrices.
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.
const FitTheta23Sterile kFitTheta23Sterile