median_contours.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Experiment/Dmsq32Constraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/FC/FCSurface.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "Utilities/rootlogon.C"
#include "OscLib/IOscCalc.h"

#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"

#include "TCanvas.h"
#include "TMarker.h"
#include "TBox.h"
#include "TLatex.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMarker.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TGaxis.h"

#include <algorithm>
#include <vector>
#include <string>
#include <fstream>

using namespace ana;

void DrawHieTag(int hie, bool th13=true, bool vert=false);

double sigma(TH1F* hist, double percentile){
     double xq[1] = {percentile};
     double yq[1];
     hist->GetQuantiles(1, yq, xq);
     return yq[0];
}

bool is_file_exist(const char *fileName)
{
    std::ifstream infile(fileName);
    std::cout<<fileName<<std::endl;
    return infile.good();
}

void median_contours(int universes=100,
                     bool corrSysts = false,
                     string folder = "2020bf",  // it will loop over all files in this folder
                     TString options="joint_realData_both_onlyIH",
                     bool dmsqSurf = false,
                     bool th13Surf = false,
                     bool zoomIn   = true,
                     bool fillContour = false)
{

  bool debug = 1;
  bool nueOnly = options.Contains("nueOnly");
  bool numuOnly = options.Contains("numuOnly");
  bool joint = options.Contains("joint");
  assert (nueOnly || numuOnly || joint);

  bool FHCOnly = options.Contains("FHCOnly");
  bool RHCOnly = options.Contains("RHCOnly");
  bool both = options.Contains("both");
  assert (FHCOnly || RHCOnly || both);

  bool fake2017 = options.Contains("fake2017");
  bool realData = options.Contains("realData");

  bool onlyNH = options.Contains("onlyNH");
  bool onlyIH = options.Contains("onlyIH");

  auto suffix = options;
  if(dmsqSurf) suffix+= "_dmsq";
  if(th13Surf) suffix+= "_th13";
  if(corrSysts) suffix+="_systs";
  else suffix+="_stats10";
  string syst = (corrSysts?"syst":"stat");

  TString outdir="./";

  //////////////////////////////////////////////////////////////////////
  // Make plots
  //////////////////////////////////////////////////////////////////////
  // load old analyses contours for comparisons 
  TFile * official =0;
  string outdir2018 = "";
  if(th13Surf) outdir2018 = "th13_delta/";
  else if (dmsqSurf) outdir2018 = "th23_dmsq/";
  else outdir2018 = "delta_th23/";
  TString ana2018 ("/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/" + outdir2018 + "stat/hist_contours_2018_" + options + "combo_"+(dmsqSurf?"dmsq_":"")+"stats");
  TFile * infile2018ana = new TFile (ana2018 + ".root", "read");
  auto mins2018 =* (TVectorD*)infile2018ana->Get("overall_minchi");
  double minchi23_2018ana = mins2018[0];

  TString ana2019 ("/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/contours/" + outdir2018 + "stat/hist_contours_2019_" + options + "_"+(dmsqSurf?"dmsq_":"")+"stats");
  TFile * infile2019ana = new TFile (ana2019 + ".root", "read");
  auto mins2019 =* (TVectorD*)infile2019ana->Get("overall_minchi");
  double minchi23_2019ana = mins2019[0];

  TString ana2020 ("/nova/ana/3flavor/Ana2020/Fits/RealData/root_contours/hist_contours_2020_" + options + "_"+(dmsqSurf?"dmsq_":"")+"stats");
  TFile * infile2020ana = new TFile (ana2020 + ".root", "read");
  auto mins2020 =* (TVectorD*)infile2020ana->Get("overall_minchi");
  double minchi23_2020ana = mins2020[0];
  
  // load universes  
  std::vector <TString> surfNames;
  if(!th13Surf && !dmsqSurf)surfNames = {"delta_th23"};
  if (dmsqSurf) surfNames.push_back("th23_dm32");
  if (th13Surf) surfNames.push_back("th13_delta");
  for(TString surfName:surfNames){
     for (int hie:{-1,+1}){
      if((onlyNH && hie<0) || (onlyIH && hie>0)) continue;
      
      auto surf_2018ana = *FrequentistSurface::LoadFrom(infile2018ana, "surf_"+surfName+(hie>0?"_NH":"_IH"));
      TH2 * surf_1Sigma_2018ana = Gaussian68Percent2D(surf_2018ana);
      TH2 * surf_2Sigma_2018ana = Gaussian2Sigma2D(surf_2018ana);
      TH2 * surf_3Sigma_2018ana = Gaussian3Sigma2D(surf_2018ana); 
      auto g1_2018ana = surf_2018ana.GetGraphs(surf_1Sigma_2018ana, minchi23_2018ana);
      auto g2_2018ana = surf_2018ana.GetGraphs(surf_2Sigma_2018ana, minchi23_2018ana);
      auto g3_2018ana = surf_2018ana.GetGraphs(surf_3Sigma_2018ana, minchi23_2018ana);

      auto surf_2019ana = *FrequentistSurface::LoadFrom(infile2019ana, "surf_"+surfName+(hie>0?"_NH":"_IH"));
      TH2 * surf_1Sigma_2019ana = Gaussian68Percent2D(surf_2019ana);
      TH2 * surf_2Sigma_2019ana = Gaussian2Sigma2D(surf_2019ana);
      TH2 * surf_3Sigma_2019ana = Gaussian3Sigma2D(surf_2019ana);
      auto g1_2019ana = surf_2019ana.GetGraphs(surf_1Sigma_2019ana, minchi23_2019ana);
      auto g2_2019ana = surf_2019ana.GetGraphs(surf_2Sigma_2019ana, minchi23_2019ana);
      auto g3_2019ana = surf_2019ana.GetGraphs(surf_3Sigma_2019ana, minchi23_2019ana);

      auto surf_2020ana = *FrequentistSurface::LoadFrom(infile2020ana, "surf_"+surfName+(hie>0?"_NH":"_IH"));
      TH2 * surf_1Sigma_2020ana = Gaussian68Percent2D(surf_2020ana);
      TH2 * surf_2Sigma_2020ana = Gaussian2Sigma2D(surf_2020ana);
      TH2 * surf_3Sigma_2020ana = Gaussian3Sigma2D(surf_2020ana);
      auto g1_2020ana = surf_2020ana.GetGraphs(surf_1Sigma_2020ana, minchi23_2020ana);
      auto g2_2020ana = surf_2020ana.GetGraphs(surf_2Sigma_2020ana, minchi23_2020ana);
      auto g3_2020ana = surf_2020ana.GetGraphs(surf_3Sigma_2020ana, minchi23_2020ana);
     
     TColor *c_coral2 = new TColor(102, 255/255.0, 169/255.0, 155/255.0);
     TColor *c_coral = new TColor(100, 255/255.0, 228/255.0, 218/255.0);
     TColor *c_wave = new TColor(101, 0/255.0, 167/255.0, 200/255.0); 
     Color_t kFillTestColor = ((hie>0)? 101:102); 
     Int_t kFillColorCont = TColor::GetColorTransparent(kFillTestColor, 0.3);
     if (onlyIH) kFillColorCont = 100;
     Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
     Int_t kFillColor;
     Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
     Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
     Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
     Int_t kDarkColor =  kGray+3; 

     int thr = 10; // submit jobs with 10 in one
     std::vector <TH2 *> hhist; // all contours
     std::vector <double> chis; // all minima

     for(int j=0; j< universes; j++){ 
       TString outfilename   ("/nova/ana/users/lkolupae/mock/contours"+ folder+ "/" + syst + "_hist_contours_"    + suffix + "_thr"+std::to_string(j) + ".root");
       if(is_file_exist(outfilename)){
         TFile * infile = new TFile (outfilename,"read");
         for (int k=0; k<10; k++){
              cout<<"univ "<<j<<endl;
              auto mins =* (TVectorD*)infile->Get(("universe_"+std::to_string(k)+"/bf_"+std::to_string(k)).c_str());
              double minchi23 = mins[0];
              cout<<minchi23<<endl;
              auto surf = *FrequentistSurface::LoadFrom(infile, "universe_"+std::to_string(k)+"/surf_"+surfName+"_univ"+std::to_string(k)+"_"+(hie>0?"NH":"IH"));
              hhist.push_back(surf.ToTH2(minchi23));
              chis.push_back(minchi23);
         }
         infile->Close();
       }
     }
     TH2D * median;
     TH2D * plsigma;
     TH2D * minsigma;
     // fill plus, minus and median histograms
     if (surfName=="delta_th23"){ 
          median = new TH2D("me", "", 30, 0-1/30.0, 2+1/30.0, 30, 0.3-0.2/30, 0.7+0.2/30);// plsigma, minsigma;
          plsigma = new TH2D("p", "", 30, 0-1/30.0, 2+1/30.0, 30, 0.3-0.2/30, 0.7+0.2/30);
          minsigma = new TH2D("mi", "", 30, 0-1/30.0, 2+1/30.0, 30, 0.3-0.2/30, 0.7+0.2/30);}
      if (surfName=="th23_dm32" && onlyNH){
          median = new TH2D("me", "", 20, 0.35-0.35/20, 0.7+0.35/20, 20, 2.1-0.6/20.0, 2.7+0.6/20.0);// plsigma, minsigma;
          plsigma = new TH2D("p", "", 20, 0.35-0.35/20, 0.7+0.35/20, 20, 2.1-0.6/20.0, 2.7+0.6/20.0);
          minsigma = new TH2D("mi", "", 20, 0.35-0.35/20, 0.7+0.35/20, 20, 2.1-0.6/20.0, 2.7+0.6/20.0);}
      if (surfName=="th23_dm32" && onlyIH){
          median = new TH2D("me", "", 20, 0.35-0.35/20, 0.7+0.35/20, 20, -2.7-0.5/20.0, -2.2+0.5/20.0);// plsigma, minsigma;
          plsigma = new TH2D("p", "", 20, 0.35-0.35/20, 0.7+0.35/20, 20, -2.7-0.5/20.0, -2.2+0.5/20.0);
          minsigma = new TH2D("mi", "", 20, 0.35-0.35/20, 0.7+0.35/20, 20, -2.7-0.5/20.0, -2.2+0.5/20.0);}
      // find plus, minus and median values for each bin
      for (int i = 0; i<=31; i++){
         for (int k = 0; k<=31; k++){
              // fill 1d histogram with all LL values for each universe
              TH1F* m = new TH1F("tmp", "tmp", 500, 0, 200);
              for(int j=0; j< hhist.size(); j++){
                    auto tmp = hhist[j]->GetBinContent(i, k);
                    m->Fill(tmp);
              }
              auto med = sigma(m, 0.5);
              median->SetBinContent(i, k, med);
              auto ps = sigma(m, 0.84);
              plsigma->SetBinContent(i, k, ps);
              auto ms = sigma(m, 0.16);
              minsigma->SetBinContent(i, k, ms);
              cout<<"m "<<med<<" ps "<<ps<<" min "<<ms<<endl;
         }
      } 

      cout<<"test "<<median->GetMinimum()<<" "<<plsigma->GetMinimum()<<" "<<minsigma->GetMinimum()<<endl;
      new TCanvas;
      median->Draw("colz");
      cout<<"minimum "<< median->GetMinimum()<<endl;
      gPad->Print(("median_"+suffix + "_" + folder + ".pdf"));

      new TCanvas;
      plsigma->Draw("colz");
      cout<<"minimum "<< plsigma->GetMinimum()<<endl;
      gPad->Print(("plsigma_"+suffix + "_" + folder + ".pdf"));

      new TCanvas;
      minsigma->Draw("colz");
      cout<<"minimum "<< minsigma->GetMinimum()<<endl;
      gPad->Print(("minsigma_"+suffix + "_" + folder +".pdf"));
      // for plus, minus and median surfaces find 1,2,3 sigma contours
      for(auto &sigma:{1,2,3}){
 
          if(surfName.Contains("delta_th23")) {
             if (zoomIn) DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
             else DrawContourCanvas(surfName, 0., 2., 0., 1.);
          }
          if(surfName.Contains("th23_dm32")) {
             if (hie>0) DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?1.8:2.05), (RHCOnly?3.3:2.85));
             else DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?-3.5:-2.95), (RHCOnly?-2.0:-2.25));
          }
          if(surfName.Contains("th13_delta")) {
             if (zoomIn) DrawContourCanvas(surfName, 0., 0.38, 0, 2);
             else DrawContourCanvas(surfName, 0., 1.0, 0, 2);
          }

          kFillColor = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);

          double value = 0;
          if (sigma == 1) value = 2.3;
          if (sigma == 2) value = 6.18;
          if (sigma == 3) value = 11.83;
         
          TVirtualPad* save = gPad; 
          TCanvas tmp;
          median->SetContour(1, &value);
          median->Draw("CONT Z LIST");
          tmp.Update();
          tmp.Paint();
          tmp.Print(("median_cont_" + folder + ".pdf").c_str());
          std::vector<TGraph*> gr_m;
          TCollection* specs = gROOT->GetListOfSpecials();
          TIter nextSpec(specs);
          while(TObject* spec = nextSpec()){
               if(!spec->InheritsFrom(TObjArray::Class())) continue;
               TObjArray* conts = (TObjArray*)spec;
               if(conts->IsEmpty()) continue;
               if(!conts->At(0)->InheritsFrom(TList::Class())) continue;
               TList* cont = (TList*)conts->At(0);
               TIter nextObj(cont);
               while(TObject* obj = nextObj()){
                     if(!obj->InheritsFrom(TGraph::Class())) continue;
                     gr_m.push_back((TGraph*)obj->Clone(UniqueName().c_str()));
               }
          }

          TCanvas tmp2;
          plsigma->SetContour(1, &value);
          plsigma->Draw("CONT Z LIST");
          tmp2.Update();
          tmp2.Paint();
          tmp2.Print(("plsigma_cont_" + folder + ".pdf").c_str());
          std::vector<TGraph*> gr_p;
          TCollection* specs_p = gROOT->GetListOfSpecials();
          TIter nextSpec_p(specs_p);
          while(TObject* spec_p = nextSpec_p()){
                if(!spec_p->InheritsFrom(TObjArray::Class())) continue;
                TObjArray* conts_p = (TObjArray*)spec_p;
                if(conts_p->IsEmpty()) continue;
                if(!conts_p->At(0)->InheritsFrom(TList::Class())) continue;
                TList* cont_p = (TList*)conts_p->At(0);
                TIter nextObj_p(cont_p);
                while(TObject* obj_p = nextObj_p()){
                      if(!obj_p->InheritsFrom(TGraph::Class())) continue;
                      gr_p.push_back((TGraph*)obj_p->Clone(UniqueName().c_str()));
                }
          }

          TCanvas tmp3;
          minsigma->SetContour(1, &value);
          minsigma->Draw("CONT Z LIST");
          tmp3.Update();
          tmp3.Paint();
          tmp3.Print(("minsigma_cont_" + folder + ".pdf").c_str());
          std::vector<TGraph*> gr_min;
          TCollection* specs_min = gROOT->GetListOfSpecials();
          TIter nextSpec_min(specs_min);
          while(TObject* spec_min = nextSpec_min()){
                if(!spec_min->InheritsFrom(TObjArray::Class())) continue;
                TObjArray* conts = (TObjArray*)spec_min;
                if(conts->IsEmpty()) continue;
                if(!conts->At(0)->InheritsFrom(TList::Class())) continue;
                TList* cont = (TList*)conts->At(0);
                TIter nextObj(cont);
                while(TObject* obj = nextObj()){
               if(!obj->InheritsFrom(TGraph::Class())) continue;
                  gr_min.push_back((TGraph*)obj->Clone(UniqueName().c_str()));
               }
          }
          gPad=save;
          cout<<"size of gr "<<gr_m.size()<<" p "<<gr_p.size()<<" min "<<gr_min.size()<<endl;

          if (folder == "2020max"){
             if(surfName.Contains("delta_th23") && onlyNH && sigma > 2) { ///contour "minus sigma"
                JoinGraphs(gr_min[0], gr_min[1], kFillColorCont)->Draw("f");
                for(TGraph* g: gr_min){
                    g->SetLineWidth(3);
                    g->SetLineStyle(kSolid);
                    g->SetLineColor(kFillColorCont);
                    g->Draw("l");
                }
             }
             else {
              for(TGraph* g: gr_min){
                  g->SetLineWidth(3);
                  g->SetLineStyle(kSolid);
                  g->SetLineColor(kFillColorCont);
                  g->SetFillColor(kFillColorCont);
                  g->Draw("f");
                  g->Draw("l");
              }
             }
             for(TGraph* g: gr_p){
                 g->SetLineWidth(3);
                 g->SetLineStyle(kSolid);
                 g->SetLineColor(kFillColorCont);
                 g->SetFillColor(kWhite);
                 g->Draw("f");
                 g->Draw("l");
             }
          }

          if (folder == "2020bf"){
              if(surfName.Contains("delta_th23") && sigma >3 && onlyIH || surfName.Contains("delta_th23") && onlyNH) { ///contour "minus sigma"
                 JoinGraphs(gr_min[0], gr_min[1], kFillColorCont)->Draw("f");
                 for(TGraph* g: gr_min){
                     g->SetLineWidth(3);
                     g->SetLineStyle(kSolid);
                     g->SetLineColor(kFillColorCont);
                     g->Draw("l");
                 }
              }
              else {
                for(TGraph* g: gr_min){
                    g->SetLineWidth(3);
                    g->SetLineStyle(kSolid);
                    g->SetLineColor(kFillColorCont);
                    g->SetFillColor(kFillColorCont);
                    g->Draw("f");
                    g->Draw("l");
                }
              }
              if (surfName.Contains("delta_th23") && sigma >2 && onlyNH){ ///contour "plus sigma"
                  if(sigma == 3)  JoinGraphs(gr_p[0], gr_p[1], kWhite)->Draw("f");
                  for(TGraph* g: gr_p){
                      g->SetLineWidth(3);
                      g->SetLineStyle(kSolid);
                      g->SetLineColor(kFillColorCont);
                      g->Draw("l");
                  }
              }
              else {
               for(TGraph* g: gr_p){
                   g->SetLineWidth(3);
                   g->SetLineStyle(kSolid);
                   g->SetLineColor(kFillColorCont);
                   g->SetFillColor(kWhite);
                   g->Draw("f");
                   g->Draw("l");
               }
              }
         }

         for(TGraph* g: gr_m){
            g->SetLineWidth(3);
            g->SetLineStyle(kSolid);
            g->SetLineColor(kFillColor);
            g->Draw("l");
         }

         auto gr2018 = g1_2018ana;
         auto gr2019 = g1_2019ana;
         auto gr2020 = g1_2020ana;
         if (sigma==1) cout<<"2018 size "<<gr2018.size()<< "2019 size "<<gr2019.size()<<" 2020 size "<<gr2020.size()<<endl;
         if (sigma==2) {gr2018 = g2_2018ana; gr2019 = g2_2019ana; gr2020 = g2_2020ana; cout<<"2018 size "<<gr2018.size()<<" 2019 size "<<gr2019.size() << " 2020 size "<<gr2020.size() <<endl;}
         if (sigma==3) {gr2018 = g3_2018ana; gr2019 = g3_2019ana; gr2020 = g3_2020ana; cout<<"2018 size "<<gr2018.size()<<" 2019 size "<<gr2019.size() << " 2020 size "<<gr2020.size() << endl; }
         
         for (auto & g:gr2018) { g->SetLineColor(kGray); g->Draw("l");} 
         for (auto & g:gr2019) { g->SetLineColor(kGray+2); g->Draw("l");}
         for (auto & g:gr2020) { g->SetLineColor(kViolet-6); g->Draw("l");}
         
         TLegend * leg = new TLegend(0.16,0.17,0.8,0.28);
         leg->SetTextSize(kBlessedLabelFontSize*0.75);
         leg->SetFillStyle(0);
         leg->SetMargin(1.3*leg->GetMargin());
         leg->SetNColumns(3);
         leg->SetMargin(1.4*leg->GetMargin());
         TGraph * dummy = new TGraph;
         dummy->SetLineWidth(3);
         dummy->SetLineColor(kFillColor);
         leg->AddEntry(dummy->Clone(),"median","l");
         dummy->SetFillColor(kFillColorCont);
         leg->AddEntry(dummy->Clone(),"68\% exp","f");
         leg->SetHeader((std::to_string(sigma) + " sigma contour").c_str());
   
         TLegend * leg2 = new TLegend(0.16,0.8,0.9,0.85);
         leg2->SetTextSize(kBlessedLabelFontSize*0.75);
         leg2->SetNColumns(3);
         dummy->SetLineColor(kGray);
         leg2->AddEntry(dummy->Clone(),"2018 stat.fit","l");
         dummy->SetLineColor(kGray+2);
         leg2->AddEntry(dummy->Clone(),"2019 stat fit","l");
         dummy->SetLineColor(kViolet-6);
         leg2->AddEntry(dummy->Clone(),"2020 stat fit","l");
         leg->Draw();
         leg2->Draw();

         Preliminary();
         DrawHieTag(hie,false);
         gPad->RedrawAxis();
        
         gPad->Print(outdir + "contour_" +  suffix + "_" + folder + "_" +
                      surfName + "_" + std::to_string(universes) + "_sigma" + std::to_string(sigma) + "_"+
                      (hie>0?"NH":"IH")  +
                      (zoomIn?"_zoom":"") +
                      ".pdf");
       }//end sigma
      }//end hie
     }//end surf
}//end all


void DrawHieTag(int hie,bool th13, bool vert){
  if(hie>0){
    TLatex * ltxNH = new TLatex(0.16,0.89,"#bf{NH}");
    ltxNH->SetTextColor(kPrimColorNH);
    ltxNH->SetNDC();
    ltxNH->SetTextAlign(22);
    ltxNH->SetTextSize(kBlessedLabelFontSize);
    if(th13)    ltxNH->Draw();
    else if(!vert) ltxNH->DrawLatex(.85,.20,"#bf{NH}");
    else ltxNH->DrawLatex(.85,.85,"#bf{NH}");
  }
  else{
    TLatex * ltxIH = new TLatex (0.16,0.46,"#bf{IH}");
    ltxIH->SetTextColor(kPrimColorIH);
    ltxIH->SetNDC();
    ltxIH->SetTextAlign(22);
    ltxIH->SetTextSize(kBlessedLabelFontSize);
    if(th13)    ltxIH->Draw();
    else if(!vert)  ltxIH->DrawLatex(.85,.20,"#bf{IH}");
    else ltxIH->DrawLatex(.85,.85,"#bf{IH}");
  }
}