caf_numu_sensitivity.C

Go to the documentation of this file.

#ifdef __CINT__
void caf_numu_sensitivity(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName){
  std::cout << "Sorry, you must run in compiled mode" << std::endl;
}
#else

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Fit.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Surface.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Vars/FitVarsNumu.h"
#include "CAFAna/Vars/Vars.h"
#include "OscLib/func/OscCalculatorPMNSOpt.h"

using namespace ana;

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TGraph.h"
#include "TList.h"
#include "TObject.h"
#include "TStyle.h"

const double pot = 6.e20;

TH1D* GetXMarginalisation(TH2* h2);
TH1D* GetYMarginalisation(TH2* h2);

const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
                         [](const caf::StandardRecord* sr)
                         {
                           if(sr->mc.nnu == 0) return false;
                           assert(sr->mc.nnu == 1);
                           return (sr->mc.nu[0].mode == caf::kMEC);
                         });

const Var kSlcMaxX({"slc.boxmax.fX"},
      [](const caf::StandardRecord *sr)
      {
        return sr->slc.boxmax.X()/100;
      });
const Var kSlcMaxY({"slc.boxmax.fY"},
      [](const caf::StandardRecord *sr)
      {
        return sr->slc.boxmax.Y()/100;
      });
const Var kSlcMaxZ({"slc.boxmax.fZ"},
      [](const caf::StandardRecord *sr)
      {
        return sr->slc.boxmax.Z()/100;
      });

void caf_numu_sensitivity(std::string kInputNonSwapFileName, std::string kInputSwapFileName, std::string kInputTauFileName, std::string kOutputFileName)
{
  std::cout << "\nrun  : ---- Running CAF NuMu sensitivity.\n";
  std::cout << "run  :     input non-swap file name:    " << kInputNonSwapFileName << std::endl;
  std::cout << "run  :         input swap file name:    " << kInputSwapFileName    << std::endl;
  std::cout << "run  :          input tau file name:    " << kInputTauFileName     << std::endl;
  std::cout << "run  :             output file name:    " << kOutputFileName       << std::endl;


  osc::OscCalculatorPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc); // A calculator set to default values

  SpectrumLoader loaderNonswap(kInputNonSwapFileName);
  SpectrumLoader loaderSwap(kInputSwapFileName);
  SpectrumLoader loaderTau(kInputTauFileName);

  loaderNonswap.SetSpillCut(kStandardSpillCuts);
  loaderSwap.SetSpillCut(kStandardSpillCuts);
  loaderTau.SetSpillCut(kStandardSpillCuts);

  // Taus are not required, you can use kNullLoader for fast tests or if they're not available
  PredictionNoExtrap* prediction = new PredictionNoExtrap(loaderNonswap, loaderSwap, loaderTau, 
      "Reconstructed neutrino Energy [GeV]", Binning::Simple(20,0,5), kCCE, kNumuFD && (!kIsDytmanMEC));

  std::cout << "\nrun  : --- run loaders.\n";
  loaderNonswap.Go();
  loaderSwap.Go();
  loaderTau.Go();
  std::cout << "\nrun  : --- done.\n";

  std::cout << "\nrun  : --- Make surfaces.\n";
  // Create fake prediction and data, normalised to the chosen POT
  Spectrum fakePred = prediction->Predict(&calc);
  Spectrum fakeData = fakePred.FakeData(pot);

  // Experiment, i.e., create set of prediction vs data
  SingleSampleExperiment fakeExpt(prediction, fakeData);

  // Fill the chi-2 (LL) 2-D histogram
  Surface* sensitivity = new Surface( &fakeExpt, &calc,
  &kFitSinSqTheta23, 40, 0.3, 0.7,
  &kFitDmSq32Scaled, 50, 2.0, 3.0 );

  std::cout << "\nrun  : --- save output.\n";
  TFile* f = new TFile(kOutputFileName.c_str(),"RECREATE");

  // Save the surface
  TH2* surf = sensitivity->ToTH2();
  surf->SetName("reco-sensitivity_surface-numu_selection");
  surf->Write();

  // Save the best fit
  TGraph* bestFit = new TGraph;
  bestFit->SetPoint(0, sensitivity->fMinX, sensitivity->fMinY);
  bestFit->SetName("reco-sensitivity_bestfit-numu_selection");
  bestFit->GetXaxis()->SetLimits(0.3, 0.7);
  bestFit->GetYaxis()->SetLimits(2.0, 3.0);
  bestFit->Write();

  // Save 1D marginalisations
  TH1D* margX = GetXMarginalisation(surf);
  TH1D* margY = GetYMarginalisation(surf);
  margX->SetName("reco-sensitivity_1DmarginX-numu_selection");
  margY->SetName("reco-sensitivity_1DmarginY-numu_selection");
  margX->Write();
  margY->Write();

  // Now different contours
  TCanvas* c = new TCanvas("c","c");
  c->SetLeftMargin(0.12);
  c->SetBottomMargin(0.15);

  TList *ObjList = c->GetListOfPrimitives();
  TObject* obj;
  TIter next(ObjList);

  TGraph* g68 = new TGraph;
  TGraph* g90 = new TGraph;

  // Draw the 68% contour
  sensitivity->DrawContour( Gaussian68Percent2D(*sensitivity), kDashed, kRed);

  while(obj=(TObject *)next())
  {
    if(!obj->InheritsFrom(TGraph::Class())) continue;
    g68 = (TGraph*) obj->Clone("reco-sensitivity_68pc_contour-numu_selection");
  }

  g68->Write();

  // Draw the 90% contour
  c->Clear();
  sensitivity->DrawContour( Gaussian90Percent2D(*sensitivity), kSolid, kRed);

  next.Begin();

  while(obj=(TObject *)next())
  {
    if(!obj->InheritsFrom(TGraph::Class())) continue;
    g90 = (TGraph*) obj->Clone("reco-sensitivity_90pc_contour-numu_selection");
  }

  g90->Write();

  TH1D* TotalPOT = new TH1D("meta-TotalPOT","TotalPOT", 1, 0, 1);
  TotalPOT->SetBinContent(1, pot);
  TotalPOT->Write();

  f->Write();
  f->Close();

  std::cout << "All plots finished!" << std::endl;

} // End of function

