nova_official_data_release.C

Go to the documentation of this file.

#include "TGraph.h"
#include "TFile.h"
#include "TH2.h"
#include "TVector.h"
#include "TCanvas.h"
#include "TLegendEntry.h"

#include "CAFAna/FC/FCSurface.h"
#include "3FlavorAna/Ana2018/nue/joint_fit_2018_tools.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "3FlavorAna/Ana2018/nue/joint_fit_2018_tools.h"
#include "Utilities/rootlogon.C"


#include <fstream>

using namespace ana;

//copied from CAFAna/nue/Ana2018/nue/FitandFC/joint_fit_2018_slices.C
TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName, 
                       bool fourierFit, std::string plot_name,
                       ofstream& summary, TString par, TString hier, TString oct,
                       bool useRoot = false, bool fccol=false);

TGraph * SmoothWithFourierFit(TGraph* gr0, 
                              bool isDelta = true,
                              bool useRoot = false,
                              int fourierFitN = 2);

void SmoothWithFourierFit(TH1* hist,
                          int nbins = 80, double minX = 0, double maxX = 2,
                          int fourierFitN=3);


TGraph* HistoToCyclicGraph(TH1* h, bool useRoot);

void AverageCyclicPoints(TGraph* g);

void AddCyclicPoints(TGraph* g);
// end from joint_fit_2018_slices.C

TH1D* ProfileAxes(TString par);

TFile* GetGaussianProfile(TString par)
{
  TString gaussDir = "/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/slices/syst/";
  TString baseName = "hist_slices_2017_joint_realData_bothcombo_systs";
  TString shortName;
  if(par.Contains("ssth23"))
    shortName = "_th23_noOct";
  else if(par.Contains("dmsq"))
    shortName = "_dmsq";
  else 
    shortName = "";
  TFile* ret = TFile::Open(gaussDir + baseName + shortName + ".root");
  return ret;
}

TString GetFCPath(TString par, TString hier)
{
  TString fcFolder = "/nova/ana/nu_e_ana/Ana2018/FC/";
  TString fcFileName = par + "_";
  fcFileName += hier.Contains("NH")? "nh":"ih";
  fcFileName += ".root";
  return fcFolder + fcFileName;
}

TString GetFCPath(TString par, TString hier, TString oct)
{
  TString fcFolder = "/nova/ana/nu_e_ana/Ana2018/FC/";
  hier.ToLower();
  if(par.Contains("ssth23"))
    return fcFolder + TString::Format("%s_%s.root", par.Data(), hier.Data());
  else
    return fcFolder + TString::Format("%s_%s_%s.root", par.Data(), oct.Data(), hier.Data());
}

TString MakeSummaryTable(FCSurface* fsurf, TTree* tfc,
                         TString par, TString hier, TString oct = "");

