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 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/LikelihoodCovMxExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Analysis/CovMxManager.h"
#include "CAFAna/Analysis/Exposures.h"

#include "NuXAna/Core/Utilities.h"
#include "NuXAna/Systs/NusAna2020Systs.h"
#include "NuXAna/macros/Nus20/Utilities.h"
#include "NuXAna/macros/Nus20/Fitting/FitUtilities.h"

#include <chrono>

using namespace ana;

using std::cout;
using std::endl;
using std::ostringstream;
using std::string;
using std::vector;

using std::chrono::high_resolution_clock;
using std::chrono::duration_cast;
using std::chrono::seconds;

using covmx::Sample;
using covmx::CovarianceMatrix;

void MakeSurface(TString opt, string systConfig="allsysts") {

  DontAddDirectory guard;

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

  FitConfig cfg = GetConfig(opt);
  if (cfg.samples.size() == 0) {
    throw std::runtime_error("No samples found! Exiting.");
  }

  for (Sample& s : cfg.samples) {
    cout << "Setting prediction for sample " << s.GetName() << "..." << endl;
    bool usePredInterp = cfg.fakeDataType == kPseudoexperiment && cfg.systType == kCovMxSysts;
    SetPrediction(s, usePredInterp);
    SetCosmic(s);
    if (usePredInterp) SetSystAlias(s);
  }

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

  cout << cfg.xAxis->var->ShortName() << "   " << cfg.yAxis->var->ShortName() << endl;

  // Oscillation calculator
  osc::IOscCalcSterile* calc = DefaultSterileCalc(4);
  SetNus20Params(calc, cfg.nullHyp);
  PrintOscParams(calc);

  // Load data
  if (cfg.fakeDataType == kAsimov) {
    cout << "Fitting with Asimov fake data with \"" << cfg.nullHyp << "\" null hypothesis:" << endl;
    for (Sample& sample: cfg.samples) {
      Spectrum* data = new Spectrum(sample.GetPrediction()->Predict(calc).FakeData(sample.GetPOT()));
      if (sample.detector == covmx::kFarDet) data->OverrideLivetime(sample.GetLivetime());
      if (sample.HasCosmic()) *data += *sample.GetCosmic();
      sample.SetData(data);
    }
  } else if (cfg.fakeDataType == kPseudoexperiment) {
    cout << "Fitting with fluctuated fake data \"" << cfg.nullHyp
      << "\" null hypothesis, pseudoexperiment " << cfg.fakeDataUniverse << endl;
    GeneratePseudoexpt(cfg.samples, calc, cfg.fakeDataUniverse, (cfg.systType == kCovMxSysts));
  }
  if (CheckOption(opt, "th24vsth34")) calc->SetDm(4, 0.5);

  // Load covariance matrix
  CovarianceMatrix* mx = nullptr;
  if (cfg.systType == kCovMxSysts) {
    TFile* f = TFile::Open((InputPath()+"/covmx/mx.root").c_str(), "read");
    ostringstream tag;
    for (size_t i = 0; i < cfg.samples.size(); ++i) {
      if (i > 0) tag << "_";
      tag << cfg.samples[i].GetID();
    }
    mx = LoadFrom<CovarianceMatrix>(f, tag.str()+"/"+systConfig).release();
    delete f;
  }

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

  // Create experiment
  IExperiment* expt(nullptr);
  if (cfg.fitType == kPoissonLL) {
    LikelihoodCovMxExperiment* ll = new LikelihoodCovMxExperiment(cfg.samples, mx, cfg.epsilon, cfg.lambda);
    expt = ll;
    if (!opt.Contains("quiet")) ll->SetVerbose(true);
  } else {
    expt = new OscCovMxExperiment(cfg.samples, mx);
    if (cfg.fitType == kCNP) dynamic_cast<OscCovMxExperiment*>(expt)->SetCNP();
  }

  // Set Gaussian constraints, make MultiExperiment
  vector<const IExperiment*> expts = { expt };
  if (!opt.Contains("numusel")) {
    cout << "No numu samples included! Adding Gaussian constraints on atmospheric parameters." << endl;
    for (const IExperiment* e : GetConstraints()) expts.push_back(e);
  } else {
    cout << "Numu samples included! Adding loose constraint on dmsq32." << endl;
    expts.push_back(GetConstraintLoose());
  }
  MultiExperiment multiExpt(expts);
  TH1D hTime("hTime", "", 1, 0, 1);
  auto start = high_resolution_clock::now();
  FrequentistSurface surf(&multiExpt, calc, *cfg.xAxis, *cfg.yAxis,
    cfg.fitVars, {}, SeedList(cfg.fitSeeds));
  auto end = high_resolution_clock::now();
  hTime.SetBinContent(1, duration_cast<seconds>(end-start).count());

  // Open output file and save surface
  string optstring = opt.Data();
  TFile* outFile = new TFile((optstring+".root").c_str(), "recreate");
  surf.SaveTo(outFile, "surface");
  outFile->WriteTObject(&hTime, "time");
  delete outFile;

} // main function