MakeExtrapSurface.C

Go to the documentation of this file.

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

// framework
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Core/IFitVar.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "NuXAna/Systs/Nus18Systs.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "NuXAna/macros/Nus18/Nus18Utilities.h"
#include "OscLib/OscCalcSterile.h"

// root
#include "TFile.h"
#include "TH1.h"

// stl
#include <sstream>

using namespace ana;

void MakeExtrapSurface(TString opt) {

  TH1::AddDirectory(0);

  // Get running options
  bool fhc = opt.Contains("FHC", TString::kIgnoreCase);
  bool rhc = opt.Contains("RHC", TString::kIgnoreCase);
  assert(fhc || rhc);

  bool oned = opt.Contains("1D", TString::kIgnoreCase);
  bool twod = opt.Contains("2D", TString::kIgnoreCase);
  assert(oned || twod);

  bool systematics = opt.Contains("syst", TString::kIgnoreCase);
  bool cos = opt.Contains("cos", TString::kIgnoreCase);
  bool data = opt.Contains("data", TString::kIgnoreCase);

  bool scale = opt.Contains("scale", TString::kIgnoreCase);

  // Get the prediction
  IPrediction* pred = LoadExtrapPrediction(opt);
  Spectrum cosmic = LoadCosmicSpectrum(rhc);

  // Get POT & livetimes
  double pot, livetime;
  if (scale) {
    double scale_pot = opt.Contains("125")? 12.5 : opt.Contains("12")? 12 : opt.Contains("35")? 35 : -1;
    assert(scale_pot > 0);
    pot = scale_pot * 1e20;
    livetime = rhc? pot*(kAna2018RHCLivetime/kAna2018RHCPOT) : pot*(kAna2018FHCLivetime/kAna2018FHCPOT);
  }
  else {
    pot = rhc? kAna2018RHCPOT : kAna2018FHCPOT;
    livetime = rhc? kAna2018RHCLivetime : kAna2018FHCLivetime;
  }

  std::cout << "POT is " << pot << ", livetime is " << livetime << std::endl;

  // Get the systs
  std::vector<const ISyst*> systs = {};
  if (systematics)
    systs = LoadISysts();
    //systs = rhc? GetNus18RHCExtrapSysts() : GetNus18FHCFDSysts("all");

  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  calc->SetNFlavors(4);
  SetAngles(calc);

  // Get the far detector data
  // -- make sure to scale appropriately
  Spectrum prediction = pred->Predict(calc);
  cosmic.Scale(livetime/cosmic.Livetime());
  cosmic.OverridePOT(pot);
  Spectrum cos_spec = cosmic.FakeData(pot);
  Spectrum fddata_spec = data ? LoadData(0, rhc) : (prediction.FakeData(pot) + cos_spec);

  // Create experiment objects
  // Gaussian constraint on th23, from NOvA Jan 2018 results (joint numu+nue)
  GaussianConstraint th23Constraint(&kFitSinSqTheta23Sterile, 0.558, 0.042);

  // Make experiment
  SingleSampleExperiment experiment_ss(pred, fddata_spec, cos_spec);
  MultiExperiment experiment({&th23Constraint, &experiment_ss});

  // Variables to profile over
  std::vector<const IFitVar*> prof_vars = {&kFitSinSqTheta23Sterile, &kFitDelta24InPiUnitsSterile};
  std::vector<const IFitVar*> prof_vars_th24 = prof_vars, prof_vars_th34 = prof_vars;
  prof_vars_th24.push_back(&kFitTheta34InDegreesSterile);
  prof_vars_th34.push_back(&kFitTheta24InDegreesSterile);

  // Binnings
  int nbins34_1d = 240, nbins34_2d = 60;
  float min34 = 0.;
  float max34 = 60.;
  float nbins24_1d = 20, nbins24_2d = 40; //160
  float min24 = 0.;
  float max24 = 40.;

  // 1D profiles
  TH1* ProfileTh24, *ProfileTh34;
  if (oned) {
    SetAngles(calc);
    ProfileTh24 = Profile(&experiment, calc,
                          &kFitTheta24InDegreesSterile, nbins24_1d, min24, max24,
                          -1, prof_vars_th24, systs);
    SetAngles(calc);
    ProfileTh34 = Profile(&experiment, calc,
                          &kFitTheta34InDegreesSterile, nbins34_1d, min34, max34,
                          -1, prof_vars_th34, systs);
    std::cout << "Profiles made!  Mean is " << ProfileTh24->GetMean() << std::endl;
    std::cout << "Profiles made!  NEntries is " << ProfileTh24->GetEntries() << std::endl;
  }

  // 2D surface
  FrequentistSurface* SurfaceTh34Th24;
  if (twod) {
    SetAngles(calc);
    std::map<const IFitVar*, std::vector<double> > seedPoints;
    //seedPoints[&kFitDelta24InPiUnitsSterile] = {TMath::Pi()/2};
    SurfaceTh34Th24 = new FrequentistSurface(&experiment, calc,
                                  &kFitTheta34InDegreesSterile, nbins34_2d, min34, max34,
                                  &kFitTheta24InDegreesSterile, nbins24_2d, min24, max24,
                                  prof_vars, systs, seedPoints);
  }

  // Open files for saving analysis output
  TFile* outFile = new TFile("ExtrapSurface.root", "RECREATE");
  if (oned) {
    ProfileTh24->Write("ProfileTh24");
    ProfileTh34->Write("ProfileTh34");
  }
  if (twod)
    SurfaceTh34Th24->SaveTo(outFile, "SurfaceTh34Th24");
  outFile->Close();
  delete outFile;

}