void nova_official_data_release(bool contours = true, bool slices = true)
{
  ofstream summary;
  summary.open("fc_summary_2018/summary_table.tex");
  
  TString outdir = "NOvA_2018_official_data_release/";
  system(TString::Format("mkdir %s", outdir.Data()).Data());
  // Let's do the contours first
  if(contours) {
    double bf_delta = 0.16576;
    double bf_th23 = 0.584509;
    double bf_dmsq32 = 2.50515;
    for(TString par : {"deltassth23", "ssth23dmsq32"}) {
      TMarker* bf = new TMarker(bf_delta, bf_th23, 43);
      if(par.Contains("ssth23dmsq32")) {
        bf->SetX(bf_th23);
        bf->SetY(bf_dmsq32);
      }
      TFile* official = new TFile(outdir + TString::Format("NOvA_2018_official_contour_%s.root", par.Data()), "recreate");                              
      bf->Write("NH_best_fit");
      
      for(TString hier: {"NH", "IH"}) {
        
        TFile* fGauss = GetGaussianSurface(par, hier);
        TFile* ffc = TFile::Open(GetFCPath(par, hier));
        TTree* tfc = (TTree*)ffc->Get("fc");
        TString rootDir = par.Contains("ssth23dmsq32")? "th23_dm32_" : "delta_th23_";
        auto surf = ISurface::LoadFrom(fGauss, "surf_" + rootDir + hier);
        TH2* hGauss = (TH2*) fGauss->Get(rootDir+hier);
        auto mins =* (TVectorD*) fGauss->Get("overall_minchi");
        double minchi23 = mins[0];
        
        int NX = hGauss->GetXaxis()->GetNbins();
        int NY = hGauss->GetYaxis()->GetNbins();
        double minX = hGauss->GetXaxis()->GetBinLowEdge(1);
        double maxX = hGauss->GetXaxis()->GetBinUpEdge(NX);
        double minY = hGauss->GetYaxis()->GetBinLowEdge(1);
        double maxY = hGauss->GetYaxis()->GetBinUpEdge(NY);
        
        std::cout << "Making FCSurface with " << std::endl;
        std::cout << "------------------------------------" << std::endl;
        std::cout << "NX: " << NX << "\tNY: " << NY << std::endl;
        std::cout << "minX: " << minX << "\tmaxX" << maxX << std::endl;
        std::cout << "minY: " << minY << "\tmaxY" << maxY << std::endl;
        
        FCSurface* fc = new FCSurface(NX, minX, maxX, NY, minY, maxY);

        auto fccol = FCCollection::LoadFromFile(GetFCPath(par, hier).Data()).release();
        fc->Add(*fccol, par.Data());

        // output TeX formatted summary table
        summary << MakeSummaryTable(fc, tfc , par, hier) << std::endl;


        TH2* fc_surf_1sig = fc->SurfaceForSignificance(0.6827); fc_surf_1sig->Smooth();
        TH2* fc_surf_2sig = fc->SurfaceForSignificance(0.9545); fc_surf_2sig->Smooth();
        TH2* fc_surf_3sig = fc->SurfaceForSignificance(0.9973);
        TH2* fc_surf_90CL = fc->SurfaceForSignificance(0.90); fc_surf_90CL->Smooth();
        if(par.Contains("deltassth23"))
          CylindricalSmoothing(fc_surf_3sig);
        else
          fc_surf_3sig->Smooth();
                
        std::vector<TGraph*> g1 = surf->GetGraphs(fc_surf_1sig, minchi23);
        std::vector<TGraph*> g2 = surf->GetGraphs(fc_surf_2sig, minchi23);
        std::vector<TGraph*> g3 = surf->GetGraphs(fc_surf_3sig, minchi23);
        std::vector<TGraph*> g90 = surf->GetGraphs(fc_surf_90CL, minchi23);
        
        // Start writing stuff out
        official->cd();

        // canvas axes are different than what is returned from the surface object
        // as a style choice. Values are pulled from
        // CAFAna/nue/Ana2018/FitandFC/joint_fit_2018_contours.C
        double axesMaxX; double axesMaxY;
        double axesMinX; double axesMinY;
        if(par.Contains("deltassth23")){
          axesMinX = 0.;
          axesMaxX = 2.;
          axesMinY = 0.275;
          axesMaxY = 0.725;
        }
        else {
          axesMinX = 0.3;
          axesMaxX = 0.7;
          if(hier.Contains("NH")) {
            axesMinY = 2.05;
            axesMaxY = 2.85;
          }
          else {
            axesMinY = -2.95;
            axesMaxY = -2.25;
          }
        }
            
        TH2F* axes = new TH2F("axes", "", 100, axesMinX, axesMaxX, 100, axesMinY, axesMaxY);
        TString xLabel = "#delta_{CP}/#pi";
        TString yLabel = "sin^{2}#theta_{23}";
        if(par.Contains("ssth23dmsq32")){
          xLabel = "sin^{2}#theta_{23}";
          yLabel = "#Deltam^{2}_{32} (10^{-3} eV^{2})";
        }
        axes->GetXaxis()->SetTitle(xLabel.Data());
        axes->GetYaxis()->SetTitle(yLabel.Data());
        CenterTitles(axes);
        
        // For consistency, only save one axis for deltassth23 as just 'axes'
        TString axesName = hier + "_axes";
        if(par.Contains("deltassth23"))
          axesName = "axes";
        if(!(par.Contains("deltassth23") && hier.Contains("IH"))){
          axes->Write(axesName.Data());
        }
        std::vector<std::pair< TString, std::vector<TGraph*> > > graphs = {{"1sigma",g1},
                                                                           {"2sigma", g2},
                                                                           {"3sigma", g3},
                                                                           {"90CL", g90}};
        for(std::pair<TString, std::vector<TGraph*>> gn : graphs){
          int i  = 0;
          for(TGraph*  g : gn.second){
            g->Write(TString::Format("%s_%s_%d", 
                                     hier.Data(),
                                     gn.first.Data(),
                                     i));
            i++;
          }
        }
      }
    }
  }

  // 1D profiles
  if(slices){
    std::map< TString, double > maxY = {{"delta", 5.0},
                                        {"ssth23", 3.0},
                                        {"dmsq", 3.0}};
    for(TString par : {"delta", "ssth23", "dmsq"}) {
      TString parShortName = par;
      if(par.Contains("ssth23")) parShortName = "th23";
      TFile* official = new TFile(outdir + TString::Format("NOvA_2018_official_profile_%s.root", par.Data()), "recreate");
      TH1D* axes = ProfileAxes(par);
      axes->SetMaximum(maxY[par]);
      axes->Write("axes");
      
      TFile* fGauss = GetGaussianProfile(par);    
      // this will hold the profiles of each parameter
      
      
      for(TString hier : {"IH", "NH"}) {
        for(TString oct : {"LO", "UO"}) {
          if(par.Contains("ssth23") && oct.Contains("LO"))
            continue; // only do ssth23 once for each hierarchy
          
          if(par.Contains("ssth23"))
            oct = "";
          TString histName = "slice_" + parShortName + "_" + hier + oct;
          std::cout << "Looking for " << histName << " in " << fGauss->GetName() << std::endl;
          TH1D* h = (TH1D*) fGauss->Get(histName);
          TGraph* graph = new TGraph(h);
          
          // only use fourier fit for period boundary conditions
          // Fourier fit doesn't like delta_UO_nh. Would need 
          // lots of harmonics to fit it correctly anyway
          bool useFourier = (par.Contains("delta") && 
                             !(hier.Contains("NH") && oct.Contains("UO")));
          
          bool smooth = true;
          graph = FCCorrectSlice(graph, GetFCPath(par, hier, oct),
                                 smooth, par.Data(),
                                 summary, par, hier, oct,
                                 !useFourier, true);
          
          if(par.Contains("dmsq")){ //HACK
            for (int i =0; i < graph->GetN(); ++i){
              graph->SetPoint(i,1000*fabs(graph->GetX()[i]),graph->GetY()[i]);
            }
          }
          official->cd();
          graph->Write(hier+oct);
        }
      }
    }
  }

  summary.close();
}

