Utilities.h
Go to the documentation of this file.
1 #include "TGraph.h"
2 #include "TStyle.h"
3 #include "TFile.h"
4 
5 enum CLNumber
6 {
13 };
14 
15 std::vector<std::string> CLName = {
16  "OneSigma",
17  "TwoSigma",
18  "ThreeSigma",
19  "NinetyPercent",
20  "NinetyFivePercent",
21  ""
22 };
23 
24 std::vector<std::string> CLNameLaTeX = {
25  "1#sigma",
26  "2#sigma",
27  "3#sigma",
28  "90\%",
29  "95\%",
30  ""
31 };
32 
33 enum CLType
34 {
39 };
40 
41 // global scope variables
42 //.............................................................................
43 
44 TFile* thisFile = nullptr;
45 std::vector<CLNumber> contourLevels;
48 bool logx;
49 bool logy;
51 
53  std::vector<int> contours){
54  thisFile = TFile::Open(f.c_str(), "read");
55  xaxis = xax;
56  yaxis = yax;
57 
58  logx = xaxis == "Th24" ? true : false;
59  logy = yaxis == "Dmsq41" ? true : false;
60 
61  for (auto &c : contours){
62  if (c == 1)
64  else if (c == 2)
66  else if (c == 3)
68  else if (c == 90)
70  else if (c == 95)
72  else {
74  "confidence level " +
75  std::to_string(c) +
76  "is not supported";
77  throw std::logic_error(errMsg);
78  }
79 
80  }
81 
82 }
83 
84 // some helpful structs and typedefs
85 //.............................................................................
86 
87 struct Contour {
88  Contour(std::vector<TGraph*> c,
89  CLNumber cn,
90  CLType ct)
91  : contours(c)
92  , cln(cn)
93  , clt(ct)
94  {
95  if (cln == CLNumber::kOneSigma){
96  clnName = "CL1";
97  clnLatex = "1#sigma";
98  }
99  if (cln == CLNumber::kTwoSigma){
100  clnName = "CL2";
101  clnLatex = "2#sigma";
102  }
103  if (cln == CLNumber::kThreeSigma){
104  clnName = "CL3";
105  clnLatex = "3#sigma";
106  }
108  clnName = "90Pct";
109  clnLatex = "90%";
110  }
112  clnName = "95Pct";
113  clnLatex = "95%";
114  }
115  if (clt == CLType::kAsimov){
116  cltName = "Asimov";
117  }
118  if (clt == CLType::kMedian){
119  cltName = "Median";
120  }
121  if (clt == CLType::kUniverse){
122  cltName = "Universe";
123  }
124  }
125 
126  std::vector<TGraph*> contours;
132 };
133 
134 typedef std::vector<Contour> ContourCollection;
135 typedef std::map<std::string, TObject*> LegendMap;
136 
137 //.............................................................................
139  // custom colour palette using Paul Tols notes:
140  // https://personal.sron.nl/~pault/
141  const Int_t Number = 3;
142  Double_t Red[Number] = { 0.00, 74.0/255.0, 152.0/255.0};
143  Double_t Green[Number] = { 0.00, 123.0/255.0, 202.0/255.0 };
144  Double_t Blue[Number] = { 0.00, 183.0/255.0, 225.0/255.0 };
145  Double_t Length[Number] = { 0.00, 0.50, 1.0};
146  Int_t nb=1000;
147  TColor::CreateGradientColorTable(Number,Length,Red,Green,Blue,nb);
148  TColor::InvertPalette();
149 
150  kPTRed = TColor::GetColor(221,61,45);
151  kPTOrange = TColor::GetColor(253,179,102);
152  kPTDarkBlue = TColor::GetColor(74, 123, 183);
153  kPTLightBlue = TColor::GetColor(152, 202, 225);
154 }
155 
156 // this is needed because we typically hadd the files together and the hists
157 // get multiplied by the number of files
158 //.............................................................................
159 int GetNFiles(){
160 
161  TDirectoryFile* dir = (TDirectoryFile*)thisFile->Get("contours/Asimov");
162  TList* lst = dir->GetListOfKeys();
163  int nFiles = 0;
164 
165  // loop everything in the directory
166  TObject* obj;
167  for (int i = 0; i < lst->GetSize(); ++i){
168  TKey* key = (TKey*)lst->At(i);
169  obj = key->ReadObj();
170 
171  TString name = obj->GetName();
172  if (name.Contains("AsimovContoursCanv"))
173  ++nFiles;
174  }
175 
176  return nFiles;
177 
178 }
179 //.............................................................................
182  int colorOverride = -1,
183  int styleOverride = -1)
184 {
185 
186  int lwidth = 2;
187  int lstyle;
188  int lcolor;
189 
190  for (auto& contourset : cc){
191  if (cln != kEmpty && cln != contourset.cln) continue;
192 
193  if (contourset.cln == CLNumber::kOneSigma) lstyle = 1;
194  if (contourset.cln == CLNumber::kTwoSigma) lstyle = 9;
195  if (contourset.cln == CLNumber::kThreeSigma) lstyle = 7;
196  if (contourset.cln == CLNumber::kNinetyPercent) lstyle = 11;
197  if (contourset.cln == CLNumber::kNinetyFivePercent) lstyle = 12;
198  if (styleOverride != -1) lstyle = styleOverride;
199 
200  if (contourset.clt == CLType::kAsimov) lcolor = kPTDarkBlue;
201  if (contourset.clt == CLType::kMedian) lcolor = kPTLightBlue;
202  if (contourset.clt == CLType::kUniverse) lcolor = kGray;
203  if (colorOverride != -1) lcolor = colorOverride;
204 
205  for (auto &gr : contourset.contours){
206  gr->SetLineColor(lcolor);
207  gr->SetLineWidth(lwidth);
208  gr->SetLineStyle(lstyle);
209  gr->Draw("same");
210  }
211 
212  }
213 
214 }
215 
216 //............................................................................
217 void Save1DContour(CLType clt, std::string var, int universeNum = -1){
218 
221  if (clt == CLType::kAsimov){
222  dirName = "contours/Asimov/";
223  histName = "Asimov"+var+"DeltaChiSqr";
224  }
225  else if (clt == CLType::kUniverse){
226  dirName = "contours/Universe_" + std::to_string(universeNum)+ "/";
227  histName = "Universe_" + std::to_string(universeNum) + var + "DeltaChiSqr";
228  }
229 
230  std::string grName = dirName+histName;
231 
232  // get histogram to set axes
233  TH1D* axes = (TH1D*)thisFile->Get((grName+"Hist").c_str());
234  TGraph* gr = (TGraph*)thisFile->Get(grName.c_str());
235  axes->GetXaxis()->SetNdivisions(505);
236  gr->SetLineWidth(2);
237  gr->SetLineColor(kPTDarkBlue);
238 
239  TLine* oneSig = new TLine(axes->GetBinLowEdge(1), 1, axes->GetBinLowEdge(axes->GetNbinsX()+1), 1);
240  TLine* twoSig = new TLine(axes->GetBinLowEdge(1), 4, axes->GetBinLowEdge(axes->GetNbinsX()+1), 4);
241  TLine* threeSig = new TLine(axes->GetBinLowEdge(1), 9, axes->GetBinLowEdge(axes->GetNbinsX()+1), 9);
242 
243  oneSig -> SetLineWidth(2);
244  oneSig -> SetLineStyle(2);
245  oneSig -> SetLineColor(kGray);
246  twoSig -> SetLineWidth(2);
247  twoSig -> SetLineStyle(2);
248  twoSig -> SetLineColor(kGray);
249  threeSig -> SetLineWidth(2);
250  threeSig -> SetLineStyle(2);
251  threeSig -> SetLineColor(kGray);
252 
253  TCanvas* canv = new TCanvas("canv", "", 600, 500);
254  axes->Draw();
255  gr->Draw("same l");
256 
257  oneSig -> Draw("same");
258  twoSig -> Draw("same");
259  threeSig -> Draw("same");
260 
261  TLatex* ltOneSig = new TLatex();
262  TLatex* ltTwoSig = new TLatex();
263  TLatex* ltThreeSig = new TLatex();
264  ltOneSig -> DrawLatexNDC(0.85, 0.2, "1#sigma");
265  ltTwoSig -> DrawLatexNDC(0.85, 0.36, "2#sigma");
266  ltThreeSig -> DrawLatexNDC(0.85, 0.62, "3#sigma");
267 
268  ltOneSig -> SetTextColor(kGray);
269  ltTwoSig -> SetTextColor(kGray);
270  ltThreeSig -> SetTextColor(kGray);
271 
272  canv->SaveAs((histName+".png").c_str());
273  canv->SaveAs((histName+".pdf").c_str());
274 
275  //c1->SetLogy();
276  axes->GetYaxis()->SetRangeUser(0,10e-2);
277  axes->Draw("");
278  gr->Draw("l same");
279 
280  canv->SaveAs((histName+"_zoom.png").c_str());
281  canv->SaveAs((histName+"_zoom.pdf").c_str());
282 
283 }
284 
285 
286 //.............................................................................
287 void getContoursFromDeltaChiSqr(TH2D* medianDeltaChiSqr,
288  std::vector<std::vector<TGraph*> > & contourGraphs)
289 {
290  std::vector<double> contourLevel({2.3, 6.18, 11.83});
291 
292  medianDeltaChiSqr->SetContour(contourLevel.size(), contourLevel.data());
293 
294  // This avoids disturbing canvases already
295  // Yes, we do need a canvas for this to work
296  // yes, we do need to call Update() on it
297  // No, ROOT is not very friendly.
298  TCanvas temp_can;
299  temp_can.cd();
300  medianDeltaChiSqr->Draw("contlist");
301  temp_can.Update();
302 
303  auto *plah = dynamic_cast<TObjArray *>(gROOT->GetListOfSpecials()->FindObject("contours"));
304  if(!plah) return;
305 
306  //cout << "TObjArray size: " << plah->GetSize() << endl;
307  contourGraphs.resize(contourLevel.size());
308 
309  for(int i = 0; i < plah->GetSize(); ++i){
310  std::vector<TGraph*> temp;
311  auto *list = dynamic_cast<TList*>(plah->At(i));
312  if(!list) break;
313  if(!(list->First())) break;
314  for(int igraph = 0; igraph < list->GetSize(); ++igraph){
315  if(list->At(igraph)) temp.emplace_back(dynamic_cast<TGraph*>(list->At(igraph)));
316  //cout << i << " " << igraph << " " << temp.back() << " " << list->At(igraph) << endl;
317  //temp.back()->SetLineColor(lineColor);
318  temp.back()->SetLineStyle(9);
319  temp.back()->SetLineWidth(3);
320  if (i == 0) temp.back()->SetLineColor(kBlue);
321  else if(i == 1) temp.back()->SetLineColor(kRed);
322  else if(i == 2) temp.back()->SetLineColor(kMagenta);
323  }
324  contourGraphs[i] = std::move(temp);
325  }
326 
327  gROOT->GetListOfSpecials()->FindObject("contours")->Clear();
328 }
329 
330 //.............................................................................
331 TH2D* makeAndFillMedianDeltaChiSqrHist(std::vector<TH2D*> const& universeDeltaChiSqrHists)
332 {
333 
334  // make the histogram we will return
335  TH2D* medianDeltaChiSqr = (TH2D*)universeDeltaChiSqrHists[0]->Clone("medianDeltaChiSqr");
336 
337  // loop over all the bins in the histograms and find the median for each bin
338  std::vector<float> binDeltaChiSqrDist;
339  for(int x = 1; x < universeDeltaChiSqrHists.front()->GetNbinsX() + 1; ++x){
340  for(int y = 1; y < universeDeltaChiSqrHists.front()->GetNbinsY() + 1; ++y){
341  binDeltaChiSqrDist.clear();
342  for(auto const* itr : universeDeltaChiSqrHists){
343  binDeltaChiSqrDist.emplace_back(itr->GetBinContent(x, y));
344  }
345 
346  // sort the values
347  std::sort(binDeltaChiSqrDist.begin(), binDeltaChiSqrDist.end());
348  medianDeltaChiSqr->SetBinContent(x, y, binDeltaChiSqrDist[binDeltaChiSqrDist.size() / 2]);
349  } // end loop over y bins
350  }// end loop over x bins
351 
352  return medianDeltaChiSqr;
353 
354 }
355 
356 //.............................................................................
358 
359  int numUniverses = ((TDirectoryFile*)thisFile->Get("contours"))->GetNkeys()-2;
360  std::vector<TH2D*> universeDeltaChiSqrHists;
363  for(size_t u = 0; u < numUniverses; ++u){
364  dirName = "contours/Universe_" + std::to_string(u) + "/";
365  histName = dirName +
366  "Universe_" +
367  std::to_string(u) +
368  xaxis + yaxis + "DeltaChiSqr";
369 
370  universeDeltaChiSqrHists.push_back((TH2D*)thisFile->Get(histName.c_str()));
371 
372  }
373 
374  // now get the median Delta chi^2 histogram
375  TH2D* medianDeltaChiSqr = makeAndFillMedianDeltaChiSqrHist(universeDeltaChiSqrHists);
376 
377  return medianDeltaChiSqr;
378 }
379 
380 //.............................................................................
381 TH2D* GetDeltaChiSqr(CLType clt, int universeNum = -1){
382 
383  int numUniverses = ((TDirectoryFile*)thisFile->Get("contours"))->GetNkeys()-2;
384 
385  TH2D* DeltaChiSqr = nullptr;
387  if (clt == CLType::kMedian){
388  DeltaChiSqr = GetMedianDeltaChiSqr();
389  }
390  else {
391  if (clt == CLType::kAsimov || clt == CLType::kNoType)
392  dirName = "contours/Asimov/Asimov";
393  else if (clt == CLType::kUniverse)
394  dirName = "contours/Universe_" + std::to_string(universeNum) +
395  "Universe_" + std::to_string(universeNum);
396 
397  std::string histName = dirName +
398  xaxis + yaxis + "DeltaChiSqr";
399 
400  DeltaChiSqr = (TH2D*)thisFile->Get(histName.c_str());
401 
402  if (clt == kAsimov)
403  DeltaChiSqr->Scale(1./GetNFiles());
404  if (clt == kNoType){
405  TH2D* clone = (TH2D*)DeltaChiSqr->Clone("clone");
406  clone->Reset();
407  return clone;
408  }
409  }
410 
411  DeltaChiSqr->GetZaxis()->SetRangeUser(0, 20);
412  DeltaChiSqr->SetContour(1000);
413 
414  return DeltaChiSqr;
415 
416 }
417 
418 //.............................................................................
420  TH2D *medianDeltaChiSqr = GetMedianDeltaChiSqr();
421 
422  std::vector<std::vector<TGraph*> > medianContours;
423  getContoursFromDeltaChiSqr(medianDeltaChiSqr, medianContours);
424 
425  Contour contour1Sigma(medianContours[0], CLNumber::kOneSigma, CLType::kMedian);
426  Contour contour2Sigma(medianContours[1], CLNumber::kTwoSigma, CLType::kMedian);
427  Contour contour3Sigma(medianContours[2], CLNumber::kThreeSigma, CLType::kMedian);
428 
429  ContourCollection cl = {contour1Sigma, contour2Sigma, contour3Sigma};
430 
431  return cl;
432 }
433 
434 
435 //.............................................................................
436 TGraph* GetBestFitPoint(CLType clt, int universeNum = -1){
437 
439  if (clt == CLType::kAsimov)
440  dirName = "contours/Asimov/Asimov";
441  else if (clt == CLType::kUniverse)
442  dirName = "contours/Universe_" + std::to_string(universeNum)
443  + "/Universe_" + std::to_string(universeNum);
444 
445  std::string grName = dirName+"FitResult";
446 
447  TGraph2D* bfp = (TGraph2D*)thisFile->Get(grName.c_str());
448  TGraph* gr = new TGraph(1);
449  gr->SetPoint(0, bfp->GetX()[0], bfp->GetY()[0]);
450 
451  return gr;
452 
453 }
454 
455 //.............................................................................
457  int numUniverses = ((TDirectoryFile*)thisFile->Get("contours"))->GetNkeys()-2;
458  TGraph* bfps = new TGraph(numUniverses);
459 
460  for (int i = 0; i < numUniverses; ++i){
461  TGraph* bfp = GetBestFitPoint(CLType::kUniverse, i);
462  bfps->SetPoint(i, bfp->GetX()[0], bfp->GetY()[0]);
463  }
464 
465  return bfps;
466 }
467 
468 //.............................................................................
469 ContourCollection GetContours(CLType clt, int universeNum = -1){
470 
471  ContourCollection Cls;
473  if (clt == CLType::kAsimov)
474  dirName = "contours/Asimov";
475  else if (clt == CLType::kUniverse)
476  dirName = "contours/Universe_" + std::to_string(universeNum);
477  else if (clt == CLType::kMedian){
478  Cls = GetMedianContours();
479  return Cls;
480  }
481 
482  TDirectoryFile* dir = (TDirectoryFile*)thisFile->Get(dirName.c_str());
483  TList* lst = dir->GetListOfKeys();
484 
485  for (auto& cln : contourLevels){
486 
487  Contour thisContour(std::vector<TGraph*>(), cln, clt);
488 
489  // loop everything in the directory
490  TObject* obj;
491  for (int i = 0; i < lst->GetSize(); ++i){
492  TKey* key = (TKey*)lst->At(i);
493  obj = key->ReadObj();
494 
495  TString name = obj->GetName();
496  if (obj->InheritsFrom("TGraph")){
497  if (name.Contains(thisContour.clnName.c_str()))
498  thisContour.contours.push_back((TGraph*)dir->Get(name));
499  }
500  }
501  Cls.push_back(thisContour);
502  }
503 
504  return Cls;
505 
506 }
507 
508 //.............................................................................
509 std::vector<ContourCollection> GetUniverseContours(){
510 
511  int numUniverses = ((TDirectoryFile*)thisFile->Get("contours"))->GetNkeys()-2;
512  std::vector<ContourCollection> universeContours;
513 
514  for (int i = 0; i < numUniverses; ++i){
515  universeContours.push_back(GetContours(CLType::kUniverse, i));
516  }
517 
518  return universeContours;
519 }
520 
521 // get the heat maps
522 //.............................................................................
523 std::vector<TH2F*> GetHeatMaps(){
524  std::vector<TH2F*> heatMaps;
526  for(size_t cl = 0; cl < 3; ++cl){
527  histName = "contours/Median/RandomUniverses" +
528  xaxis + yaxis + std::to_string(cl + 1) +
529  "sigmaHeatMap";
530  heatMaps.push_back((TH2F*)thisFile->Get(histName.c_str()));
531  histName = heatMaps.back()->GetTitle();
532  histName += " Heat Map";
533  heatMaps.back()->SetTitle(histName.c_str());
534  }
535  return heatMaps;
536 }
537 
538 //............................................................................
539 TH2D* GetHiddenParameter(std::string hiddenParName)
540 {
541  TH2D* hidPar = (TH2D*)thisFile->Get(("contours/Asimov/Asimov"+hiddenParName).c_str());
542 
543  std::string title("#theta_{13}");
544  if (hiddenParName == "Th23" ) title = "sin^{2}#theta_{23}";
545  else if(hiddenParName == "Th34" ) title = "sin^{2}#theta_{34}";
546  else if(hiddenParName == "dCP" ) title = "#delta_{CP}";
547  else if(hiddenParName == "Dmsq32" ) title = "#Delta m^{2}_{32}";
548  else if(hiddenParName == "Dmsq41" ) title = "#Delta m^{2}_{41}";
549 
550  hidPar->SetTitle(title.c_str());
551  hidPar->Scale(1./GetNFiles());
552  return hidPar;
553 }
554 
CLType clt
Definition: Utilities.h:128
bool logx
Definition: Utilities.h:48
TFile * thisFile
Definition: Utilities.h:44
const XML_Char * name
Definition: expat.h:151
tree Draw("slc.nhit")
enum BeamMode kRed
std::string cltName
Definition: Utilities.h:129
TH2D * GetMedianDeltaChiSqr()
Definition: Utilities.h:357
Int_t kPTLightBlue
Definition: Utilities.h:50
Contour(std::vector< TGraph * > c, CLNumber cn, CLType ct)
Definition: Utilities.h:88
std::vector< Contour > ContourCollection
Definition: Utilities.h:134
hmean SetLineStyle(2)
TH2D * GetDeltaChiSqr(CLType clt, int universeNum=-1)
Definition: Utilities.h:381
void SetCMFColorPalette()
Definition: Utilities.h:138
void getContoursFromDeltaChiSqr(TH2D *medianDeltaChiSqr, std::vector< std::vector< TGraph * > > &contourGraphs)
Definition: Utilities.h:287
TH2D * GetHiddenParameter(std::string hiddenParName)
Definition: Utilities.h:539
TCanvas * canv
CLType
Definition: Utilities.h:33
std::vector< TGraph * > contours
Definition: Utilities.h:126
std::map< std::string, TObject * > LegendMap
Definition: Utilities.h:135
std::vector< ContourCollection > GetUniverseContours()
Definition: Utilities.h:509
std::vector< TH2F * > GetHeatMaps()
Definition: Utilities.h:523
TGraph * GetBestFitPoint(CLType clt, int universeNum=-1)
Definition: Utilities.h:436
std::vector< std::string > CLNameLaTeX
Definition: Utilities.h:24
TGraph * GetUniverseBestFitPoints()
Definition: Utilities.h:456
static const char *const errMsg[]
Definition: Error.h:69
double lst
Definition: runWimpSim.h:113
hmean SetLineWidth(2)
Int_t kPTDarkBlue
Definition: Utilities.h:50
hmean SetLineColor(4)
int GetNFiles()
Definition: Utilities.h:159
std::vector< std::string > CLName
Definition: Utilities.h:15
ContourCollection GetMedianContours()
Definition: Utilities.h:419
std::string yaxis
Definition: Utilities.h:47
Int_t kPTRed
Definition: Utilities.h:50
std::string dirName
Definition: PlotSpectra.h:47
TDirectory * dir
Definition: macro.C:5
void cc()
Definition: test_ana.C:28
std::string clnLatex
Definition: Utilities.h:131
void DrawContour(ContourCollection cc, CLNumber cln=CLNumber::kEmpty, int colorOverride=-1, int styleOverride=-1)
Definition: Utilities.h:180
std::string clnName
Definition: Utilities.h:130
bool logy
Definition: Utilities.h:49
void SetGlobalVariables(std::string f, std::string xax, std::string yax, std::vector< int > contours)
Definition: Utilities.h:52
std::string xaxis
Definition: Utilities.h:46
CLNumber cln
Definition: Utilities.h:127
size_t numUniverses
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
prelim SetTextColor(kBlue+1)
static const double nb
Definition: Units.h:89
Float_t e
Definition: plot.C:35
enum BeamMode kBlue
void Save1DContour(CLType clt, std::string var, int universeNum=-1)
Definition: Utilities.h:217
TH2D * makeAndFillMedianDeltaChiSqrHist(std::vector< TH2D * > const &universeDeltaChiSqrHists)
Definition: Utilities.h:331
ContourCollection GetContours(CLType clt, int universeNum=-1)
Definition: Utilities.h:469
Int_t kPTOrange
Definition: Utilities.h:50
std::vector< CLNumber > contourLevels
Definition: Utilities.h:45
enum BeamMode string