MakeSurfaceJoint.C

Go to the documentation of this file.

#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Sample.h"
#include "CAFAna/Prediction/CovarianceMatrix.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/OscCovMxExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "NuXAna/Systs/Nus18Systs.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/CovMxSurface.h"

#include "CAFAna/covmx/Utilities.h"

using namespace ana;

void MakeSurfaceJoint()
{
  DontAddDirectory guard;

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

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

  std::vector<covmx::Sample> samples;
  std::vector<bool> unoscillated(samples.size(), false);
  for (covmx::Polarity pol : { covmx::Polarity::kFHC, covmx::Polarity::kRHC }) {
    for (covmx::Detector det : { covmx::Detector::kNearDet, covmx::Detector::kFarDet }) {
      samples.push_back(covmx::Sample(covmx::kNC, pol, det));
      unoscillated.push_back(det == covmx::Detector::kNearDet);
      samples.push_back(covmx::Sample(covmx::Analysis::kCCNue, pol, det));
      unoscillated.push_back(det == covmx::Detector::kNearDet);
      for (covmx::Quantile q : { covmx::Quantile::kQ1, covmx::Quantile::kQ2, covmx::Quantile::kQ3, covmx::Quantile::kQ4 }) {
        samples.push_back(covmx::Sample(covmx::Analysis::kCCNumu, pol, det, q));
        unoscillated.push_back(det == covmx::Detector::kNearDet);
      }
    }
  }
  for (size_t i = 0; i < samples.size(); ++i) {
    std::cout << "Sample: " << samples[i].GetTag() << ", disabling oscillations? " << unoscillated[i] << std::endl;
  }

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

  // Create sterile oscillation calculator
  osc::OscCalcSterile* calc = DefaultSterileCalc(4);
  SetNus20Params(calc);

  // Load predictions
  std::vector<IPrediction*> preds;
  for (covmx::Sample s : samples) preds.push_back(LoadPrediction(s, "systs", kUseCVMFS));

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

  // Load data
  std::vector<Spectrum*> data;
  // if (fakedata_dir) {
  //   for (auto sample : samples)
  //     data.push_back(Spectrum::LoadFrom(fakedata_dir->GetDirectory(sample.GetName().c_str())).release());
  // }
  // else {
    for (size_t i = 0; i < preds.size(); ++i) {
      double pot = (samples[i].detector == covmx::Detector::kNearDet) ? potND : potFD;
      data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(pot)));
      if (samples[i].detector == covmx::kFarDet)
        data.back()->OverrideLivetime(livetimeFD);
      if (cosmic[i]) *data[i] += *cosmic[i];
    }
  // }

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

  // Create experiment
  OscCovMxExperiment* expt = nullptr;
  // if (mx) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
  // else expt = new OscCovMxExperiment(preds, data, cosmic);
  expt = new OscCovMxExperiment(preds, data, cosmic);
  auto noosccalc = DefaultSterileCalc(4);
  for (int i = 2; i <= 4; ++i) calc->SetDm(i, 0);
  expt->DisableOscillations(unoscillated, noosccalc);

  std::vector<const IChiSExperiment*> expts;
  expts.push_back(new const GaussianConstraint(&kFitSinSq2Theta13Sterile, 0.082, 0.004));
  expts.push_back(expt);
  MultiExperiment multiExpt(expts);

  // Set seed values and fit vars
  std::vector<const IFitVar*> fitVars;
  std::map<const IFitVar*, std::vector<double> > seedValues;
  std::vector<double> fit_seeds;

  fitVars.push_back(&kFitDelta13InPiUnitsSterile);
  seedValues[&kFitDelta13InPiUnitsSterile] = { 1.5 };

  fitVars.push_back(&kFitSinSq2Theta13Sterile);
  seedValues[&kFitSinSq2Theta13Sterile] = { 0.082 };

  // Get max length of fit var name, for cout padding
  TH2D* varHists[fitVars.size()];
  unsigned int padWidth = (*std::max_element(fitVars.begin(), fitVars.end(),
    [] (const IFitVar* lhs, const IFitVar* rhs) {
      return lhs->ShortName().length() < rhs->ShortName().length();
  }))->ShortName().length();

  Binning xbins(Binning::Simple(50, 0.3, 0.7));
  Binning ybins(Binning::Simple(50, 2e-3, 3.3e-3));

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

  CovMxSurface surf(xbins, ybins, &multiExpt, calc,
    &kFitSinSqTheta23Sterile, &kFitDmSq32Sterile,
    fitVars, seedValues, {}, SystShifts::Nominal(),
    true, nullptr, preds, data, cosmic, part, njobs);

  // Open output file
  TFile* f = TFile::Open("surface_joint.root", "recreate");
  surf.SaveTo(f->mkdir("surface"));
  delete f;

}