MakeSurface.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/Fit/Fit.h"

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

#include <ctime>

using namespace ana;

void MakeSurface(TString opt) {

  std::cout << "MakeSurface: Begin." << std::endl;

  DontAddDirectory guard;

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

  std::vector<std::string> polarity = { "fhc", "rhc" };
  std::cout << "MakeSurface: Getting Polarity Option." << std::endl;
  std::vector<bool> use_polarity = GetOptionConfig(opt, polarity);
  std::vector<std::string> detector = { "neardet", "fardet" };
  std::cout << "MakeSurface: Getting Detector Option." << std::endl;
  std::vector<bool> use_detector = GetOptionConfig(opt, detector);
  std::cout << "MakeSurface: ...done" << std::endl;

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

  // Get options for a particular systematic treatment.
  // Expect:
  //  Single  cm  syst : _covmxsyst=name_
  //  Single pull syst : _pullsyst=name_
  //  All  cm  systs   : _covmxsysts_
  //  All pull systs   : _pullsyst=all_
  //  All cm, single pull swapped : _covmxsysts_pullsyst=name_

  // Single pull syst
  std::string pull_syst = "none";
  size_t pspos = optstring.find("pullsyst=");
  if (pspos != std::string::npos) {
    if (optstring.substr(pspos+8, 1) != "=")
      assert(false and "The pullsyst option must take the form \"pullsyst=xxx\".");
    size_t start = pspos+9;
    size_t end = optstring.find("_", pspos);
    if (end != std::string::npos)
      pull_syst = optstring.substr(start, end-start);
    else
      pull_syst = std::stod(optstring.substr(start));
    std::cout << "Treating syst " << pull_syst << " as a pull term." << std::endl;
  }

  // All CM systs or All Pull systs
  bool covmx_systs = CheckOption(opt, "covmxsysts");
  if (pull_syst == "all" and covmx_systs) assert (false and
    "Treating systematics as both pull terms and covariance matrix is not permitted.");

  // Single CM syst
  std::string cm_syst = "none";
  size_t cmspos = optstring.find("covmxsyst=");
  if (cmspos != std::string::npos) {
    if (optstring.substr(cmspos+9, 1) != "=")
      assert(false and "The covmxsyst option must take the form \"covmxsyst=xxx\".");
    if (covmx_systs)
      assert(false and "All covmxsysts already, \"covmxsyst=xxx\" is redundant");
    if (pull_syst!="none") 
      assert(false and "Are you sure you want to use a single CM syst with a single pull term?");
    size_t start = cmspos+10;
    size_t end = optstring.find("_", cmspos);
    if (end != std::string::npos)
      cm_syst = optstring.substr(start, end-start);
    else
      cm_syst = std::stod(optstring.substr(start));
    std::cout << "Treating syst " << cm_syst << " as a cm term." << std::endl;
  }

  // Load saved fake data, if requested
  TDirectory* fakedata_dir = nullptr;
  size_t upos = optstring.find("universe");
  if (upos != std::string::npos) {
    size_t universe;
    size_t start = upos + 8;
    size_t end = optstring.find("_", upos);
    if (end != std::string::npos)
      universe = std::stoi(optstring.substr(start, end - start));
    else
      universe = std::stoi(optstring.substr(start));

    std::cout << "Using saved fake data universe " << universe << std::endl;

    std::string fakedata_filename = FindCAFAnaDir() + "data/joint/fakedata.root";
    TFile* fakedata_file = TFile::Open(fakedata_filename.c_str(), "read");
    std::string dirname = "universe_" + std::to_string(universe);
    fakedata_dir = (TDirectory*)fakedata_file->Get(dirname.c_str());
    if (!fakedata_dir)
      assert(false and "Unable to load fake data directory.");
  }

  // 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 prof_delta24 = not CheckOption(opt, "noprofdelta24");

  // 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);
  SetAngles(calc);

  // Load predictions
  std::vector<covmx::Sample> samples;
  std::vector<IPrediction*> preds;
  std::vector<const ISyst*> isysts;
  for (size_t d = 0; d < 2; ++d) {   // for detector
    for (size_t p = 0; p < 2; ++p) { // for beam polarity
      if (!use_polarity[p] || !use_detector[d]) continue;
      samples.push_back(covmx::Sample(covmx::Analysis::kNC,
        (covmx::Polarity)p, (covmx::Detector)d));
      if(samples.size()==1) 
          isysts = GetSystsFromFile(samples.back());
      std::cout << "MakeSurface: Loading Prediction..." << std::endl;
      preds.push_back(LoadPrediction(samples.back(), pull_syst == "none" ? "nosysts" : "systs"));
      std::cout << "MakeSurface: ...Loaded Prediction." << std::endl;
    } // for beam polarity
  } // for detector

  // Load cosmics
  std::vector<Spectrum*> cosmic;
  for (auto sample : samples){
    std::cout << "MakeSurface: Loading Cosmic..." << std::endl;
    cosmic.push_back(LoadCosmic(sample).release());
    std::cout << "MakeSurface: ...Loaded Cosmic." << std::endl;
  }

  // Load data
  std::vector<Spectrum*> data;
  if (fakedata_dir) {
    for (auto sample : samples){
      std::cout << "MakeSurface: Loading Fake Data..." << std::endl;
      data.push_back(Spectrum::LoadFrom(fakedata_dir, sample.GetName().c_str())).release();
      std::cout << "MakeSurface: ...Loaded Fake Data." << std::endl;
    }
  }
  else {
    for (size_t i = 0; i < preds.size(); ++i) {
      data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(12e20)));
      if(samples[i].detector == covmx::Detector::kFarDet)
        data.back()->OverrideLivetime(kNus19SensLivetime);
      if (cosmic[i]) *data[i] += *cosmic[i];
    }
  }

  // Load systematics
  //std::vector<const ISyst*> isysts = GetSystsFromFile(samples[0]);

  // Pull term systematics
  std::vector<const ISyst*> pull_systs = {};
  SystConcat* sc = nullptr;
  if (pull_syst == "all") {
    // All pull systs, single detector
    if (samples.size() == 1)
      pull_systs = isysts;
    // All pull systs, both detectors
    else {
      sc = new SystConcat(GetSystConcat(samples, isysts, preds));
      pull_systs = sc->GetActiveSysts();
    }
  }
  else if (pull_syst != "none") {
    // Single pull syst, single detector
    if (samples.size() == 1) {
      for (const ISyst* syst : isysts) {
        if (GetSystBaseName(syst->ShortName()) == pull_syst) {
          pull_systs = { syst };
          break;
        }
      }
      if (pull_systs.empty()) {
        std::ostringstream err;
        err << "Couldn't find requested pull systematic " << pull_syst;
        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()) == pull_syst) {
          sc = new SystConcat(it_sc);
          pull_systs = sc->GetActiveSysts();
          break;
        }
      }
    }
  }

  // Load covariance matrix
  CovarianceMatrix* mx = nullptr;
  if (covmx_systs || cm_syst!="none") {
    std::vector<CovarianceMatrix*> CMs;
    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)) {
      std::string sysname = it_mx->GetName();
      // if single CM syst configured, ignore all others
      if (cm_syst!="none" && cm_syst!=sysname) continue;
      // otherwise add the syst as long as it's not a pull term
      if (GetSystBaseName(sysname) != pull_syst) {
        std::cout << "Adding " << sysname
          << " to covariance matrix." << std::endl;
        mx->AddMatrix(it_mx);
      }
      else
        std::cout << "Omitting pull term " << pull_syst
          << " from matrix." << std::endl;
    }
  }

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

  // Create experiment
  OscCovMxExperiment* expt = nullptr;
  if (covmx_systs||cm_syst!="none") expt = new OscCovMxExperiment(preds, mx, data, cosmic);
  else expt = new OscCovMxExperiment(preds, data, cosmic);

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

  // Set fixed dm41 for th24 vs th34 case
  if (th24vsth34){
    calc->SetDm(4, dm41);
    std::cout << "\n\n Set the dm41 \n\n" << std::endl;
  }

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

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

  std::vector<const IFitVar*> xy_vars;

  if (th24vsth34) {
    xbins = new Binning(Binning::Simple(50, 0, 50));
    ybins = new Binning(Binning::Simple(50, 0, 45));
    xy_vars.push_back(&kFitTheta34InDegreesSterile);
    xy_vars.push_back(&kFitTheta24InDegreesSterile);
  }
  else if (th24vsdm41 || th34vsdm41) {
    xbins = new Binning(Binning::LogUniform(50, 1e-4, 1));
    ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
    if (th24vsdm41) {
      xy_vars.push_back(&kFitSinSqTheta24Sterile);
      xy_vars.push_back(&kFitDmSq41Sterile);
      fit_vars.push_back(&kFitSinSqTheta34Sterile);
      seed_values[&kFitSinSqTheta34Sterile] = { 0.01, 0.5, 0.99 };
    }
    if (th34vsdm41) {
      xy_vars.push_back(&kFitSinSqTheta34Sterile);
      xy_vars.push_back(&kFitDmSq41Sterile);
      fit_vars.push_back(&kFitSinSqTheta24Sterile);
      seed_values[&kFitSinSqTheta24Sterile] = { 0.01, 0.5, 0.99 };
    }
  }

  // Create surface
  size_t part = 1275;
  size_t njobs = 2500;
  if (RunningOnGrid()){
    part=JobNumber();
    njobs=25;
  }

  CovMxSurface surf(*xbins, *ybins, &multi_expt, calc, xy_vars[0], xy_vars[1],
    fit_vars, seed_values, pull_systs, 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, optstring.c_str());
  delete f;

}