MakeSurfaceNoNDOsc.C

Go to the documentation of this file.

/*
 * MakeSurface.C
 *
 *  Created on: May 7, 2018
 *  
 *      Author: J. Hewes <jhewes15@fnal.gov>
 *
 *  Create contours for 2019 covmx analysis
 */

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Sample.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Experiment/OscCovMxExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Analysis/CovMxSurface.h"
#include "CAFAna/Analysis/CovMxManager.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Fit.h"

#include "NuXAna/macros/Nus20/Utilities.h"

#include <ctime>

using namespace ana;

void MakeSurfaceNoNDOsc(TString opt) {

  DontAddDirectory guard;

  // Define exposures
  double potND = kAna2018SensitivityFHCNDPOT;
  double potFD = kAna2018FHCPOT;
  double livetimeFD = kAna2018FHCLivetime;

  /////////////////////////////////////////////////////////////////////
  //                                                                 //
  //                      ARGUMENT PARSING                           //
  //                                                                 //
  /////////////////////////////////////////////////////////////////////

  std::vector<std::string> polarity = { "fhc", "rhc" };
  std::vector<bool> usePolarity = GetOptionConfig(opt, polarity);
  std::vector<std::string> detector = { "neardet", "fardet" };
  std::vector<bool> useDetector = GetOptionConfig(opt, detector);

  std::string optstring(opt.Data());

  std::vector<const ISyst*> isysts;
  if (opt.Contains("covmxsysts") or opt.Contains("pullsyst")) {
    if (useDetector[1]) isysts = LoadSysts(kNusFHCFarDet, kUseCVMFS);
    if (useDetector[0]) isysts = LoadSysts(kNusFHCNearDet, kUseCVMFS);
  }

  // Get pull and covariance matrix syst options
  bool covmxSysts = CheckOption(opt, "covmxsysts");
  std::string pullSyst = "none";
  size_t pspos = optstring.find("pullsyst=");
  if (pspos != std::string::npos) {
    if (optstring.substr(pspos+8, 1) != "=")
      assert(false and "The dm41 option must take the form \"dm41=xxx\".");
    size_t start = pspos+9;
    size_t end = optstring.find("_", pspos);
    pullSyst = optstring.substr(start, end-start);
    std::cout << "Treating syst " << pullSyst << " as a pull term." << std::endl;
  }
  if (pullSyst == "all" and covmxSysts) assert (false and
    "Treating systematics as both pull terms and covariance matrix is not permitted.");

  // Load saved fake data, if requested
  std::string fakedata;
  size_t universe = 999999;
  size_t fdpos = optstring.find("fakedata=");
  if (optstring.substr(fdpos+8, 1) != "=")
    assert(false and "The data option must take the form \"fakedata=xxx\".");
  if (fdpos != std::string::npos) {
    size_t start = fdpos + 9;
    if (optstring.substr(start, 6) == "asimov") {
      fakedata = "asimov";
    }
    else if (optstring.substr(start, 5) == "stats") {
      fakedata = "stats";
      size_t end = optstring.find("_", fdpos);
      universe = std::stoi(optstring.substr(start+5, end-(start+5)));
    }
    else if (optstring.substr(start, 5) == "systs") {
      fakedata = "systs";
      size_t end = optstring.find("_", fdpos);
      universe = std::stoi(optstring.substr(start+5, end-(start+5)));
    }
    else {
      assert(false and "Fake data type not recognised!");
    }
  }

  // Set delta m, if requested
  double dm41 = -1;
  size_t dmpos = optstring.find("dm41=");
  if (dmpos != std::string::npos) {
    if (optstring.substr(dmpos+4, 1) != "=")
      throw std::runtime_error("The dm41 option must take the form \"dm41=xxx\".");
    size_t start = dmpos+5;
    size_t end = optstring.find("_", dmpos);
    if (end != std::string::npos)
      dm41 = std::stod(optstring.substr(start, end-start));
    else
      dm41 = std::stod(optstring.substr(start));
    std::cout << "Using dm41 value " << dm41 << std::endl;
  }

  // Don't profile delta24
  bool profDelta24 = not CheckOption(opt, "noprofdelta24");
  bool profVars = not CheckOption(opt, "noprofvars");

  // Check what type of plot we're making
  bool th24vsth34 = false;
  bool th24vsdm41 = false;
  bool th34vsdm41 = false;
  if (CheckOption(opt, "th24vsth34")) th24vsth34 = true;
  else if (CheckOption(opt, "th24vsdm41")) th24vsdm41 = true;
  else if (CheckOption(opt, "th34vsdm41")) th34vsdm41 = true;
  else throw std::runtime_error("The options string must contain plot type: \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");

  // Make sure dm41 is set for th24 vs th34 case
  if (th24vsth34 && dm41 == -1)
    throw std::runtime_error("If you want to run a th24 vs th34 surface, you must specify a value for dm41.");

  /////////////////////////////////////////////////////////////////////
  //                                                                 //
  //                        PREPARE INPUTS                           //
  //                                                                 //
  /////////////////////////////////////////////////////////////////////

  // Oscillation calculator
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  SetNus20Params(calc);

  // Load predictions
  std::vector<covmx::Sample> samples;
  std::vector<IPrediction*> preds;
  for (size_t d = 0; d < 2; ++d) {   // for detector
    for (size_t p = 0; p < 2; ++p) { // for beam polarity
      if (!usePolarity[p] || !useDetector[d]) continue;
      samples.push_back(covmx::Sample(covmx::Analysis::kNC,
        (covmx::Polarity)p, (covmx::Detector)d));
      preds.push_back(LoadPrediction(samples.back(), pullSyst == "none" ? "nosysts" : "systs", kUseCVMFS));
    } // for beam polarity
  } // for detector

  // Load cosmics
  std::vector<Spectrum*> cosmic;
  for (auto sample : samples)
    cosmic.push_back(LoadCosmic(sample, kUseCVMFS).release());

  // Load data
  std::vector<Spectrum*> data;
  if (fakedata == "asimov") {
    std::cout << "Using Asimov fake data." << std::endl;
    for (size_t i = 0; i < preds.size(); ++i) {
      double pot = (samples[i].detector == covmx::Detector::kNearDet) ? potND : potFDl
      data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(pot)));
      if(samples[i].detector == covmx::Detector::kFarDet)
        data.back()->OverrideLivetime(livetimeFD);
      if (cosmic[i]) *data[i] += *cosmic[i];
    }
  }
  else {
    std::cout << "Loading universe " << universe << " of " << fakedata
      << " fake data." << std::endl;
    data = LoadFakeData(samples, universe, (fakedata=="systs"), kUseCVMFS);
  }

  // Pull term systematics
  std::vector<const ISyst*> pullSysts = {};
  SystConcat* sc = nullptr;
  if (pullSyst == "all") {
    // All pull systs, single detector
    if (samples.size() == 1)
      pullSysts = isysts;
    // All pull systs, both detectors
    else {
      sc = new SystConcat(GetSystConcat(samples, isysts, preds));
      pullSysts = sc->GetActiveSysts();
    }
  }
  else if (pullSyst != "none") {
    // Single pull syst, single detector
    if (samples.size() == 1) {
      for (const ISyst* syst : isysts) {
        if (GetSystBaseName(syst->ShortName()) == pullSyst) {
          pullSysts = { syst };
          break;
        }
      }
      if (pullSysts.empty()) {
        std::ostringstream err;
        err << "Couldn't find requested pull systematic " << pullSyst;
        assert (false and err.str().c_str());
      }
    }
    // Single pull syst, both detectors
    else {
      std::vector<SystConcat> all_scs = GetSystConcats(samples, preds, isysts);
      for (auto it_sc : all_scs) {
        if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pullSyst) {
          sc = new SystConcat(it_sc);
          pullSysts = sc->GetActiveSysts();
          break;
        }
      }
    }
  }

  // Load covariance matrix
  CovarianceMatrix* mx = nullptr;
  if (covmxSysts) {
    std::vector<Binning> bins;
    for (auto d : data) bins.push_back(d->GetBinnings()[0]);
    mx = new CovarianceMatrix(bins);
    std::cout << std::endl;
    for (CovarianceMatrix* it_mx : GetMatrices(samples, kUseCVMFS)) {
      if (GetSystBaseName(it_mx->GetName()) != pullSyst) {
        mx->AddMatrix(it_mx);
      }
      else
        std::cout << "Omitting pull term " << pullSyst
          << " from matrix." << std::endl;
    }
  }

  /////////////////////////////////////////////////////////////////////
  //                                                                 //
  //                          MAKE SURFACE                           //
  //                                                                 //
  /////////////////////////////////////////////////////////////////////

  // Create experiment
  OscCovMxExperiment* expt = nullptr;
  if (covmxSysts) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
  else expt = new OscCovMxExperiment(preds, data, cosmic);
  if (samples.size() > 1 and sc) expt->ConfigureSysts(sc);
  auto noosccalc = DefaultSterileCalc(4);
  for (int i = 1; i <= 4; ++i) {
    noosccalc->SetDm(i, 0);
  }
  expt->DisableOscillations({true, false}, noosccalc);

  // Set Gaussian constraints, make MultiExperiment
  std::vector<const IChiSqExperiment*> expts = GetConstraints();
  expts.push_back(expt);
  MultiExperiment multiExpt(expts);

  // Set fixed dm41 for th24 vs th34 case
  if (th24vsth34) calc->SetDm(4, dm41);

  // Binning, fit vars and seeds
  std::vector<const IFitVar*> fitVars = { &kFitTheta23Sterile,
    &kFitDmSq32Sterile };
  if (profDelta24) fitVars.push_back(&kFitDelta24InPiUnitsSterile);
  else calc->SetDelta(2, 4, 0.5);
  std::map<const IFitVar*, std::vector<double> > seedValues;
  seedValues[&kFitTheta23Sterile] = { kNus19Th23 };
  seedValues[&kFitDmSq32Sterile] = { kNus19Dm32 };
  if (profDelta24) seedValues[&kFitDelta24InPiUnitsSterile] = { 0.5 };

  Binning* xbins = nullptr;
  Binning* ybins = nullptr;

  std::vector<const IFitVar*> xyVars;

  if (th24vsth34) {
    xbins = new Binning(Binning::Simple(50, 0, 50));
    ybins = new Binning(Binning::Simple(50, 0, 45));
    xyVars.push_back(&kFitTheta34InDegreesSterile);
    xyVars.push_back(&kFitTheta24InDegreesSterile);
  }
  else if (th24vsdm41) {
    xbins = new Binning(Binning::LogUniform(50, 1e-4, 1));
    ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
    xyVars.push_back(&kFitSinSqTheta24Sterile);
    xyVars.push_back(&kFitDmSq41Sterile);
    fitVars.push_back(&kFitSinSqTheta34Sterile);
    seedValues[&kFitSinSqTheta34Sterile] = { 1. };
  }
  else if (th34vsdm41) {
    xbins = new Binning(Binning::Simple(50, 0, 1));
    ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
    xyVars.push_back(&kFitSinSqTheta34Sterile);
    xyVars.push_back(&kFitDmSq41Sterile);
    fitVars.push_back(&kFitSinSqTheta24Sterile);
    seedValues[&kFitSinSqTheta24Sterile] = { 1. };
  }

  // Create surface
  size_t part = 127;
  size_t njobs = 250;
  if (RunningOnGrid()){
    part=JobNumber();
    njobs=10;
  }

  // Override fit variables if requested
  if (not profVars) {
    fitVars.clear();
    seedValues.clear();
  }

  CovMxSurface surf(*xbins, *ybins, &multiExpt, calc, xyVars[0], xyVars[1],
    fitVars, seedValues, pullSysts, SystShifts::Nominal(), true, nullptr,
    preds, data, cosmic);//, part, njobs);

  // Open output file and save surface
  std::string outfile = optstring + ".root";
  TFile* f = new TFile(outfile.c_str(), "recreate");
  surf.SaveTo(f->mkdir(optstring.c_str()));
  delete f;

}