TString MakeSummaryTable(FCSurface* fcsurf,  TTree* tfc,
                         TString par, TString hier, TString oct)
{
  int maxNPts = 0;
  int minNPts = 100000;
  int totNPts = tfc->GetEntries();
  int NBins = 0;
  int empty;
  double avgNPts;
  for(int i = 1; i < fcsurf->Binning()->GetXaxis()->GetNbins()+1; i++){
    for(int j = 1; j < fcsurf->Binning()->GetYaxis()->GetNbins()+1; j++){
      if(fcsurf->GetBin(i, j)->Empty()) {
        empty++;
        continue;
      }
      int NPts =  fcsurf->GetBin(i, j)->NPts();
      maxNPts = std::max(maxNPts, NPts);
      minNPts = std::min(minNPts, NPts);
      NBins++;
    }
  }
  avgNPts = totNPts / double(NBins);

  TString label = "";
  if(par == "ssth23dmsq32")
    label = "$\\sin^2\\theta_{23} vs \\Deltam^2_{32}$";
  if(par == "deltassth23")
    label = "$\\delta_{CP} vs \\sin^2\\theta_{23}$";
  if(par == "delta")
    label = "$\\delta_{CP}$";
  if(par == "ssth23")
    label = "$\\sin^2\\theta_{23}$";
  if(par == "dmsq")
    label = "$\\Deltam^2_{32}$";
  label += " ";
  label += oct;
  label += hier;
  
  TString row(TString::Format("%s &",label.Data()));
  row += TString::Format("%d & ", NBins);
  row += TString::Format("%d & ", totNPts);
  row += TString::Format("%d & ", maxNPts);
  row += TString::Format("%d & ", minNPts);
  row += TString::Format("%.2f \\\\ ", avgNPts);

  std::cout << row << std::endl;

  return row;

}

