makeMedianContours.C

Go to the documentation of this file.

#include "TFile.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TDirectory.h"
#include "TKey.h"
#include "TMarker.h"
#include "TH1.h"
#include "TH2.h"
#include "TLegend.h"
#include "TMath.h"
#include "TGraph.h"
#include "TPad.h"
#include "TROOT.h"

#include <string>
#include <map>
#include <iostream>

//------------------------------------------------------------------------------
// User should change the following values
//
size_t numUniverses = 500;
std::string inputFileName("cmf_contour_3F_Th23Dmsq32_stat_only_500_uni.root");
std::string pars("Th23Dmsq32");
//
// User intervention stops here
//------------------------------------------------------------------------------


//------------------------------------------------------------------------------
TH2F* makeAndFillMedianDeltaChiSqrHist(std::vector<TH2F*> const& universeDeltaChiSqrHists)
{
  // make the histogram we will return
  TH2F* medianDeltaChiSqr = dynamic_cast<TH2F*>(universeDeltaChiSqrHists.front()->Clone("medianDeltaChiSqr"));
  medianDeltaChiSqr->Reset();

  // loop over all the bins in the histograms and find the median for each bin
  std::vector<float> binDeltaChiSqrDist;
  for(int x = 1; x < universeDeltaChiSqrHists.front()->GetNbinsX() + 1; ++x){
    for(int y = 1; y < universeDeltaChiSqrHists.front()->GetNbinsY() + 1; ++y){
      binDeltaChiSqrDist.clear();
      for(auto const* itr : universeDeltaChiSqrHists){
        binDeltaChiSqrDist.emplace_back(itr->GetBinContent(x, y));
      }

      // sort the values
      std::sort(binDeltaChiSqrDist.begin(), binDeltaChiSqrDist.end());
      medianDeltaChiSqr->SetBinContent(x, y, binDeltaChiSqrDist[binDeltaChiSqrDist.size() / 2]);
    } // end loop over y bins
  }// end loop over x bins

  return medianDeltaChiSqr;

} // end makeAndFillMedianDeltaChiSqrHist

//------------------------------------------------------------------------------
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)
{
  // make a TCanvas and divide it into 2 x 2
  TCanvas *heatMapCanv = new TCanvas("heatMapCanv", "heatMapCanv", 1200, 800);

  heatMapCanv->Divide(2, 2);

  for(size_t c = 0; c < heatMaps.size(); ++c){
    heatMapCanv->cd(c + 1);
    heatMaps[c]->Draw("colz");

    if(c == 0){
      for(auto& itr : median1SigmaContours)
        itr->Draw("lsame");

      for(auto& itr : asimov1SigmaContours)
        itr->Draw("lsame");
    }
    if(c == 1){
      for(auto& itr : median2SigmaContours)
        itr->Draw("lsame");

      for(auto& itr : asimov2SigmaContours)
        itr->Draw("lsame");
    }
    if(c == 2){
      for(auto& itr : median3SigmaContours)
        itr->Draw("lsame");

      for(auto& itr : asimov3SigmaContours)
        itr->Draw("lsame");
    }

  }

  heatMapCanv->cd(4);

  TLegend *leg = new TLegend(0.1, 0.1, 0.9, 0.9);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->AddEntry(asimov1SigmaContours.front(), "Asimov - 1#sigma", "l");
  leg->AddEntry(asimov2SigmaContours.front(), "Asimov - 2#sigma", "l");
  leg->AddEntry(asimov3SigmaContours.front(), "Asimov - 3#sigma", "l");
  leg->AddEntry(median1SigmaContours.front(), "Median - 1#sigma", "l");
  leg->AddEntry(median2SigmaContours.front(), "Median - 2#sigma", "l");
  leg->AddEntry(median3SigmaContours.front(), "Median - 3#sigma", "l");
  leg->Draw();

  heatMapCanv->Update();
}

