MakeSurfaceBinningStudy.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/CovMxSurface.h"
#include "CAFAna/Analysis/CovMxManager.h"
#include "CAFAna/Analysis/Exposures.h"

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

#include <ctime>

using namespace ana;

using std::back_inserter;
using std::begin;
using std::endl;
using std::end;
using std::endl;
using std::exception;
using std::ostringstream;
using std::runtime_error;
using std::string;
using std::transform;
using std::vector;

using osc::OscCalcSterile;

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

void MakeSurfaceBinningStudy(TString opt, size_t ndPred, size_t fdPred) {

  try {

    DontAddDirectory guard;

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

    FitConfig cfg = GetConfig(opt);

    if (ndPred < 1 || ndPred > 6) {
      ostringstream err;
      err << "Near detector binning option must be an integer between 1 and 6!"
          << " Value provided was " << ndPred;
      throw runtime_error(err.str());
    }
    if (fdPred < 1 || fdPred > 6) {
      ostringstream err;
      err << "Far detector binning option must be an integer between 1 and 6!"
          << " Value provided was " << fdPred;
      throw runtime_error(err.str());
    }


    //string path = "/cvmfs/nova.osgstorage.org/analysis/nux/binning_studies/";
    string path = pnfs2xrootd("/pnfs/nova/persistent/analysis/nux/binning_studies/");
    vector<string> suffix = { "", "_res10", "_res20", "_res30", "_res40", "_res50" };
    for (Sample& s : cfg.samples) {
      IPrediction* pred;
      if (s.detector == covmx::kNearDet) {
        string predFileName = path+"pred/pred_"+s.GetTag()+suffix[ndPred-1]+".root";
        cout << "Loading prediction from " << predFileName << endl;
        TFile* f = TFile::Open(predFileName.c_str(), "read");
        pred = LoadFromFile<NDPredictionSterile>(f, "pred_"+s.GetTag()).release();
        delete f;
        SetAxis(s, 10*(ndPred-1));
      } else {
        string predFileName = path+"pred/pred_"+s.GetTag()+suffix[fdPred-1]+".root";
        cout << "Loading prediction from " << predFileName << endl;
        TFile* f = TFile::Open(predFileName.c_str(), "read");
        pred = LoadFromFile<FDPredictionSterile>(f, "pred_"+s.GetTag()).release();
        delete f;
        SetAxis(s, 10*(fdPred-1));
      }
      s.SetPrediction(pred);
    }

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

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

    // Load cosmics
    for (Sample& s : cfg.samples) {
      if (s.detector == covmx::kFarDet) {
        string cosmicFileName = path+"cosmics/"+s.GetTag()+suffix[fdPred-1]+".root";
        cout << "Loading cosmic spectrum from " << cosmicFileName << endl;
        TFile* f = TFile::Open(cosmicFileName.c_str(), "read");
        s.SetCosmic(LoadFromFile<Spectrum>(f, s.GetTag()).release());
        delete f;
      }
    }

    cout << cfg.xVar->ShortName() << "   " << cfg.yVar->ShortName() << endl;

    // Load data
    if (cfg.fakeDataType == kAsimov || cfg.fakeDataType == kFakeSignal) {
      if (cfg.fakeDataType == kFakeSignal) {
        cout << "Using fake signal oscillation parameter set " << cfg.fakeDataSet << endl;
        SetFakeSignalParams(calc, cfg.fakeDataSet);
      } else {
        cout << "Using Asimov fake data." << endl;
      }
      cout << "Oscillation parameters used for fake data:" << 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());
        TH1D* tmp = data->ToTH1(9e20);
        if (sample.HasCosmic()) *data += *sample.GetCosmic();
        tmp = data->ToTH1(9e20);
        sample.SetData(data);
      }
      SetNus20Params(calc);
    }

    else if (cfg.fakeDataType == kFakeSignalFluctuated) {

      cout << "Loading fake signal data, type " << cfg.fakeDataSet << ", pseudoexperiment "
        << cfg.fakeDataUniverse << endl;

      ostringstream filepath;
      filepath << kNus20Path << "/fakedata/fhc_fake_signal_" << 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->GetDirectory(s.GetTag().c_str())).release());
      }
      delete f; // Close fake data file
    }

    // Load covariance matrix
    CovarianceMatrix* mx = nullptr;
    if (cfg.systType == kCovMxSysts) {
      mx = new CovarianceMatrix(cfg.samples);
      for (CovarianceMatrix* itMx : GetMatrices(cfg.samples, kUseCVMFS)) {
        bool use = false;
        if (cfg.systConfig == kAllSysts) use = true;
        else if (cfg.systConfig == kDetectorSysts) {
          for (string name : kDetectorSystNames) {
            if (itMx->GetName() == name) {
              use = true;
              break;
            }
          } // for detector syst name
        } else if (cfg.systConfig == kBeamSysts) {
          for (string name : kBeamSystNames) {
            if (itMx->GetName() == name) {
              use = true;
              break;
            }
          } // for beam syst name
        } else if (cfg.systConfig == kXsecSysts) {
          for (const ISyst* syst : getAllXsecSysts_2018_RPAFix()) {
            if (itMx->GetName() == syst->ShortName()) {
              use = true;
              break;
            }
          } // for xsec isyst
        }
        if (use) mx->AddMatrix(itMx);
        delete itMx;
      }
    } else if (cfg.systType == kCovMxSingleSyst) {
      for (CovarianceMatrix* itMx : GetMatrices(cfg.samples, kUseCVMFS)) {
        if (itMx->GetName() == cfg.systName) {
          mx = itMx;
        } else {
          delete itMx;
        }
      }
      if (!mx) { // Throw an error if we didn't find the right matrix
        ostringstream err;
        err << "Couldn't find matrix with name \"" << cfg.systName << "\"!";
        throw std::runtime_error(err.str());
      }
    } else if (cfg.systType == kCovMxWithSinglePullSyst) {
      mx = new covmx::CovarianceMatrix(cfg.samples);
      for (CovarianceMatrix* itMx : GetMatrices(cfg.samples, kUseCVMFS)) {
        if (itMx->GetName() != cfg.systName) {
          mx->AddMatrix(itMx);
        }
        delete itMx;
      }
    }
    /////////////////////////////////////////////////////////////////////
    //                                                                 //
    //                          MAKE SURFACE                           //
    //                                                                 //
    /////////////////////////////////////////////////////////////////////

    // Create experiment
    IChiSqExperiment* expt(nullptr);
    if (cfg.fitType == kPoissonLL) {
      expt = new LikelihoodCovMxExperiment(cfg.samples, mx, cfg.epsilon);
    } else {
      expt = new OscCovMxExperiment(cfg.samples, mx);
      if (cfg.fitType == kCNP) dynamic_cast<OscCovMxExperiment*>(expt)->SetCNP();
    }

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

    CovMxSurface surf(*cfg.xBins, *cfg.yBins, &multiExpt, calc, cfg.xVar, cfg.yVar,
      cfg.fitVars, cfg.fitSeeds, {} /*FIX THIS!!*/, SystShifts::Nominal(),
      cfg.jobNumber, cfg.nJobs);

    // Open output file and save surface
    string optstring = opt.Data();
    string fitName = optstring.substr(0, optstring.find("_njobs"));
    string outFileName = fitName + ".root";
    TFile* outFile = TFile::Open(outFileName.c_str(), "recreate");
    surf.SaveTo(outFile);

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

} // main function