make_prediction_noextrap.C

Go to the documentation of this file.

#pragma once

#include "3FlavorAna/Ana2018/numu/includes.h"
#include "3FlavorAna/Ana2018/numu/set_loaders_ana2018.h"

using namespace ana;


void make_prediction_noextrap(std::string polarity = "rhc",
                              std::string period   = "full",
                              bool nh = true,
                              bool systematics = false,
                              bool calc3a = true)
{

  std::cout << "\n=================================" << std::endl;
  std::cout << "    make predictions no extrap     " << std::endl;
  std::cout << "=================================\n" << std::endl;

  std::string extrapName = "noextrap";


  // production files
  std::string production = "";
  if(polarity == "fhc"){
    production += ".d";
  }
  if(polarity == "rhc"){
    production += ".e";
  }

  // systematics
  std::string systName = "";
  auto systs = getAllAna2018Systs(kNumuAna2018);
  if(systematics){
    systName = "systs";
  }
  if(!systematics){
    systName = "stats";
    systs = {};
  }
  for (auto & sys:systs) std::cout << sys->ShortName() << std::endl;

  // calculator
  std::string calcName = "2018calc";
  std::string hierarchy = "";
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();
  if(nh){
    hierarchy = "nh";
    calc->SetTh23(asin(sqrt(0.559455)));
    calc->SetDmsq32(0.00244);
  }
  if(!nh){
    hierarchy = "ih";
    calc->SetTh23(asin(sqrt(0.559455)));
    calc->SetDmsq32(-0.00244);
  }


  // selection and weight
  const Cut kNumuNDCut = kNumuCutND2018;
  const Cut kNumuFDCut = kNumuCutFD2018;
  const Var kXSecWeight = kXSecCVWgt2018;

  std::string OutDir = "";
  std::cout << "polarity: " << polarity << ", period: " << period << ", files: " << production << std::endl;
  std::cout << "systematics: " << systName << " , calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;


  // quantiles
  std::string InDirQuant = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles";
  std::string InFileQuant = InDirQuant + "/quantiles__" + polarity + "_full__numu2018.root";
  std::cout << "\n\n --- get quantile cuts from file ---" << std::endl;
  TFile* inFile = TFile::Open( pnfs2xrootd(InFileQuant).c_str() );
  TH2* FDSpec2D = (TH2*)inFile->FindObjectAny("FDSpec2D");
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
  inFile->Close(); 
  delete inFile;


  //////////////////// make loader for prediction ////////////////////                                                                                                             
  std::cout << "loaders" << std::endl;
  Loaders loaders;
  loaders.SetSpillCut(kStandardSpillCuts);

  std::cout << "set loaders" << std::endl;
  set_loaders_ana2018(loaders, polarity, period, production);

  std::cout << "declare prediction" << std::endl;
  std::vector<NoExtrapGenerator*> NoExtrapVec;
  std::vector<PredictionInterp*> predictionVec;

  std::cout << "loop over quantiles for prediction" << std::endl;
  for(auto thisCut : HadEFracQuantCuts){
    NoExtrapVec.push_back(new NoExtrapGenerator(kNumuCCOptimisedAxis, thisCut && kNumuFDCut, kPPFXFluxCVWgt*kXSecWeight));
    predictionVec.push_back(new PredictionInterp(systs, calc, *(NoExtrapVec.back()), loaders));
  }

  std::cout << "loaders.Go()" << std::endl;
  loaders.Go();
  
  
  std::cerr << "Saving Predictions" << std::endl;
  for(size_t quant=1; quant <= predictionVec.size();quant++){
    auto prediction = predictionVec[quant-1];
    std::string OutFileName = Form("prediction_quant%d__", quant) + extrapName + "_" + systName + "__" + polarity + "_" + calcName + "_" + hierarchy + "__" + period + "__numu2018.root";
    std::cerr << "Saving Quantile " << quant << " to Prediction file: " << OutFileName << std::endl;
    TFile *f = TFile::Open(OutFileName.c_str(), "RECREATE");
    prediction -> SaveTo(f, "prediction");
    f->Close();
  }


}//make_predictions_noextrap