makeMedianContours.C
Go to the documentation of this file.
1 #include "TFile.h"
2 #include "TGraph.h"
3 #include "TCanvas.h"
4 #include "TDirectory.h"
5 #include "TKey.h"
6 #include "TMarker.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TLegend.h"
10 #include "TMath.h"
11 #include "TGraph.h"
12 #include "TPad.h"
13 #include "TROOT.h"
14 
15 #include <string>
16 #include <map>
17 #include <iostream>
18 
19 //------------------------------------------------------------------------------
20 // User should change the following values
21 //
22 size_t numUniverses = 500;
23 std::string inputFileName("cmf_contour_3F_Th23Dmsq32_stat_only_500_uni.root");
24 std::string pars("Th23Dmsq32");
25 //
26 // User intervention stops here
27 //------------------------------------------------------------------------------
28 
29 
30 //------------------------------------------------------------------------------
31 TH2F* makeAndFillMedianDeltaChiSqrHist(std::vector<TH2F*> const& universeDeltaChiSqrHists)
32 {
33  // make the histogram we will return
34  TH2F* medianDeltaChiSqr = dynamic_cast<TH2F*>(universeDeltaChiSqrHists.front()->Clone("medianDeltaChiSqr"));
35  medianDeltaChiSqr->Reset();
36 
37  // loop over all the bins in the histograms and find the median for each bin
38  std::vector<float> binDeltaChiSqrDist;
39  for(int x = 1; x < universeDeltaChiSqrHists.front()->GetNbinsX() + 1; ++x){
40  for(int y = 1; y < universeDeltaChiSqrHists.front()->GetNbinsY() + 1; ++y){
41  binDeltaChiSqrDist.clear();
42  for(auto const* itr : universeDeltaChiSqrHists){
43  binDeltaChiSqrDist.emplace_back(itr->GetBinContent(x, y));
44  }
45 
46  // sort the values
47  std::sort(binDeltaChiSqrDist.begin(), binDeltaChiSqrDist.end());
48  medianDeltaChiSqr->SetBinContent(x, y, binDeltaChiSqrDist[binDeltaChiSqrDist.size() / 2]);
49  } // end loop over y bins
50  }// end loop over x bins
51 
52  return medianDeltaChiSqr;
53 
54 } // end makeAndFillMedianDeltaChiSqrHist
55 
56 //------------------------------------------------------------------------------
57 void makeHeatMapPlots(std::vector<TGraph*> const& median1SigmaContours,
58  std::vector<TGraph*> const& median2SigmaContours,
59  std::vector<TGraph*> const& median3SigmaContours,
60  std::vector<TGraph*> const& asimov1SigmaContours,
61  std::vector<TGraph*> const& asimov2SigmaContours,
62  std::vector<TGraph*> const& asimov3SigmaContours,
63  std::vector<TH2F*> const& heatMaps)
64 {
65  // make a TCanvas and divide it into 2 x 2
66  TCanvas *heatMapCanv = new TCanvas("heatMapCanv", "heatMapCanv", 1200, 800);
67 
68  heatMapCanv->Divide(2, 2);
69 
70  for(size_t c = 0; c < heatMaps.size(); ++c){
71  heatMapCanv->cd(c + 1);
72  heatMaps[c]->Draw("colz");
73 
74  if(c == 0){
75  for(auto& itr : median1SigmaContours)
76  itr->Draw("lsame");
77 
78  for(auto& itr : asimov1SigmaContours)
79  itr->Draw("lsame");
80  }
81  if(c == 1){
82  for(auto& itr : median2SigmaContours)
83  itr->Draw("lsame");
84 
85  for(auto& itr : asimov2SigmaContours)
86  itr->Draw("lsame");
87  }
88  if(c == 2){
89  for(auto& itr : median3SigmaContours)
90  itr->Draw("lsame");
91 
92  for(auto& itr : asimov3SigmaContours)
93  itr->Draw("lsame");
94  }
95 
96  }
97 
98  heatMapCanv->cd(4);
99 
100  TLegend *leg = new TLegend(0.1, 0.1, 0.9, 0.9);
101  leg->SetBorderSize(0);
102  leg->SetFillStyle(0);
103  leg->AddEntry(asimov1SigmaContours.front(), "Asimov - 1#sigma", "l");
104  leg->AddEntry(asimov2SigmaContours.front(), "Asimov - 2#sigma", "l");
105  leg->AddEntry(asimov3SigmaContours.front(), "Asimov - 3#sigma", "l");
106  leg->AddEntry(median1SigmaContours.front(), "Median - 1#sigma", "l");
107  leg->AddEntry(median2SigmaContours.front(), "Median - 2#sigma", "l");
108  leg->AddEntry(median3SigmaContours.front(), "Median - 3#sigma", "l");
109  leg->Draw();
110 
111  heatMapCanv->Update();
112 }
113 
114 //------------------------------------------------------------------------------
115 void makeUniverseContourPlots(TH2F* backdropHist,
116  std::vector<TGraph*> const& median1SigmaContours,
117  std::vector<TGraph*> const& median2SigmaContours,
118  std::vector<TGraph*> const& median3SigmaContours,
119  std::vector<TGraph*> const& asimov1SigmaContours,
120  std::vector<TGraph*> const& asimov2SigmaContours,
121  std::vector<TGraph*> const& asimov3SigmaContours,
122  std::vector<TGraph*> const& universe1SigmaContours,
123  std::vector<TGraph*> const& universe2SigmaContours,
124  std::vector<TGraph*> const& universe3SigmaContours)
125 {
126  // make a TCanvas and divide it into 2 x 2
127  TCanvas *contoursCanv = new TCanvas("contoursCanv", "contoursCanv", 1200, 800);
128 
129  contoursCanv->Divide(2, 2);
130 
131  for(size_t c = 0; c < 3; ++c){
132  contoursCanv->cd(c + 1);
133  backdropHist->SetTitle((std::to_string(c + 1) + "#sigma").c_str());
134  backdropHist->Draw();
135 
136  if(c == 0){
137  for(auto & itr : universe1SigmaContours){
138  itr->SetLineColor(kGray);
139  itr->SetLineWidth(1);
140  itr->Draw("lsame");
141  }
142 
143  for(auto & itr : median1SigmaContours)
144  itr->Draw("lsame");
145 
146  for(auto & itr : asimov1SigmaContours)
147  itr->Draw("lsame");
148  }
149  if(c == 1){
150  for(auto & itr : universe2SigmaContours){
151  itr->SetLineColor(kGray);
152  itr->SetLineWidth(1);
153  itr->Draw("lsame");
154  }
155 
156  for(auto & itr : median2SigmaContours)
157  itr->Draw("lsame");
158 
159  for(auto & itr : asimov2SigmaContours)
160  itr->Draw("lsame");
161  }
162  if(c == 2){
163  for(auto & itr : universe3SigmaContours){
164  itr->SetLineColor(kGray);
165  itr->SetLineWidth(1);
166  itr->Draw("lsame");
167  }
168 
169  for(auto & itr : median3SigmaContours)
170  itr->Draw("lsame");
171 
172  for(auto & itr : asimov3SigmaContours)
173  itr->Draw("lsame");
174  }
175  contoursCanv->Update();
176  }
177 
178  contoursCanv->cd(4);
179 
180  TLegend *leg = new TLegend(0.1, 0.1, 0.9, 0.9);
181  leg->SetBorderSize(0);
182  leg->SetFillStyle(0);
183  leg->AddEntry(asimov1SigmaContours.front(), "Asimov - 1#sigma", "l");
184  leg->AddEntry(asimov2SigmaContours.front(), "Asimov - 2#sigma", "l");
185  leg->AddEntry(asimov3SigmaContours.front(), "Asimov - 3#sigma", "l");
186  leg->AddEntry(median1SigmaContours.front(), "Median - 1#sigma", "l");
187  leg->AddEntry(median2SigmaContours.front(), "Median - 2#sigma", "l");
188  leg->AddEntry(median3SigmaContours.front(), "Median - 3#sigma", "l");
189  leg->Draw();
190 
191  contoursCanv->Update();
192 }
193 
194 //------------------------------------------------------------------------------
195 void FindContours(std::vector<TGraph*> & oneSigmaGrs,
196  std::vector<TGraph*> & twoSigmaGrs,
197  std::vector<TGraph*> & threeSigmaGrs,
198  TDirectory* directory)
199 {
200  //cout << "FindUniverseContours from directory " << directory << endl;
201 
202  std::string keyName;
203  TKey *key;
204  TIter next(directory->GetListOfKeys());
205  while((key = dynamic_cast<TKey*>(next()))){
206  keyName = key->GetName();
207 
208  if(keyName.find("1#sigma") != std::string::npos){
209  oneSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
210  oneSigmaGrs.back()->SetLineColor(kBlue);
211  oneSigmaGrs.back()->SetLineWidth(3);
212  }
213  else if(keyName.find("2#sigma") != std::string::npos){
214  twoSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
215  twoSigmaGrs.back()->SetLineColor(kRed);
216  twoSigmaGrs.back()->SetLineWidth(3);
217  }
218  else if(keyName.find("3#sigma") != std::string::npos){
219  threeSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
220  threeSigmaGrs.back()->SetLineColor(kMagenta);
221  threeSigmaGrs.back()->SetLineWidth(3);
222  }
223  } // end loop to find the Xsigma graphs and Delta chi^2 surfaces
224 
225  //cout << "there are " << 1SigmaGrs.size() << " 1 sigma contours "
226  // << 2SigmaGrs.size() << " 2 sigma contours "
227  // << 3SigmaGrs.size() << " 3 sigma contours." << endl;
228 }
229 
230 
231 //------------------------------------------------------------------------------
232 void getContoursFromDeltaChiSqr(TH2F* medianDeltaChiSqr,
233  std::vector<std::vector<TGraph*> > & contourGraphs)
234 {
235  std::vector<double> contourLevel({2.3, 6.18, 11.83});
236 
237  medianDeltaChiSqr->SetContour(contourLevel.size(), contourLevel.data());
238 
239  // This avoids disturbing canvases already
240  // Yes, we do need a canvas for this to work
241  // yes, we do need to call Update() on it
242  // No, ROOT is not very friendly.
243  TCanvas temp_can;
244  temp_can.cd();
245  medianDeltaChiSqr->Draw("contlist");
246  temp_can.Update();
247 
248  auto *plah = dynamic_cast<TObjArray *>(gROOT->GetListOfSpecials()->FindObject("contours"));
249  if(!plah) return;
250 
251  //cout << "TObjArray size: " << plah->GetSize() << endl;
252  contourGraphs.resize(contourLevel.size());
253 
254  for(int i = 0; i < plah->GetSize(); ++i){
255  std::vector<TGraph*> temp;
256  auto *list = dynamic_cast<TList*>(plah->At(i));
257  if(!list) break;
258  if(!(list->First())) break;
259  for(int igraph = 0; igraph < list->GetSize(); ++igraph){
260  if(list->At(igraph)) temp.emplace_back(dynamic_cast<TGraph*>(list->At(igraph)));
261  //cout << i << " " << igraph << " " << temp.back() << " " << list->At(igraph) << endl;
262  //temp.back()->SetLineColor(lineColor);
263  temp.back()->SetLineStyle(9);
264  temp.back()->SetLineWidth(3);
265  if (i == 0) temp.back()->SetLineColor(kBlue);
266  else if(i == 1) temp.back()->SetLineColor(kRed);
267  else if(i == 2) temp.back()->SetLineColor(kMagenta);
268  }
269  contourGraphs[i] = std::move(temp);
270  }
271 
272  gROOT->GetListOfSpecials()->FindObject("contours")->Clear();
273 }
274 
275 //------------------------------------------------------------------------------
276 void makeMedianCanvas(TH2F* backdropHist,
277  std::vector<TGraph*> const& median1SigmaGrs,
278  std::vector<TGraph*> const& median2SigmaGrs,
279  std::vector<TGraph*> const& median3SigmaGrs,
280  std::vector<TGraph*> const& asimov1SigmaGrs,
281  std::vector<TGraph*> const& asimov2SigmaGrs,
282  std::vector<TGraph*> const& asimov3SigmaGrs)
283 {
284  TCanvas *medianCanv = new TCanvas("medianCanv", "medianCanv", 1200, 800);
285 
286  // draw a histogram defining the boundaries of the space
287  backdropHist->Draw();
288 
289  for(auto & itr : median1SigmaGrs) itr->Draw("lsame");
290  for(auto & itr : median2SigmaGrs) itr->Draw("lsame");
291  for(auto & itr : median3SigmaGrs) itr->Draw("lsame");
292 
293  for(auto & itr : asimov1SigmaGrs) itr->Draw("lsame");
294  for(auto & itr : asimov2SigmaGrs) itr->Draw("lsame");
295  for(auto & itr : asimov3SigmaGrs) itr->Draw("lsame");
296 
297  TLegend *leg = new TLegend(0.6, 0.6, 0.9, 0.9);
298  leg->SetBorderSize(0);
299  leg->SetFillStyle(0);
300 
301  leg->AddEntry(asimov1SigmaGrs[0], "Asimov - 1#sigma", "l");
302  leg->AddEntry(asimov2SigmaGrs[0], "Asimov - 2#sigma", "l");
303  leg->AddEntry(asimov3SigmaGrs[0], "Asimov - 3#sigma", "l");
304  leg->AddEntry(median1SigmaGrs[0], "Median - 1#sigma", "l");
305  leg->AddEntry(median2SigmaGrs[0], "Median - 2#sigma", "l");
306  leg->AddEntry(median3SigmaGrs[0], "Median - 3#sigma", "l");
307  leg->Draw();
308 
309  medianCanv->Update();
310 }
311 
312 //------------------------------------------------------------------------------
314 {
315  // open the TFile containing the contours for the
316  // different universes
317  TFile *tf = new TFile(inputFileName.c_str(), "READ");
318 
319  // get the asimov contours
320  std::vector<TGraph*> asimov1SigmaContours;
321  std::vector<TGraph*> asimov2SigmaContours;
322  std::vector<TGraph*> asimov3SigmaContours;
323 
324  FindContours(asimov1SigmaContours,
325  asimov2SigmaContours,
326  asimov3SigmaContours,
327  tf->GetDirectory("contours/Asimov"));
328 
329  // fill the vector of Delta chi^2 histograms for the different universes
330  std::vector<TH2F*> universeDeltaChiSqrHists;
331  std::vector<TGraph*> universe1SigmaGrs;
332  std::vector<TGraph*> universe2SigmaGrs;
333  std::vector<TGraph*> universe3SigmaGrs;
334 
336  std::string grName;
338  for(size_t u = 0; u < numUniverses; ++u){
339  dirName = "contours/Universe_" + std::to_string(u) + "/";
340  histName = dirName + "Universe_" + std::to_string(u) + pars + "DeltaChiSqr";
341  universeDeltaChiSqrHists.emplace_back(dynamic_cast<TH2F*>(tf->Get(histName.c_str())));
342 
343  FindContours(universe1SigmaGrs,
344  universe2SigmaGrs,
345  universe3SigmaGrs,
346  tf->GetDirectory(dirName.c_str()));
347  }
348 
349  // now get the median Delta chi^2 histogram
350  auto *medianDeltaChiSqr = makeAndFillMedianDeltaChiSqrHist(universeDeltaChiSqrHists);
351 
352  std::vector<std::vector<TGraph*> > medianContours;
353  getContoursFromDeltaChiSqr(medianDeltaChiSqr, medianContours);
354 
355  // make a generic back drop histogram
356  TH2F *backdropHist = dynamic_cast<TH2F*>(universeDeltaChiSqrHists.front()->Clone("BackDropHist"));
357  backdropHist->Reset();
358 
359  makeMedianCanvas(backdropHist,
360  medianContours[0],
361  medianContours[1],
362  medianContours[2],
363  asimov1SigmaContours,
364  asimov2SigmaContours,
365  asimov3SigmaContours);
366 
367  makeUniverseContourPlots(backdropHist,
368  medianContours[0],
369  medianContours[1],
370  medianContours[2],
371  asimov1SigmaContours,
372  asimov2SigmaContours,
373  asimov3SigmaContours,
374  universe1SigmaGrs,
375  universe2SigmaGrs,
376  universe3SigmaGrs);
377 
378 
379  // get the heat maps
380  std::vector<TH2F*> heatMaps;
381  for(size_t cl = 0; cl < 3; ++cl){
382  histName = "contours/Median/RandomUniverses" + pars + std::to_string(cl + 1) + "sigmaHeatMap";
383  heatMaps.emplace_back(dynamic_cast<TH2F*>(tf->Get(histName.c_str())));
384  }
385 
386  makeHeatMapPlots(medianContours[0],
387  medianContours[1],
388  medianContours[2],
389  asimov1SigmaContours,
390  asimov2SigmaContours,
391  asimov3SigmaContours,
392  heatMaps);
393 
394 }
size_t numUniverses
enum BeamMode kRed
void makeMedianCanvas(TH2F *backdropHist, std::vector< TGraph * > const &median1SigmaGrs, std::vector< TGraph * > const &median2SigmaGrs, std::vector< TGraph * > const &median3SigmaGrs, std::vector< TGraph * > const &asimov1SigmaGrs, std::vector< TGraph * > const &asimov2SigmaGrs, std::vector< TGraph * > const &asimov3SigmaGrs)
std::string inputFileName("cmf_contour_3F_Th23Dmsq32_stat_only_500_uni.root")
string directory
projection from multiple chains
Timing fit.
void FindContours(std::vector< TGraph * > &oneSigmaGrs, std::vector< TGraph * > &twoSigmaGrs, std::vector< TGraph * > &threeSigmaGrs, TDirectory *directory)
void makeHeatMapPlots(std::vector< TGraph * > const &median1SigmaContours, std::vector< TGraph * > const &median2SigmaContours, std::vector< TGraph * > const &median3SigmaContours, std::vector< TGraph * > const &asimov1SigmaContours, std::vector< TGraph * > const &asimov2SigmaContours, std::vector< TGraph * > const &asimov3SigmaContours, std::vector< TH2F * > const &heatMaps)
void makeUniverseContourPlots(TH2F *backdropHist, std::vector< TGraph * > const &median1SigmaContours, std::vector< TGraph * > const &median2SigmaContours, std::vector< TGraph * > const &median3SigmaContours, std::vector< TGraph * > const &asimov1SigmaContours, std::vector< TGraph * > const &asimov2SigmaContours, std::vector< TGraph * > const &asimov3SigmaContours, std::vector< TGraph * > const &universe1SigmaContours, std::vector< TGraph * > const &universe2SigmaContours, std::vector< TGraph * > const &universe3SigmaContours)
TH2F * makeAndFillMedianDeltaChiSqrHist(std::vector< TH2F * > const &universeDeltaChiSqrHists)
std::string dirName
Definition: PlotSpectra.h:47
std::string pars("Th23Dmsq32")
void getContoursFromDeltaChiSqr(TH2F *medianDeltaChiSqr, std::vector< std::vector< TGraph * > > &contourGraphs)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
enum BeamMode kBlue
void next()
Definition: show_event.C:84
void makeMedianContours()
enum BeamMode string