//------------------------------------------------------------------------------
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)
{
  // make a TCanvas and divide it into 2 x 2
  TCanvas *contoursCanv = new TCanvas("contoursCanv", "contoursCanv", 1200, 800);

  contoursCanv->Divide(2, 2);

  for(size_t c = 0; c < 3; ++c){
    contoursCanv->cd(c + 1);
    backdropHist->SetTitle((std::to_string(c + 1) + "#sigma").c_str());
    backdropHist->Draw();

    if(c == 0){
      for(auto & itr : universe1SigmaContours){
        itr->SetLineColor(kGray);
        itr->SetLineWidth(1);
        itr->Draw("lsame");
      }

      for(auto & itr : median1SigmaContours)
        itr->Draw("lsame");

      for(auto & itr : asimov1SigmaContours)
        itr->Draw("lsame");
    }
    if(c == 1){
      for(auto & itr : universe2SigmaContours){
        itr->SetLineColor(kGray);
        itr->SetLineWidth(1);
        itr->Draw("lsame");
      }

      for(auto & itr : median2SigmaContours)
        itr->Draw("lsame");

      for(auto & itr : asimov2SigmaContours)
        itr->Draw("lsame");
    }
    if(c == 2){
      for(auto & itr : universe3SigmaContours){
        itr->SetLineColor(kGray);
        itr->SetLineWidth(1);
        itr->Draw("lsame");
      }

      for(auto & itr : median3SigmaContours)
        itr->Draw("lsame");

      for(auto & itr : asimov3SigmaContours)
        itr->Draw("lsame");
    }
    contoursCanv->Update();
  }

  contoursCanv->cd(4);

  TLegend *leg = new TLegend(0.1, 0.1, 0.9, 0.9);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->AddEntry(asimov1SigmaContours.front(), "Asimov - 1#sigma", "l");
  leg->AddEntry(asimov2SigmaContours.front(), "Asimov - 2#sigma", "l");
  leg->AddEntry(asimov3SigmaContours.front(), "Asimov - 3#sigma", "l");
  leg->AddEntry(median1SigmaContours.front(), "Median - 1#sigma", "l");
  leg->AddEntry(median2SigmaContours.front(), "Median - 2#sigma", "l");
  leg->AddEntry(median3SigmaContours.front(), "Median - 3#sigma", "l");
  leg->Draw();

  contoursCanv->Update();
}

//------------------------------------------------------------------------------
void FindContours(std::vector<TGraph*> & oneSigmaGrs,
                  std::vector<TGraph*> & twoSigmaGrs,
                  std::vector<TGraph*> & threeSigmaGrs,
                  TDirectory*            directory)
{
  //cout << "FindUniverseContours from directory " << directory << endl;

  std::string keyName;
  TKey       *key;
  TIter next(directory->GetListOfKeys());
  while((key = dynamic_cast<TKey*>(next()))){
    keyName = key->GetName();

    if(keyName.find("1#sigma") != std::string::npos){
      oneSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
      oneSigmaGrs.back()->SetLineColor(kBlue);
      oneSigmaGrs.back()->SetLineWidth(3);
    }
    else if(keyName.find("2#sigma") != std::string::npos){
      twoSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
      twoSigmaGrs.back()->SetLineColor(kRed);
      twoSigmaGrs.back()->SetLineWidth(3);
    }
    else if(keyName.find("3#sigma") != std::string::npos){
      threeSigmaGrs.emplace_back(dynamic_cast<TGraph*>(key->ReadObj()));
      threeSigmaGrs.back()->SetLineColor(kMagenta);
      threeSigmaGrs.back()->SetLineWidth(3);
    }
  } // end loop to find the Xsigma graphs and Delta chi^2 surfaces

  //cout << "there are " << 1SigmaGrs.size() << " 1 sigma contours "
  //     << 2SigmaGrs.size() << " 2 sigma contours "
  //     << 3SigmaGrs.size() << " 3 sigma contours." << endl;
}


//------------------------------------------------------------------------------
void getContoursFromDeltaChiSqr(TH2F*                                medianDeltaChiSqr,
                                std::vector<std::vector<TGraph*> > & contourGraphs)
{
  std::vector<double> contourLevel({2.3, 6.18, 11.83});

  medianDeltaChiSqr->SetContour(contourLevel.size(), contourLevel.data());

  // This avoids disturbing canvases already
  // Yes, we do need a canvas for this to work
  // yes, we do need to call Update() on it
  // No, ROOT is not very friendly.
  TCanvas temp_can;
  temp_can.cd();
  medianDeltaChiSqr->Draw("contlist");
  temp_can.Update();

  auto *plah = dynamic_cast<TObjArray *>(gROOT->GetListOfSpecials()->FindObject("contours"));
  if(!plah) return;

  //cout << "TObjArray size: " << plah->GetSize() << endl;
  contourGraphs.resize(contourLevel.size());

  for(int i = 0; i < plah->GetSize(); ++i){
    std::vector<TGraph*> temp;
    auto *list = dynamic_cast<TList*>(plah->At(i));
    if(!list) break;
    if(!(list->First())) break;
    for(int igraph = 0; igraph < list->GetSize(); ++igraph){
      if(list->At(igraph)) temp.emplace_back(dynamic_cast<TGraph*>(list->At(igraph)));
      //cout << i << " " << igraph << " " << temp.back() << " " << list->At(igraph) << endl;
      //temp.back()->SetLineColor(lineColor);
      temp.back()->SetLineStyle(9);
      temp.back()->SetLineWidth(3);
      if     (i == 0) temp.back()->SetLineColor(kBlue);
      else if(i == 1) temp.back()->SetLineColor(kRed);
      else if(i == 2) temp.back()->SetLineColor(kMagenta);
    }
    contourGraphs[i] = std::move(temp);
  }

  gROOT->GetListOfSpecials()->FindObject("contours")->Clear();
}

