makeContourPlotsOld.C
Go to the documentation of this file.
1 // script to make some contour plots for CMF sterile analysis.
2 //
3 // run : root -l -b 'makeContourPlots.C("/path/to/input/file", {1,2,3}, "xpar", "ypar")'
4 // where:
5 // -- arg1: a histogram file from running a cmf_combineresultsjob.fcl job
6 // -- arg2: a vector of the contour levels to run, 1 sigma, 2 sigma, etc.
7 // -- arg3/4: parameter names - must match those in the file _exactly_
8 //
9 // author: Adam Lister
10 // email : alister1@fnal.gov
11 
12 // global scope variables
13 TFile* thisFile = nullptr;
14 int nContours = 0;
15 bool is1Sigma = false;
16 bool is2Sigma = false;
17 bool is3Sigma = false;
20 bool logx;
21 bool logy;
22 
23 // enum for confidence limits
24 enum CLNumber
25 {
29 };
30 
31 // function for plotting a vector of confidence limits
32 // isMono and isMedian are for drawing comparisons of asimov and median
33 // contours
34 //.............................................................................
35 void PlotCLs(std::vector<TGraph*> CLs, CLNumber cln,
36  bool isMono = false, bool isMedian = false, int ColorOverride = -1){
37 
38  int lcolor;
39  int lwidth;
40  int lstyle;
41 
42  if (cln == CLNumber::kOneSigma){
43  lcolor = kWhite;
44  lwidth = 2;
45  lstyle = 1;
46  }
47  if (cln == CLNumber::kTwoSigma){
48  lcolor = kWhite;
49  lwidth = 2;
50  lstyle = 9;
51  }
52  if (cln == CLNumber::kThreeSigma){
53  lcolor = kWhite;
54  lwidth = 2;
55  lstyle = 7;
56  }
57 
58  if (isMono){
59  if (isMedian) lcolor = TColor::GetColor(152, 202, 225);
60  else lcolor = TColor::GetColor(74, 123, 183);
61  }
62 
63  if (ColorOverride != -1)
64  lcolor = ColorOverride;
65 
66  for (auto &graph : CLs){
67  graph->SetLineColor(lcolor);
68  graph->SetLineWidth(lwidth);
69  graph->SetLineStyle(lstyle);
70  graph->Draw("same");
71  }
72 
73 }
74 
75 // make and fill median contours
76 //.............................................................................
77 TH2F* makeAndFillMedianDeltaChiSqrHist(std::vector<TH2F*> const& universeDeltaChiSqrHists)
78 {
79  // make the histogram we will return
80  TH2F* medianDeltaChiSqr = dynamic_cast<TH2F*>(universeDeltaChiSqrHists.front()->Clone("medianDeltaChiSqr"));
81  medianDeltaChiSqr->Reset();
82 
83  // loop over all the bins in the histograms and find the median for each bin
84  std::vector<float> binDeltaChiSqrDist;
85  for(int x = 1; x < universeDeltaChiSqrHists.front()->GetNbinsX() + 1; ++x){
86  for(int y = 1; y < universeDeltaChiSqrHists.front()->GetNbinsY() + 1; ++y){
87  binDeltaChiSqrDist.clear();
88  for(auto const* itr : universeDeltaChiSqrHists){
89  binDeltaChiSqrDist.emplace_back(itr->GetBinContent(x, y));
90  }
91 
92  // sort the values
93  std::sort(binDeltaChiSqrDist.begin(), binDeltaChiSqrDist.end());
94  medianDeltaChiSqr->SetBinContent(x, y, binDeltaChiSqrDist[binDeltaChiSqrDist.size() / 2]);
95  } // end loop over y bins
96  }// end loop over x bins
97 
98  return medianDeltaChiSqr;
99 
100 } // end makeAndFillMedianDeltaChiSqrHist
101 
102 // make a contour from a dela chisqr surface (mostly for median universes)
103 //.............................................................................
104 void getContoursFromDeltaChiSqr(TH2F* medianDeltaChiSqr,
105  std::vector<std::vector<TGraph*> > & contourGraphs)
106 {
107  std::vector<double> contourLevel({2.3, 6.18, 11.83});
108 
109  medianDeltaChiSqr->SetContour(contourLevel.size(), contourLevel.data());
110 
111  // This avoids disturbing canvases already
112  // Yes, we do need a canvas for this to work
113  // yes, we do need to call Update() on it
114  // No, ROOT is not very friendly.
115  TCanvas temp_can;
116  temp_can.cd();
117  medianDeltaChiSqr->Draw("contlist");
118  temp_can.Update();
119 
120  auto *plah = dynamic_cast<TObjArray *>(gROOT->GetListOfSpecials()->FindObject("contours"));
121  if(!plah) return;
122 
123  //cout << "TObjArray size: " << plah->GetSize() << endl;
124  contourGraphs.resize(contourLevel.size());
125 
126  for(int i = 0; i < plah->GetSize(); ++i){
127  std::vector<TGraph*> temp;
128  auto *list = dynamic_cast<TList*>(plah->At(i));
129  if(!list) break;
130  if(!(list->First())) break;
131  for(int igraph = 0; igraph < list->GetSize(); ++igraph){
132  if(list->At(igraph)) temp.emplace_back(dynamic_cast<TGraph*>(list->At(igraph)));
133  //cout << i << " " << igraph << " " << temp.back() << " " << list->At(igraph) << endl;
134  //temp.back()->SetLineColor(lineColor);
135  temp.back()->SetLineStyle(9);
136  temp.back()->SetLineWidth(3);
137  if (i == 0) temp.back()->SetLineColor(kBlue);
138  else if(i == 1) temp.back()->SetLineColor(kRed);
139  else if(i == 2) temp.back()->SetLineColor(kMagenta);
140  }
141  contourGraphs[i] = std::move(temp);
142  }
143 
144  gROOT->GetListOfSpecials()->FindObject("contours")->Clear();
145 }
146 
147 
148 // draw delta chi^2 map with contours on top for a given contour type
149 // allowed contour types are:
150 // -- Asimov
151 // -- Median
152 // -- Universe_<n>_
153 //.............................................................................
155 
156  std::string dir = "contours/" + contourtype;
157  TDirectoryFile* ContourDir = (TDirectoryFile*)thisFile->Get(dir.c_str());
158 
159  // make sure there's stuff in the directory
160  TList* lst = ContourDir->GetListOfKeys();
161  if (!lst) {
162  std::cout << "No keys found in file\n";
163  exit(1);
164  }
165 
166  // interested in getting:
167  // -- TH2 of deltachi^2
168  // -- all 2D contours
169  // -- 1D contours
170 
171  TH2F* deltaChiSqr = nullptr;
172  TGraph2D* bestFit = nullptr;
173  std::vector<TGraph*> CL1s;
174  std::vector<TGraph*> CL2s;
175  std::vector<TGraph*> CL3s;
176 
177  // loop everything in the directory
178  TObject* obj;
179  for (int i = 0; i < lst->GetSize(); ++i){
180  TKey* key = (TKey*)lst->At(i);
181  obj = key->ReadObj();
182 
183  // ignore anything that's not a TGraph or a TH2
184  if (!obj->InheritsFrom("TH2") &&
185  !obj->InheritsFrom("TGraph") &&
186  !obj->InheritsFrom("TGraph2D")){
187  continue;
188  }
189 
190  TString name = obj->GetName();
191 
192  if (obj->InheritsFrom("TH2")){
193  if (name.Contains("Prob") ||
194  (name.Contains("CLContours")) ||
195  (!name.Contains("DeltaChiSqr"))) continue;
196  std::cout << obj->GetName() << std::endl;
197  deltaChiSqr = (TH2F*)ContourDir->Get(obj->GetName());
198  }
199 
200  if (obj->InheritsFrom("TGraph")){
201  if (name.Contains("CL1")) CL1s.push_back((TGraph*)ContourDir->Get(obj->GetName()));
202  if (name.Contains("CL2")) CL2s.push_back((TGraph*)ContourDir->Get(obj->GetName()));
203  if (name.Contains("CL3")) CL3s.push_back((TGraph*)ContourDir->Get(obj->GetName()));
204  }
205 
206  if (obj->InheritsFrom("TGraph2D") && contourtype != "Median"){
207  bestFit = (TGraph2D*)ContourDir->Get(obj->GetName());
208  }
209 
210  }
211 
212  TCanvas* c1 = new TCanvas("c1", "", 700, 500);
213  c1->SetLeftMargin(0.18);
214  c1->SetRightMargin(0.18);
215  c1->SetBottomMargin(0.18);
216  if (logx)
217  c1->SetLogx();
218  if (logy)
219  c1->SetLogy();
220  //c1->SetLogz();
221  c1->cd();
222  deltaChiSqr->GetXaxis()->SetTitleOffset(1.2);
223  deltaChiSqr->GetZaxis()->SetRangeUser(10e-2,20);
224  deltaChiSqr->SetContour(1000);
225  deltaChiSqr->Draw("colz");
229 
230  TGraph *bf = nullptr;
231  if (contourtype != "Median"){
232  Double_t x[1] = {bestFit->GetX()[0]};
233  Double_t y[1] = {bestFit->GetY()[0]};
234  bf = new TGraph(1, x, y);
235  bf->SetMarkerSize(1.5);
236  bf->SetMarkerStyle(29);
237  bf->SetMarkerColor(TColor::GetColor(221,61,45));
238  bf->Draw("same p");
239  }
240 
241  std::string label = contourtype;
242  if (label.find("Contours") != std::string::npos)
243  label.erase(label.find("Contours"), label.find("Contours"));
244  label.erase(remove(label.begin(), label.end(), '_'), label.end());
245 
246  TPaveText* pt = new TPaveText(0.42, 0.85, 0.85, 0.90, "NDC");
247  pt->SetFillStyle(0);
248  pt->SetLineWidth(0);
249  pt->SetTextColor(kWhite);
250  pt->AddText((label+" Sensitivity").c_str());
251  pt->Draw("same");
252 
253  float legendSize = nContours * 0.05;
254  TLegend* leg = new TLegend(0.6, 0.85-legendSize, 0.8, 0.85);
255  if (is1Sigma) leg->AddEntry(CL1s.at(0), "1#sigma", "l");
256  if (is2Sigma) leg->AddEntry(CL2s.at(0), "2#sigma", "l");
257  if (is3Sigma) leg->AddEntry(CL3s.at(0), "3#sigma", "l");
258  if (contourtype != "Median"){
259  leg->SetTextAlign(22);
260  leg->SetY1(leg->GetY1()-0.05);
261  leg->AddEntry(bf, "Best Fit", "p");
262  }
263  leg->SetTextColor(kWhite);
264  leg->SetLineWidth(0);
265  leg->SetFillStyle(0);
266  leg->Draw("same");
267 
268  gPad->Modified();
269  gPad->Update();
270  gPad->RedrawAxis();
271 
272  c1->SaveAs((std::string(deltaChiSqr->GetName())+".png").c_str());
273  c1->SaveAs((std::string(deltaChiSqr->GetName())+".pdf").c_str());
274 
275 }
276 
277 // compare asimov and median contours for a given confidence level
278 //.............................................................................
280 
281  std::string clnName;
282  std::string clnLatex;
283  if (cln == CLNumber::kOneSigma)
284  {
285  clnName = "CL1";
286  clnLatex = "1#sigma";
287  }
288  if (cln == CLNumber::kTwoSigma)
289  {
290  clnName = "CL2";
291  clnLatex = "2#sigma";
292  }
293  if (cln == CLNumber::kThreeSigma)
294  {
295  clnName = "CL3";
296  clnLatex = "3#sigma";
297  }
298 
299  std::vector<TGraph*> CLsAsimov;
300  std::vector<TGraph*> CLsMedian;
301 
302  // Asimov
303  std::string asimovDirName = "contours/Asimov";
304  TDirectoryFile* AsimovDir = (TDirectoryFile*)thisFile->Get(asimovDirName.c_str());
305  TList* lstAsimov = AsimovDir->GetListOfKeys();
306 
307  // loop everything in the directory
308  TObject* obj;
309  for (int i = 0; i < lstAsimov->GetSize(); ++i){
310  TKey* key = (TKey*)lstAsimov->At(i);
311  obj = key->ReadObj();
312 
313  // ignore anything that's not a TGraph or a TH2
314  if (!obj->InheritsFrom("TGraph")) continue;
315 
316  TString name = obj->GetName();
317 
318  if (obj->InheritsFrom("TGraph")){
319  if (name.Contains(clnName.c_str())) CLsAsimov.push_back((TGraph*)AsimovDir->Get(obj->GetName()));
320  }
321 
322  }
323 
324  // Median
325  std::string medianDirName = "contours/Median";
326  TDirectoryFile* MedianDir = (TDirectoryFile*)thisFile->Get(medianDirName.c_str());
327  TList* lstMedian = MedianDir->GetListOfKeys();
328 
329  // loop everything in the directory
330  for (int i = 0; i < lstMedian->GetSize(); ++i){
331  TKey* key = (TKey*)lstMedian->At(i);
332  obj = key->ReadObj();
333 
334  // ignore anything that's not a TGraph or a TH2
335  if (!obj->InheritsFrom("TGraph")) continue;
336 
337  TString name = obj->GetName();
338 
339  if (obj->InheritsFrom("TGraph")){
340  if (name.Contains(clnName.c_str())) CLsMedian.push_back((TGraph*)MedianDir->Get(obj->GetName()));
341  }
342 
343  }
344 
345 
346  TCanvas* c1 = new TCanvas("c1", "", 700, 500);
347  c1->SetLeftMargin(0.18);
348  c1->SetRightMargin(0.18);
349  c1->SetBottomMargin(0.18);
350  if (logx)
351  c1->SetLogx();
352  if (logy)
353  c1->SetLogy();
354  c1->SetLogz();
355  c1->cd();
356  // sets axis range
357  std::string axesName = "contours/Asimov/Asimov"+xaxis+yaxis+"CLContours";
358  TH2F* axes = (TH2F*)thisFile->Get(axesName.c_str());
359  axes->GetXaxis()->SetTitleOffset(1.2);
360  axes->Draw();
361  PlotCLs(CLsAsimov, cln, true, false);
362  PlotCLs(CLsMedian, cln, true, true);
363 
364  TLegend* leg = new TLegend(0.2, 0.2, 0.5, 0.3);
365  leg->AddEntry(CLsAsimov.at(0), ("Asimov "+clnLatex).c_str(), "l");
366  leg->AddEntry(CLsMedian.at(0), ("Median "+clnLatex).c_str(), "l");
367  leg->SetLineWidth(0);
368  leg->SetFillStyle(0);
369  leg->Draw("same");
370 
371  gPad->Modified();
372  gPad->Update();
373 
374  std::string filename = "AsimovVersusMedian_"+clnName;
375  c1->SaveAs((filename+".png").c_str());
376  c1->SaveAs((filename+".pdf").c_str());
377 
378 }
379 
381  std::string clnName;
382  std::string clnLatex;
383  if (cln == CLNumber::kOneSigma)
384  {
385  clnName = "CL1";
386  clnLatex = "1#sigma";
387  }
388  if (cln == CLNumber::kTwoSigma)
389  {
390  clnName = "CL2";
391  clnLatex = "2#sigma";
392  }
393  if (cln == CLNumber::kThreeSigma)
394  {
395  clnName = "CL3";
396  clnLatex = "3#sigma";
397  }
398 
399  std::vector<TGraph*> CLsMedian;
400  std::vector<TGraph*> CLsAsimov;
401  std::vector<TGraph*> CLsUniverses;
402 
403  // get best fit points straight up
404  //TGraph* bfps = (TGraph*)thisFile->Get("contours/UniverseBestFitPoints");
405  //bfps->SetMarkerSize(1.5);
406  //bfps->SetMarkerStyle(29);
407  //bfps->SetMarkerColor(TColor::GetColor(221,61,45));
408 
409  TDirectoryFile* contoursDir = (TDirectoryFile*)thisFile->Get("contours");
410  int nUniverses = contoursDir->GetNkeys()-2; // there's Asimov + Median
411  TGraph* bfps = new TGraph(nUniverses);
412  bfps->SetMarkerSize(1.5);
413  bfps->SetMarkerStyle(29);
414  bfps->SetMarkerColor(TColor::GetColor(221,61,45));
415 
416  // Asimov
417  std::string asimovDirName = "contours/Asimov";
418  TDirectoryFile* AsimovDir = (TDirectoryFile*)thisFile->Get(asimovDirName.c_str());
419  TList* lstAsimov = AsimovDir->GetListOfKeys();
420 
421  // loop everything in the directory
422  TObject* obj;
423  for (int i = 0; i < lstAsimov->GetSize(); ++i){
424  TKey* key = (TKey*)lstAsimov->At(i);
425  obj = key->ReadObj();
426 
427  // ignore anything that's not a TGraph or a TH2
428  if (!obj->InheritsFrom("TGraph")) continue;
429 
430  TString name = obj->GetName();
431 
432  if (obj->InheritsFrom("TGraph")){
433  if (name.Contains(clnName.c_str())) CLsAsimov.push_back((TGraph*)AsimovDir->Get(obj->GetName()));
434  }
435 
436  }
437 
438 
439  // Median contour
440  std::string medianDirName = "contours/Median";
441  TDirectoryFile* MedianDir = (TDirectoryFile*)thisFile->Get(medianDirName.c_str());
442  TList* lstMedian = MedianDir->GetListOfKeys();
443 
444  // loop everything in the directory
445  for (int i = 0; i < lstMedian->GetSize(); ++i){
446  TKey* key = (TKey*)lstMedian->At(i);
447  obj = key->ReadObj();
448 
449  TString name = obj->GetName();
450 
451  if (obj->InheritsFrom("TGraph") && name.Contains()){
452  if (name.Contains(clnName.c_str())) CLsMedian.push_back((TGraph*)MedianDir->Get(obj->GetName()));
453  }
454 
455  }
456 
457  // Universe contours
458  int minUniverse = 0;
459  float minUniverseVal = 1e10;
460  int maxUniverse = 0;
461  float maxUniverseVal = 0;
462  for (int i = 0; i < nUniverses; ++i){
463  std::string universeDirName = "contours/Universe_"+std::to_string(i);
464  TDirectoryFile* UniverseDir = (TDirectoryFile*)thisFile->Get(universeDirName.c_str());
465  TList* lstUniverse = UniverseDir->GetListOfKeys();
466 
467  // loop everything in the directory
468  for (int i = 0; i < lstUniverse->GetSize(); ++i){
469  TKey* key = (TKey*)lstUniverse->At(i);
470  obj = key->ReadObj();
471 
472  // ignore anything that's not a TGraph
473  if (!obj->InheritsFrom("TGraph")) continue;
474 
475  TString name = obj->GetName();
476 
477  if (obj->InheritsFrom("TGraph2D")){
478  if (name.Contains("FitResult")){
479  TGraph2D* bfp = (TGraph2D*)UniverseDir->Get(obj->GetName());
480  bfps->SetPoint(i, bfp->GetX()[0], bfp->GetY()[0]);
481  }
482  }
483  else if (obj->InheritsFrom("TGraph")){
484  if (name.Contains(clnName.c_str())) CLsUniverses.push_back((TGraph*)UniverseDir->Get(obj->GetName()));
485  double xmin = TMath::MinElement(((TGraph*)UniverseDir->Get(obj->GetName()))->GetN(),
486  ((TGraph*)UniverseDir->Get(obj->GetName()))->GetX());
487  if (xmin < minUniverseVal){
488  minUniverseVal = xmin;
489  minUniverse = i;
490  }
491  double xmax = TMath::MaxElement(((TGraph*)UniverseDir->Get(obj->GetName()))->GetN(),
492  ((TGraph*)UniverseDir->Get(obj->GetName()))->GetX());
493  if (xmax > maxUniverseVal){
494  maxUniverseVal = xmax;
495  maxUniverse = i;
496  }
497 
498  }
499  }
500  }
501 
502  std::cout
503  << "Outlying Universes: "
504  << "\n -- Minimum: " << minUniverse << " at " << minUniverseVal
505  << "\n -- Maximum: " << maxUniverse << " at " << maxUniverseVal
506  << std::endl;
507 
508  TCanvas* c1 = new TCanvas("c1", "", 700, 500);
509  c1->SetLeftMargin(0.18);
510  c1->SetRightMargin(0.18);
511  c1->SetBottomMargin(0.18);
512  if (logx)
513  c1->SetLogx();
514  if (logy)
515  c1->SetLogy();
516  c1->SetLogz();
517  c1->cd();
518  // sets axis range
519  std::string axesName = "contours/Asimov/Asimov"+xaxis+yaxis+"CLContours";
520  TH2F* axes = (TH2F*)thisFile->Get(axesName.c_str());
521  axes->GetXaxis()->SetTitleOffset(1.2);
522  axes->Draw();
523  for (auto& uni : CLsUniverses){
524 
525  uni->SetLineColor(kGray);
526  uni->SetLineWidth(2);
527  uni->Draw("same");
528 
529  }
530  PlotCLs(CLsMedian, cln, true, true);
531  PlotCLs(CLsAsimov, cln, true, false);
532 
533  TLegend* leg = new TLegend(0.2, 0.2, 0.5, 0.3);
534  leg->AddEntry(CLsUniverses.at(0), ("Universes "+clnLatex).c_str(), "l");
535  leg->AddEntry(CLsMedian.at(0) , ("Median "+clnLatex).c_str(), "l");
536  leg->AddEntry(CLsAsimov.at(0) , ("Asimov "+clnLatex).c_str(), "l");
537  leg->SetLineWidth(0);
538  leg->SetFillStyle(0);
539  leg->Draw("same");
540 
541  gPad->Modified();
542  gPad->Update();
543 
544  std::string filename = "MedianAndUniverses"+clnName;
545  c1->SaveAs((filename+".png").c_str());
546  c1->SaveAs((filename+".pdf").c_str());
547 
548  bfps->Draw("p same");
549  leg->Draw("same");
550  c1->SaveAs((filename+"_withBestFits.png").c_str());
551  c1->SaveAs((filename+"_withBestFits.pdf").c_str());
552 
553 }
554 
555 //............................................................................
557 
558  std::string dir = "contours/" + contourtype;
559  TDirectoryFile* ContourDir = (TDirectoryFile*)thisFile->Get(dir.c_str());
560 
561  std::string name = "Asimov"+var+"DeltaChiSqr";
562 
563  // get histogram to set axes
564  TH1D* axes = (TH1D*)ContourDir->Get((name+"Hist").c_str());
565  TGraph* gr = (TGraph*)ContourDir->Get(name.c_str());
566 
567  gr->SetLineWidth(2);
568  gr->SetLineColor(TColor::GetColor(74, 123, 183));
569 
570  TCanvas *c1 = new TCanvas("c1", "c1", 700, 500);
571  c1->cd();
572  if (var == "Th24" || var == "Dmsq41")
573  c1->SetLogx();
574 
575  axes->Draw();
576  axes->GetXaxis()->SetTitleOffset(1.15);
577  gr->Draw("same l");
578 
579  TLine* oneSig = new TLine(axes->GetBinLowEdge(1), 1, axes->GetBinLowEdge(axes->GetNbinsX()+1), 1);
580  TLine* twoSig = new TLine(axes->GetBinLowEdge(1), 4, axes->GetBinLowEdge(axes->GetNbinsX()+1), 4);
581  TLine* threeSig = new TLine(axes->GetBinLowEdge(1), 9, axes->GetBinLowEdge(axes->GetNbinsX()+1), 9);
582 
583  oneSig->SetLineWidth(2);
584  oneSig->SetLineStyle(2);
585  oneSig->SetLineColor(kGray);
586  twoSig->SetLineWidth(2);
587  twoSig->SetLineStyle(2);
588  twoSig->SetLineColor(kGray);
589  threeSig->SetLineWidth(2);
590  threeSig->SetLineStyle(2);
591  threeSig->SetLineColor(kGray);
592 
593  oneSig->Draw("same");
594  twoSig->Draw("same");
595  threeSig->Draw("same");
596 
597  TLatex* ltOneSig = new TLatex();
598  TLatex* ltTwoSig = new TLatex();
599  TLatex* ltThreeSig = new TLatex();
600  ltOneSig->DrawLatexNDC(0.85, 0.2, "1#sigma");
601  ltTwoSig->DrawLatexNDC(0.85, 0.36, "2#sigma");
602  ltThreeSig->DrawLatexNDC(0.85, 0.62, "3#sigma");
603 
604  ltOneSig->SetTextColor(kGray);
605  ltTwoSig->SetTextColor(kGray);
606  ltThreeSig->SetTextColor(kGray);
607 
608  c1->SaveAs((name+".png").c_str());
609  c1->SaveAs((name+".pdf").c_str());
610 
611  //c1->SetLogy();
612  axes->GetYaxis()->SetRangeUser(0,10e-2);
613  axes->Draw("");
614  gr->Draw("l same");
615 
616  c1->SaveAs((name+"_zoom.png").c_str());
617  c1->SaveAs((name+"_zoom.pdf").c_str());
618 
619 }
620 
621 //............................................................................
622 void DrawHiddenParameter(std::string hiddenParName)
623 {
624 
625  gStyle->SetPalette(kCandy);
626 
627  TH2D* hidPar = (TH2D*)thisFile->Get(("contours/Asimov/Asimov"+hiddenParName).c_str());
628 
629  TCanvas* c1 = new TCanvas("c1", "", 700, 500);
630  c1->SetLeftMargin(0.18);
631  c1->SetRightMargin(0.18);
632  c1->SetBottomMargin(0.18);
633  if (logx)
634  c1->SetLogx();
635  if (logy)
636  c1->SetLogy();
637  c1->cd();
638 
639  float min = 10e10;
640  float max = 0;
641 
642  for (int i = 1; i < hidPar->GetNbinsX(); ++i){
643  for (int j = 1; j < hidPar->GetNbinsY(); ++j){
644  if (hidPar->GetBinContent(i,j) < min) min = hidPar->GetBinContent(i,j);
645  if (hidPar->GetBinContent(i,j) > max) max = hidPar->GetBinContent(i,j);
646  }
647  }
648  if (min == 0) min = -0.001;
649 
650  hidPar->GetZaxis()->SetRangeUser(min, max);
651  hidPar->SetContour(1000);
652 
653  hidPar->Draw("colz");
654 
655  c1->SaveAs(("HiddenParameterSpectra"+hiddenParName+".png").c_str());
656  c1->SaveAs(("HiddenParameterSpectra"+hiddenParName+".pdf").c_str());
657 
658 }
659 
660 // main function
661 //.............................................................................
662 void makeContourPlots(std::string file, std::vector<int> contours, std::string xlabel, std::string ylabel){
663 
664  thisFile = new TFile(file.c_str(), "read");
665  nContours = (int)contours.size();
666 
667  xaxis = xlabel;
668  yaxis = ylabel;
669 
670  logx = xaxis == "Th24" ? true : false;
671  logy = yaxis == "Dmsq41" ? true : false;
672 
673  std::cout << "Configured to make the following contours" << std::endl;
674  for(auto const& contour : contours){
675  if (contour < 11)
676  std::cout << " -- " << contour << " sigma" << std::endl;
677  else
678  std::cout << " -- " << contour << " percent" << std::endl;
679 
680  if (contour == 1)
681  is1Sigma = true;
682  if (contour == 2)
683  is2Sigma = true;
684  if (contour == 3)
685  is3Sigma = true;
686  }
687  std::cout << "Note that if these contours do not exist in the input file "
688  << "this macro will fail."
689  << std::endl;
690 
691  // custom colour palette using Paul Tols notes:
692  // https://personal.sron.nl/~pault/
693  const Int_t Number = 3;
694  Double_t Red[Number] = { 0.00, 74.0/255.0, 152.0/255.0};
695  Double_t Green[Number] = { 0.00, 123.0/255.0, 202.0/255.0 };
696  Double_t Blue[Number] = { 0.00, 183.0/255.0, 225.0/255.0 };
697  Double_t Length[Number] = { 0.00, 0.50, 1.0};
698  Int_t nb=1000;
699  TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
700  TColor::InvertPalette();
701 
702  DrawDeltaChiSqrWithContours("Asimov");
703  DrawDeltaChiSqrWithContours("Median");
704  DrawDeltaChiSqrWithContours("Universe_1");
711  Draw1DContour("Asimov", xaxis);
712  Draw1DContour("Asimov", yaxis);
713  //DrawHiddenParameter("Th34");
714  DrawHiddenParameter("Th13");
715 
716 
717 }
TH2F * makeAndFillMedianDeltaChiSqrHist(std::vector< TH2F * > const &universeDeltaChiSqrHists)
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
std::map< std::string, double > xmax
void makeContourPlots(std::string file, std::vector< int > contours, std::string xlabel, std::string ylabel)
void PlotCLs(std::vector< TGraph * > CLs, CLNumber cln, bool isMono=false, bool isMedian=false, int ColorOverride=-1)
std::string yaxis
string filename
Definition: shutoffs.py:106
int nContours
void Draw1DContour(std::string contourtype, std::string var)
std::string xaxis
bool logy
const char * label
void DrawMedianAndUniverses(CLNumber cln)
double lst
Definition: runWimpSim.h:113
bool is3Sigma
const double j
Definition: BetheBloch.cxx:29
bool is2Sigma
OStream cout
Definition: OStream.cxx:6
TFile * thisFile
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
bool logx
void getContoursFromDeltaChiSqr(TH2F *medianDeltaChiSqr, std::vector< std::vector< TGraph * > > &contourGraphs)
void DrawDeltaChiSqrWithContours(std::string contourtype)
void DrawHiddenParameter(std::string hiddenParName)
TDirectory * dir
Definition: macro.C:5
exit(0)
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24
void CompareAsimovAndMedian(CLNumber cln)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
static const double nb
Definition: Units.h:89
Float_t e
Definition: plot.C:35
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kBlue
cons_index_list< index_uni, nil_index_list > uni
bool is1Sigma
enum BeamMode string