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/Fit/CovMxSurface.h"
#include "CAFAna/Analysis/CovMxManager.h"
#include "CAFAna/Analysis/Exposures.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::string;
using std::vector;
using std::ostringstream;
using std::cout;
using std::endl;
using std::exception;

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") {

  try {

    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;
      SetPrediction(s, false);
    }

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

    // Oscillation calculator
    osc::OscCalcSterile* calc = DefaultSterileCalc(4);
    SetNus20Params(calc, cfg.nullHyp);
    if (CheckOption(opt, "th24vsth34")) calc->SetDm(4, 0.5);

    // Load cosmics
    for (Sample& sample : cfg.samples) SetCosmic(sample);

    cout << cfg.xAxis->var->ShortName() << "   " << cfg.yAxis->var->ShortName() << endl;
 
    // Load data
    if (cfg.fakeDataType == kAsimov) {
      cout << "Fitting with Asimov fake data with \"" << cfg.nullHyp << "\" null hypothesis:" << endl;
      PrintOscParams(calc);
      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);
      }
      SetNus20Params(calc, cfg.nullHyp);
      if (CheckOption(opt, "th24vsth34")) calc->SetDm(4, 0.5);
    } else if (cfg.fakeDataType == kPseudoexperiment) {
      cout << "Loading fake data, type " << cfg.fakeDataSet << ", pseudoexperiment "
        << cfg.fakeDataUniverse << endl;
      ostringstream filepath;
      filepath << InputPath() << "/fakedata/fakedata_" << cfg.fakeDataSet << ".root";
      TFile* f = TFile::Open(filepath.str().c_str(), "read");
      // Get universe and spectra
      TDirectory* dir = f->GetDirectory(Form("universe_%u", cfg.fakeDataUniverse));
      for (Sample& s : cfg.samples) {
        s.SetData(Spectrum::LoadFrom(dir, s.GetTag()).release());
      }
      delete f; // Close fake data file
    }

    // 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;
      ll->SaveHists();
      if (!opt.Contains("quiet")) ll->SetVerbose(true);
      ll->SetResetShifts(false);
    } 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! Omitting Gaussian constraints on atmospheric parameters." << endl;
    }
    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;

  } catch(const exception& e) {
    cout << endl << "Exception occurred! " << e.what() << endl;
  }

} // main function