compareAsimovUniverseCL.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 size_t numUniverses = 500;
20 double inputX = 0.5;
21 double inputY = 2.51;
22 
23 //------------------------------------------------------------------------------
24 void FindUniverseContours(std::vector<TGraph*> & oneSigmaGrs,
25  std::vector<TGraph*> & twoSigmaGrs,
26  std::vector<TGraph*> & threeSigmaGrs,
27  TDirectory* directory)
28 {
29  //cout << "FindUniverseContours from directory " << directory << endl;
30 
31  std::string keyName;
32  TKey *key;
33  TIter next(directory->GetListOfKeys());
34  while((key = dynamic_cast<TKey*>(next()))){
35  keyName = key->GetName();
36 
37  if(keyName.find("1#sigma") != std::string::npos){
38  oneSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
39  }
40  else if(keyName.find("2#sigma") != std::string::npos){
41  twoSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
42  }
43  else if(keyName.find("3#sigma") != std::string::npos){
44  threeSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
45  }
46  } // end loop to find the Xsigma graphs and Delta chi^2 surfaces
47 
48  //cout << "there are " << oneSigmaGrs.size() << " 1 sigma contours "
49  // << twoSigmaGrs.size() << " 2 sigma contours "
50  // << threeSigmaGrs.size() << " 3 sigma contours." << endl;
51 }
52 
53 //------------------------------------------------------------------------------
55  std::map<size_t, TH2F*> & deltaChiSqrMap,
56  TDirectory* directory,
57  std::string const& keyName)
58 {
59  deltaChiSqrMap.emplace(u, dynamic_cast<TH2F*>(directory->Get(keyName.c_str())));
60  deltaChiSqrMap[u]->SetDirectory(0);
61 
62  ////cout << "deltaChiSqr from directory " << deltaChiSqrMap[u] << endl;
63 }
64 
65 //------------------------------------------------------------------------------
67  size_t u,
68  std::map<size_t, TH1D*> & universeAsimovSpectrum,
69  std::map<size_t, TH1D*> & universeFHCSpectrum,
70  std::map<size_t, TH1D*> & universeRHCSpectrum)
71 {
72  // get objects from the random universe file
73  auto *directory = dynamic_cast<TDirectory*>(inputFile->Get(("fake/Universe_" + std::to_string(u)).c_str()));
74 
75  //cout << "loop over directory " << directory << endl;
76 
77  std::string keyName;
78  TKey *key;
79  TIter next(directory->GetListOfKeys());
80  while((key = dynamic_cast<TKey*>(next()))){
81  keyName = key->GetName();
82 
83  if (keyName.find("Asimov") != std::string::npos) universeAsimovSpectrum.emplace(u, dynamic_cast<TH1D*>(key->ReadObj()));
84  else if(keyName.find("Poisson") != std::string::npos){
85  universeFHCSpectrum.emplace(u, dynamic_cast<TH1D*>(key->ReadObj()->Clone(("FHCUniverse_" + std::to_string(u)).c_str())));
86  universeRHCSpectrum.emplace(u, dynamic_cast<TH1D*>(key->ReadObj()->Clone(("RHCUniverse_" + std::to_string(u)).c_str())));
87  }
88  }
89 
90  ////cout << "got spectra " << universeAsimovSpectrum[u] << " " << universeFHCSpectrum[u] << " " << universeRHCSpectrum[u] << endl;
91 }
92 
93 //------------------------------------------------------------------------------
94 void FillUniverseContainers(TFile *contourFile,
95  TFile *randomFile,
96  size_t numUniverses,
97  std::map<size_t, std::vector<TGraph*>> & universe1sigmaGrs,
98  std::map<size_t, std::vector<TGraph*>> & universe2sigmaGrs,
99  std::map<size_t, std::vector<TGraph*>> & universe3sigmaGrs,
100  std::map<size_t, TH2F*> & universeDeltaChiSqr,
101  std::map<size_t, TH1D*> & universeAsimovSpectrum,
102  std::map<size_t, TH1D*> & universeFHCSpectrum,
103  std::map<size_t, TH1D*> & universeRHCSpectrum)
104 {
105  TDirectory *td = nullptr;
106 
107  std::vector<TGraph*> grVec1;
108  std::vector<TGraph*> grVec2;
109  std::vector<TGraph*> grVec3;
110 
111  for(size_t u = 0; u < numUniverses; ++u){
112 
113  // get the objects from the contour file
114  td = dynamic_cast<TDirectory*>(contourFile->Get(("contours/UniverseContours_" + std::to_string(u) + "_").c_str()));
115 
116  grVec1.clear();
117  grVec2.clear();
118  grVec3.clear();
119 
120  FindUniverseContours(grVec1,
121  grVec2,
122  grVec3,
123  td);
124 
125  universe1sigmaGrs.emplace(u, grVec1);
126  universe2sigmaGrs.emplace(u, grVec2);
127  universe3sigmaGrs.emplace(u, grVec3);
128 
130  universeDeltaChiSqr,
131  td,
132  ("UniverseContours_" + std::to_string(u) + "_Th23Dmsq32DeltaChiSqr").c_str());
133 
134  //cout << u << " " << universeDeltaChiSqr[u] << endl;
135 
136  FindUniverseSpectra(randomFile,
137  u,
138  universeAsimovSpectrum,
139  universeFHCSpectrum,
140  universeRHCSpectrum);
141 
142  //cout << "spectra in containers" << endl;
143 
144  } // end loop over universes
145 
146 }
147 
148 //------------------------------------------------------------------------------
149 void DrawBestFitsOverContours(TH2F* limitsHist,
150  TGraph *bestFitType1Gr,
151  TGraph *bestFitType2Gr,
152  std::string const& pdfFileName,
153  std::string const& type1Legend,
154  std::string const& type2Legend,
155  std::vector<TGraph*> const& asimovType11sigmaGrs,
156  std::vector<TGraph*> const& asimovType12sigmaGrs,
157  std::vector<TGraph*> const& asimovType13sigmaGrs,
158  std::vector<TGraph*> const& medianType11sigmaGrs,
159  std::vector<TGraph*> const& medianType12sigmaGrs,
160  std::vector<TGraph*> const& medianType13sigmaGrs,
161  std::vector<TGraph*> const& asimovType21sigmaGrs,
162  std::vector<TGraph*> const& asimovType22sigmaGrs,
163  std::vector<TGraph*> const& asimovType23sigmaGrs,
164  std::vector<TGraph*> const& medianType21sigmaGrs,
165  std::vector<TGraph*> const& medianType22sigmaGrs,
166  std::vector<TGraph*> const& medianType23sigmaGrs)
167 {
168  TMarker *marker = new TMarker(inputX, inputY, 30);
169  marker->SetMarkerColor(kBlack);
170  marker->SetMarkerSize(3);
171 
172  TCanvas *canv = new TCanvas("bestFitCanv");
173 
174  limitsHist->Draw();
175  marker->Draw("same");
176 
177  TLegend *leg1 = new TLegend(0.25, 0.67, 0.35, 0.9);
178  leg1->SetFillStyle(0);
179  leg1->SetBorderSize(0);
180  leg1->AddEntry(asimovType21sigmaGrs.front(), "1#sigma", "l");
181  leg1->AddEntry(asimovType22sigmaGrs.front(), "2#sigma", "l");
182  leg1->AddEntry(asimovType23sigmaGrs.front(), "3#sigma", "l");
183  leg1->Draw();
184 
185  TLegend *leg2 = new TLegend(0.3, 0.67, 0.65, 0.9);
186  leg2->SetFillStyle(0);
187  leg2->SetBorderSize(0);
188  leg2->AddEntry(asimovType11sigmaGrs.front(), (type1Legend + " Asimov").c_str(), "l");
189  leg2->AddEntry(asimovType21sigmaGrs.front(), (type2Legend + " Asimov").c_str(), "l");
190  if(medianType13sigmaGrs.size() > 0) leg2->AddEntry(medianType13sigmaGrs.front(), (type1Legend + " Median").c_str(), "l");
191  if(medianType23sigmaGrs.size() > 0) leg2->AddEntry(medianType23sigmaGrs.front(), (type2Legend + " Median").c_str(), "l");
192  leg2->Draw();
193 
194  for(auto const& itr : asimovType11sigmaGrs){
195  itr->SetLineColor(kBlue);
196  itr->SetLineStyle(3);
197  itr->Draw("same");
198  }
199  for(auto const& itr : asimovType12sigmaGrs){
200  itr->SetLineColor(kRed);
201  itr->SetLineStyle(3);
202  itr->Draw("same");
203  }
204  for(auto const& itr : asimovType13sigmaGrs){
205  itr->SetLineColor(kMagenta);
206  itr->SetLineStyle(3);
207  itr->Draw("same");
208  }
209 
210  for(auto const& itr : medianType11sigmaGrs){
211  itr->SetLineColor(kBlue);
212  itr->SetLineStyle(7);
213  itr->Draw("same");
214  }
215  for(auto const& itr : medianType12sigmaGrs){
216  itr->SetLineColor(kRed);
217  itr->SetLineStyle(7);
218  itr->Draw("same");
219  }
220  for(auto const& itr : medianType13sigmaGrs){
221  itr->SetLineColor(kMagenta);
222  itr->SetLineStyle(7);
223  itr->Draw("same");
224  }
225 
226  for(auto const& itr : asimovType21sigmaGrs){
227  itr->SetLineColor(kBlue);
228  itr->Draw("same");
229  }
230  for(auto const& itr : asimovType22sigmaGrs){
231  itr->SetLineColor(kRed);
232  itr->Draw("same");
233  }
234  for(auto const& itr : asimovType23sigmaGrs){
235  itr->SetLineColor(kMagenta);
236  itr->Draw("same");
237  }
238 
239  for(auto const& itr : medianType21sigmaGrs){
240  itr->SetLineColor(kBlue);
241  itr->SetLineStyle(5);
242  itr->Draw("same");
243  }
244  for(auto const& itr : medianType22sigmaGrs){
245  itr->SetLineColor(kRed);
246  itr->SetLineStyle(5);
247  itr->Draw("same");
248  }
249  for(auto const& itr : medianType23sigmaGrs){
250  itr->SetLineColor(kMagenta);
251  itr->SetLineStyle(5);
252  itr->Draw("same");
253  }
254 
255  canv->Print(pdfFileName.c_str(), "pdf");
256 
257  std::string modifiedFileName(pdfFileName);
258  modifiedFileName.erase(modifiedFileName.find("("), 1);
259 
260  // now draw the Asimov contours for the Poisson LL and the covmat only
261  TCanvas *canv1 = new TCanvas("asimovCanv");
262 
263  limitsHist->Draw();
264  marker->Draw("same");
265 
266  TLegend *leg3 = new TLegend(0.25, 0.67, 0.35, 0.9);
267  leg3->SetFillStyle(0);
268  leg3->SetBorderSize(0);
269  leg3->AddEntry(asimovType21sigmaGrs.front(), "1#sigma", "l");
270  leg3->AddEntry(asimovType22sigmaGrs.front(), "2#sigma", "l");
271  leg3->AddEntry(asimovType23sigmaGrs.front(), "3#sigma", "l");
272  leg3->Draw();
273 
274  TLegend *leg4 = new TLegend(0.3, 0.67, 0.65, 0.9);
275  leg4->SetFillStyle(0);
276  leg4->SetBorderSize(0);
277  leg4->AddEntry(asimovType11sigmaGrs.front(), (type1Legend + " Asimov").c_str(), "l");
278  leg4->AddEntry(asimovType21sigmaGrs.front(), (type2Legend + " Asimov").c_str(), "l");
279  leg4->Draw();
280 
281  for(auto const& itr : asimovType11sigmaGrs){
282  itr->SetLineColor(kBlue);
283  itr->SetLineStyle(3);
284  itr->Draw("same");
285  }
286  for(auto const& itr : asimovType12sigmaGrs){
287  itr->SetLineColor(kRed);
288  itr->SetLineStyle(3);
289  itr->Draw("same");
290  }
291  for(auto const& itr : asimovType13sigmaGrs){
292  itr->SetLineColor(kMagenta);
293  itr->SetLineStyle(3);
294  itr->Draw("same");
295  }
296 
297  for(auto const& itr : asimovType21sigmaGrs){
298  itr->SetLineColor(kBlue);
299  itr->Draw("same");
300  }
301  for(auto const& itr : asimovType22sigmaGrs){
302  itr->SetLineColor(kRed);
303  itr->Draw("same");
304  }
305  for(auto const& itr : asimovType23sigmaGrs){
306  itr->SetLineColor(kMagenta);
307  itr->Draw("same");
308  }
309 
310  canv1->Print(modifiedFileName.c_str(), "pdf");
311 
312  // now draw the median contours for the Poisson LL and the covmat only
313  TCanvas *canv2 = new TCanvas("medianCanv");
314 
315  limitsHist->Draw();
316  marker->Draw("same");
317 
318  TLegend *leg5 = new TLegend(0.25, 0.67, 0.35, 0.9);
319  leg5->SetFillStyle(0);
320  leg5->SetBorderSize(0);
321  if(medianType21sigmaGrs.size() > 0) leg5->AddEntry(medianType21sigmaGrs.front(), "1#sigma", "l");
322  if(medianType22sigmaGrs.size() > 0) leg5->AddEntry(medianType22sigmaGrs.front(), "2#sigma", "l");
323  if(medianType23sigmaGrs.size() > 0) leg5->AddEntry(medianType23sigmaGrs.front(), "3#sigma", "l");
324  leg5->Draw();
325 
326  TLegend *leg6 = new TLegend(0.3, 0.67, 0.65, 0.9);
327  leg6->SetFillStyle(0);
328  leg6->SetBorderSize(0);
329  if(medianType13sigmaGrs.size() > 0) leg6->AddEntry(medianType13sigmaGrs.front(), (type1Legend + " Median").c_str(), "l");
330  if(medianType23sigmaGrs.size() > 0) leg6->AddEntry(medianType23sigmaGrs.front(), (type2Legend + " Median").c_str(), "l");
331  leg6->Draw();
332 
333  if(medianType11sigmaGrs.size() > 0) for(auto const& itr : medianType11sigmaGrs) itr->Draw("same");
334  if(medianType12sigmaGrs.size() > 0) for(auto const& itr : medianType12sigmaGrs) itr->Draw("same");
335  if(medianType13sigmaGrs.size() > 0) for(auto const& itr : medianType13sigmaGrs) itr->Draw("same");
336 
337  if(medianType21sigmaGrs.size() > 0) for(auto const& itr : medianType21sigmaGrs) itr->Draw("same");
338  if(medianType22sigmaGrs.size() > 0) for(auto const& itr : medianType22sigmaGrs) itr->Draw("same");
339  if(medianType23sigmaGrs.size() > 0) for(auto const& itr : medianType23sigmaGrs) itr->Draw("same");
340 
341  canv2->Print(modifiedFileName.c_str(), "pdf");
342 
343  // now draw the Asimov and median contours for the Poisson LL
344  TCanvas *canv3 = new TCanvas((type2Legend + "Canv").c_str());
345 
346  limitsHist->Draw();
347  marker->Draw("same");
348 
349  TLegend *leg7 = new TLegend(0.25, 0.67, 0.35, 0.9);
350  leg7->SetFillStyle(0);
351  leg7->SetBorderSize(0);
352  leg7->AddEntry(asimovType21sigmaGrs.front(), "1#sigma", "l");
353  leg7->AddEntry(asimovType22sigmaGrs.front(), "2#sigma", "l");
354  leg7->AddEntry(asimovType23sigmaGrs.front(), "3#sigma", "l");
355  leg7->Draw();
356 
357  TLegend *leg8 = new TLegend(0.3, 0.67, 0.65, 0.9);
358  leg8->SetFillStyle(0);
359  leg8->SetBorderSize(0);
360  leg8->AddEntry(asimovType23sigmaGrs.front(), (type2Legend + " Asimov").c_str(), "l");
361  if(medianType23sigmaGrs.size() > 0) leg8->AddEntry(medianType23sigmaGrs.front(), (type2Legend + " Median").c_str(), "l");
362  leg8->Draw();
363 
364  for(auto const& itr : asimovType21sigmaGrs) itr->Draw("same");
365  for(auto const& itr : asimovType22sigmaGrs) itr->Draw("same");
366  for(auto const& itr : asimovType23sigmaGrs) itr->Draw("same");
367 
368  if(medianType21sigmaGrs.size() > 0) for(auto const& itr : medianType21sigmaGrs) itr->Draw("same");
369  if(medianType22sigmaGrs.size() > 0) for(auto const& itr : medianType22sigmaGrs) itr->Draw("same");
370  if(medianType23sigmaGrs.size() > 0) for(auto const& itr : medianType23sigmaGrs) itr->Draw("same");
371  bestFitType2Gr->Draw("psame");
372 
373  canv3->Print(modifiedFileName.c_str(), "pdf");
374 
375  // now draw the Asimov and median contours for the covmat
376  TCanvas *canv4 = new TCanvas((type1Legend + "Canv").c_str());
377 
378  limitsHist->Draw();
379  marker->Draw("same");
380 
381  TLegend *leg9 = new TLegend(0.25, 0.67, 0.35, 0.9);
382  leg9->SetFillStyle(0);
383  leg9->SetBorderSize(0);
384  leg9->AddEntry(asimovType11sigmaGrs.front(), "1#sigma", "l");
385  leg9->AddEntry(asimovType12sigmaGrs.front(), "2#sigma", "l");
386  leg9->AddEntry(asimovType13sigmaGrs.front(), "3#sigma", "l");
387  leg9->Draw();
388 
389  TLegend *leg10 = new TLegend(0.3, 0.67, 0.65, 0.9);
390  leg10->SetFillStyle(0);
391  leg10->SetBorderSize(0);
392  leg10->AddEntry(asimovType11sigmaGrs.front(), (type1Legend + " Asimov").c_str(), "l");
393  if(medianType13sigmaGrs.size() > 0) leg10->AddEntry(medianType13sigmaGrs.front(), (type1Legend + " Median").c_str(), "l");
394  leg10->Draw();
395 
396  for(auto const& itr : asimovType11sigmaGrs) itr->Draw("same");
397  for(auto const& itr : asimovType12sigmaGrs) itr->Draw("same");
398  for(auto const& itr : asimovType13sigmaGrs) itr->Draw("same");
399 
400  if(medianType11sigmaGrs.size() > 0) for(auto const& itr : medianType11sigmaGrs) itr->Draw("same");
401  if(medianType12sigmaGrs.size() > 0) for(auto const& itr : medianType12sigmaGrs) itr->Draw("same");
402  if(medianType13sigmaGrs.size() > 0) for(auto const& itr : medianType13sigmaGrs) itr->Draw("same");
403 
404  bestFitType1Gr->Draw("psame");
405  canv4->Print(modifiedFileName.c_str(), "pdf");
406 
407  delete canv;
408  delete canv1;
409  delete canv2;
410  delete canv3;
411  delete canv4;
412 }
413 
414 //------------------------------------------------------------------------------
416  std::string const& pdfFileName,
417  std::string const& type1Legend,
418  std::string const& type2Legend,
419  std::vector<TGraph*> & asimov3sigmaGrs,
420  std::map<size_t, std::vector<TGraph*>> & universe1sigmaGrs,
421  std::map<size_t, std::vector<TGraph*>> & universe2sigmaGrs,
422  std::map<size_t, std::vector<TGraph*>> & universe3sigmaGrs,
423  std::map<size_t, std::vector<TGraph*>> & universeType21sigmaGrs,
424  std::map<size_t, std::vector<TGraph*>> & universeType22sigmaGrs,
425  std::map<size_t, std::vector<TGraph*>> & universeType23sigmaGrs,
426  std::map<size_t, TH2F*> & universeDeltaChiSqr,
427  std::map<size_t, TH1D*> & universeAsimovSpectrum,
428  std::map<size_t, TH1D*> & universeFHCSpectrum,
429  std::map<size_t, TH1D*> & universeRHCSpectrum)
430 {
431  TCanvas *canv = new TCanvas();
432  canv->cd();
433 
434  TPad *deltaChiSqrPad = new TPad("deltaChiSqrPad", "#Delta#chi^{2} Pad", 0.00, 0.60, 1.00, 1.00);
435  TPad *legendPad = new TPad("legendPad", "Legend Pad", 0.75, 0.00, 1.00, 0.60);
436  TPad *fhcSpectrumPad = new TPad("fhcSpectrumPad", "FHC Spectrum", 0.00, 0.30, 0.75, 0.60);
437  TPad *rhcSpectrumPad = new TPad("rhcSpectrumPad", "RHC Spectrum", 0.00, 0.00, 0.75, 0.30);
438 
439  deltaChiSqrPad->Draw();
440  legendPad ->Draw();
441  fhcSpectrumPad->Draw();
442  rhcSpectrumPad->Draw();
443 
444  // keep track of the number of times the input oscillation parameters
445  // have Delta chi^2 < 2.3
446  int cnt = 0;
447 
448  int bin = universeDeltaChiSqr[0]->FindBin(0.58, 2.51);
449 
450  for(size_t u = 0; u < numUniverses; ++u){
451 
452  cout << "Draw Universe " << u << endl;
453 
454  deltaChiSqrPad->cd();
455 
456  ////cout << universeDeltaChiSqr.find(u)->second << endl;
457 
458  universeDeltaChiSqr[u]->SetMinimum(0.);
459  universeDeltaChiSqr[u]->SetMaximum(15.);
460  universeDeltaChiSqr[u]->Draw("colz");
461  for(auto const& itr : asimov3sigmaGrs){
462  itr->SetLineColor(kMagenta);
463  itr->SetLineStyle(3);
464  itr->Draw("lsame");
465  }
466 
467  ////cout << "asimov contours drawn on universe " << u << endl;
468 
469  if(universeDeltaChiSqr[u]->GetBinContent(bin) < 2.3) ++cnt;
470 
471  ////cout << u << " " << universeDeltaChiSqr[u]->GetBinContent(universeDeltaChiSqr[u]->FindBin(0.58, 2.51)) << endl;
472 
473  for(auto const& itr : universe1sigmaGrs[u]){
474  itr->SetLineColor(kBlue);
475  itr->SetLineWidth(2);
476  itr->Draw("same");
477  }
478 
479  for(auto const& itr : universe2sigmaGrs[u]){
480  itr->SetLineColor(kRed);
481  itr->SetLineWidth(2);
482  itr->Draw("same");
483  }
484 
485  for(auto const& itr : universe3sigmaGrs[u]){
486  itr->SetLineColor(kMagenta);
487  itr->SetLineWidth(2);
488  itr->Draw("same");
489  }
490 
491  for(auto const& itr : universeType21sigmaGrs[u]){
492  itr->SetLineColor(kBlue);
493  itr->SetLineWidth(2);
494  itr->SetLineStyle(2);
495  itr->Draw("same");
496  }
497 
498  for(auto const& itr : universeType22sigmaGrs[u]){
499  itr->SetLineColor(kRed);
500  itr->SetLineWidth(2);
501  itr->SetLineStyle(2);
502  itr->Draw("same");
503  }
504 
505  for(auto const& itr : universeType23sigmaGrs[u]){
506  itr->SetLineColor(kMagenta);
507  itr->SetLineWidth(2);
508  itr->SetLineStyle(2);
509  itr->Draw("same");
510  }
511 
512  TMarker *marker = new TMarker(inputX, inputY, 30);
513  marker->SetMarkerColor(kBlack);
514  marker->SetMarkerSize(3);
515  marker->Draw();
516 
517  legendPad->cd();
518  TLegend *leg1 = new TLegend(0.1, 0.1, 0.9, 0.9);
519  leg1->SetBorderSize(0);
520  leg1->SetFillStyle(0);
521  leg1->AddEntry(asimov3sigmaGrs.front(), (type1Legend + " Asimov 3#sigma").c_str(), "l");
522  leg1->AddEntry(universe1sigmaGrs[u].front(), (type1Legend + " 1#sigma").c_str(), "l");
523  leg1->AddEntry(universe2sigmaGrs[u].front(), (type1Legend + " 2#sigma").c_str(), "l");
524  leg1->AddEntry(universe3sigmaGrs[u].front(), (type1Legend + " 3#sigma").c_str(), "l");
525  leg1->AddEntry(universeType21sigmaGrs[u].front(), (type2Legend + " 1#sigma").c_str(), "l");
526  leg1->AddEntry(universeType22sigmaGrs[u].front(), (type2Legend + " 2#sigma").c_str(), "l");
527  leg1->AddEntry(universeType23sigmaGrs[u].front(), (type2Legend + " 3#sigma").c_str(), "l");
528  leg1->AddEntry(marker, "Input Point", "p");
529  leg1->Draw();
530 
531  fhcSpectrumPad->cd();
532 
533  universeFHCSpectrum[u]->GetXaxis()->SetRangeUser(93, 180);
534  universeFHCSpectrum[u]->SetLineColor(kRed);
535  universeFHCSpectrum[u]->SetMaximum(25);
536  universeFHCSpectrum[u]->Draw("hist");
537  universeAsimovSpectrum[u]->Draw("histsame");
538 
539  TLegend *leg2 = new TLegend(0.5, 0.5, 0.8, 0.8);
540  leg2->SetBorderSize(0);
541  leg2->SetFillStyle(0);
542  leg2->AddEntry(universeAsimovSpectrum[u], "Asimov Spectrum", "l");
543  leg2->AddEntry(universeFHCSpectrum[u], "Fluctuated Spectrum", "l");
544  leg2->Draw();
545 
546  rhcSpectrumPad->cd();
547 
548  universeRHCSpectrum[u]->GetXaxis()->SetRangeUser(274, 361);
549  universeRHCSpectrum[u]->SetLineColor(kRed);
550  universeRHCSpectrum[u]->SetMaximum(15);
551  universeRHCSpectrum[u]->Draw("hist");
552  universeAsimovSpectrum[u]->Draw("histsame");
553 
554  if(u == numUniverses-1) canv->Print((pdfFileName + ")").c_str(), "pdf");
555  else canv->Print(pdfFileName.c_str(), "pdf");
556 
557 
558  delete marker;
559  delete leg1;
560  delete leg2;
561 
562  deltaChiSqrPad->Clear();
563  fhcSpectrumPad->Clear();
564  rhcSpectrumPad->Clear();
565  legendPad ->Clear();
566  }
567 
568  delete deltaChiSqrPad;
569  delete fhcSpectrumPad;
570  delete rhcSpectrumPad;
571  delete legendPad;
572  delete canv;
573 
574  cout << cnt << endl;
575 }
576 
577 //------------------------------------------------------------------------------
578 void DrawUniverseNSigmaContours(TH2F* limitsHist,
579  std::vector<TGraph*> const& median1sigmaGrs,
580  std::vector<TGraph*> const& median2sigmaGrs,
581  std::vector<TGraph*> const& median3sigmaGrs,
582  std::string const& pdfFileName,
583  std::map<size_t, std::vector<TGraph*>> const& universeGrs)
584 {
585  TCanvas *canv = new TCanvas("nsigmaCanv");
586  limitsHist->Draw();
587  for(auto const& uitr : universeGrs){
588  for(auto const& gritr : uitr.second){
589  gritr->Draw("same");
590  }
591  } // end loop over universes
592 
593  for(auto const& itr : median1sigmaGrs){
594  itr->SetLineColor(kBlue);
595  itr->SetLineStyle(7);
596  itr->SetLineWidth(3);
597  itr->Draw("same");
598  }
599  for(auto const& itr : median2sigmaGrs){
600  itr->SetLineColor(kRed);
601  itr->SetLineStyle(7);
602  itr->SetLineWidth(3);
603  itr->Draw("same");
604  }
605  for(auto const& itr : median3sigmaGrs){
606  itr->SetLineColor(kMagenta);
607  itr->SetLineStyle(7);
608  itr->SetLineWidth(3);
609  itr->Draw("same");
610  }
611 
612  canv->Print(pdfFileName.c_str(), "pdf");
613  delete canv;
614 }
615 
616 //------------------------------------------------------------------------------
617 void FindMedianCLContours(TH2F* clHeatMap,
618  std::vector<std::vector<TGraph*>> & contourGraphs,
619  int lineColor)
620 {
621  std::vector<double> contourLevel({1.});
622 
623  clHeatMap->SetContour(contourLevel.size(), contourLevel.data());
624 
625  // This avoids disturbing canvases already
626  // Yes, we do need a canvas for this to work
627  // yes, we do need to call Update() on it
628  // No, ROOT is not very friendly.
629  TCanvas temp_can;
630  temp_can.cd();
631  clHeatMap->Draw("contlist");
632  temp_can.Update();
633 
634  auto *plah = dynamic_cast<TObjArray *>(gROOT->GetListOfSpecials()->FindObject("contours"));
635  if(!plah) return;
636 
637  //cout << "TObjArray size: " << plah->GetSize() << endl;
638  contourGraphs.resize(contourLevel.size());
639 
640  for(int i = 0; i < plah->GetSize(); ++i){
641  std::vector<TGraph*> temp;
642  auto *list = dynamic_cast<TList*>(plah->At(i));
643  if(!list) break;
644  if(!(list->First())) break;
645  for(int igraph = 0; igraph < list->GetSize(); ++igraph){
646  if(list->At(igraph)) temp.emplace_back(dynamic_cast<TGraph*>(list->At(igraph)));
647  //cout << i << " " << igraph << " " << temp.back() << " " << list->At(igraph) << endl;
648  temp.back()->SetLineColor(lineColor);
649  if(i == 0) temp.back()->SetLineStyle(2);
650  if(i == 2) temp.back()->SetLineStyle(5);
651  }
652  contourGraphs[i] = std::move(temp);
653  }
654 
655  gROOT->GetListOfSpecials()->FindObject("contours")->Clear();
656 }
657 
658 //------------------------------------------------------------------------------
659 void DrawMedianContours(std::map<size_t, TH2F*> & universeDeltaChiSqr,
660  std::string const& pdfFileName,
661  std::vector<TGraph*> const& median1sigmaGrs,
662  std::vector<TGraph*> const& median2sigmaGrs,
663  std::vector<TGraph*> const& median3sigmaGrs,
664  std::string const& legendBase)
665 {
666  std::vector<TH2F*> sigmaMaps({dynamic_cast<TH2F*>(universeDeltaChiSqr[0]->Clone("oneSigmaMap" )),
667  dynamic_cast<TH2F*>(universeDeltaChiSqr[0]->Clone("twoSigmaMap" )),
668  dynamic_cast<TH2F*>(universeDeltaChiSqr[0]->Clone("threeSigmaMap"))});
669  sigmaMaps[0]->Reset();
670  sigmaMaps[1]->Reset();
671  sigmaMaps[2]->Reset();
672  sigmaMaps[0]->SetTitle("1#sigma Map");
673  sigmaMaps[1]->SetTitle("2#sigma Map");
674  sigmaMaps[2]->SetTitle("3#sigma Map");
675 
676  double xVal;
677  double yVal;
678  TH2F *hist = nullptr;
679  for(auto const& itr : universeDeltaChiSqr){
680  hist = itr.second;
681  for(int x = 0; x < hist->GetNbinsX(); ++x){
682  xVal = hist->GetXaxis()->GetBinCenter(x);
683  for(int y = 0; y < hist->GetNbinsY(); ++y){
684  yVal = hist->GetYaxis()->GetBinCenter(y);
685 
686  if(hist->GetBinContent(x + 1, y + 1) < 2.30) sigmaMaps[0]->Fill(xVal, yVal);
687  if(hist->GetBinContent(x + 1, y + 1) < 6.18) sigmaMaps[1]->Fill(xVal, yVal);
688  if(hist->GetBinContent(x + 1, y + 1) < 11.83) sigmaMaps[2]->Fill(xVal, yVal);
689  }
690  }
691  } // end loop over universes
692 
693  std::string canvName;
694  for(auto itr : sigmaMaps){
695  itr->SetMaximum(100);
696  itr->SetMinimum(0);
697 
698  canvName = itr->GetName();
699  canvName += "HeatMapCanv";
700  TCanvas *canv = new TCanvas(canvName.c_str());
701  canv->SetGridx();
702  canv->SetGridy();
703  itr->Draw("colz");
704  canv->Print(pdfFileName.c_str(), "pdf");
705  delete canv;
706  }
707 
708  //cout << "heat maps drawn" << endl;
709 
710  std::vector< std::vector<TGraph*> > oneSigmaGr;
711  std::vector< std::vector<TGraph*> > twoSigmaGr;
712  std::vector< std::vector<TGraph*> > threeSigmaGr;
713 
714  FindMedianCLContours(sigmaMaps[0],
715  oneSigmaGr,
716  kBlue);
717  FindMedianCLContours(sigmaMaps[1],
718  twoSigmaGr,
719  kRed);
720  FindMedianCLContours(sigmaMaps[2],
721  threeSigmaGr,
722  kMagenta);
723 
724  //cout << " found all median CL contours" << endl;
725 
726  // make a histogram to draw the boundaries of the space
727  TH2F *boundingBox = dynamic_cast<TH2F*>(universeDeltaChiSqr[0]->Clone("boundingBox"));
728  boundingBox->Reset();
729 
730  TCanvas *canv = new TCanvas("medianContourCanv");
731  canv->SetGridy();
732  canv->SetGridx();
733  canv->cd();
734  boundingBox->Draw();
735 
736  TMarker *marker = new TMarker(inputX, inputY, 30);
737  marker->SetMarkerColor(kBlack);
738  marker->SetMarkerSize(3);
739  marker->Draw();
740 
741  for(auto & itr : median1sigmaGrs){
742  itr->SetLineWidth(3);
743  itr->SetLineStyle(1);
744  itr->Draw("same");
745  }
746  for(auto & itr : median2sigmaGrs){
747  itr->SetLineWidth(3);
748  itr->SetLineStyle(1);
749  itr->Draw("same");
750  }
751  for(auto & itr : median3sigmaGrs){
752  itr->SetLineWidth(3);
753  itr->SetLineStyle(1);
754  itr->Draw("same");
755  }
756 
757  //cout << "median sigma grs drawn" << endl;
758 
759  for(auto & itr : oneSigmaGr){
760  for(auto & gritr : itr) gritr->Draw("same");
761  }
762  //cout << " 1 sigma gr drawn" << endl;
763  for(auto & itr : twoSigmaGr){
764  for(auto & gritr : itr) gritr->Draw("same");
765  }
766  //cout << " 2 sigma gr drawn" << endl;
767 
768  for(auto & itr : threeSigmaGr){
769  for(auto & gritr : itr) gritr->Draw("same");
770  }
771 
772  //cout << " 3 sigma gr drawn" << endl;
773 
774  TLegend *leg1 = new TLegend(0.25, 0.67, 0.35, 0.9);
775  leg1->SetFillStyle(0);
776  leg1->SetBorderSize(0);
777  if(median1sigmaGrs.size() > 0) leg1->AddEntry(median1sigmaGrs.front(), "1#sigma", "l");
778  if(median2sigmaGrs.size() > 0) leg1->AddEntry(median2sigmaGrs.front(), "2#sigma", "l");
779  if(median3sigmaGrs.size() > 0) leg1->AddEntry(median3sigmaGrs.front(), "3#sigma", "l");
780  leg1->Draw();
781 
782  //cout << "leg1 drawn" << endl;
783 
784  TLegend *leg2 = new TLegend(0.3, 0.67, 0.65, 0.9);
785  leg2->SetFillStyle(0);
786  leg2->SetBorderSize(0);
787  if(median1sigmaGrs.size() > 0) leg2->AddEntry(median1sigmaGrs.front(), (legendBase + " Median").c_str(), "l");
788  leg2->AddEntry(oneSigmaGr.front().front(), (legendBase + " Heat Map > 1 Universe").c_str(), "l");
789  //leg2->AddEntry(oneSigmaGr[1].front(), (legendBase + " Heat Map 50%").c_str(), "l");
790  //leg2->AddEntry(oneSigmaGr.back().front(), (legendBase + " Heat Map 57%").c_str(), "l");
791  leg2->Draw();
792 
793  canv->Print(pdfFileName.c_str(), "pdf");
794  delete canv;
795 
796  //cout << "done with raw median contours for " << legendBase << endl;
797 }
798 
799 
800 //------------------------------------------------------------------------------
801 // Need to fix this function up
802 /* void DrawPoissonDistributionsByBin() */
803 /* { */
804 /* // loop over the bins and draw the projection of the poisson draws */
805 /* // along with the poisson distribution for the mean given by the Asimov */
806 /* // distribution */
807 /* TCanvas *poissonCanv = new TCanvas("poissonCanv"); */
808 /* TH1D *histBin; */
809 /* std::vector<double> counts(300); */
810 /* std::vector<double> prob(300); */
811 /* double max; */
812 /* int nBins = universeAsimovSpectrum[0]->GetNbinsX(); */
813 
814 /* TH1D *chiSqrPerBinByBin = new TH1D("chiSqrPerBinByBin", */
815 /* ";Logical Bin;#chi^{2}/ndf", */
816 /* 480, */
817 /* 0., */
818 /* 480.); */
819 
820 /* TH2D *poissonDrawsByBin = dynamic_cast<TH2D*>(tfRandom->Get("fake/PoissonDraws")); */
821 
822 /* for(int b = 0; b < nBins; ++b){ */
823 
824 /* if((b + 1 > 130 && b + 1 < 220) || */
825 /* (b + 1 > 387 && b + 1 < 477)){ */
826 /* max = 0.; */
827 
828 /* histBin = poissonDrawsByBin->ProjectionY(("poissonProjBin" + std::to_string(b+1)).c_str(), b+1, b+1); */
829 
830 /* TH1D binPoissonByBin("binPoissonByBin", */
831 /* ("Poisson Distribution For Bin " + std::to_string(b)+";Events in Bin;Universes").c_str(), */
832 /* 30, */
833 /* 0., */
834 /* 30.); */
835 
836 /* /\* //cout << b + 1 *\/ */
837 /* /\* << " " *\/ */
838 /* /\* << universeAsimovSpectrum[0]->GetBinContent(b + 1) *\/ */
839 /* /\* << " " *\/ */
840 /* /\* << 100. * TMath::Poisson(1, universeAsimovSpectrum[0]->GetBinContent(b + 1)) *\/ */
841 /* /\* << endl; *\/ */
842 
843 /* for(int i = 0; i < histBin->GetNbinsX() + 1; ++i){ */
844 /* binPoissonByBin.Fill(histBin->GetBinContent(i + 1)); */
845 /* } */
846 
847 /* double scale = binPoissonByBin.Integral(); */
848 /* double chiSqrByHist = 0.; */
849 /* double chiSqrByBin = 0.; */
850 /* double byHistContent; */
851 /* double byBinContent; */
852 
853 /* for(size_t i = 0; i < counts.size(); ++i){ */
854 /* counts[i] = i * 0.1; */
855 /* prob[i] = scale * TMath::Poisson(counts[i], universeAsimovSpectrum[0]->GetBinContent(b + 1)); */
856 /* max = std::max(max, prob[i]); */
857 
858 /* if(i%10 == 0){ */
859 /* byBinContent = binPoissonByBin.GetBinContent(binPoissonByBin.FindBin(counts[i] + 0.5)); */
860 /* chiSqrByBin += std::pow(prob[i] - byBinContent, 2.) / 29.; */
861 
862 /* ////cout << counts[i] << " " << prob[i] << " " << binPoissonByHist.FindBin(counts[i]) << " " << byHistContent << " " << chiSqrByHist << " " << byBinContent << " " << chiSqrByBin << endl; */
863 /* } */
864 /* } */
865 
866 
867 /* chiSqrPerBinByBin->Fill(b, chiSqrByBin); */
868 
869 /* ////cout << "bin " << b + 1 << " " << binPoisson.Integral() << endl; */
870 
871 /* max = std::max(max, binPoissonByBin.GetMaximum()); */
872 
873 /* binPoissonByBin.SetMaximum(1.2 * max); */
874 /* binPoissonByBin.Draw("hist"); */
875 
876 /* TGraph gr(counts.size(), counts.data(), prob.data()); */
877 /* gr.SetLineColor(kRed); */
878 /* gr.Draw("lsame"); */
879 
880 /* TLegend leg(0.5, 0.5, 0.9, 0.9); */
881 /* leg.SetBorderSize(0); */
882 /* leg.SetFillStyle(0); */
883 /* leg.AddEntry(&binPoissonByBin, ("Poisson Draw Bin By Bin #chi^{2} = " + std::to_string(chiSqrByBin) ).c_str(), "l"); */
884 /* leg.Draw(); */
885 
886 /* if (b == 130) poissonCanv->Print("poissonDistributions.pdf(", "pdf"); */
887 /* else poissonCanv->Print("poissonDistributions.pdf" , "pdf"); */
888 /* } */
889 /* } */
890 
891 
892 /* chiSqrPerBinByBin->GetXaxis()->SetRangeUser(130, 220); */
893 /* chiSqrPerBinByBin->Draw("hist"); */
894 
895 /* TLegend leg(0.5, 0.5, 0.9, 0.9); */
896 /* leg.SetBorderSize(0); */
897 /* leg.SetFillStyle(0); */
898 /* leg.AddEntry(chiSqrPerBinByBin, "Poisson Draw Bin By Bin" , "l"); */
899 /* leg.Draw(); */
900 /* poissonCanv->Print("poissonDistributions.pdf", "pdf"); */
901 
902 /* chiSqrPerBinByBin->GetXaxis()->SetRangeUser(386, 480); */
903 /* chiSqrPerBinByBin->Draw("hist"); */
904 
905 /* leg.Draw(); */
906 /* poissonCanv->Print("poissonDistributions.pdf)", "pdf"); */
907 
908 /* } */
909 
910 //------------------------------------------------------------------------------
911 void MakeBiasPlots(TGraph* bestFitGr,
912  std::vector<double> const& inputPoint,
913  std::string const& legend,
914  std::string const& pdfFileName)
915 {
916 
917  // loop over the points in the graph and fill histograms of the difference between
918  // the point and the input x and y values
919 
920  TH1D* xDiffHist = new TH1D((legend + "xDiffHist").c_str(),
921  (legend + ";#Delta(sin^{2}#theta_{23});Universes").c_str(),
922  20,
923  -0.2,
924  0.2);
925 
926  TH1D* yDiffHist = new TH1D((legend + "yDiffHist").c_str(),
927  (legend + ";#Delta(#Delta m^{2}_{32}) #times 10^{-3}eV^{2};Universes").c_str(),
928  15,
929  -0.3,
930  0.3);
931 
932  TH2D* diffHist = new TH2D((legend + "pointDiffHist").c_str(),
933  (legend + ";#Delta(sin^{2}#theta_{23});#Delta(#Delta m^{2}_{32}) #times 10^{-3}eV^{2}").c_str(),
934  20,
935  -0.2,
936  0.2,
937  15,
938  -0.3,
939  0.3);
940 
941  double x, xInput, y;
942  for(int i = 0; i < bestFitGr->GetN(); ++i){
943  bestFitGr->GetPoint(i, x, y);
944  xInput = (x < 0.5) ? 0.5 - (inputPoint[0] - 0.5) : inputPoint[0];
945  xDiffHist->Fill(x - xInput);
946  yDiffHist->Fill(y - inputPoint[1]);
947  diffHist ->Fill(x - xInput,
948  y - inputPoint[1]);
949  }
950 
951  TCanvas *canv = new TCanvas((legend + "diffsCanv").c_str());
952  canv->Divide(2, 2);
953 
954  canv->cd(1)->SetGridy();
955  canv->cd(1)->SetGridx();
956  xDiffHist->Draw("hist");
957 
958  canv->cd(2)->SetGridy();
959  canv->cd(2)->SetGridx();
960  yDiffHist->Draw("hist");
961 
962  canv->cd(3)->SetGridy();
963  canv->cd(3)->SetGridx();
964  diffHist->Draw("colz");
965 
966  canv->Print(pdfFileName.c_str(), "pdf");
967 }
968 
969 //------------------------------------------------------------------------------
971 {
972  // set some variables we need
973  std::string randomUType1FileName = "cmf_random_universes_no_syst_poisson_fluctuations_1x_exposure_max_mix_hist.root";
974  std::string randomUType2FileName = "cmf_random_universes_no_syst_poisson_fluctuations_1x_exposure_max_mix_hist.root";
975  std::string contourType1FileName = "cmf_3flavor_test_cont_no_systs_fluctuations_1x_exposure_statonly_covmat_max_mix_hist.root";
976  std::string contourType2FileName = "cmf_3flavor_test_cont_no_systs_fluctuations_1x_exposure_statonly_poisson_max_mix_hist.root";
977 
978  std::string type1Legend("CovMat Stats Only");
979  std::string type2Legend("Poisson LL Stats Only");
980 
981  auto *tfType1Contour = new TFile(contourType1FileName.c_str(), "READ");
982  auto *tfType2Contour = new TFile(contourType2FileName.c_str(), "READ");
983  auto *tfType1Random = new TFile(randomUType1FileName.c_str(), "READ");
984  auto *tfType2Random = new TFile(randomUType2FileName.c_str(), "READ");
985 
986  ////cout << tfContour << " " << tfType1Random << endl;
987 
988  // get a histogram that defines the axis ranges
989  auto *limitsHist = dynamic_cast<TH2F*>(tfType1Contour->Get("contours/Asimov/AsimovTh23Dmsq32CLContours"));
990 
991  //cout << "got limits hist " << limitsHist << endl;
992 
993  std::vector<TGraph*> asimov1sigmaGrs;
994  std::vector<TGraph*> asimov2sigmaGrs;
995  std::vector<TGraph*> asimov3sigmaGrs;
996 
997  FindUniverseContours(asimov1sigmaGrs,
998  asimov2sigmaGrs,
999  asimov3sigmaGrs,
1000  dynamic_cast<TDirectory*>(tfType1Contour->Get("contours/Asimov")));
1001 
1002  std::vector<TGraph*> asimovType21sigmaGrs;
1003  std::vector<TGraph*> asimovType22sigmaGrs;
1004  std::vector<TGraph*> asimovType23sigmaGrs;
1005 
1006  FindUniverseContours(asimovType21sigmaGrs,
1007  asimovType22sigmaGrs,
1008  asimovType23sigmaGrs,
1009  dynamic_cast<TDirectory*>(tfType2Contour->Get("contours/Asimov")));
1010 
1011 
1012  //cout << "got asimov contours" << endl;
1013 
1014  std::vector<TGraph*> median1sigmaGrs;
1015  std::vector<TGraph*> median2sigmaGrs;
1016  std::vector<TGraph*> median3sigmaGrs;
1017 
1018  FindUniverseContours(median1sigmaGrs,
1019  median2sigmaGrs,
1020  median3sigmaGrs,
1021  dynamic_cast<TDirectory*>(tfType1Contour->Get("contours/Median")));
1022 
1023  std::vector<TGraph*> medianType21sigmaGrs;
1024  std::vector<TGraph*> medianType22sigmaGrs;
1025  std::vector<TGraph*> medianType23sigmaGrs;
1026 
1027  FindUniverseContours(medianType21sigmaGrs,
1028  medianType22sigmaGrs,
1029  medianType23sigmaGrs,
1030  dynamic_cast<TDirectory*>(tfType2Contour->Get("contours/Median")));
1031 
1032  //cout << "got median contours" << endl;
1033 
1034  // store the collection of contours for each universe
1035  std::map<size_t, std::vector<TGraph*>> universe1sigmaGrs;
1036  std::map<size_t, std::vector<TGraph*>> universe2sigmaGrs;
1037  std::map<size_t, std::vector<TGraph*>> universe3sigmaGrs;
1038  std::map<size_t, TH2F*> universeDeltaChiSqr;
1039  std::map<size_t, TH1D*> universeAsimovSpectrum;
1040  std::map<size_t, TH1D*> universeFHCSpectrum;
1041  std::map<size_t, TH1D*> universeRHCSpectrum;
1042 
1043  FillUniverseContainers(tfType1Contour,
1044  tfType1Random,
1045  numUniverses,
1046  universe1sigmaGrs,
1047  universe2sigmaGrs,
1048  universe3sigmaGrs,
1049  universeDeltaChiSqr,
1050  universeAsimovSpectrum,
1051  universeFHCSpectrum,
1052  universeRHCSpectrum);
1053 
1054  //cout << "filled type1 containers" << endl;
1055 
1056  std::map<size_t, std::vector<TGraph*>> universeType21sigmaGrs;
1057  std::map<size_t, std::vector<TGraph*>> universeType22sigmaGrs;
1058  std::map<size_t, std::vector<TGraph*>> universeType23sigmaGrs;
1059  std::map<size_t, TH2F*> universeType2DeltaChiSqr;
1060  std::map<size_t, TH1D*> universeType2AsimovSpectrum;
1061  std::map<size_t, TH1D*> universeType2FHCSpectrum;
1062  std::map<size_t, TH1D*> universeType2RHCSpectrum;
1063 
1064  FillUniverseContainers(tfType2Contour,
1065  tfType1Random,
1066  numUniverses,
1067  universeType21sigmaGrs,
1068  universeType22sigmaGrs,
1069  universeType23sigmaGrs,
1070  universeType2DeltaChiSqr,
1071  universeType2AsimovSpectrum,
1072  universeType2FHCSpectrum,
1073  universeType2RHCSpectrum);
1074 
1075  //cout << "filled type2 containers" << endl;
1076 
1077  // get the median canvas and plot the best fits from the universes on it
1078  auto *bestFitType1Gr = dynamic_cast<TGraph* >(tfType1Contour->Get("contours/UniverseBestFitPoints"));
1079  auto *bestFitType2Gr = dynamic_cast<TGraph* >(tfType2Contour->Get("contours/UniverseBestFitPoints"));
1080 
1081  //cout << "best fit graph " << bestFitGr << endl;
1082 
1083  DrawBestFitsOverContours(limitsHist,
1084  bestFitType1Gr,
1085  bestFitType2Gr,
1086  "universes.pdf(",
1087  type1Legend,
1088  type2Legend,
1089  asimov1sigmaGrs,
1090  asimov2sigmaGrs,
1091  asimov3sigmaGrs,
1092  median1sigmaGrs,
1093  median2sigmaGrs,
1094  median3sigmaGrs,
1095  asimovType21sigmaGrs,
1096  asimovType22sigmaGrs,
1097  asimovType23sigmaGrs,
1098  medianType21sigmaGrs,
1099  medianType22sigmaGrs,
1100  medianType23sigmaGrs);
1101 
1102  //cout << "best fits drawn over graphs" << endl;
1103 
1104  MakeBiasPlots(bestFitType1Gr,
1105  std::vector<double>({inputX, inputY}),
1106  type1Legend,
1107  "universes.pdf");
1108 
1109  MakeBiasPlots(bestFitType2Gr,
1110  std::vector<double>({inputX, inputY}),
1111  type2Legend,
1112  "universes.pdf");
1113 
1114  //cout << "best fits drawn" << endl;
1115 
1116  /* // draw the contour ranges around the median */
1117  /* medianCanv->Draw(); */
1118  /* medianCanv->cd(); */
1119 
1120  /* asimov1sigma1->Draw("lsame"); */
1121  /* asimov1sigma2->Draw("lsame"); */
1122  /* asimov2sigma->Draw("lsame"); */
1123  /* asimov3sigma->Draw("lsame"); */
1124  /* bestFitGr->Draw("psame"); */
1125  /* medianCanv->Print("universes.pdf", "pdf"); */
1126 
1127  limitsHist->SetTitle((type1Legend + " 1#sigma Contours").c_str());
1128  DrawUniverseNSigmaContours(limitsHist,
1129  median1sigmaGrs,
1130  median2sigmaGrs,
1131  median3sigmaGrs,
1132  "universes.pdf",
1133  universe1sigmaGrs);
1134 
1135  //cout << "1sigma contours drawn" << endl;
1136 
1137  limitsHist->SetTitle((type1Legend + " 2#sigma Contours").c_str());
1138  DrawUniverseNSigmaContours(limitsHist,
1139  median1sigmaGrs,
1140  median2sigmaGrs,
1141  median3sigmaGrs,
1142  "universes.pdf",
1143  universe2sigmaGrs);
1144 
1145  //cout << "2sigma contours drawn" << endl;
1146 
1147  limitsHist->SetTitle((type1Legend + " 3#sigma Contours").c_str());
1148  DrawUniverseNSigmaContours(limitsHist,
1149  median1sigmaGrs,
1150  median2sigmaGrs,
1151  median3sigmaGrs,
1152  "universes.pdf",
1153  universe3sigmaGrs);
1154 
1155  //cout << "3sigma contours drawn" << endl;
1156 
1157  DrawMedianContours(universeDeltaChiSqr,
1158  "universes.pdf",
1159  median1sigmaGrs,
1160  median2sigmaGrs,
1161  median3sigmaGrs,
1162  type1Legend);
1163 
1164  //cout << "type 1 median contours drawn" << endl;
1165 
1166  DrawMedianContours(universeType2DeltaChiSqr,
1167  "universes.pdf",
1168  medianType21sigmaGrs,
1169  medianType22sigmaGrs,
1170  medianType23sigmaGrs,
1171  type2Legend);
1172 
1173  //cout << "type 2 median contours drawn" << endl;
1174 
1176  "universes.pdf",
1177  type1Legend,
1178  type2Legend,
1179  asimov3sigmaGrs,
1180  universe1sigmaGrs,
1181  universe2sigmaGrs,
1182  universe3sigmaGrs,
1183  universeType21sigmaGrs,
1184  universeType22sigmaGrs,
1185  universeType23sigmaGrs,
1186  universeDeltaChiSqr,
1187  universeAsimovSpectrum,
1188  universeFHCSpectrum,
1189  universeRHCSpectrum);
1190 
1191  //cout << "universes drawn" << endl;
1192 
1193 
1194  /* // get the asimov canvas */
1195  /* auto *asimovCanv = dynamic_cast<TCanvas*>(tf->Get((topLevelDirName + asimovDirName + asimovCanvNam).c_str())); */
1196 
1197  /* asimovCanv->Draw(); */
1198  /* asimovCanv->cd(); */
1199 
1200  /* // loop over the contours for the universes and draw those on the same canv */
1201  /* int color = kOrange; */
1202  /* for(auto const& uitr : universe3sigmaGrs){ */
1203  /* for(auto const& gritr : uitr.second){ */
1204  /* gritr->SetLineWidth(4); */
1205  /* gritr->SetLineColor(color); */
1206  /* gritr->Draw("lsame"); */
1207  /* } */
1208  /* if (uitr.first > 9 && uitr.first < 19) color = kGreen; */
1209  /* else if(uitr.first > 19 && uitr.first < 29) color = kCyan; */
1210  /* else if(uitr.first > 29 && uitr.first < 39) color = kViolet; */
1211  /* if(uitr.first > 39) break; */
1212  /* } */
1213 
1214  /* asimovCanv->Update(); */
1215 
1216 } // fin
TFile * inputFile
Definition: PlotXSec.C:134
enum BeamMode kRed
void DrawBestFitsOverContours(TH2F *limitsHist, TGraph *bestFitType1Gr, TGraph *bestFitType2Gr, std::string const &pdfFileName, std::string const &type1Legend, std::string const &type2Legend, std::vector< TGraph * > const &asimovType11sigmaGrs, std::vector< TGraph * > const &asimovType12sigmaGrs, std::vector< TGraph * > const &asimovType13sigmaGrs, std::vector< TGraph * > const &medianType11sigmaGrs, std::vector< TGraph * > const &medianType12sigmaGrs, std::vector< TGraph * > const &medianType13sigmaGrs, std::vector< TGraph * > const &asimovType21sigmaGrs, std::vector< TGraph * > const &asimovType22sigmaGrs, std::vector< TGraph * > const &asimovType23sigmaGrs, std::vector< TGraph * > const &medianType21sigmaGrs, std::vector< TGraph * > const &medianType22sigmaGrs, std::vector< TGraph * > const &medianType23sigmaGrs)
void FillUniverseContainers(TFile *contourFile, TFile *randomFile, size_t numUniverses, std::map< size_t, std::vector< TGraph * >> &universe1sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universe2sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universe3sigmaGrs, std::map< size_t, TH2F * > &universeDeltaChiSqr, std::map< size_t, TH1D * > &universeAsimovSpectrum, std::map< size_t, TH1D * > &universeFHCSpectrum, std::map< size_t, TH1D * > &universeRHCSpectrum)
TCanvas * canv
string directory
projection from multiple chains
TLegend * leg1
Definition: plot_hist.C:105
void DrawUniverses(size_t numUniverses, std::string const &pdfFileName, std::string const &type1Legend, std::string const &type2Legend, std::vector< TGraph * > &asimov3sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universe1sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universe2sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universe3sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universeType21sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universeType22sigmaGrs, std::map< size_t, std::vector< TGraph * >> &universeType23sigmaGrs, std::map< size_t, TH2F * > &universeDeltaChiSqr, std::map< size_t, TH1D * > &universeAsimovSpectrum, std::map< size_t, TH1D * > &universeFHCSpectrum, std::map< size_t, TH1D * > &universeRHCSpectrum)
void FindUniverseSpectra(TFile *inputFile, size_t u, std::map< size_t, TH1D * > &universeAsimovSpectrum, std::map< size_t, TH1D * > &universeFHCSpectrum, std::map< size_t, TH1D * > &universeRHCSpectrum)
double inputX
void compareAsimovUniverseCL()
void MakeBiasPlots(TGraph *bestFitGr, std::vector< double > const &inputPoint, std::string const &legend, std::string const &pdfFileName)
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
void FindUniverseDeltaChiSqr(size_t u, std::map< size_t, TH2F * > &deltaChiSqrMap, TDirectory *directory, std::string const &keyName)
void DrawUniverseNSigmaContours(TH2F *limitsHist, std::vector< TGraph * > const &median1sigmaGrs, std::vector< TGraph * > const &median2sigmaGrs, std::vector< TGraph * > const &median3sigmaGrs, std::string const &pdfFileName, std::map< size_t, std::vector< TGraph * >> const &universeGrs)
double inputY
void FindMedianCLContours(TH2F *clHeatMap, std::vector< std::vector< TGraph * >> &contourGraphs, int lineColor)
size_t numUniverses
void FindUniverseContours(std::vector< TGraph * > &oneSigmaGrs, std::vector< TGraph * > &twoSigmaGrs, std::vector< TGraph * > &threeSigmaGrs, TDirectory *directory)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void DrawMedianContours(std::map< size_t, TH2F * > &universeDeltaChiSqr, std::string const &pdfFileName, std::vector< TGraph * > const &median1sigmaGrs, std::vector< TGraph * > const &median2sigmaGrs, std::vector< TGraph * > const &median3sigmaGrs, std::string const &legendBase)
enum BeamMode kBlue
void next()
Definition: show_event.C:84
enum BeamMode string