TH1D* ProfileAxes(TString par)
{
  double minX = 0;
  double maxX = 2;
  TString xLabel = "#delta_{CP}/#pi";
  TString yLabel = "Significance (#sigma)";
  // min and max values of X axis are copied from 
  // CAFAna/nue/Ana2018/FitandFC/joint_fit_2018_slices.C
  if(par.Contains("dmsq")){
    minX = 2.15;
    maxX = 2.75;
    xLabel = "|#Deltam^{2}_{32}| (10^{-3} eV^{2})";
  }
  else if(par.Contains("th23")){
    minX = 0.32;
    maxX = 0.72;
    xLabel = "sin^{2}#theta_{23}";
  } 
  
  TH1D* ret = new TH1D("axes","", 100, minX, maxX);
  ret->GetXaxis()->SetTitle(xLabel);
  ret->GetYaxis()->SetTitle(yLabel);
  CenterTitles(ret);

  return ret;
}

TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName, 
                       bool fourierFit, std::string plot_name,
                       ofstream& summary,
                       TString par, TString hier, TString oct,
                       bool useRoot, bool fccol)
{
  bool th23plot = fcFileName.Contains("th23");
  bool dmsqplot = fcFileName.Contains("dmsq");
  TString hiername = fcFileName.Contains("nh")? "NH":"IH";
  double minX = 0;
  double maxX = 2;
  if (th23plot) {minX = 0.3; maxX = 0.7;}
  if (dmsqplot) {
    minX = hiername=="NH" ? 2e-3 : -3e-3; 
    maxX = hiername=="NH" ? 3e-3 : -2e-3; 
  }
  //We want slice, not sqrt
  auto n= sqrtslice->GetN();
  auto sliceXH = new TGraph(n);
  for (int i = 0; i < n;i++ ){
    sliceXH->SetPoint(i,sqrtslice->GetX()[i],util::sqr(sqrtslice->GetY()[i]));
  }
  FCSurface *fcXH = 0;
  if(!fccol)  fcXH = FCSurface::LoadFromFile(fcFileName.Data()).release();
  else {
    auto fccol = FCCollection::LoadFromFile(fcFileName.Data()).release(); 
    //Only one Y bin required for slices.
    // All slices this year had 60 bins. Need to find a better way
    // next year
    fcXH = new FCSurface(60,minX,maxX,1,0,1);
    fcXH->Add(*fccol, plot_name);
  }
  TFile* ffc = TFile::Open(fcFileName.Data());
  TTree* tfc = (TTree*) ffc->Get("fc");

  // output TeX formatted summary table
  summary << MakeSummaryTable(fcXH, tfc , par, hier, oct) << std::endl;

  TGraph* sigXH = new TGraph;
  TGraph* pcXH = new TGraph;
  for(int i = 0; i < n; ++i){
    const double delta = sliceXH->GetX()[i];
    const double upXH = sliceXH->GetY()[i];
    // Fall back to gaussian up if dchisq is too large or FCBin is empty.
    // Harder to believe these dchisq are covered with our
    // pseudoexperiments anyway
    if(upXH > 14 || fcXH->GetBin(i+1,1)->Empty()){
      sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
    }
    else {
      double pXH;
      if(upXH > 0) pXH = fcXH->GetBin(i+1, 1)->SignificanceForUpValue(upXH);
      else pXH = 0; //HACK HACK
      sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
    }
  }
  if(fourierFit) {
    sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
  }
  return sigXH; 
}


