SystsBirksAna.C
Go to the documentation of this file.
4 
5 using namespace ana;
6 
7 void PlotSysts(TH1* nom, TH1* hsB, TH1* hsC,
8  TDirectory* rootOUT, FILE* textOFS,
9  strings strs);
10 
12 {
13  TH1::AddDirectory(0);
14 
15  // A map of systematic labels
16  std::map<std::string, std::string> shiftlabels;
17  shiftlabels["BirksB"] = "BirksB";
18  shiftlabels["BirksC"] = "BirksC";
19 
20  // Define a map of samples and cuts/axes to set them up
21  // Each sample, essentially, is a way to set up the analysis
22  std::string labelRecoE = "Calorimetric Energy (GeV)";
23  std::vector<std::string> cut_samples;
24  cut_samples.push_back("Ana01");
25 
26  // Set up maps to the decompositions/predictions (each flavor component)
27  // Nominal maps are indexed purely by sample label
28  // Shifted maps are indexed by sample label, systematic label, then sigma of the shift
29  std::map<std::string, IDecomp*> decomps_nominal;
30  std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted;
31  std::map<std::string, PredictionNoExtrap*> predsNE_nominal;
32  std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted;
33  std::map<std::string, PredictionSterile*> predsSt_nominal;
34  std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted;
35 
36  std::string folder = "/nova/ana/steriles/Ana01/Systematics/";
37  std::string filenm = "SystsBirks";
38 
39  std::string loadLocation = folder + filenm + ".root";
40  std::string saveLocation = folder + filenm + "Ana.root";
41  std::string textLocation = folder + filenm + "Ana.txt";
42 
43  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
44  LoadMaps(rootL, &decomps_nominal, &decomps_shifted);
45  LoadMaps(rootL, &predsNE_nominal, &predsNE_shifted);
46  LoadMaps(rootL, &predsSt_nominal, &predsSt_shifted);
47  rootL->Close();
48 
49  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
50  FILE* textF;
51  textF = fopen(textLocation.c_str(), "w");
52 
53  for(const auto& sample : cut_samples) {
54  strings strs;
55  strs.fComponent = "NC";
56  strs.fDet = "ND";
58  strs.fSample = sample;
59  strs.fSystType = "Birks";
60  strs.fXLabel = labelRecoE;
61 
62  InitializeSystText(textF, strs);
63 
64  TH1* hNomNC_decomp = GetNC(decomps_nominal[sample], NCSCALE);
65  TH1* hNomBG_decomp = GetBG(decomps_nominal[sample], NCSCALE);
66 
67  TH1* hNomNC_predNE = GetNC(predsNE_nominal[sample], NCSCALE);
68  TH1* hNomBG_predNE = GetBG(predsNE_nominal[sample], NCSCALE);
69 
70  TH1* hNomNC_predSt = GetNC(predsSt_nominal[sample], NCSCALE);
71  TH1* hNomBG_predSt = GetBG(predsSt_nominal[sample], NCSCALE);
72 
73  for(const auto& shifts : shiftlabels) {
74  strs.fSystS = shifts.first;
75  strs.fSystL = shifts.second;
76 
77  strs.fComponent = "NC";
78  strs.fDet = "ND";
79  PlotSyst(hNomNC_decomp,
80  decomps_shifted[sample][shifts.first][1],
81  rootF, textF, strs, NCSCALE, true, true);
82 
83  strs.fComponent = "BG";
84  PlotSyst(hNomBG_decomp,
85  decomps_shifted[sample][shifts.first][1],
86  rootF, textF, strs, NCSCALE, false, true);
87 
88  strs.fComponent = "NC";
89  strs.fDet = "FD";
90  PlotSyst(hNomNC_predNE,
91  predsNE_shifted[sample][shifts.first][1],
92  rootF, textF, strs, NCSCALE, true, true);
93 
94  strs.fComponent = "BG";
95  PlotSyst(hNomBG_predNE,
96  predsNE_shifted[sample][shifts.first][1],
97  rootF, textF, strs, NCSCALE, false, true);
98 
99  strs.fComponent = "NC";
100  strs.fDet = "EX";
101  PlotSyst(hNomNC_predSt,
102  predsSt_shifted[sample][shifts.first][1],
103  rootF, textF, strs, NCSCALE, true, true);
104 
105  strs.fComponent = "BG";
106  PlotSyst(hNomBG_predSt,
107  predsSt_shifted[sample][shifts.first][1],
108  rootF, textF, strs, NCSCALE, false, true);
109  }
110 
111  TH1* hSftNC_B_decomp = GetNC(decomps_shifted[sample]["BirksB"][1], NCSCALE);
112  TH1* hSftBG_B_decomp = GetBG(decomps_shifted[sample]["BirksB"][1], NCSCALE);
113  TH1* hSftNC_C_decomp = GetNC(decomps_shifted[sample]["BirksC"][1], NCSCALE);
114  TH1* hSftBG_C_decomp = GetBG(decomps_shifted[sample]["BirksC"][1], NCSCALE);
115 
116  TH1* hSftNC_B_predNE = GetNC(predsNE_shifted[sample]["BirksB"][1], NCSCALE);
117  TH1* hSftBG_B_predNE = GetBG(predsNE_shifted[sample]["BirksB"][1], NCSCALE);
118  TH1* hSftNC_C_predNE = GetNC(predsNE_shifted[sample]["BirksC"][1], NCSCALE);
119  TH1* hSftBG_C_predNE = GetBG(predsNE_shifted[sample]["BirksC"][1], NCSCALE);
120 
121  TH1* hSftNC_B_predSt = GetNC(predsSt_shifted[sample]["BirksB"][1], NCSCALE);
122  TH1* hSftBG_B_predSt = GetBG(predsSt_shifted[sample]["BirksB"][1], NCSCALE);
123  TH1* hSftNC_C_predSt = GetNC(predsSt_shifted[sample]["BirksC"][1], NCSCALE);
124  TH1* hSftBG_C_predSt = GetBG(predsSt_shifted[sample]["BirksC"][1], NCSCALE);
125 
126  strs.fComponent = "NC";
127  strs.fDet = "ND";
128  PlotSysts(hNomNC_decomp, hSftNC_B_decomp, hSftNC_C_decomp,
129  rootF, textF, strs);
130 
131  strs.fComponent = "BG";
132  PlotSysts(hNomBG_decomp, hSftBG_B_decomp, hSftBG_C_decomp,
133  rootF, textF, strs);
134 
135  strs.fComponent = "NC";
136  strs.fDet = "FD";
137  PlotSysts(hNomNC_predNE, hSftNC_B_predNE, hSftNC_C_predNE,
138  rootF, textF, strs);
139 
140  strs.fComponent = "BG";
141  PlotSysts(hNomBG_predNE, hSftBG_B_predNE, hSftBG_C_predNE,
142  rootF, textF, strs);
143 
144  strs.fComponent = "NC";
145  strs.fDet = "EX";
146  PlotSysts(hNomNC_predSt, hSftNC_B_predSt, hSftNC_C_predSt,
147  rootF, textF, strs);
148 
149  strs.fComponent = "BG";
150  PlotSysts(hNomBG_predSt, hSftBG_B_predSt, hSftBG_C_predSt,
151  rootF, textF, strs);
152  }
153 
154  fclose(textF);
155  rootF->Close();
156 
157  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
158 }
159 
160 void PlotSysts(TH1* nom, TH1* hsB, TH1* hsC,
161  TDirectory* rootOUT, FILE* textOFS,
162  strings strs)
163 {
164  double nNom = nom->Integral();
165  double nSBB = hsB->Integral();
166  double nSBC = hsC->Integral();
167 
168  ZeroError(nom);
169  ZeroError(hsB);
170  ZeroError(hsC);
171 
172  fprintf(textOFS, "\nEvent counts for %.2s %.2s %s Systs:\n",
173  strs.fDet.c_str(), strs.fComponent.c_str(), strs.fSystType.c_str());
174  fprintf(textOFS, "NC: %6.2f, BirksB: %6.2f, BirksC: %6.2f\n",
175  nNom, nSBB, nSBC);
176  fprintf(textOFS, "BirksB Diff: %6.2f, BirksC Diff: %6.2f",
177  100.*std::abs(nSBB/nNom - 1.), 100.*std::abs(nSBC/nNom - 1.));
178 
179  if(std::abs(nSBB - nNom) > std::abs(nSBC - nNom)) {
180  AddErrorInQuadrature(nom, hsB, true);
181  }
182  else {
183  AddErrorInQuadrature(nom, hsC, true);
184  }
185 
186  std::string xlabel = strs.fXLabel;
187  std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
188  std::string title = strs.fDet + " " + strs.fComponent +
189  " Errors for All " + strs.fSystType + " Systematics";
190 // std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
191  std::string fullTitle = ";" + xlabel + ";" + ylabel;
192  std::string name = "c" + strs.fComponent + strs.fDet +
193  strs.fSystType + "SystsAll" + strs.fSample;
194  std::string ylabelRat = "Shifted / Nominal";
195  std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
196 
197  TH1* hRat = (TH1*)nom->Clone();
198  TH1* hSBBRat = (TH1*)hsB->Clone();
199  TH1* hSBCRat = (TH1*)hsC->Clone();
200 
201  int nBins = nom->GetNbinsX();
202  for(int i = 0; i <= nBins+1; ++i) {
203  hRat->SetBinContent(i, 1.);
204  }
205  hSBBRat->Divide(nom);
206  hSBCRat->Divide(nom);
207 
208  double maxval = nom->GetBinContent(nom->GetMaximumBin());
209  maxval = std::max(maxval, hsB->GetBinContent(hsB->GetMaximumBin()));
210  maxval = std::max(maxval, hsC->GetBinContent(hsC->GetMaximumBin()));
211 
212  double minval = nom->GetBinContent(nom->GetMinimumBin());
213  minval = std::min(minval, hsB->GetBinContent(hsB->GetMinimumBin()));
214  minval = std::min(minval, hsC->GetBinContent(hsC->GetMinimumBin()));
215 
216  if(maxval < minval) { std::swap(maxval, minval); }
217 
218  double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
219  maxvalRat = std::max(maxvalRat, hSBBRat->GetBinContent(hSBCRat->GetMaximumBin()));
220  maxvalRat = std::max(maxvalRat, hSBCRat->GetBinContent(hSBCRat->GetMaximumBin()));
221 
222  double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
223  for(int i = 1; i <= nBins; ++i) {
224  if(hSBBRat->GetBinContent(i) > 0.) {
225  minvalRat = std::min(minvalRat, hSBBRat->GetBinContent(i));
226  }
227  if(hSBCRat->GetBinContent(i) > 0.) {
228  minvalRat = std::min(minvalRat, hSBCRat->GetBinContent(i));
229  }
230  }
231 
232  if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
233  minvalRat = 1 - std::abs(maxvalRat - 1.);
234  maxvalRat = 1 + std::abs(maxvalRat - 1.);
235  }
236  else {
237  minvalRat = 1 - std::abs(1. - minvalRat);
238  maxvalRat = 1 + std::abs(1. - minvalRat);
239  }
240 
241  CenterTitles(nom);
242  nom->SetLineColor(kBlack);
243  nom->SetMaximum(1.1*maxval);
244  nom->SetMinimum(0);
245  nom->SetTitle(fullTitle.c_str());
246 
247  CenterTitles(hsB);
248  hsB->SetLineColor(kRed);
249  hsB->SetLineStyle(2);
250  hsB->SetMaximum(1.1*maxval);
251  hsB->SetMinimum(0);
252  hsB->SetTitle(fullTitle.c_str());
253 
254  CenterTitles(hsC);
255  hsC->SetLineColor(kBlue);
256  hsC->SetLineStyle(2);
257  hsC->SetMaximum(1.1*maxval);
258  hsC->SetMinimum(0);
259  hsC->SetTitle(fullTitle.c_str());
260 
261  CenterTitles(hRat);
262  hRat->SetLineColor(kBlack);
263  hRat->SetMaximum(1.1*maxvalRat);
264  hRat->SetMinimum(0.9*minvalRat);
265  hRat->SetTitle(fullTitleRat.c_str());
266 
267  CenterTitles(hSBBRat);
268  hSBBRat->SetLineColor(kRed);
269  hSBBRat->SetLineStyle(2);
270  hSBBRat->SetMaximum(1.1*maxvalRat);
271  hSBBRat->SetMinimum(0.9*minvalRat);
272  hSBBRat->SetTitle(fullTitleRat.c_str());
273 
274  CenterTitles(hSBCRat);
275  hSBCRat->SetLineColor(kBlue);
276  hSBCRat->SetLineStyle(2);
277  hSBCRat->SetMaximum(1.1*maxvalRat);
278  hSBCRat->SetMinimum(0.9*minvalRat);
279  hSBCRat->SetTitle(fullTitleRat.c_str());
280 
281  double xL = 0.6, xR = 0.85;
282  double yB = 0.6, yT = 0.85;
283  TLegend* leg = new TLegend(xL, yB, xR, yT);
284  leg->AddEntry(nom, "Nominal", "l");
285  leg->AddEntry(hsB, "BirksB", "l");
286  leg->AddEntry(hsC, "BirksC", "l");
287 
288  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
289  TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
290  TPad* pRatio = new TPad("pRatio", "", 0., 0., 1., 0.375);
291  pSpecs->Draw();
292  pRatio->Draw();
293 
294  pSpecs->cd();
295  nom->Draw("hist");
296  hsB->Draw("hist same");
297  hsC->Draw("hist same");
298  leg->Draw();
299  Simulation();
300 
301  pRatio->cd();
302  hRat->Draw("hist");
303  hSBBRat->Draw("hist same");
304  hSBCRat->Draw("hist same");
305 
306  c->Update();
307  rootOUT->WriteTObject(c);
308  return;
309 }
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
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string fComponent
Definition: PPFXHelper.h:103
void PlotSyst(TH1 *h, TH1 *hp1, TH1 *hm1, TH1 *hp2, TH1 *hm2, TDirectory *out, FILE *text, strings strs)
Definition: PPFXHelper.h:662
void SystsBirksAna()
Definition: SystsBirksAna.C:11
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
float abs(float number)
Definition: d0nt_math.hpp:39
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TH1 * GetNC(IDecomp *specs, double POT)
Definition: PPFXHelper.h:211
fclose(fg1)
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 PlotSysts(TH1 *nom, TH1 *hsB, TH1 *hsC, TDirectory *rootOUT, FILE *textOFS, strings strs)
OStream cout
Definition: OStream.cxx:6
std::string fXLabel
Definition: PPFXHelper.h:110
std::string fSystType
Definition: PPFXHelper.h:107
std::string fSystS
Definition: PPFXHelper.h:108
void LoadMaps(TDirectory *dir, std::map< std::string, IDecomp * > *nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > *shifted)
Definition: PPFXHelper.h:273
static double NCSCALE
T min(const caf::Proxy< T > &a, T b)
TH1 * GetBG(IDecomp *specs, double POT)
Definition: PPFXHelper.h:236
enum BeamMode kBlue
std::string StringFromDouble(double pot)
Definition: PPFXHelper.h:128
A helper structure to contain a group of string for plotting.
Definition: PPFXHelper.h:101
std::string fSample
Definition: PPFXHelper.h:106
enum BeamMode string