34 void SetErrors(TH2D* hFlat_avg_energy, TH2* hFlat);
37 Double_t
PolN(Double_t *
x, Double_t *
par);
41 PolyDef(Int_t _nparam, TString _name, TF2* _func)
42 : nparam(_nparam),
name(_name),
func(_func)
44 PolyDef(TString _name,
int _degree,
bool _cross_terms)
48 cross_terms(_cross_terms)
56 func->FixParameter(4, 0);
57 func->FixParameter(5, 0);
58 func->FixParameter(7, 0);
59 func->FixParameter(8, 0);
61 func->FixParameter(6, 0);
66 func->FixParameter(6, 0);
67 func->FixParameter(7, 0);
68 func->FixParameter(8, 0);
88 TFile* fFlat = TFile::Open((directory +
"flat_energy_spectra.root").c_str());
98 double POT = sTrueE_flat->
POT();
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);
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");
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";
136 fout.open(
"FitResults/all_fit_results.txt");
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" );
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);
162 "(Reco - True) / True #bar{#nu} Energy " + pol->
GetName(),
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) {
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));
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]");
185 TDirectory*
temp = gDirectory;
186 TCanvas*
c =
new TCanvas();
189 c->SaveAs(name.c_str());
196 TF1* fun = hist->GetFunction(pol->
name);
197 int dof = fun->GetNDF();
199 std::cout <<
"-----------------------------------------\n";
206 std::cout <<
"-----------------------------------------\n";
209 fout <<
"-----------------------------------------\n";
211 fout <<
"chi2 = " << fun->GetChisquare() <<
std::endl;
212 fout <<
"chi2/dof = " << fun->GetChisquare()/dof <<
std::endl;
214 fout <<
"double p" <<
i <<
" = " << fun->GetParameter(
i) <<
";" <<
std::endl;
216 fout <<
"-----------------------------------------\n";
219 std::ofstream polyout;
220 polyout.open(
"FitResults/" + pol->
name +
"_fit_results.txt");
222 polyout <<
"p" <<
i <<
"=" << fun->GetParameter(
i) <<
std::endl;
224 polyout <<
"chi2=" << fun->GetChisquare() <<
std::endl;
225 polyout <<
"chi2/dof= " << fun->GetChisquare()/dof <<
std::endl;
232 Double_t
result = par[0]*x[0] + par[1]*x[1];
233 result += par[2]*x[0]*x[0] + par[3]*x[1]*x[1];
234 result += par[4]*x[0]*x[0]*x[0] + par[5]*x[1]*x[1]*x[1];
235 result += par[6]*x[0]*x[1];
236 result += par[7]*x[0]*x[0]*x[1] + par[8]*x[0]*x[1]*x[1];
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);
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++) {
255 double avg_energy = hFlat_energy->GetBinContent(
i,
j) /
n_events;
256 double error = avg_energy /
sqrt(n_events);
258 ret->SetBinContent(
i,
j, avg_energy);
259 ret->SetBinError(
i,
j, error);
268 TDirectory*
tmp = gDirectory;
270 TCanvas*
c =
new TCanvas();
272 h->Draw(option.c_str());
273 c->SaveAs(name.c_str());
281 TDirectory*
tmp = gDirectory;
283 TCanvas*
c =
new TCanvas();
285 h->Draw(option.c_str());
286 c->SaveAs(name.c_str());
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++) {
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);
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Cuts and Vars for the 2020 FD DiF Study.
Float_t y1[n_points_granero]
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.
Float_t x1[n_points_granero]
PolyDef(Int_t _nparam, TString _name, TF2 *_func)
Double_t PolN(Double_t *x, Double_t *par)
string directory
projection from multiple chains
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
void PrintErrorSurface(TH2D *hist, PolyDef *pol, std::string name)
Representation of a spectrum in any variable, with associated POT.
void PrintFitResults(TH2D *hist, std::ofstream &fout, PolyDef *pol)
PolyDef(TString _name, int _degree, bool _cross_terms)
void SetPaletteBlueRedWhite()
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
void PrintPlot(TH1D *h, std::string name, std::string option)
void SetErrors(TH2D *hFlat_avg_energy, TH2 *hFlat)
double func(double x, double y)
std::vector< double > POT
void makeEnergyEstimator()
TH2D * MakeAverageTrueEnergySurface(TH2D *hFlat_energy, TH2D *hFlat)
static constexpr Double_t degree