makeEnergyEstimator.C
Go to the documentation of this file.
1 // I found that the energy estimator used by the nuebar oscillation analysis
2 // was not sufficient for this analysis due to training on a sample of
3 // neutrinos with energy less than 4.5 GeV.
4 //
5 // This script marks the first piece of training a new estimator.
6 // Following docdb-26696, we need to normalize the true energy distribution
7 // to avoid biasing the energy estimation toward our peak beam energy.
8 //
9 // Using the spectra obtained from getFlatWeightedFittingSpectra.C
10 // Fit the energy estimator
11 
12 // dddoyle@colostate.edu
13 
14 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Analysis/Style.h"
16 
17 #include "TH1D.h"
18 #include "TH2D.h"
19 #include "TCanvas.h"
20 #include "TVectorD.h"
21 #include "TF2.h"
22 #include "TF1.h"
23 
27 
28 #include <fstream>
29 
30 using namespace ana;
31 
32 void PrintPlot(TH1D* h, std::string name, std::string option);
33 void PrintPlot(TH2D* h, std::string name, std::string option);
34 void SetErrors(TH2D* hFlat_avg_energy, TH2* hFlat);
35 TH2D* MakeAverageTrueEnergySurface(TH2D* hFlat_energy, TH2D* hFlat);
36 
37 Double_t PolN(Double_t *x, Double_t *par);
38 
39 class PolyDef {
40 public:
41  PolyDef(Int_t _nparam, TString _name, TF2* _func)
42  : nparam(_nparam), name(_name), func(_func)
43  {}
44  PolyDef(TString _name, int _degree, bool _cross_terms)
45  : nparam(9),
46  name(_name),
47  degree(_degree),
48  cross_terms(_cross_terms)
49  {
50  func = new TF2(name, PolN, 0, 10, 0, 10, 9);
51  FixTerms();
52  }
53  void FixTerms()
54  {
55  if(degree == 2) {
56  func->FixParameter(4, 0);
57  func->FixParameter(5, 0);
58  func->FixParameter(7, 0);
59  func->FixParameter(8, 0);
60  if(!cross_terms) {
61  func->FixParameter(6, 0);
62  }
63  }
64  else {
65  if(!cross_terms) {
66  func->FixParameter(6, 0);
67  func->FixParameter(7, 0);
68  func->FixParameter(8, 0);
69  }
70  }
71  }
72  TF2* GetFunction() {return func;}
73  TString GetName() {return name;}
74  const Int_t nparam;
75  TString name;
76  TF2* func;
78  int degree;
79 
80 };
81 
82 void PrintFitResults(TH2D* hist, std::ofstream& fout, PolyDef* pol);
83 void PrintErrorSurface(TH2D* hist, PolyDef* pol, std::string name);
84 
86 {
87  std::string directory = "/nova/ana/users/ddoyle/NuebarEnergyEstimatorFiles/";
88  TFile* fFlat = TFile::Open((directory + "flat_energy_spectra.root").c_str());
89 
90  Spectrum* sEM_vs_HAD_nom = Spectrum::LoadFrom(fFlat, "em_vs_had_nominal")).release(;
91  Spectrum* sEM_vs_HAD_flat = Spectrum::LoadFrom(fFlat, "em_vs_had_flat")).release(;
92  Spectrum* sEM_vs_HAD_flat_energy_sum = Spectrum::LoadFrom(fFlat, "em_vs_had_flat_energy_sum")).release(;
93 
94  Spectrum* sTrueE_nom = Spectrum::LoadFrom(fFlat, "trueE_nominal")).release(;
95  Spectrum* sTrueE_flat = Spectrum::LoadFrom(fFlat, "trueE_flat")).release(;
96  fFlat->Close();
97 
98  double POT = sTrueE_flat->POT();
99 
100  TH2D* hEM_vs_HAD_nom = (TH2D*) sEM_vs_HAD_nom->ToTH2(POT);
101  TH2D* hEM_vs_HAD_flat = (TH2D*) sEM_vs_HAD_flat->ToTH2(POT);
102  TH2D* hEM_vs_HAD_flat_energy_sum = (TH2D*) sEM_vs_HAD_flat_energy_sum->ToTH2(POT);
103  TH1D* hTrueE_nom = (TH1D*)sTrueE_nom->ToTH1(POT);
104  TH1D* hTrueE_flat = (TH1D*) sTrueE_flat->ToTH1(POT);
105 
106  // Calculate average true nuetrino energy surface
107  // Set errors on the average energy as avg_energy / sqrt(N) (error on the mean)
108  TH2D* avg_energy_surface = MakeAverageTrueEnergySurface(hEM_vs_HAD_flat_energy_sum, hEM_vs_HAD_flat);
109 
110  // Make some plots
111  PrintPlot(hEM_vs_HAD_nom, directory + "Plots/em_vs_had_nom.pdf", "colz");
112  PrintPlot(hEM_vs_HAD_flat, directory + "Plots/em_vs_had_flat.pdf", "colz");
113  PrintPlot(hEM_vs_HAD_flat_energy_sum, directory + "Plots/em_vs_had_flat_energy_sum.pdf", "colz");
114  PrintPlot(avg_energy_surface, directory + "Plots/avg_energy_surface.pdf", "colz");
115  PrintPlot(hTrueE_nom, directory + "Plots/trueE_nom.pdf", "hist");
116  PrintPlot(hTrueE_flat, directory + "Plots/trueE_flat.pdf", "hist");
117 
118  // Do a fit to a polynomial
119  // Try a few different function
120  PolyDef* pol2 = new PolyDef("pol2", 2, false);
121  PolyDef* pol2X = new PolyDef("pol2x", 2, true);
122  PolyDef* pol3 = new PolyDef("pol3", 3, false);
123  PolyDef* pol3X = new PolyDef("pol3x", 3, true);
124 
125  // Do the fits
126  std::cout << "Starting fits...." << std::endl;
127  std::cout << "-----------------------------------------\n";
128  avg_energy_surface->Fit(pol2->name, "+");
129  avg_energy_surface->Fit(pol2X->name, "+");
130  avg_energy_surface->Fit(pol3->name, "+");
131  avg_energy_surface->Fit(pol3X->name, "+");
132  std::cout << "-----------------------------------------\n";
133 
134  // Print the fit parameters so they can be put into a header
135  std::ofstream fout;
136  fout.open("FitResults/all_fit_results.txt");
137  PrintFitResults(avg_energy_surface, fout, pol2 );
138  PrintFitResults(avg_energy_surface, fout, pol2X );
139  PrintFitResults(avg_energy_surface, fout, pol3 );
140  PrintFitResults(avg_energy_surface, fout, pol3X );
141  fout.close();
142 
143  PrintErrorSurface(avg_energy_surface, pol2 , directory + "Plots/error_surface_pol2.pdf" );
144  PrintErrorSurface(avg_energy_surface, pol2X, directory + "Plots/error_surface_pol2x.pdf" );
145  PrintErrorSurface(avg_energy_surface, pol3 , directory + "Plots/error_surface_pol3.pdf" );
146  PrintErrorSurface(avg_energy_surface, pol3X, directory + "Plots/error_surface_pol3x.pdf" );
147 }
148 
150 {
151  Double_t* pars = pol->GetFunction()->GetParameters();
152  int dof = pol->GetFunction()->GetNDF();
153 
154  int NX = hist->GetNbinsX();
155  int NY = hist->GetNbinsY();
156  double x1 = hist->GetXaxis()->GetBinLowEdge(1);
157  double x2 = hist->GetXaxis()->GetBinUpEdge(NX);
158  double y1 = hist->GetYaxis()->GetBinLowEdge(1);
159  double y2 = hist->GetYaxis()->GetBinUpEdge(NY);
160 
161  TH2D* ratio = new TH2D("ratio_" + pol->GetName(),
162  "(Reco - True) / True #bar{#nu} Energy " + pol->GetName(),
163  NX, x1, x2, NY, y1, y2);
164  int nonempty = 0;
165  for(int i = 1; i < NX + 1; i++) {
166  for(int j = 1; j < NY + 1; j++) {
167  if(hist->GetBinContent(i, j) > 0) {
168  nonempty++;
169  double x = hist->GetXaxis()->GetBinCenter(i);
170  double y = hist->GetYaxis()->GetBinCenter(j);
171  double pcoord[2] = {x, y};
172  Double_t* coord = pcoord;
173  double fval = PolN(coord, pars);
174  ratio->SetBinContent(i, j, (fval - hist->GetBinContent(i, j)) / hist->GetBinContent(i, j));
175  }
176  }
177  }
178  ratio->SetMaximum(0.5);
179  ratio->SetMinimum(-0.5);
180  ratio->GetXaxis()->SetTitle("EM Shower Energy [GeV]");
181  ratio->GetYaxis()->SetTitle("Hadronic Energy [GeV]");
182 
184 
185  TDirectory* temp = gDirectory;
186  TCanvas* c = new TCanvas();
187  c->cd();
188  ratio->Draw("colz");
189  c->SaveAs(name.c_str());
190  temp->cd();
191  delete c;
192 }
193 
194 void PrintFitResults(TH2D* hist, std::ofstream& fout, PolyDef* pol)
195 {
196  TF1* fun = hist->GetFunction(pol->name);
197  int dof = fun->GetNDF();
198  // First display on screen
199  std::cout << "-----------------------------------------\n";
200  std::cout << "Fit Results: " << pol->name << std::endl;
201  std::cout << "chi2 = " << fun->GetChisquare() << std::endl;
202  std::cout << "chi2/dof = " << fun->GetChisquare()/dof << std::endl;
203  for(int i = 0; i < pol->nparam; i++) {
204  std::cout << "double p" << i << " = " << fun->GetParameter(i) << ";" << std::endl;
205  }
206  std::cout << "-----------------------------------------\n";
207 
208  // Now write it to file
209  fout << "-----------------------------------------\n";
210  fout << "Fit Results: " << pol->name << std::endl;
211  fout << "chi2 = " << fun->GetChisquare() << std::endl;
212  fout << "chi2/dof = " << fun->GetChisquare()/dof << std::endl;
213  for(int i = 0; i < pol->nparam; i++) {
214  fout << "double p" << i << " = " << fun->GetParameter(i) << ";" << std::endl;
215  }
216  fout << "-----------------------------------------\n";
217 
218  // now output each individual poly for parsing later
219  std::ofstream polyout;
220  polyout.open("FitResults/" + pol->name + "_fit_results.txt");
221  for(int i = 0; i < pol->nparam; i++) {
222  polyout << "p" << i << "=" << fun->GetParameter(i) << std::endl;
223  }
224  polyout << "chi2=" << fun->GetChisquare() << std::endl;
225  polyout << "chi2/dof= " << fun->GetChisquare()/dof << std::endl;
226  polyout.close();
227 
228 }
229 
230 Double_t PolN(Double_t *x, Double_t *par)
231 {
232  Double_t result = par[0]*x[0] + par[1]*x[1]; // 1 order
233  result += par[2]*x[0]*x[0] + par[3]*x[1]*x[1]; // 2 order
234  result += par[4]*x[0]*x[0]*x[0] + par[5]*x[1]*x[1]*x[1]; // 3 order
235  result += par[6]*x[0]*x[1]; // 2 order cross terms
236  result += par[7]*x[0]*x[0]*x[1] + par[8]*x[0]*x[1]*x[1]; // 3 order cross terms
237  return result;
238 }
239 
240 TH2D* MakeAverageTrueEnergySurface(TH2D* hFlat_energy, TH2D* hFlat)
241 {
242  int NX = hFlat->GetNbinsX();
243  int NY = hFlat->GetNbinsY();
244  double x1 = hFlat->GetXaxis()->GetBinLowEdge(1);
245  double x2 = hFlat->GetXaxis()->GetBinUpEdge(NX);
246  double y1 = hFlat->GetYaxis()->GetBinLowEdge(1);
247  double y2 = hFlat->GetYaxis()->GetBinUpEdge(NY);
248 
249  TH2D* ret = new TH2D("avg_energy", "Average True #bar{#nu} Energy [GeV]", NX, x1, x2, NY, y1, y2);
250  ret->GetXaxis()->SetTitle("EM Shower Energy [GeV]");
251  ret->GetYaxis()->SetTitle("Hadronic Energy [GeV]");
252  for(int i = 1; i < NX + 1; i++) {
253  for(int j = 1; j < NY + 1; j++) {
254  double n_events = hFlat->GetBinContent(i, j);
255  double avg_energy = hFlat_energy->GetBinContent(i, j) / n_events;
256  double error = avg_energy / sqrt(n_events);
257  if(avg_energy > 0) {
258  ret->SetBinContent(i, j, avg_energy);
259  ret->SetBinError(i, j, error);
260  }
261  }
262  }
263  return ret;
264 }
265 
266 void PrintPlot(TH1D* h, std::string name, std::string option)
267 {
268  TDirectory* tmp = gDirectory;
269 
270  TCanvas* c = new TCanvas();
271  c->cd();
272  h->Draw(option.c_str());
273  c->SaveAs(name.c_str());
274  delete c;
275 
276  tmp->cd();
277 }
278 
279 void PrintPlot(TH2D* h, std::string name, std::string option)
280 {
281  TDirectory* tmp = gDirectory;
282 
283  TCanvas* c = new TCanvas();
284  c->cd();
285  h->Draw(option.c_str());
286  c->SaveAs(name.c_str());
287  delete c;
288 
289  tmp->cd();
290 }
291 
292 void SetErrors(TH2D* hFlat_avg_energy, TH2* hFlat)
293 {
294  int NX = hFlat->GetNbinsX();
295  int NY = hFlat->GetNbinsY();
296  for(int i = 1; i < NX + 1; i++) {
297  for(int j = 1; j < NY + 1; j++) {
298  double n_events = hFlat->GetBinContent(i, j);
299  double avg_energy = hFlat_avg_energy->GetBinContent(i, j);
300  double error = avg_energy / sqrt(n_events);
301  hFlat_avg_energy->SetBinError(i, j, error);
302  }
303  }
304 }
const XML_Char * name
Definition: expat.h:151
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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
Float_t x1[n_points_granero]
Definition: compare.C:5
PolyDef(Int_t _nparam, TString _name, TF2 *_func)
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
Int_t par
Definition: SimpleIterate.C:24
Double_t PolN(Double_t *x, Double_t *par)
string directory
projection from multiple chains
const Int_t nparam
Float_t tmp
Definition: plot.C:36
TF2 * GetFunction()
void PrintErrorSurface(TH2D *hist, PolyDef *pol, std::string name)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void PrintFitResults(TH2D *hist, std::ofstream &fout, PolyDef *pol)
PolyDef(TString _name, int _degree, bool _cross_terms)
void SetPaletteBlueRedWhite()
Definition: Style.cxx:7
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
void PrintPlot(TH1D *h, std::string name, std::string option)
void SetErrors(TH2D *hFlat_avg_energy, TH2 *hFlat)
double func(double x, double y)
const double j
Definition: BetheBloch.cxx:29
double POT() const
Definition: Spectrum.h:219
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void makeEnergyEstimator()
std::string pars("Th23Dmsq32")
TString GetName()
TH2D * MakeAverageTrueEnergySurface(TH2D *hFlat_energy, TH2D *hFlat)
static constexpr Double_t degree
Definition: Munits.h:157
int n_events
Definition: test_Tier.py:22