compareCMFCAFResults.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TH1.h"
3 #include "TCanvas.h"
4 #include "TLegend.h"
5 #include "TGraph.h"
6 #include <string>
7 #include <iostream>
8 #include <vector>
9 
10 //------------------------------------------------------------------------------
11 void grabNuMuHistograms(std::vector<TH1D*> & cmfFHCFDNuMuQuants,
12  std::vector<TH1D*> & cafFHCFDNuMuQuants,
13  std::vector<TH1D*> & ratFHCFDNuMuQuants,
14  std::vector<TH1D*> & cmfRHCFDNuMuQuants,
15  std::vector<TH1D*> & cafRHCFDNuMuQuants,
16  std::vector<TH1D*> & ratRHCFDNuMuQuants,
17  TFile * cmfFile,
18  TFile * cafFile,
19  bool useOscillated)
20 {
21  std::string quantStr;
23  for(int i = 1; i < 5; ++i){
24  quantStr = std::to_string(i);
25  name = "fit/FitResults/DataEnergySpectrumFHCFarDetNuMuSelQ" + quantStr;
26  cmfFHCFDNuMuQuants.push_back(dynamic_cast<TH1D*>(cmfFile->Get(name.c_str())));
27  cmfFHCFDNuMuQuants.back()->SetMarkerStyle(20);
28 
29  name = "fit/FitResults/DataEnergySpectrumRHCFarDetNuMuSelQ" + quantStr;
30  cmfRHCFDNuMuQuants.push_back(dynamic_cast<TH1D*>(cmfFile->Get(name.c_str())));
31  cmfRHCFDNuMuQuants.back()->SetMarkerStyle(20);
32 
33  ratFHCFDNuMuQuants.push_back(new TH1D(("ratioFHCNuMuQ" + quantStr).c_str(),
34  ";Energy (GeV);CMF/CAF Ratio",
35  cmfFHCFDNuMuQuants.back()->GetXaxis()->GetNbins(),
36  cmfFHCFDNuMuQuants.back()->GetXaxis()->GetXbins()->GetArray()));
37  ratFHCFDNuMuQuants.back()->GetXaxis()->CenterTitle();
38  ratFHCFDNuMuQuants.back()->GetYaxis()->CenterTitle();
39  ratFHCFDNuMuQuants.back()->GetXaxis()->SetDecimals();
40  ratFHCFDNuMuQuants.back()->GetYaxis()->SetDecimals();
41  ratFHCFDNuMuQuants.back()->SetMarkerStyle(22);
42 
43  ratRHCFDNuMuQuants.push_back(new TH1D(("ratioRHCNuMuQ" + quantStr).c_str(),
44  ";Energy (GeV);CMF/CAF Ratio",
45  cmfRHCFDNuMuQuants.back()->GetXaxis()->GetNbins(),
46  cmfRHCFDNuMuQuants.back()->GetXaxis()->GetXbins()->GetArray()));
47  ratRHCFDNuMuQuants.back()->GetXaxis()->CenterTitle();
48  ratRHCFDNuMuQuants.back()->GetYaxis()->CenterTitle();
49  ratRHCFDNuMuQuants.back()->GetXaxis()->SetDecimals();
50  ratRHCFDNuMuQuants.back()->GetYaxis()->SetDecimals();
51  ratRHCFDNuMuQuants.back()->SetMarkerStyle(22);
52 
53  name = (useOscillated) ? "NuMuQ" + quantStr + "FHCOscillated" : "NuMuQ" + quantStr + "FHCUnoscillated";
54  cafFHCFDNuMuQuants.push_back(dynamic_cast<TH1D*>(cafFile->Get(name.c_str())));
55  cafFHCFDNuMuQuants.back()->SetLineColor(2);
56 
57  name = (useOscillated) ? "NuMuQ" + quantStr + "RHCOscillated" : "NuMuQ" + quantStr + "RHCUnoscillated";
58  cafRHCFDNuMuQuants.push_back(dynamic_cast<TH1D*>(cafFile->Get(name.c_str())));
59  cafRHCFDNuMuQuants.back()->SetLineColor(2);
60 
61  // make the ratios of CMF to CAF counts
62  ratFHCFDNuMuQuants.back()->Divide(cmfFHCFDNuMuQuants.back(), cafFHCFDNuMuQuants.back());
63  ratRHCFDNuMuQuants.back()->Divide(cmfRHCFDNuMuQuants.back(), cafRHCFDNuMuQuants.back());
64 
65  } // end loop to fill NuMu selections
66 }
67 
68 //------------------------------------------------------------------------------
69 void makeNuMuCanvas(std::vector<TH1D*> & cmfFHCFDNuMuQuants,
70  std::vector<TH1D*> & cafFHCFDNuMuQuants,
71  std::vector<TH1D*> & ratFHCFDNuMuQuants,
72  std::vector<TH1D*> & cmfRHCFDNuMuQuants,
73  std::vector<TH1D*> & cafRHCFDNuMuQuants,
74  std::vector<TH1D*> & ratRHCFDNuMuQuants,
75  TCanvas * canv)
76 {
77  // divide the canvas into 16 parts, 1 each for the FHC/RHC spectra in quantiles and
78  // some for the ratios too
79  canv->Divide(4,4);
80 
81  TVirtualPad* pad;
82  for(size_t i = 0; i < cmfFHCFDNuMuQuants.size(); ++i){
83 
84  pad = canv->cd(i + 1);
85  cafFHCFDNuMuQuants[i]->SetMaximum(cafFHCFDNuMuQuants[i]->GetMaximum() * 1.8);
86  cafFHCFDNuMuQuants[i]->Draw();
87  cmfFHCFDNuMuQuants[i]->Draw("psame");
88  if(i == 3){
89  TLegend *leg = new TLegend(0.6, 0.6, 0.9, 0.9);
90  leg->AddEntry(cmfFHCFDNuMuQuants[i], "CMF", "p");
91  leg->AddEntry(cafFHCFDNuMuQuants[i], "CAF", "l");
92  leg->SetFillColor(0);
93  leg->SetFillStyle(0);
94  leg->Draw();
95  pad->Update();
96  }
97 
98  pad = canv->cd(i + 5);
99  pad->SetGridy();
100  ratFHCFDNuMuQuants[i]->GetYaxis()->SetRangeUser(0.98, 1.02);
101  ratFHCFDNuMuQuants[i]->Draw("phist");
102 
103  canv->cd(i + 9);
104  cafRHCFDNuMuQuants[i]->SetMaximum(cafRHCFDNuMuQuants[i]->GetMaximum() * 1.8);
105  cafRHCFDNuMuQuants[i]->Draw();
106  cmfRHCFDNuMuQuants[i]->Draw("psame");
107 
108  pad = canv->cd(i + 13);
109  pad->SetGridy();
110  ratRHCFDNuMuQuants[i]->GetYaxis()->SetRangeUser(0.98, 1.02);
111  ratRHCFDNuMuQuants[i]->Draw("phist");
112  }
113 
114  canv->Update();
115 }
116 
117 //------------------------------------------------------------------------------
118 void grabNuEHistograms(std::vector<TH1D*> & cmfFHCFDNuESels,
119  std::vector<TH1D*> & cafFHCFDNuESels,
120  std::vector<TH1D*> & ratFHCFDNuESels,
121  std::vector<TH1D*> & cmfRHCFDNuESels,
122  std::vector<TH1D*> & cafRHCFDNuESels,
123  std::vector<TH1D*> & ratRHCFDNuESels,
124  TFile * cmfFile,
125  TFile * cafFile,
126  bool useOscillated)
127 {
128  // the caf histogram has the NuE selections all in a single histogram,
129  // cmf has 3 different histograms, so we have to fix that.
130  std::vector<std::string> selections({"NuESel_LowPID", "NuESel_HighPID", "NuESel_Peripheral"});
131 
133  for(auto const& itr : selections){
134  name = "fit/FitResults/DataEnergySpectrumFHCFarDet" + itr;
135  cmfFHCFDNuESels.push_back(dynamic_cast<TH1D*>(cmfFile->Get(name.c_str())));
136  cmfFHCFDNuESels.back()->SetMarkerStyle(20);
137 
138  name = "fit/FitResults/DataEnergySpectrumRHCFarDet" + itr;
139  cmfRHCFDNuESels.push_back(dynamic_cast<TH1D*>(cmfFile->Get(name.c_str())));
140  cmfRHCFDNuESels.back()->SetMarkerStyle(20);
141 
142  std::cout << cmfFHCFDNuESels.back()->GetName() << " "
143  << cmfFHCFDNuESels.back()->GetXaxis()->GetNbins() << " "
144  << cmfFHCFDNuESels.back()->GetXaxis()->GetXmin() << " "
145  << cmfFHCFDNuESels.back()->GetXaxis()->GetXmax() << std::endl;
146 
147  name = "cafFHCFD" + itr;
148  cafFHCFDNuESels.push_back(new TH1D(name.c_str(),
149  ";Energy (GeV);Events",
150  cmfFHCFDNuESels.back()->GetXaxis()->GetNbins(),
151  cmfFHCFDNuESels.back()->GetXaxis()->GetXmin(),
152  cmfFHCFDNuESels.back()->GetXaxis()->GetXmax()));
153  cafFHCFDNuESels.back()->SetLineColor(2);
154 
155  std::cout << cafFHCFDNuESels.back()->GetName() << " " << cafFHCFDNuESels.size() << std::endl;
156 
157  name = "cafRHCFD" + itr;
158  cafRHCFDNuESels.push_back(new TH1D(name.c_str(),
159  ";Energy (GeV);Events",
160  cmfRHCFDNuESels.back()->GetXaxis()->GetNbins(),
161  cmfRHCFDNuESels.back()->GetXaxis()->GetXmin(),
162  cmfRHCFDNuESels.back()->GetXaxis()->GetXmax()));
163  cafRHCFDNuESels.back()->SetLineColor(2);
164 
165  name = "ratioFDFHC" + itr;
166  ratFHCFDNuESels.push_back(new TH1D(name.c_str(),
167  ";Energy (GeV);CMF/CAF Ratio",
168  cmfFHCFDNuESels.back()->GetXaxis()->GetNbins(),
169  cmfFHCFDNuESels.back()->GetXaxis()->GetXmin(),
170  cmfFHCFDNuESels.back()->GetXaxis()->GetXmax()));
171  ratFHCFDNuESels.back()->GetXaxis()->CenterTitle();
172  ratFHCFDNuESels.back()->GetYaxis()->CenterTitle();
173  ratFHCFDNuESels.back()->GetXaxis()->SetDecimals();
174  ratFHCFDNuESels.back()->GetYaxis()->SetDecimals();
175  ratFHCFDNuESels.back()->SetMarkerStyle(22);
176 
177  name = "ratioFDRHC" + itr;
178  ratRHCFDNuESels.push_back(new TH1D(name.c_str(),
179  ";Energy (GeV);CMF/CAF Ratio",
180  cmfFHCFDNuESels.back()->GetXaxis()->GetNbins(),
181  cmfFHCFDNuESels.back()->GetXaxis()->GetXmin(),
182  cmfFHCFDNuESels.back()->GetXaxis()->GetXmax()));
183 
184  ratRHCFDNuESels.back()->GetXaxis()->CenterTitle();
185  ratRHCFDNuESels.back()->GetYaxis()->CenterTitle();
186  ratRHCFDNuESels.back()->GetXaxis()->SetDecimals();
187  ratRHCFDNuESels.back()->GetYaxis()->SetDecimals();
188  ratRHCFDNuESels.back()->SetMarkerStyle(22);
189 
190  }
191 
192  // now fill the caf histograms
193  name = (useOscillated) ? "Oscillated" : "Unoscillated";
194  TH1D* cafFHC = dynamic_cast<TH1D*>(cafFile->Get(("NuEFHC" + name).c_str()));
195  TH1D* cafRHC = dynamic_cast<TH1D*>(cafFile->Get(("NuERHC" + name).c_str()));
196 
197  // in the CAF histograms the bins for the Low PID sample are
198  // 3 - 8, in 0.5 GeV steps starting at 1 GeV
199  for(int i = 3; i < 9; ++i){
200  std::cout << i << " " << (i - 3) * 0.5 + 1.25 << " " << cafFHC->GetBinContent(i) << " " << cafRHC->GetBinContent(i) << std::endl;
201  cafFHCFDNuESels[0]->Fill((i - 3) * 0.5 + 1.25, cafFHC->GetBinContent(i));
202  cafRHCFDNuESels[0]->Fill((i - 3) * 0.5 + 1.25, cafRHC->GetBinContent(i));
203  }
204 
205  // in the CAF histograms the bins for the High PID sample are
206  // 12 - 17, in 0.5 GeV steps starting at 1 GeV
207  for(int i = 12; i < 18; ++i){
208  std::cout << i << " " << (i - 12) * 0.5 + 1.25 << " " << cafFHC->GetBinContent(i) << " " << cafRHC->GetBinContent(i) << std::endl;
209  cafFHCFDNuESels[1]->Fill((i - 12) * 0.5 + 1.25, cafFHC->GetBinContent(i));
210  cafRHCFDNuESels[1]->Fill((i - 12) * 0.5 + 1.25, cafRHC->GetBinContent(i));
211  }
212 
213  // the peripheral samples are easy - last ones in the FD
214  // vector and we want just 1 number, but have to take the
215  // integral of several bins in the CAF prediction histograms
216  cafFHCFDNuESels.back()->Fill(5., cafFHC->Integral(20, 27));
217  cafRHCFDNuESels.back()->Fill(5., cafRHC->Integral(20, 27));
218 
219  // get the ratios
220  for(size_t h = 0; h < cmfFHCFDNuESels.size(); ++h){
221  cout << cmfFHCFDNuESels[h]->GetName() << " " << cafFHCFDNuESels[h]->GetName() << std::endl;
222  ratFHCFDNuESels[h]->Divide(cmfFHCFDNuESels[h], cafFHCFDNuESels[h]);
223  cout << cmfRHCFDNuESels[h]->GetName() << " " << cafRHCFDNuESels[h]->GetName() << std::endl;
224  ratRHCFDNuESels[h]->Divide(cmfRHCFDNuESels[h], cafRHCFDNuESels[h]);
225  }
226 }
227 
228 //------------------------------------------------------------------------------
229 void makeNuECanvas(std::vector<TH1D*> & cmfFHCFDNuESels,
230  std::vector<TH1D*> & cafFHCFDNuESels,
231  std::vector<TH1D*> & ratFHCFDNuESels,
232  std::vector<TH1D*> & cmfRHCFDNuESels,
233  std::vector<TH1D*> & cafRHCFDNuESels,
234  std::vector<TH1D*> & ratRHCFDNuESels,
235  TCanvas * canv)
236 {
237 
238  canv->Divide(3,4);
239 
240  TVirtualPad* pad;
241  for(size_t i = 0; i < cmfFHCFDNuESels.size(); ++i){
242 
243  pad = canv->cd(i + 1);
244  cafFHCFDNuESels[i]->SetMaximum(cafFHCFDNuESels[i]->GetMaximum() * 1.8);
245  cafFHCFDNuESels[i]->Draw();
246  cmfFHCFDNuESels[i]->Draw("psame");
247  if(i == 2){
248  TLegend *leg = new TLegend(0.15, 0.15, 0.45, 0.45);
249  leg->AddEntry(cmfFHCFDNuESels[i], "CMF", "p");
250  leg->AddEntry(cafFHCFDNuESels[i], "CAF", "l");
251  leg->SetFillColor(0);
252  leg->SetFillStyle(0);
253  leg->Draw();
254  pad->Update();
255  }
256 
257  pad = canv->cd(i + 4);
258  pad->SetGridy();
259  ratFHCFDNuESels[i]->GetYaxis()->SetRangeUser(0.97, 1.03);
260  ratFHCFDNuESels[i]->Draw("phist");
261 
262  pad = canv->cd(i + 7);
263  cafRHCFDNuESels[i]->SetMaximum(cafRHCFDNuESels[i]->GetMaximum() * 1.8);
264  cafRHCFDNuESels[i]->Draw();
265  cmfRHCFDNuESels[i]->Draw("psame");
266 
267  pad = canv->cd(i + 10);
268  pad->SetGridy();
269  ratRHCFDNuESels[i]->GetYaxis()->SetRangeUser(0.97, 1.03);
270  ratRHCFDNuESels[i]->Draw("phist");
271  }
272 
273  canv->Update();
274 }
275 
276 //------------------------------------------------------------------------------
277 void compareCMFCAFResults(bool useOscillated=false)
278 {
279  // open the files and pull the NuMu histograms for the FD
280  std::string cmfFileName = (useOscillated) ? "cmf_covariance_matrix_bestfit_hist.root" : "cmf_covariance_matrix_noosc_hist.root";
281 
282  TFile* cmfFile = new TFile(cmfFileName.c_str(), "READ");
283  TFile* cafFile = new TFile("caf_fd_predictions.root", "READ");
284 
285  std::vector<TH1D*> cmfFHCFDNuMuQuants;
286  std::vector<TH1D*> cafFHCFDNuMuQuants;
287  std::vector<TH1D*> ratFHCFDNuMuQuants;
288  std::vector<TH1D*> cmfFHCFDNuESels;
289  std::vector<TH1D*> cafFHCFDNuESels;
290  std::vector<TH1D*> ratFHCFDNuESels;
291  std::vector<TH1D*> cmfRHCFDNuMuQuants;
292  std::vector<TH1D*> cafRHCFDNuMuQuants;
293  std::vector<TH1D*> ratRHCFDNuMuQuants;
294  std::vector<TH1D*> cmfRHCFDNuESels;
295  std::vector<TH1D*> cafRHCFDNuESels;
296  std::vector<TH1D*> ratRHCFDNuESels;
297 
298  // get the NuMu histograms from the files and make ratios too
299  grabNuMuHistograms(cmfFHCFDNuMuQuants,
300  cafFHCFDNuMuQuants,
301  ratFHCFDNuMuQuants,
302  cmfRHCFDNuMuQuants,
303  cafRHCFDNuMuQuants,
304  ratRHCFDNuMuQuants,
305  cmfFile,
306  cafFile,
307  useOscillated);
308 
309  // Grab the NuE histograms from the files
310  grabNuEHistograms(cmfFHCFDNuESels,
311  cafFHCFDNuESels,
312  ratFHCFDNuESels,
313  cmfRHCFDNuESels,
314  cafRHCFDNuESels,
315  ratRHCFDNuESels,
316  cmfFile,
317  cafFile,
318  useOscillated);
319 
320  // Make the canvases
321 
322  TCanvas *numuCanv = new TCanvas("numuCanv", "#nu_{#mu} Comparison", 800, 800);
323  TCanvas *nueCanv = new TCanvas("nueCanv", "#nu_{e} Comparison", 800, 800);
324 
325  makeNuMuCanvas(cmfFHCFDNuMuQuants,
326  cafFHCFDNuMuQuants,
327  ratFHCFDNuMuQuants,
328  cmfRHCFDNuMuQuants,
329  cafRHCFDNuMuQuants,
330  ratRHCFDNuMuQuants,
331  numuCanv);
332 
333  makeNuECanvas(cmfFHCFDNuESels,
334  cafFHCFDNuESels,
335  ratFHCFDNuESels,
336  cmfRHCFDNuESels,
337  cafRHCFDNuESels,
338  ratRHCFDNuESels,
339  nueCanv);
340 
341 }
const XML_Char * name
Definition: expat.h:151
TCanvas * canv
void grabNuEHistograms(std::vector< TH1D * > &cmfFHCFDNuESels, std::vector< TH1D * > &cafFHCFDNuESels, std::vector< TH1D * > &ratFHCFDNuESels, std::vector< TH1D * > &cmfRHCFDNuESels, std::vector< TH1D * > &cafRHCFDNuESels, std::vector< TH1D * > &ratRHCFDNuESels, TFile *cmfFile, TFile *cafFile, bool useOscillated)
void compareCMFCAFResults(bool useOscillated=false)
void makeNuMuCanvas(std::vector< TH1D * > &cmfFHCFDNuMuQuants, std::vector< TH1D * > &cafFHCFDNuMuQuants, std::vector< TH1D * > &ratFHCFDNuMuQuants, std::vector< TH1D * > &cmfRHCFDNuMuQuants, std::vector< TH1D * > &cafRHCFDNuMuQuants, std::vector< TH1D * > &ratRHCFDNuMuQuants, TCanvas *canv)
void makeNuECanvas(std::vector< TH1D * > &cmfFHCFDNuESels, std::vector< TH1D * > &cafFHCFDNuESels, std::vector< TH1D * > &ratFHCFDNuESels, std::vector< TH1D * > &cmfRHCFDNuESels, std::vector< TH1D * > &cafRHCFDNuESels, std::vector< TH1D * > &ratRHCFDNuESels, TCanvas *canv)
void grabNuMuHistograms(std::vector< TH1D * > &cmfFHCFDNuMuQuants, std::vector< TH1D * > &cafFHCFDNuMuQuants, std::vector< TH1D * > &ratFHCFDNuMuQuants, std::vector< TH1D * > &cmfRHCFDNuMuQuants, std::vector< TH1D * > &cafRHCFDNuMuQuants, std::vector< TH1D * > &ratRHCFDNuMuQuants, TFile *cmfFile, TFile *cafFile, bool useOscillated)
OStream cout
Definition: OStream.cxx:6
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
enum BeamMode string