//------------------------------------------------------------------------------
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)
{
  TCanvas *medianCanv = new TCanvas("medianCanv", "medianCanv", 1200, 800);

  // draw a histogram defining the boundaries of the space
  backdropHist->Draw();

  for(auto & itr : median1SigmaGrs) itr->Draw("lsame");
  for(auto & itr : median2SigmaGrs) itr->Draw("lsame");
  for(auto & itr : median3SigmaGrs) itr->Draw("lsame");

  for(auto & itr : asimov1SigmaGrs) itr->Draw("lsame");
  for(auto & itr : asimov2SigmaGrs) itr->Draw("lsame");
  for(auto & itr : asimov3SigmaGrs) itr->Draw("lsame");

  TLegend *leg = new TLegend(0.6, 0.6, 0.9, 0.9);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);

  leg->AddEntry(asimov1SigmaGrs[0], "Asimov - 1#sigma", "l");
  leg->AddEntry(asimov2SigmaGrs[0], "Asimov - 2#sigma", "l");
  leg->AddEntry(asimov3SigmaGrs[0], "Asimov - 3#sigma", "l");
  leg->AddEntry(median1SigmaGrs[0], "Median - 1#sigma", "l");
  leg->AddEntry(median2SigmaGrs[0], "Median - 2#sigma", "l");
  leg->AddEntry(median3SigmaGrs[0], "Median - 3#sigma", "l");
  leg->Draw();

  medianCanv->Update();
}

//------------------------------------------------------------------------------
void makeMedianContours()
{
  // open the TFile containing the contours for the
  // different universes
  TFile *tf = new TFile(inputFileName.c_str(), "READ");

  // get the asimov contours
  std::vector<TGraph*> asimov1SigmaContours;
  std::vector<TGraph*> asimov2SigmaContours;
  std::vector<TGraph*> asimov3SigmaContours;

  FindContours(asimov1SigmaContours,
               asimov2SigmaContours,
               asimov3SigmaContours,
               tf->GetDirectory("contours/Asimov"));

  // fill the vector of Delta chi^2 histograms for the different universes
  std::vector<TH2F*> universeDeltaChiSqrHists;
  std::vector<TGraph*> universe1SigmaGrs;
  std::vector<TGraph*> universe2SigmaGrs;
  std::vector<TGraph*> universe3SigmaGrs;

  std::string histName;
  std::string grName;
  std::string dirName;
  for(size_t u = 0; u < numUniverses; ++u){
    dirName  = "contours/Universe_" + std::to_string(u) + "/";
    histName = dirName + "Universe_" + std::to_string(u) + pars + "DeltaChiSqr";
    universeDeltaChiSqrHists.emplace_back(dynamic_cast<TH2F*>(tf->Get(histName.c_str())));

    FindContours(universe1SigmaGrs,
                 universe2SigmaGrs,
                 universe3SigmaGrs,
                 tf->GetDirectory(dirName.c_str()));
  }

  // now get the median Delta chi^2 histogram
  auto *medianDeltaChiSqr = makeAndFillMedianDeltaChiSqrHist(universeDeltaChiSqrHists);

  std::vector<std::vector<TGraph*> > medianContours;
  getContoursFromDeltaChiSqr(medianDeltaChiSqr, medianContours);

  // make a generic back drop histogram
  TH2F *backdropHist = dynamic_cast<TH2F*>(universeDeltaChiSqrHists.front()->Clone("BackDropHist"));
  backdropHist->Reset();

  makeMedianCanvas(backdropHist,
                   medianContours[0],
                   medianContours[1],
                   medianContours[2],
                   asimov1SigmaContours,
                   asimov2SigmaContours,
                   asimov3SigmaContours);

  makeUniverseContourPlots(backdropHist,
                           medianContours[0],
                           medianContours[1],
                           medianContours[2],
                           asimov1SigmaContours,
                           asimov2SigmaContours,
                           asimov3SigmaContours,
                           universe1SigmaGrs,
                           universe2SigmaGrs,
                           universe3SigmaGrs);


  // get the heat maps
  std::vector<TH2F*> heatMaps;
  for(size_t cl = 0; cl < 3; ++cl){
    histName = "contours/Median/RandomUniverses" + pars + std::to_string(cl + 1) + "sigmaHeatMap";
    heatMaps.emplace_back(dynamic_cast<TH2F*>(tf->Get(histName.c_str())));
  }

  makeHeatMapPlots(medianContours[0],
                   medianContours[1],
                   medianContours[2],
                   asimov1SigmaContours,
                   asimov2SigmaContours,
                   asimov3SigmaContours,
                   heatMaps);

}