DrawSystShifts.C
Go to the documentation of this file.
1 // DrawSystShifts by J. Hewes (jhewes15@fnal.gov)
2 /// Draw systematic shifts in each detector bin and compare to covariance matrix
3 
4 #include "CAFAna/Core/ISyst.h"
9 
11 
12 #include "TSystem.h"
13 #include "TCanvas.h"
14 #include "TGraph.h"
15 #include "TRandom3.h"
16 #include "TText.h"
17 
18 #include <fstream>
19 #include <algorithm>
20 
21 using namespace ana;
22 
23 float Gaus(float mean, float var, float x) {
24  float norm = 1.0/sqrt(2.0*M_PI*var);
25  return norm * exp(-1*pow(x-mean,2)/(2*var));
26 }
27 
28 void DrawBins(IPrediction* pred, const ISyst* syst, TMatrixD mx,
29  TVectorD vec, size_t nbins, std::ofstream& f, std::string det);
30 
31 void DrawPulls(IPrediction* pred, const ISyst* syst, size_t nbins,
32  std::ofstream& f, std::string det);
33 
34 size_t nPoints = 1001;
35 std::vector<std::string> one_d = { "AbsCalibShape", "Cherenkov", "XSecTune" };
36 
38 
39  DontAddDirectory guard;
40 
41  // Define exposures
42  double potND = kAna2018SensitivityFHCNDPOT;
43  double potFD = kAna2018FHCPOT;
44 
45  // Set up latex output
46  gSystem->MakeDirectory("tex_pulls");
47  gSystem->MakeDirectory("tex_pulls/plots");
48  gSystem->CopyFile(Form("%s/nus/Nus20/Plots/header_nus.tex",
49  FindCAFAnaDir().c_str()), "tex_pulls/header_nus.tex");
50 
51  // Write header
52  std::ofstream maintex;
53  maintex.open("tex_pulls/main.tex");
54  maintex << "\n\\input{header_nus.tex}\n\n";
55  maintex << "\\newpage\n\n";
56 
57  // Load predictions
58  IPrediction* predNDSysts = LoadPrediction(kNusFHCNearDet, "systs", kUseCVMFS);
59  IPrediction* predFDSysts = LoadPrediction(kNusFHCFarDet, "systs", kUseCVMFS);
60  IPrediction* predNDGENIE = LoadPrediction(kNusFHCNearDet, "genie", kUseCVMFS);
61  IPrediction* predFDGENIE = LoadPrediction(kNusFHCFarDet, "genie", kUseCVMFS);
62  auto calc = DefaultSterileCalc(4);
64 
65  // Spectrum vectors
68  TH1D* hND = predNDSysts->Predict(calc).ToTH1(potND);
69  TH1D* hFD = predFDSysts->Predict(calc).ToTH1(potND);
70  for (size_t i = 0; i < size_t(kNCNDBinning.NBins()); ++i)
71  vND(i) = hND->GetBinContent(i+1);
72  for (size_t i = 0; i < size_t(kNCFDBinning.NBins()); ++i)
73  vFD(i) = hFD->GetBinContent(i+1);
74  HistCache::Delete(hND);
75  HistCache::Delete(hFD);
76 
77  // Load covariance matrices
78  std::map<std::string, CovarianceMatrix*> mxND;
79  for (auto itMx : GetMatrices({kNusFHCNearDet}, kUseCVMFS))
80  mxND[itMx->GetName()] = itMx;
81  std::map<std::string, CovarianceMatrix*> mxFD;
82  for (auto itMx : GetMatrices({kNusFHCFarDet}, kUseCVMFS))
83  mxFD[itMx->GetName()] = itMx;
84 
85  // Near detector systematics
86  maintex << "\\section{Near Detector}\n\n";
87  for (auto syst : LoadSysts(kNusFHCNearDet, kUseCVMFS)) {
88  maintex << "\\subsection{" << GetSystBaseName(syst->ShortName()) << "}\n\n";
89  TMatrixD mx(mxND[GetSystBaseName(syst->ShortName())]->GetCovMxRelative(vND));
90  DrawBins(predNDSysts, potND, syst, mx, vND, kNCNDBinning.NBins(), maintex, "nd");
91  DrawPulls(predNDSysts, potND, syst, kNCNDBinning.NBins(), maintex, "nd");
92  }
93  for (auto syst : getAllXsecSysts_2018()) {
94  std::string systname = syst->ShortName();
95  std::replace(systname.begin(), systname.end(), '_', ' ');
96  maintex << "\\subsection{" << systname << "}\n\n";
97  TMatrixD mx(mxND[syst->ShortName()]->GetCovMxRelative(vND));
98  DrawBins(predNDGENIE, potND, syst, mx, vND, kNCNDBinning.NBins(), maintex, "nd");
99  DrawPulls(predNDGENIE, potND, syst, kNCNDBinning.NBins(), maintex, "nd");
100  }
101 
102  // Far detector systematics
103  maintex << "\\section{Far Detector}\n\n";
104  for (auto syst : LoadSysts(kNusFHCFarDet, kUseCVMFS)) {
105  maintex << "\\subsection{" << GetSystBaseName(syst->ShortName()) << "}\n\n";
106  TMatrixD mx(mxFD[GetSystBaseName(syst->ShortName())]->GetCovMxRelative(vFD));
107  DrawBins(predFDSysts, potFD, syst, mx, vFD, kNCFDBinning.NBins(), maintex, "fd");
108  DrawPulls(predFDSysts, potFD, syst, kNCFDBinning.NBins(), maintex, "fd");
109  }
110  for (auto syst : getAllXsecSysts_2018()) {
111  std::string systname = syst->ShortName();
112  std::replace(systname.begin(), systname.end(), '_', ' ');
113  maintex << "\\subsection{" << systname << "}\n\n";
114  TMatrixD mx(mxFD[syst->ShortName()]->GetCovMxRelative(vFD));
115  DrawBins(predFDGENIE, potFD, syst, mx, vFD, kNCFDBinning.NBins(), maintex, "fd");
116  DrawPulls(predFDGENIE, potFD, syst, kNCFDBinning.NBins(), maintex, "fd");
117  }
118 
119  maintex << "\\end{document}\n";
120  maintex.close();
121 
122 }
123 
124 void DrawBins(IPrediction* pred, double pot, const ISyst* syst, TMatrixD mx,
125  TVectorD vec, size_t nbins, std::ofstream& f, std::string det) {
126 
127  std::cout << "Processing syst " << syst->ShortName() << std::endl;
128  f << "\\begin{figure}[h!]\n \\centering\n";
129  auto calc = DefaultSterileCalc(4);
131 
132  // Make output dirs
133  std::string systname = syst->ShortName().substr(0,7) == "nus_fhc" ?
134  GetSystBaseName(syst->ShortName()) : syst->ShortName();
135 
136  std::vector<TH1D*> bins(nbins);
137  bool allNoEffect = true;
138  bool allOneDirectional = true;
139  std::vector<bool> oneDirectional(nbins, false);
140  for (size_t i = 0; i < nbins; ++i)
141  bins[i] = new TH1D("", "", 60, 1-(3*sqrt(mx(i,i))), 1+(3*sqrt(mx(i,i))));
142 
143  TH1D* hNom = pred->Predict(calc).ToTH1();
144 
145  // Loop over bins
146  TRandom3 rnd(0);
147  for (size_t i = 0; i <= 1e5; ++i) {
148  SystShifts shifts;
149  if (std::any_of(one_d.begin(), one_d.end(), [systname](std::string name)
150  { return name == systname; })) shifts.SetShift(syst, fabs(rnd.Gaus()));
151  else shifts.SetShift(syst, rnd.Gaus());
152  TH1D* h = pred->PredictSyst(calc, shifts).ToTH1(pot);
153  for (size_t j = 0; j < nbins; ++j)
154  bins[j]->Fill(h->GetBinContent(j+1)/vec(j));
155  HistCache::Delete(h);
156  }
157 
158  // Here we want to classify each systematic in each bin.
159  // Here's how this is going to go. First we want to classify
160  // as one-directional, or not. This should be super easy.
161  // In each bin, we add up bins below and above centre, and if
162  // either integral is zero, it's one-directional. That's
163  // all there is to it.
164  for (size_t i = 0; i < nbins; ++i) {
165  if (bins[i]->Integral(1, 29) != 0 or bins[i]->Integral(32, 60) != 0)
166  allNoEffect = false;
167  if (bins[i]->Integral(1, 29) == 0 or bins[i]->Integral(32, 60) == 0)
168  oneDirectional[i] = true;
169  else allOneDirectional = false;
170  }
171 
172  // Calculate area fraction. For this we literally just loop over a bunch
173  // of slices in x, and take the minimum height at each one. then we add
174  // up a little area slice.
175  std::vector<float> nonGaussianScore(nbins, 0);
176  for (size_t i = 0; i < nbins; ++i) {
177  float xrange = 1 - bins[i]->GetXaxis()->GetBinLowEdge(1);
178  float norm = 0.1*bins[i]->Integral()*sqrt(mx(i,i));
179  for (int j = -999; j <= +999; ++j) {
180  float x = (0.001 * float(j) * xrange)+1;
181  float yHist = bins[i]->GetBinContent(bins[i]->GetXaxis()->FindBin(x));
182  float yGaus = norm * Gaus(1, mx(i,i), x);
183  if (oneDirectional[i]) yGaus *= 2;
184  float yMin = yHist < yGaus ? yHist : yGaus;
185  nonGaussianScore[i] += (yMin * 0.1 * xrange) / norm;
186  }
187  }
188 
189  // Save histograms to file
190  TCanvas c("", "", 900, 1600);
191  std::vector<TPad*> pads;
192  std::vector<TGraph*> graphs;
193  pads.push_back(new TPad("classification", "", 0, 0.95, 1, 0.999));
194  pads.back()->cd();
195  std::string classification;
196  if (allNoEffect) classification = "No effect";
197  else if (allOneDirectional) classification = "One-directional";
198  else classification = "Two-directional";
199  float meanScore = 0.;
200  for (size_t i = 0; i < nbins; ++i) meanScore += nonGaussianScore[i] / float(nbins);
201  TText* t = new TText(0.5, 0.5, Form("%s, score = %.1f%%", classification.c_str(), meanScore));
202  t->SetTextSize(0.6);
203  t->SetTextAlign(22);
204  t->Draw();
205  for (size_t i = 0; i < nbins; ++i) {
206  float x1 = 0.25*(float(i%4));
207  float y1 = 0.95-(0.19*(1.+std::floor(0.25*float(i))));
208  pads.push_back(new TPad(Form("bin%zu", i), "", x1, y1, x1+0.25, y1+0.18));
209  pads.back()->cd();
210  std::vector<float> x(nPoints);
211  std::vector<float> y(nPoints);
212  float binWidth = 6.0*sqrt(mx(i,i))/float(nPoints-1);
213  float norm = 0.1*bins[i]->Integral()*sqrt(mx(i,i));
214  for (size_t j = 0; j < nPoints; ++j) {
215  x[j] = 1 + (-3.0*sqrt(mx(i,i))) + (float(j)*binWidth);
216  y[j] = norm * Gaus(1,mx(i,i),x[j]);
217  if (oneDirectional[i]) y[j] *= 2;
218  }
219  graphs.push_back(new TGraph(nPoints, &x[0], &y[0]));
220  bins[i]->Draw("hist");
221  bins[i]->SetTitle(Form("Bin %zu - %.2f < E (GeV) < %.2f", i+1,
222  hNom->GetXaxis()->GetBinLowEdge(i+1), hNom->GetXaxis()->GetBinUpEdge(i+1)));
223  graphs.back()->Draw("c same");
224  graphs.back()->SetLineColor(kRed-7);
225  graphs.back()->SetLineWidth(3);
226  TText* tmp = new TText(0.85, 0.85, Form("%.1f%%", nonGaussianScore[i]));
227  tmp->SetTextAlign(33);
228  tmp->SetTextSize(0.1);
229  tmp->SetNDC(true);
230  tmp->Draw();
231  // pads.back()->Update();
232  }
233 
234  c.cd();
235  for (auto pad : pads) pad->Draw();
236 
237  c.SaveAs(Form("tex_pulls/plots/gaus_%s_%s.pdf",
238  det.c_str(), systname.c_str()));
239  f << " \\includegraphics[height=0.8\\textheight]{"
240  << Form("plots/gaus_%s_%s.pdf", det.c_str(), systname.c_str())
241  << "}\n\\end{figure}\n\n\\clearpage\n\n";
242 
243  HistCache::Delete(hNom);
244  delete calc;
245 }
246 
247 void DrawPulls(IPrediction* pred, double pot, const ISyst* syst, size_t nbins,
248  std::ofstream& f, std::string det) {
249 
250  f << "\\begin{figure}[h!]\n \\centering\n";
251  auto calc = DefaultSterileCalc(4);
253 
254  // Make output dirs
255  std::string systname = syst->ShortName().substr(0,7) == "nus_fhc" ?
256  GetSystBaseName(syst->ShortName()) : syst->ShortName();
257 
258  std::vector<float> x(nPoints);
259  std::vector<std::vector<float>> y(nbins, std::vector<float>(nPoints));
260 
261  TH1D* hNom = pred->Predict(calc).ToTH1(pot);
262 
263  for (size_t i = 0; i < nPoints; ++i) {
264  SystShifts shifts;
265  x[i] = (6./float(nPoints-1))*(float(i)-(0.5*(nPoints-1)));
266  shifts.SetShift(syst, x[i]);
267  TH1D* h = pred->PredictSyst(calc, shifts).ToTH1(pot);
268  for (size_t j = 0; j < nbins; ++j)
269  y[j][i] = h->GetBinContent(j+1);
270  HistCache::Delete(h);
271  }
272 
273  // Save histograms to file
274  TCanvas c("", "", 900, 1600);
275  std::vector<TPad*> pads;
276  std::vector<TGraph*> graphs;
277  for (size_t i = 0; i < nbins; ++i) {
278  float x1 = 0.25*(float(i%4));
279  float y1 = 1-(0.2*(1.+std::floor(0.25*float(i))));
280  pads.push_back(new TPad(Form("bin%zu", i), "", x1, y1, x1+0.25, y1+0.1999));
281  pads.back()->cd();
282  graphs.push_back(new TGraph(nPoints, &x[0], &y[i][0]));
283  graphs.back()->Draw("al");
284  graphs.back()->SetTitle(Form("Bin %zu - %.2f < E (GeV) < %.2f", i+1,
285  hNom->GetXaxis()->GetBinLowEdge(i+1), hNom->GetXaxis()->GetBinUpEdge(i+1)));
286  graphs.back()->SetLineWidth(3);
287  }
288 
289  c.cd();
290  for (auto pad : pads) pad->Draw();
291 
292  c.SaveAs(Form("tex_pulls/plots/pull_%s_%s.pdf",
293  det.c_str(), systname.c_str()));
294  f << " \\includegraphics[height=0.8\\textheight]{"
295  << Form("plots/pull_%s_%s.pdf", det.c_str(), systname.c_str())
296  << "}\n\\end{figure}\n\n\\clearpage\n\n";
297 
298  delete calc;
299 }
300 
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string GetSystBaseName(std::string name)
Get the last part of a systematic&#39;s ShortName.
Definition: Utilities.h:232
size_t nPoints
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Float_t y1[n_points_granero]
Definition: compare.C:5
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:148
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
double Integral(const Spectrum &s, const double pot, int cvnregion)
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
const bool kUseCVMFS
Float_t tmp
Definition: plot.C:36
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
osc::OscCalcDumb calc
std::string FindCAFAnaDir()
Definition: Utilities.cxx:208
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
#define M_PI
Definition: SbMath.h:34
void DrawBins(IPrediction *pred, const ISyst *syst, TMatrixD mx, TVectorD vec, size_t nbins, std::ofstream &f, std::string det)
const int nbins
Definition: cellShifts.C:15
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
std::vector< std::string > one_d
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
correl_xv GetXaxis() -> SetDecimals()
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
#define pot
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:49
const double j
Definition: BetheBloch.cxx:29
Eigen::VectorXd vec
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
float Gaus(float mean, float var, float x)
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Binning kNCNDBinning
NC custom binning.
Definition: Binning.cxx:104
Float_t norm
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
int NBins() const
Definition: Binning.h:29
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
const Binning kNCFDBinning
Definition: Binning.cxx:105
std::vector< const ISyst * > LoadSysts(covmx::Sample sample)
Get systematics for a given sample.
Definition: Utilities.h:524
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void DrawPulls(IPrediction *pred, const ISyst *syst, size_t nbins, std::ofstream &f, std::string det)
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
void DrawSystShifts()