TGraph * SmoothWithFourierFit(TGraph* gr0, bool isDelta, bool useRoot, int fourierFitN)
{
  int nbins = gr0->GetN();
  double minX = gr0->GetX()[0];
  double maxX = gr0->GetX()[nbins-1];
  //TGraph to center bin
  minX = minX  - (gr0->GetX()[1] - gr0->GetX()[0])/2;
  maxX = maxX  + (gr0->GetX()[1] - gr0->GetX()[0])/2;
  //Fourier fit  relies on periodicity, ignore extra points
  if(isDelta && !useRoot){
    nbins = 60;
    minX = 0;
    maxX = 2;
  }
  std::cout << "Smooth " << nbins << " nbins " << minX << ", " << maxX << "\n";
  TH1D * hist = new TH1D(UniqueName().c_str(),"",nbins,minX,maxX);
  for (int i = 0; i < nbins; ++i){
    double val = gr0->Eval(hist->GetBinCenter(i+1));
    hist->SetBinContent(i + 1, val * val); //use dchisq not chi2 to smooth
  }
  if(useRoot) hist->Smooth(2);
  else  SmoothWithFourierFit(hist, nbins, minX, maxX, fourierFitN);
  //now return sqrt
  for (int i = 1; i <= nbins; ++i){
    hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
  }
  if(isDelta) return HistoToCyclicGraph(hist, useRoot);
  else return new TGraph(hist);
}


TGraph* HistoToCyclicGraph(TH1* h, bool useRoot = false)
{
  TGraph* g = new TGraph;
  g->SetPoint(0, 0, 0); // temporary
  for(int i = 1; i <= h->GetNbinsX(); ++i)
    g->SetPoint(i, h->GetXaxis()->GetBinCenter(i), h->GetBinContent(i));
  g->SetLineColor(h->GetLineColor());
  g->SetLineWidth(2);
  
  if(useRoot)
    AverageCyclicPoints(g);
  else
    AddCyclicPoints(g);

  return g;
}
void SmoothWithFourierFit(TH1* hist, int nbins, double minX, double maxX, int fourierFitN)
{
  FitToFourier fitF (hist, minX, maxX, fourierFitN);
  auto f1= fitF.Fit();
  hist->Reset();
  hist->Add(f1);
}

void AddCyclicPoints(TGraph* g)
{
  double x, y;
  g->GetPoint(g->GetN()-1, x, y);
  g->SetPoint(0, x-2, y);
  g->GetPoint(1, x, y);
  g->SetPoint(g->GetN(), x+2, y);
  std::cout << "\n Added cyclic points " << 0 << " " << g->GetX()[0]<< " " <<  g->GetY()[0]
            << ", " << g->GetN() <<  x+2 <<  y << "\n\n";
}
void AverageCyclicPoints(TGraph* g)
{
  double x0, y0;
  double x2pi, y2pi;
  g->GetPoint(g->GetN()-1, x2pi, y2pi);
  g->GetPoint(1, x0, y0);
  std::cout << "(x0, y0) = (" << x0 << ", " << y0 << ")\n";
  std::cout << "(x2pi, y2pi) = (" << x2pi << ", " << y2pi << ")\n";
  
  double avgY = ( y0 + y2pi ) / 2;
  std::cout << "Average = " << avgY << std::endl;
  g->SetPoint(0, 0, avgY);
  g->SetPoint(g->GetN(), 2, avgY);
  
  std::cout << "\nFound average cyclic point: (0, " << avgY << ")\n";
  std::cout << "\nAdded average cyclic points " << 0 << ": " << g->GetX()[0] << " " << g->GetY()[0];
  std::cout << ", " << g->GetN() << ": " << g->GetX()[g->GetN()-1] << " " << g->GetY()[g->GetN()-1] << "\n\n";
}