TH1D* GetXMarginalisation(TH2* h2)
{
  int nY = h2->GetYaxis()->GetNbins();
  int nX = h2->GetXaxis()->GetNbins();

  double Xmin = h2->GetXaxis()->GetXmin();
  double Xmax = h2->GetXaxis()->GetXmax();
  double Ymin = h2->GetYaxis()->GetXmin();
  double Ymax = h2->GetYaxis()->GetXmax();

  TH1D *h1Dim = new TH1D("","",nX,Xmin,Xmax);

  for(int binx = 1; binx <= nX; binx++) {
    double minchi = 999999.;
    for(int biny = 1; biny <= nY; biny++) {
      double val = h2->GetBinContent(binx,biny);
      if(val < minchi) minchi = val;
    }
    h1Dim->SetBinContent(binx,minchi);
  }

  return h1Dim;
}

TH1D* GetYMarginalisation(TH2* h2)
{
  int nY = h2->GetYaxis()->GetNbins();
  int nX = h2->GetXaxis()->GetNbins();

  double Xmin = h2->GetXaxis()->GetXmin();
  double Xmax = h2->GetXaxis()->GetXmax();
  double Ymin = h2->GetYaxis()->GetXmin();
  double Ymax = h2->GetYaxis()->GetXmax();

  TH1D *h1Dim = new TH1D("","",nY,Ymin,Ymax);

  for(int biny = 1; biny <= nY; biny++) {
    double minchi = 999999.;
    for(int binx = 1; binx <= nX; binx++) {
      double val = h2->GetBinContent(binx,biny);
      if(val < minchi) minchi = val;
    }
    h1Dim->SetBinContent(biny,minchi);
  }

  return h1Dim;
}

#endif