compareChiSqrCalcs.C
Go to the documentation of this file.
1 #include <vector>
2 #include <string>
3 #include <iostream>
4 #include "TFile.h"
5 #include "TGraph2D.h"
6 #include "TGraph.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TCanvas.h"
10 #include "TLegend.h"
11 
12 struct BestFitInfo{
13 
14  BestFitInfo(double x,
15  double y,
16  double z)
17  : xVal(x)
18  , yVal(y)
19  , minChiSqr(z)
20  {}
21 
22  double xVal;
23  double yVal;
24  double minChiSqr;
25 };
26 
27 //------------------------------------------------------------------------------
29  size_t numUniverses,
30  std::vector<BestFitInfo> & bestFitInfoVec)
31 {
32  // loop over the directories for the different universes and get the best
33  // fit points and chi^2 values
35  for(size_t u = 0; u < numUniverses; ++u){
36  name = "contours/UniverseContours_" + std::to_string(u) + "_/UniverseContours_" + std::to_string(u) + "_FitResult";
37  auto *fitGr = dynamic_cast<TGraph2D*>(file->Get(name.c_str()));
38 
39  // get the point location and chi^2
40  if(fitGr != nullptr){
41  double x, y, z;
42 
43  fitGr->GetPoint(0, x, y, z);
44 
45  bestFitInfoVec.emplace_back(x, y, z);
46 
47  std::cout << "adding point ("
48  << x
49  << ", "
50  << y
51  << ", "
52  << z
53  << ")"
54  << std::endl;
55  }
56  else
57  std::cout << "could not find graph "
58  << name
59  << std::endl;
60  }
61 
62 }
63 
64 //------------------------------------------------------------------------------
65 void DrawBestFitValues(std::string const& pdfName,
66  std::string const& calc1Label,
67  std::string const& calc2Label,
68  std::vector<BestFitInfo> const& calc1Vec,
69  std::vector<BestFitInfo> const& calc2Vec)
70 {
71  std::cout << "Drawing best fit plots" << std::endl;
72 
73  // now make some histograms showing the differences in the best fit outcomes
74  // for the two methods
75  TGraph *calc1XVals = new TGraph(calc1Vec.size());
76  TGraph *calc2XVals = new TGraph(calc1Vec.size());
77  TGraph *calc1YVals = new TGraph(calc1Vec.size());
78  TGraph *calc2YVals = new TGraph(calc1Vec.size());
79  TGraph *calc1ChiSqrVals = new TGraph(calc1Vec.size());
80  TGraph *calc2ChiSqrVals = new TGraph(calc1Vec.size());
81 
82  /* calc1XVals ->SetMarkerSize(4); */
83  /* calc1YVals ->SetMarkerSize(4); */
84  /* calc1ChiSqrVals->SetMarkerSize(4); */
85  calc1XVals ->SetMarkerStyle(kFullCircle);
86  calc1YVals ->SetMarkerStyle(kFullCircle);
87  calc1ChiSqrVals->SetMarkerStyle(kFullCircle);
88  calc1XVals ->SetMarkerColor(kRed);
89  calc1YVals ->SetMarkerColor(kRed);
90  calc1ChiSqrVals->SetMarkerColor(kRed);
91  calc1XVals ->SetLineColor(kRed);
92  calc1YVals ->SetLineColor(kRed);
93  calc1ChiSqrVals->SetLineColor(kRed);
94  calc2XVals ->SetMarkerStyle(kFullCircle);
95  calc2YVals ->SetMarkerStyle(kFullCircle);
96  calc2ChiSqrVals->SetMarkerStyle(kFullCircle);
97 
98  for(size_t u = 0; u < calc1Vec.size(); ++u){
99  calc1XVals ->SetPoint(u, u, calc1Vec[u].xVal );
100  calc2XVals ->SetPoint(u, u, calc2Vec[u].xVal );
101  calc1YVals ->SetPoint(u, u, calc1Vec[u].yVal );
102  calc2YVals ->SetPoint(u, u, calc2Vec[u].yVal );
103  calc1ChiSqrVals->SetPoint(u, u, calc1Vec[u].minChiSqr);
104  calc2ChiSqrVals->SetPoint(u, u, calc2Vec[u].minChiSqr);
105  }
106 
107  TH2F* xVals = new TH2F("xVals",
108  ";Universe;Best Fit - sin^{2}#theta_{23}",
109  100,
110  0,
111  100,
112  100,
113  0.35,
114  0.75);
115 
116  TH2F* yVals = new TH2F("yVals",
117  ";Universe;Best Fit - #Delta m^{2}_{32}",
118  100,
119  0,
120  100,
121  100,
122  2.35,
123  2.75);
124 
125  TH2F* chiSqrVals = new TH2F("chiSqrVals",
126  ";Universe;Best Fit - #chi^{2}",
127  100,
128  0,
129  100,
130  100,
131  0,
132  300);
133 
134 
135  TCanvas *valsCanv = new TCanvas("valsCanv");
136  valsCanv->Divide(2, 2);
137  valsCanv->cd(1);
138 
139  xVals->Draw();
140  calc1XVals->Draw("psame");
141  calc2XVals->Draw("psame");
142 
143  valsCanv->cd(2);
144  yVals->Draw();
145  calc1YVals->Draw("psame");
146  calc2YVals->Draw("psame");
147 
148  valsCanv->cd(3);
149  chiSqrVals->Draw();
150  calc1ChiSqrVals->Draw("psame");
151  calc2ChiSqrVals->Draw("psame");
152 
153  valsCanv->cd(4);
154  TLegend *leg = new TLegend(0.25, 0.25, 0.75, 0.75);
155  leg->AddEntry(calc1XVals, calc1Label.c_str(), "p");
156  leg->AddEntry(calc2XVals, calc2Label.c_str(), "p");
157  leg->Draw();
158 
159  valsCanv->Print((pdfName + "(").c_str(), "pdf");
160 }
161 
162 //------------------------------------------------------------------------------
164  std::vector<BestFitInfo> const& calc1Vec,
165  std::vector<BestFitInfo> const& calc2Vec)
166 {
167  std::cout << "Drawing best fit difference plots" << std::endl;
168 
169  // create some TGraphs showing the differences between the calculations
170 
171  TGraph* xDiffs = new TGraph(calc1Vec.size());
172  TGraph* yDiffs = new TGraph(calc1Vec.size());
173  TGraph* chiSqrDiffs = new TGraph(calc1Vec.size());
174 
175  xDiffs ->SetMarkerStyle(kFullCircle);
176  yDiffs ->SetMarkerStyle(kFullCircle);
177  chiSqrDiffs->SetMarkerStyle(kFullCircle);
178  /* xDiffs ->SetMarkerSize(4); */
179  /* yDiffs ->SetMarkerSize(4); */
180  /* chiSqrDiffs->SetMarkerSize(4); */
181 
182  for(size_t i = 0; i < calc1Vec.size(); ++i){
183  xDiffs ->SetPoint(i, i, calc1Vec[i].xVal - calc2Vec[i].xVal );
184  yDiffs ->SetPoint(i, i, calc1Vec[i].yVal - calc2Vec[i].yVal );
185  chiSqrDiffs->SetPoint(i, i, calc1Vec[i].minChiSqr - calc2Vec[i].minChiSqr);
186  }
187 
188  TH2F* xDiffHist = new TH2F("xDiffs",
189  ";Universe;#Delta Best Fit - sin^{2}#theta_{23}",
190  100,
191  0,
192  100,
193  100,
194  -1.,
195  1.);
196 
197  TH2F* yDiffHist = new TH2F("yDiffs",
198  ";Universe;#Delta Best Fit - #Delta m^{2}_{32}",
199  100,
200  0,
201  100,
202  100,
203  -1,
204  1);
205 
206  TH2F* chiSqrDiffHist = new TH2F("chiSqrDiffs",
207  ";Universe;#Delta Best Fit - #chi^{2}",
208  100,
209  0,
210  100,
211  100,
212  -150,
213  150);
214 
215 
216  TCanvas *diffsCanv = new TCanvas("diffsCanv");
217  diffsCanv->Divide(2, 2);
218  diffsCanv->cd(1);
219 
220  xDiffHist->Draw();
221  xDiffs->Draw("psame");
222 
223  diffsCanv->cd(2);
224  yDiffHist->Draw();
225  yDiffs->Draw("psame");
226 
227  diffsCanv->cd(3);
228  chiSqrDiffHist->Draw();
229  chiSqrDiffs->Draw("psame");
230 
231  diffsCanv->Print((pdfName + ")").c_str(), "pdf");
232 }
233 
234 //------------------------------------------------------------------------------
236 {
237  std::string fileName1("cmf_3flavor_test_cont_no_systs_FD_Poission_by_bin_statonly_covmat_hist.root");
238  std::string calc1Label("Stats Only Covariance Matrix");
239 
240  std::string fileName2("cmf_3flavor_test_cont_no_systs_FD_Poission_by_bin_poisson_calc_hist.root");
241  std::string calc2Label("Poisson LL");
242 
243  TFile *file1 = new TFile(fileName1.c_str(), "READ");
244  TFile *file2 = new TFile(fileName2.c_str(), "READ");
245 
246  std::cout << "files open" << std::endl;
247 
248  size_t numUniverses = 100;
249 
250  std::vector<BestFitInfo> calc1BestFitInfo;
251  std::vector<BestFitInfo> calc2BestFitInfo;
252 
253  extractBestFitInfo(file1, numUniverses, calc1BestFitInfo);
254  extractBestFitInfo(file2, numUniverses, calc2BestFitInfo);
255 
256  file1->Close();
257  file2->Close();
258 
259  std::cout << "files closed" << endl;
260 
261  std::string pdfFileName("compareChiSqrCalcs_statsonly.pdf");
262 
263  DrawBestFitValues(pdfFileName,
264  calc1Label,
265  calc2Label,
266  calc1BestFitInfo,
267  calc2BestFitInfo);
268 
269  DrawBestFitDifferences(pdfFileName,
270  calc1BestFitInfo,
271  calc2BestFitInfo);
272 }
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
void DrawBestFitValues(std::string const &pdfName, std::string const &calc1Label, std::string const &calc2Label, std::vector< BestFitInfo > const &calc1Vec, std::vector< BestFitInfo > const &calc2Vec)
void DrawBestFitDifferences(std::string const &pdfName, std::vector< BestFitInfo > const &calc1Vec, std::vector< BestFitInfo > const &calc2Vec)
BestFitInfo(double x, double y, double z)
void compareChiSqrCalcs()
z
Definition: test.py:28
void extractBestFitInfo(TFile *file, size_t numUniverses, std::vector< BestFitInfo > &bestFitInfoVec)
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
size_t numUniverses
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
enum BeamMode string