DrawExtrapSurface.C

Go to the documentation of this file.

//////////////////////////////////////////////////////////////
// DrawExtrapSurface.C
// Mike Wallbank (wallbank@fnal.gov)
//
// Make 1D profiles and 2D surface in th24-th34 parameter
// space, using extrapolated predictions.
//////////////////////////////////////////////////////////////

// framwork
#include "CAFAna/Fit/FrequentistSurface.h"
#include "NuXAna/macros/Nus18/Nus18Utilities.h"

// root
#include "TFile.h"
#include "TCanvas.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TH2.h"
#include "TLegend.h"
#include "TText.h"
#include "TLine.h"

// stl
#include <iomanip>

using namespace ana;

void Draw1DProfile(TH1* hist, std::string name, TFile* outFile);
void Draw2DSurface(FrequentistSurface* surface, std::string name, TFile* outFile);

// -------------------------------------------------------------------------------------------
void DrawExtrapSurface(std::string inFileName, TString opt) {

  TH1::AddDirectory(0);

  // Get running options
  bool oned = opt.Contains("1D", TString::kIgnoreCase);
  bool twod = opt.Contains("2D", TString::kIgnoreCase);
  assert(oned or twod);

  // Get the input and output file names
  std::string fileName = inFileName.substr(0, inFileName.find(".root"));
  std::string saveName = fileName + "Ana.root";

  // Load the analysis objects from the in file
  TFile* inFile = new TFile(inFileName.c_str(), "READ");
  TH1* ProfileTh24 = (TH1*)inFile->Get("ProfileTh24");
  TH1* ProfileTh34 = (TH1*)inFile->Get("ProfileTh34");
  FrequentistSurface* Surface = FrequentistSurface::LoadFrom(inFile, TString("SurfaceTh34Th24")).release();

  // Drawing
  TFile* outFile = new TFile(saveName.c_str(), "RECREATE");

  // Draw 1D profiles
  if (oned) {
    Draw1DProfile(ProfileTh24, "1DProfileTheta24", outFile);
    Draw1DProfile(ProfileTh34, "1DProfileTheta34", outFile);
  }

  // Draw 2D surface
  if (twod)
    Draw2DSurface(Surface, "2DSurfaceTheta24Theta34", outFile);

  // Clean up
  outFile->Close();
  delete outFile;

}

// -------------------------------------------------------------------------------------------
void Draw1DProfile(TH1* hist, std::string name, TFile* outFile) {

  TCanvas* canv = new TCanvas("canv", "", 800, 600);
  TLine line;
  TText text;

  float xmin = hist->GetBinLowEdge(1);
  float xmax = hist->GetBinLowEdge(hist->GetXaxis()->GetLast()+1);
  float ymin = hist->GetMinimum();
  float ymax = hist->GetMaximum();

  canv->cd();
  canv->Clear();
  hist->Draw("hist");
  line.SetLineStyle(kDashed);
  line.SetLineColor(kBlack);
  line.SetLineWidth(2);
  text.SetTextColor(kBlack);
  text.SetTextFont(42);
  text.SetTextSize(0.037);
  text.SetTextAlign(11);
  line.DrawLine(xmin, 1, xmax, 1);
  text.DrawText(0+(xmax-xmin)/10., 1+(ymax-ymin)/100., "68% C.L.");
  line.DrawLine(xmin, 9, xmax, 9);
  text.DrawText(0+(xmax-xmin)/10., 9+(ymax-ymin)/100., "90% C.L.");
  Preliminary();
  outFile->cd();
  canv->Write(name.c_str());

  delete canv;

  return;

}

// -------------------------------------------------------------------------------------------
void Draw2DSurface(FrequentistSurface* surface, std::string name, TFile* outFile) {

  TCanvas* canv = new TCanvas("canv", "", 800, 600);
  TLine line;
  TText text;

  surface->DrawContour(Gaussian68Percent1D(*surface), kDashed, kOfficialNus18RHCColor);
  surface->DrawContour(Gaussian90Percent1D(*surface), kSolid,  kOfficialNus18RHCColor);

  // Draw legend
  //TLegend* leg = new TLegend(0.6,0.6,0.88,0.85);
  TLegend* leg = new TLegend(0.645,0.7,0.82,0.88);
  TH1* dummy = new TH1F("","",1,0,1);
  dummy->SetLineColor(kOfficialNus18RHCColor);
  dummy->SetLineWidth(2);
  dummy->SetLineStyle(kDashed);
  leg->AddEntry(dummy->Clone(), "68% C.L.", "l");
  dummy->SetLineStyle(kSolid);
  leg->AddEntry(dummy->Clone(), "90% C.L.", "l");
  leg->Draw();

  // Annotate
  Simulation();
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  calc->SetNFlavors(4);
  SetAngles(calc);
  std::stringstream description, params4f, params3f, paramsDelta;
  description << "#bf{NOvA Sensitivity}, 12 #times 10^{20} POT";
  params4f << "#Deltam^{2}_{41} = " << setprecision(3) << calc->GetDm(4) << " eV^{2}";
  params3f << "sin^{2}#theta_{23} = " << setprecision(3) << TMath::Power(TMath::Sin(calc->GetAngle(2,3)),2) << ", "
           << "#Deltam^{2}_{32} = " << setprecision(3) << (calc->GetDm(3)-calc->GetDm(2))/1e-3 << " #times 10^{-3} eV^{2}";
  paramsDelta << "#delta_{13} = " << setprecision(3) << calc->GetDelta(1,3)/TMath::Pi() << "#pi";
    //<< "#delta_{24} = " << setprecision(3) << calc->GetDelta(2,4)/TMath::Pi() << "#pi";
  TLatex* tex = new TLatex();
  tex->SetNDC();
  tex->SetTextFont(42);
  tex->SetTextSize(0.037);
  tex->SetTextAlign(11);
  tex->DrawLatex(0.135, 0.82, description.str().c_str());
  tex->DrawLatex(0.135, 0.82 - 1*0.055, params3f.str().c_str());
  tex->DrawLatex(0.135, 0.82 - 2*0.055, params4f.str().c_str());
  tex->DrawLatex(0.135, 0.82 - 3*0.055, paramsDelta.str().c_str());
  canv->RedrawAxis();

  outFile->cd();
  canv->Write(name.c_str());

  return;

}