SystsOscParamAna.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
13 
14 #include "OscLib/OscCalcSterile.h"
15 
16 #include "TArrow.h"
17 #include "TCanvas.h"
18 #include "TDirectory.h"
19 #include "TFile.h"
20 #include "TH1.h"
21 #include "TKey.h"
22 #include "TLatex.h"
23 #include "TLegend.h"
24 #include "TLine.h"
25 #include "TPad.h"
26 
27 using namespace ana;
28 
31 
34 
37 
38 void PlotSyst(TH1* nom, TH1* hp1, TH1* hm1,
39  TDirectory* out, FILE* text, strings strs, bool use_max);
40 void PlotSyst(TH1* h, TH1* hp1, TH1* hm1,
41  TDirectory* out, FILE* text, strings strs);
42 
43 // After running a fit, reset the values in a calculator to preset defaults
45 
47 {
48  TH1::AddDirectory(0);
49 
50  // A map of systematic labels
51  std::map<std::string, std::string> shiftlabels;
52  shiftlabels["theta13"] = "#theta_{13}";
53  shiftlabels["delta13"] = "#delta_{13}";
54  shiftlabels["Deltam232"] = "#Deltam^2_{32}";
55 
56  std::map<std::string, PredictionCombinePeriods*> preds;
57 
58  // Set up path to file with filled CAFAna objects and for output of analysis
59  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
60  std::string filenm = "SystsOscParam";
61 
62  std::string loadLocation = "/nova/ana/steriles/Ana01/FitSystEffects.root";
63  std::string saveLocation = filenm + "Ana.root";
64  std::string textLocation = filenm + "Ana.txt";
65 
66  // Load filled objects from file
67  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
68  TDirectory* tmpL = gDirectory;
69  TDirectory* loadDir = gDirectory;
70  loadDir->cd((loadLocation + ":/pred").c_str());
72  rootL->Close();
73 
74  // Set up oscillation calculator that uses default 3 flavor parameters
76  calc4fnom->SetNFlavors(4);
77  ResetAngles(calc4fnom);
78 
79  osc::OscCalcSterile* calc4fplus = DefaultSterileCalc(4);
80  calc4fplus->SetNFlavors(4);
81  ResetAngles(calc4fplus);
82 
83  osc::OscCalcSterile* calc4fminus = DefaultSterileCalc(4);
84  calc4fminus->SetNFlavors(4);
85  ResetAngles(calc4fminus);
86 
87 
88 
89  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
90  FILE* textF;
91  textF = fopen(textLocation.c_str(), "w");
92 
93  strings strs;
94  strs.fComponent = "NC";
95  strs.fDet = "EX";
96  strs.fPOT = StringFromDouble(NCSCALE);
97  strs.fSample = "Ana01";
98  strs.fSystType = "OscParam";
99  strs.fXLabel = "Calorimetric Energy (GeV)";
100 
101  InitializeSystText(textF, strs);
102 
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};
107 
108  TH1* hNomNC_predSt = GetNC(&pred, NCSCALE, calc4fnom);
109  TH1* hNomBG_predSt = GetBG(&pred, NCSCALE, calc4fnom);
110 
111  for(const auto& shifts : shiftlabels) {
112  ResetAngles(calc4fplus);
113 
114  strs.fSystS = shifts.first;
115  strs.fSystL = shifts.second;
116 
117  if(shifts.first.compare("theta13") == 0) {
118  calc4fplus->SetAngle(1, 3, calc4fplus->GetAngle(1, 3) + parammap[shifts.first].second);
119  }
120  else if(shifts.first.compare("delta13") == 0) {
121  calc4fplus->SetDelta(1, 3, calc4fplus->GetDelta(1,3) + parammap[shifts.first].second);
122  }
123  else {
124  calc4fplus->SetDm(3, calc4fplus->GetDm(3) + parammap[shifts.first].second);
125  }
126 
127  TH1* hSp1NC_predSt = GetNC(&pred, NCSCALE, calc4fplus);
128  TH1* hSp1BG_predSt = GetBG(&pred, NCSCALE, calc4fplus);
129 
130  ResetAngles(calc4fminus);
131 
132  if(shifts.first.compare("theta13") == 0) {
133  calc4fminus->SetAngle(1, 3, calc4fminus->GetAngle(1, 3) - parammap[shifts.first].second);
134  }
135  else if(shifts.first.compare("delta13") == 0) {
136  calc4fminus->SetDelta(1, 3, calc4fminus->GetDelta(1,3) - parammap[shifts.first].second);
137  }
138  else {
139  calc4fminus->SetDm(3, calc4fminus->GetDm(3) - parammap[shifts.first].second);
140  }
141 
142  TH1* hSm1NC_predSt = GetNC(&pred, NCSCALE, calc4fminus);
143  TH1* hSm1BG_predSt = GetBG(&pred, NCSCALE, calc4fminus);
144 
145  strs.fComponent = "NC";
146  strs.fDet = "EX";
147  std::cout << "Plotting signal" << std::endl;
148  PlotSyst(hNomNC_predSt, hSp1NC_predSt, hSm1NC_predSt,
149  rootF, textF, strs, false);
150 
151  strs.fComponent = "BG";
152  PlotSyst(hNomBG_predSt, hSp1BG_predSt, hSm1BG_predSt,
153  rootF, textF, strs, false);
154  }
155 
156  strs.fComponent = "NC";
157  PlotSystBand(hNomNC_predSt, rootF, textF, strs);
158 
159  strs.fComponent = "BG";
160  PlotSystBand(hNomBG_predSt, rootF, textF, strs);
161 
162  fclose(textF); // Text file is done now, close it
163  rootF->Close(); // Close the output root file
164 }
165 
166 //------------------------------------------------------------------------------
168 {
169  Spectrum s = specs->PredictComponent(
171  );
172  TH1* ret = s.ToTH1(POT);
173  std::cout << "NC: " << ret->Integral(0,-1) << std::endl;
174  ZeroError(ret);
175  return ret;
176 }
177 
178 //------------------------------------------------------------------------------
180 {
181  Spectrum s = specs->PredictComponent(
183  );
184  TH1* ret = s.ToTH1(POT);
185  std::cout << " Bkg: " << ret->Integral(0,-1) << std::endl;
186  ZeroError(ret);
187  return ret;
188 }
189 
190 //------------------------------------------------------------------------------
193 {
194  // This pattern must match that in NusSystFromHist::LoadHists
195  std::string ret = "h" + chan + "_" + det + "_" + syst + "_" + sigma;
196  return ret;
197 }
198 
199 //------------------------------------------------------------------------------
202 {
203  std::string ret = chan + " " + det + " " + syst + " Systematic, "
204  + sigma + "#sigma";
205  return ret;
206 }
207 
208 //----------------------------------------------------------------------------
209 /// Plot FD syst shifted and nominal spectra
210 /// Should be able to change the input to IPrediction and make this redundant
211 /// This is for plotting an on/off style shift
212 void PlotSyst(TH1* nom, TH1* hp1, TH1* hm1,
213  TDirectory* out, FILE* text, strings strs, bool use_max)
214 {
215  AddErrorInQuadrature(nom, hp1, hm1, use_max);
216  PlotSyst(nom, hp1, hm1, out, text, strs);
217  return;
218 }
219 
220 //----------------------------------------------------------------------------
221 /// Plot a nominal spectrum and a shifted spectra
222 /// This is for a systematic with multiple sigma levels
223 /// \param h The nominal spectrum
224 /// \param hp1 The shifted spectrum, +1 sigma
225 /// \param hm1 The shifted spectrum, -1 sigma
226 /// \param out The directory/file for root output
227 /// \param text The file for text output
228 /// \param strs Contains strings for pretty plots
229 void PlotSyst(TH1* h, TH1* hp1, TH1* hm1,
230  TDirectory* out, FILE* text, strings strs)
231 {
232  std::string xlabel = strs.fXLabel;
233  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
234  std::string title = strs.fDet + " " + strs.fComponent +
235  " Error for " + strs.fSystL + " Systematic";
236  std::string fullTitle = ";" + xlabel + ";" + ylabel;
237 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
238  std::string name = "c" + strs.fComponent + strs.fDet +
239  strs.fSystS + strs.fSample;
240  std::string ylabelRat = "Shifted / Nominal";
241  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
242 
243  TH1* hRat = (TH1*)h->Clone();
244  TH1* hp1Rat = (TH1*)hp1->Clone();
245  TH1* hm1Rat = (TH1*)hm1->Clone();
246 
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.);
255 
256  if(i != 0 && i != nBins+1) {
257  evts_diff += std::max(std::abs(hp1->GetBinContent(i) - h->GetBinContent(i)),
258  std::abs(hm1->GetBinContent(i) - h->GetBinContent(i)));
259  }
260  }
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);
264 
265  fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
266  strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
267  percent_diff, percent_diff_p, percent_diff_m);
268 
269  hp1Rat->Divide(h);
270  hm1Rat->Divide(h);
271 
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()));
275 
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()));
279 
280  if(maxval < minval) { std::swap(maxval, minval); }
281 
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()));
285 
286  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
287  for(int i = 1; i <= nBins; ++i) {
288  if(hp1Rat->GetBinContent(i) > 0.) {
289  minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
290  }
291  if(hm1Rat->GetBinContent(i) > 0.) {
292  minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
293  }
294  }
295 
296  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
297  minvalRat = 1 - std::abs(maxvalRat - 1.);
298  maxvalRat = 1 + std::abs(maxvalRat - 1.);
299  }
300  else {
301  minvalRat = 1 - std::abs(1. - minvalRat);
302  maxvalRat = 1 + std::abs(1. - minvalRat);
303  }
304 
305  CenterTitles(h);
306  h->SetLineColor(kBlack);
307  h->SetMaximum(1.1*maxval);
308  h->SetMinimum(0);
309  h->SetTitle(fullTitle.c_str());
310 
311  CenterTitles(hp1);
312  hp1->SetLineColor(kRed);
313  hp1->SetLineStyle(2);
314  hp1->SetMaximum(1.1*maxval);
315  hp1->SetMinimum(0);
316  hp1->SetTitle(fullTitle.c_str());
317 
318  CenterTitles(hm1);
319  hm1->SetLineColor(kRed);
320  hm1->SetLineStyle(2);
321  hm1->SetMaximum(1.1*maxval);
322  hm1->SetMinimum(0);
323  hm1->SetTitle(fullTitle.c_str());
324 
325  CenterTitles(hRat);
326  hRat->SetLineColor(kBlack);
327  hRat->SetMaximum(1.1*maxvalRat);
328  hRat->SetMinimum(0.9*minvalRat);
329  hRat->SetTitle(fullTitleRat.c_str());
330 
331  CenterTitles(hp1Rat);
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());
337 
338  CenterTitles(hm1Rat);
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());
344 
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");
350 
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);
354  pSpecs->Draw();
355  pRatio->Draw();
356 
357  pSpecs->cd();
358  h->Draw("hist");
359  hp1->Draw("hist same");
360  hm1->Draw("hist same");
361  leg->Draw();
362  Simulation();
363 
364  pRatio->cd();
365  hRat->Draw("hist");
366  hp1Rat->Draw("hist same");
367  hm1Rat->Draw("hist same");
368 
369  c->Update();
370  out->WriteTObject(c);
371 
372  hRat->SetName( MakeName(strs.fComponent, strs.fDet, strs.fSystType, "0").c_str() );
373  hRat->SetTitle(MakeTitle(strs.fComponent, strs.fDet, strs.fSystType, "0").c_str() );
374  hp1Rat->SetName( MakeName (strs.fComponent, strs.fDet, strs.fSystType, "+1").c_str() );
375  hp1Rat->SetTitle(MakeTitle(strs.fComponent, strs.fDet, strs.fSystType, "+1").c_str() );
376  hm1Rat->SetName( MakeName (strs.fComponent, strs.fDet, strs.fSystType, "-1").c_str() );
377  hm1Rat->SetTitle(MakeTitle(strs.fComponent, strs.fDet, strs.fSystType, "-1").c_str() );
378  ZeroError(hRat);
379  ZeroError(hp1Rat);
380  ZeroError(hm1Rat);
381  out->WriteTObject(hRat);
382  out->WriteTObject(hp1Rat);
383  out->WriteTObject(hm1Rat);
384  return;
385 }
386 
387 //------------------------------------------------------------------------------
389 {
391  calc->SetDm(4, 0.5);
392  calc->SetAngle(2, 4, 0.17453293); // 10 degrees
393  calc->SetAngle(3, 4, 0.61086524); // 35 degrees
394 
395  return;
396 }
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
int nBins
Definition: plotROC.py:16
std::string fDet
Definition: PPFXHelper.h:104
double GetDelta(int i, int j) const
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
std::string fComponent
Definition: PPFXHelper.h:103
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.
Definition: Spectrum.cxx:209
void SystsOscParamAna()
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)
Definition: PPFXHelper.h:662
void ZeroError(TH1 *nom)
Set the errors of the input histogram to 0.
Definition: PPFXHelper.h:149
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
osc::OscCalcDumb calc
std::string MakeName(std::string chan, std::string det, std::string syst, std::string sigma)
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
fclose(fg1)
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
void AddErrorInQuadrature(TH1 *nom, TH1 *up, bool use_max)
Definition: PPFXHelper.h:163
void InitializeSystText(FILE *text, strings strs)
Print some initial text about a systematic – the systematic type and sample.
Definition: PPFXHelper.h:263
std::string fSystL
Definition: PPFXHelper.h:109
std::string fPOT
Definition: PPFXHelper.h:105
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
void PlotSystBand(TH1 *h, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:968
std::string fXLabel
Definition: PPFXHelper.h:110
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetDm(int i, double dm)
void ResetAngles(osc::OscCalcSterile *calc)
std::string fSystType
Definition: PPFXHelper.h:107
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
double GetDm(int i) const
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::string fSystS
Definition: PPFXHelper.h:108
Sum MC predictions from different periods scaled according to data POT targets.
static double NCSCALE
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
std::string StringFromDouble(double pot)
Definition: PPFXHelper.h:128
double GetAngle(int i, int j) const
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
std::string fSample
Definition: PPFXHelper.h:106