18 #include "TDirectory.h" 38 void PlotSyst(TH1* nom, TH1* hp1, TH1* hm1,
51 std::map<std::string, std::string> shiftlabels;
52 shiftlabels[
"theta13"] =
"#theta_{13}";
53 shiftlabels[
"delta13"] =
"#delta_{13}";
54 shiftlabels[
"Deltam232"] =
"#Deltam^2_{32}";
56 std::map<std::string, PredictionCombinePeriods*> preds;
62 std::string loadLocation =
"/nova/ana/steriles/Ana01/FitSystEffects.root";
67 TFile* rootL =
new TFile(loadLocation.c_str(),
"READ");
68 TDirectory* tmpL = gDirectory;
69 TDirectory* loadDir = gDirectory;
70 loadDir->cd((loadLocation +
":/pred").c_str());
89 TFile* rootF =
new TFile(saveLocation.c_str(),
"RECREATE");
91 textF = fopen(textLocation.c_str(),
"w");
97 strs.fSample =
"Ana01";
98 strs.fSystType =
"OscParam";
99 strs.fXLabel =
"Calorimetric Energy (GeV)";
103 std::map<std::string, std::pair<double, double> > parammap;
104 parammap[
"theta13"] = {8.475, 0.0353};
105 parammap[
"delta13"] = {1.39, 0.325};
106 parammap[
"Deltam232"] = {0.00244, 0.00006};
111 for(
const auto& shifts : shiftlabels) {
114 strs.fSystS = shifts.first;
115 strs.fSystL = shifts.second;
117 if(shifts.first.compare(
"theta13") == 0) {
118 calc4fplus->
SetAngle(1, 3, calc4fplus->
GetAngle(1, 3) + parammap[shifts.first].second);
120 else if(shifts.first.compare(
"delta13") == 0) {
121 calc4fplus->
SetDelta(1, 3, calc4fplus->
GetDelta(1,3) + parammap[shifts.first].second);
124 calc4fplus->
SetDm(3, calc4fplus->
GetDm(3) + parammap[shifts.first].second);
132 if(shifts.first.compare(
"theta13") == 0) {
133 calc4fminus->
SetAngle(1, 3, calc4fminus->
GetAngle(1, 3) - parammap[shifts.first].second);
135 else if(shifts.first.compare(
"delta13") == 0) {
136 calc4fminus->
SetDelta(1, 3, calc4fminus->
GetDelta(1,3) - parammap[shifts.first].second);
139 calc4fminus->
SetDm(3, calc4fminus->
GetDm(3) - parammap[shifts.first].second);
142 TH1* hSm1NC_predSt =
GetNC(&pred,
NCSCALE, calc4fminus);
143 TH1* hSm1BG_predSt =
GetBG(&pred,
NCSCALE, calc4fminus);
145 strs.fComponent =
"NC";
148 PlotSyst(hNomNC_predSt, hSp1NC_predSt, hSm1NC_predSt,
149 rootF, textF, strs,
false);
151 strs.fComponent =
"BG";
152 PlotSyst(hNomBG_predSt, hSp1BG_predSt, hSm1BG_predSt,
153 rootF, textF, strs,
false);
156 strs.fComponent =
"NC";
159 strs.fComponent =
"BG";
216 PlotSyst(nom, hp1, hm1, out, text, strs);
235 " Error for " + strs.
fSystL +
" Systematic";
236 std::string fullTitle =
";" + xlabel +
";" + ylabel;
241 std::string fullTitleRat =
";" + xlabel +
";" + ylabelRat;
243 TH1* hRat = (TH1*)h->Clone();
244 TH1* hp1Rat = (TH1*)hp1->Clone();
245 TH1* hm1Rat = (TH1*)hm1->Clone();
247 int nBins = h->GetNbinsX();
248 double nom_evts = h->Integral(1, nBins);
249 double evts_diff = 0.;
250 double percent_diff = 100.;
251 double percent_diff_p = 100.;
252 double percent_diff_m = 100.;
253 for(
int i = 0;
i <= nBins+1; ++
i) {
254 hRat->SetBinContent(
i, 1.);
256 if(
i != 0 &&
i != nBins+1) {
258 std::abs(hm1->GetBinContent(
i) - h->GetBinContent(
i)));
261 percent_diff *= (evts_diff/nom_evts);
262 percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
263 percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);
265 fprintf(text,
"%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
267 percent_diff, percent_diff_p, percent_diff_m);
272 double maxval = h->GetBinContent(h->GetMaximumBin());
273 maxval =
std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
274 maxval =
std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
276 double minval = h->GetBinContent(h->GetMinimumBin());
277 minval =
std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
278 minval =
std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
280 if(maxval < minval) {
std::swap(maxval, minval); }
282 double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
283 maxvalRat =
std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
284 maxvalRat =
std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
286 double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
288 if(hp1Rat->GetBinContent(
i) > 0.) {
289 minvalRat =
std::min(minvalRat, hp1Rat->GetBinContent(
i));
291 if(hm1Rat->GetBinContent(
i) > 0.) {
292 minvalRat =
std::min(minvalRat, hm1Rat->GetBinContent(
i));
297 minvalRat = 1 -
std::abs(maxvalRat - 1.);
298 maxvalRat = 1 +
std::abs(maxvalRat - 1.);
301 minvalRat = 1 -
std::abs(1. - minvalRat);
302 maxvalRat = 1 +
std::abs(1. - minvalRat);
306 h->SetLineColor(kBlack);
307 h->SetMaximum(1.1*maxval);
309 h->SetTitle(fullTitle.c_str());
312 hp1->SetLineColor(
kRed);
313 hp1->SetLineStyle(2);
314 hp1->SetMaximum(1.1*maxval);
316 hp1->SetTitle(fullTitle.c_str());
319 hm1->SetLineColor(
kRed);
320 hm1->SetLineStyle(2);
321 hm1->SetMaximum(1.1*maxval);
323 hm1->SetTitle(fullTitle.c_str());
326 hRat->SetLineColor(kBlack);
327 hRat->SetMaximum(1.1*maxvalRat);
328 hRat->SetMinimum(0.9*minvalRat);
329 hRat->SetTitle(fullTitleRat.c_str());
332 hp1Rat->SetLineColor(
kRed);
333 hp1Rat->SetLineStyle(2);
334 hp1Rat->SetMaximum(1.1*maxvalRat);
335 hp1Rat->SetMinimum(0.9*minvalRat);
336 hp1Rat->SetTitle(fullTitleRat.c_str());
339 hm1Rat->SetLineColor(
kRed);
340 hm1Rat->SetLineStyle(2);
341 hm1Rat->SetMaximum(1.1*maxvalRat);
342 hm1Rat->SetMinimum(0.9*minvalRat);
343 hm1Rat->SetTitle(fullTitleRat.c_str());
345 double xL = 0.6, xR = 0.85;
346 double yB = 0.6, yT = 0.85;
347 TLegend*
leg =
new TLegend(xL, yB, xR, yT);
348 leg->AddEntry(h,
"Nominal",
"l");
349 leg->AddEntry(hp1,
"+/- 1 #sigma",
"l");
351 TCanvas*
c =
new TCanvas(name.c_str(), title.c_str(), 800, 800);
352 TPad* pSpecs =
new TPad(
"pSpecs",
"", 0., 0.375, 1., 1.);
353 TPad* pRatio =
new TPad(
"pRatio",
"", 0., 0., 1., 0.375);
359 hp1->Draw(
"hist same");
360 hm1->Draw(
"hist same");
366 hp1Rat->Draw(
"hist same");
367 hm1Rat->Draw(
"hist same");
370 out->WriteTObject(c);
381 out->WriteTObject(hRat);
382 out->WriteTObject(hp1Rat);
383 out->WriteTObject(hm1Rat);
T max(const caf::Proxy< T > &a, T b)
double GetDelta(int i, int j) const
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
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
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.
void SetDelta(int i, int j, double delta)
std::string MakeTitle(std::string chan, std::string det, std::string syst, std::string sigma)
Adapt the PMNS_Sterile calculator to standard interface.
void PlotSyst(TH1 *h, TH1 *hp1, TH1 *hm1, TH1 *hp2, TH1 *hm2, TDirectory *out, FILE *text, strings strs)
void ZeroError(TH1 *nom)
Set the errors of the input histogram to 0.
void CenterTitles(TH1 *histo)
std::string MakeName(std::string chan, std::string det, std::string syst, std::string sigma)
Representation of a spectrum in any variable, with associated POT.
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TH1 * GetNC(IDecomp *specs, double POT)
Charged-current interactions.
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
void AddErrorInQuadrature(TH1 *nom, TH1 *up, bool use_max)
void InitializeSystText(FILE *text, strings strs)
Print some initial text about a systematic – the systematic type and sample.
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
std::vector< double > POT
void PlotSystBand(TH1 *h, TDirectory *out, FILE *text, strings strs)
void SetAngle(int i, int j, double th)
void SetDm(int i, double dm)
void ResetAngles(osc::OscCalcSterile *calc)
Neutral-current interactions.
Both neutrinos and antineutrinos.
double GetDm(int i) const
Standard interface to all prediction techniques.
Sum MC predictions from different periods scaled according to data POT targets.
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
TH1 * GetBG(IDecomp *specs, double POT)
std::string StringFromDouble(double pot)
double GetAngle(int i, int j) const
A helper structure to contain a group of string for plotting.