template_basic.C

Go to the documentation of this file.

// most basic template.  FHC only, no cosmics file, no systematics, 2 choices of theta_23

#include "OscLib/OscCalcGeneral.h"
#include "OscLib/OscCalcPMNS.h"
#include "OscLib/OscCalc.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Core/Utilities.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Binning.h"


#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TMath.h"
#include "TGaxis.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TStyle.h"
#include "THStack.h"
#include "TPaveText.h"

#include <cmath>
#include <iostream>
#include <sstream>
#include <string>
#include <fstream>
#include <iomanip>

#include "Utilities/func/MathUtil.h"

// Use with Nova style ROOT logon script

#define location 0 // 0 for FNAL 1 for caltech
#define mcset 1 // 0 = ideal MC goes with pot like 18e20, or 1 = realistic conditions mc for 1st analysis POT
#define pot 18e20

using namespace ana;

void Simulation()  // for some reason doesn't always remember definition in Style script
{
  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Simulation");
  prelim->SetTextColor(kGray+1);
  prelim->SetNDC();
  prelim->SetTextSize(2/30.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

void template_basic()
{

  TGaxis::SetMaxDigits(2);

  gStyle->SetTitleOffset(.85, "y");
  gStyle->SetTitleOffset(.84, "x");
  gStyle->SetTitleSize(0.05,"x");

  std::string fname, fnameSwap;

  if(mcset==1) {

    if(location==1) {
      fname = "$NOVA_DATA/mc/forFA/fdgenie.nonswap.decaf.root";
      fnameSwap = "$NOVA_DATA/mc/forFA/fdgenie.swap.decaf.root";
    }
    else if(location==0) {
      fname = "prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_numu_contain";
      fnameSwap = "prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_numu_contain";
    }
  }

  else if(mcset==0) {

    if(location==1) {
      fname = "prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_ideal_numu_contain";
      fnameSwap = "prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_ideal_numu_contain";
    }
    else if(location==0) {
      fname = "prod_decaf_S15-05-22a_fd_genie_fhc_nonswap_ideal_numu_contain";
      fnameSwap = "prod_decaf_S15-05-22a_fd_genie_fhc_fluxswap_ideal_numu_contain";
    }
  }

  osc::OscCalcPMNS calc;

  calc.SetL(810);
  calc.SetRho(0); // No matter effects
  calc.SetDmsq21(7.59e-5);
  calc.SetDmsq32(2.4e-3);
  calc.SetTh12(.601);
  calc.SetTh13(.1567);
  //  calc.SetdCP(TMath::Pi()/2);
  calc.SetdCP(0);
  calc.SetTh23(TMath::Pi()/4);  // sin^2(2theta) = 1

  const std::string filename = fname;
  const std::string filenameSwap = fnameSwap;

  SpectrumLoader loader(filename);
  SpectrumLoader loaderSwap(filenameSwap);

  // Load up the necessary histograms, seperated by flavour etc
  PredictionNoExtrap predCC(loader, loaderSwap,
                            "CC Neutrino Energy (GeV)", kNumuEnergyBinning,
                            kCCE, kNumuFD);

  loader.Go();
  loaderSwap.Go();

  Spectrum obsCC = predCC.Predict(&calc);  // the prediction
  Spectrum fakeCC = obsCC.FakeData(pot);  // the data, faked
  SingleSampleExperiment exptCC(&predCC, fakeCC);  // made into an 'experiment'

  TLatex text;
  text.SetNDC(true);
  text.SetTextSize(0.04);

  // now, the same but with a different truth value of theta_23
  calc.SetTh23(.67264);  // sin^2(2theta) = 0.95

  Spectrum obsCC2 = predCC.Predict(&calc);
  Spectrum fakeCC2 = obsCC2.FakeData(pot);
  SingleSampleExperiment exptCC2(&predCC, fakeCC2);

  new TCanvas("comp1b","CC spectra .95");  // plot spectra for each sample
  DataMCComparison(fakeCC2, obsCC2);
  TH1* ccu2 = predCC.PredictUnoscillated().ToTH1(pot);
  ccu2->SetLineStyle(7);
  ccu2->Draw("hist same");
  gPad->Print("plots/comp1b.ps");

  calc.SetTh23(1); // Seed value

  TH1F *hcc = new TH1F("cc", "cc", 10, 0, 1);
  TH1F *hqe = new TH1F("qe", "qe", 10, 0, 1);
  TH1F *huc = new TH1F("uc", "uc", 10, 0, 1);
  TH1F *hall = new TH1F("all", "all", 10, 0, 1);

  hcc->SetLineWidth(2);
  hqe->SetLineWidth(2);
  huc->SetLineWidth(2);
  hall->SetLineWidth(2);
  hcc->SetLineColor(2);
  hqe->SetLineColor(4);
  huc->SetLineColor(kMagenta);
  hall->SetLineColor(1);

  //The log-likelihood surface: FHC
  new TCanvas("contoursfhc1","contoursfhc1") ;

  FrequentistSurface surfCC(&exptCC, &calc,
                 &kFitSinSq2Theta23, 30, .9, 1,
                 &kFitDmSq32Scaled, 30, 2.2, 2.7);
  surfCC.DrawBestFit(kBlack);
  surfCC.DrawContour(Gaussian90Percent2D(surfCC), kSolid, kBlack);

  Simulation();

  TLegend *leg = new TLegend(0.11,0.11,0.42,0.35);
  leg->SetFillColor(10);
  leg->SetFillStyle(0);
  leg->SetBorderSize(0);
  if(pot==18e20) leg->SetHeader("18e20 POT, FHC, 14 kton: 90% CL");
  else if(pot==6e20) leg->SetHeader("6e20 POT, FHC, 14 kton: 90% CL");
  else if(pot==2e20) leg->SetHeader("2e20 POT, FHC, 14 kton: 90% CL");
  else if(pot==1e20) leg->SetHeader("1e20 POT, FHC, 14 kton: 90% CL");

  leg->AddEntry(hall,"Contained CC","L");
  leg->Draw("same");
  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");

  gPad->Print("plots/basic.1.ps");
  gPad->Print("plots/basic.1.C");

  //2

  new TCanvas("contoursfhc2","contoursfhc2") ;

  FrequentistSurface surfCC2(&exptCC2, &calc,
                  &kFitSinSq2Theta23, 30, .87, 1,
                  &kFitDmSq32Scaled, 30, 2.2, 2.7);
  surfCC2.DrawBestFit(kRed);
  surfCC2.DrawContour(Gaussian90Percent2D(surfCC2), kSolid, kRed);

  leg->Draw("same");
  text.DrawLatex(0.3,0.85,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");

  Simulation();

  gPad->Print("plots/basic.fhc.2.ps");
  gPad->Print("plots/basic.fhc.2.C");

  new TCanvas("contourscombine","contours90percent") ;
  surfCC2.DrawBestFit(kRed);
  surfCC2.DrawContour(Gaussian2Sigma2D(surfCC2), kSolid, kRed);
  surfCC2.DrawContour(Gaussian68Percent2D(surfCC2), 7, kRed);

  surfCC.DrawBestFit(kBlack);
  surfCC.DrawContour(Gaussian2Sigma2D(surfCC), kSolid, kBlack);
  surfCC.DrawContour(Gaussian68Percent2D(surfCC), 7, kBlack);

  TH1F *solid  = new TH1F("90% CL", "90% CL", 10, 0, 1);
  TH1F *dashed = new TH1F("68% CL", "68% CL", 10, 0, 1);

  solid->SetLineWidth(2);
  dashed->SetLineWidth(2);
  dashed->SetLineStyle(7);
  solid->SetLineColor(1);
  dashed->SetLineColor(1);

  TLegend *leg4 = new TLegend(0.12,0.12,0.5,0.4);
  leg4->SetFillColor(10);
  leg4->SetFillStyle(0);
  leg4->SetHeader("18e20 POT FHC, 14 kton");

  leg4->AddEntry(solid, "2 #sigma","L");
  leg4->AddEntry(dashed,"1 #sigma","L");
  leg4->SetBorderSize(0);
  leg4->Draw("");
  text.DrawLatex(0.3,0.92,"NO#nuA #nu_{#mu} Sensitivity (14kton, 700kW)");

  gPad->Print("plots/basic.overlay.ps");
  gPad->Print("plots/basic.overlay.C");
}