make_nus17_fc_surfs.C

Go to the documentation of this file.

#include "TCanvas.h"
#include "TFile.h"
#include "TH2.h"
#include "TLegend.h"
#include "TMarker.h"
#include "TRandom3.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TLatex.h"
#include "TLine.h"
#include "TPad.h"

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/FC/FCCollection.h"
#include "CAFAna/FC/FCPoint.h"
#include "CAFAna/FC/FCSurface.h"
#include "CAFAna/Vars/FitVars.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/IFitVar.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Analysis/Plots.h"
#include "NuXAna/Systs/Nus17Systs.h"
#include "OscLib/OscCalcSterile.h"

#include "NuXAna/macros/NuSLoadMacroDefs.h"

using namespace ana;

// Global variables
bool kLowerOctant = true;

// Reset calculator values after each fit to preset defaults
void ResetAngles(osc::OscCalcSterile* calc);

void make_nus17_fc_surfs(int nTrials = 100, int bin = 0, int Idx = 0, bool lowerOctant = true) {
  kLowerOctant = lowerOctant;
  TH1::AddDirectory(0);

  // Set up path to file with filled CAFAna objects and for output of analysis
  std::string folder = "/nova/ana/steriles/Nus17/";
  std::string filenm = "predInterpSysts_nus17_v1";

  std::string loadLocation = "/nova/ana/steriles/Nus17/predInterpSysts_nus17_v1.root";
  std::string saveLocation = "plot_" + filenm + "_shape_2D_lower.root";
  std::string textLocation = "count_" + filenm + "_shape_2D_lower.txt";

  // Load in the NuMI data spectrum
  Spectrum* intime = LoadFromFile<Spectrum>("/nova/ana/steriles/Nus17/fout_nus17_box_opening_restricted.root","spec_nus_vise_numi").release();

  TH1 *hSpectrum = intime->ToTH1(intime->POT());
  Spectrum data_spec(hSpectrum, intime->POT(), intime->Livetime());
  std::cout << "Data spectrum: " << data_spec.POT()
      <<" POT and "
      << data_spec.ToTH1(intime->POT())->GetSumOfWeights()
      << " events" << std::endl;
  delete hSpectrum;

  TFile *cos = new TFile("/nova/ana/steriles/Nus17/nus17_cosmic_background.root","read");
  TH1   *hcosmic = (TH1*)cos->Get("cosall");
  Spectrum cosmic(hcosmic, kAna2017POT, kAna2017Livetime);
  std::cout << "Cosmic Data spectrum: " << cosmic.POT()
      <<" POT and "
      << cosmic.ToTH1(cosmic.POT())->GetSumOfWeights()
      << " events"
      << ", Livetime: "
      << cosmic.Livetime() << std::endl;
  delete hcosmic;
  cos->Close();

  // Load filled objects from file
  TFile* rootL = new TFile(loadLocation.c_str(), "read");
  TDirectory* tmpL = gDirectory;
  TDirectory* loadDir = gDirectory;

  loadDir->cd((loadLocation + ":/nus17_predI").c_str());
  std::unique_ptr<PredictionInterp> pred = PredictionInterp::LoadFrom(gDirectory);

  tmpL->cd();
  rootL->Close(); // Don't forget to close the file!

  // OscCalc for Sterile 3+1 Model
  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  ResetAngles(calc4f);

  // theta_34
  double minX = 0.;
  double maxX = 50.;
  int    NX   = 50;
  double widthX = (maxX-minX)/(NX-1);

  // theta_24
  double minY = 0.;
  double maxY = 45.;
  double NY   = 45.;
  double widthY = (maxY-minY)/(NY-1);

  assert(bin < NX*NY && "Invalid bin number");
  int binY = bin/NX;
  int binX = bin - NX*binY;
  double centerX = minX + binX*widthX;
  double centerY = minY + binY*widthY;

  // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}, #deltam_{32}
  std::vector<const IFitVar*> fitvars;
  fitvars.push_back(&kFitTheta34InDegreesSterile);
  fitvars.push_back(&kFitTheta24InDegreesSterile);
  fitvars.push_back(&kFitSinSqTheta23Sterile); // to be constrained to PDG
  fitvars.push_back(&kFitDmSq32Sterile); // to be constrained to PDG
  fitvars.push_back(&kFitDelta24InPiUnitsSterile); // profile delta_24

  FCCollection fccol;
  TRandom3 rnd(0); // zero seeds randomly
  
  for(int iTrial = 0; iTrial < nTrials; ++iTrial){
    ResetSterileCalcToDefault(calc4f);
    ResetAngles(calc4f);

    std::cout << "On trial: " << iTrial << std::endl;

    // Throw a true X,Y value, randomly accross current bin
    double trueX = -1;
    // Edge bins centered on a physical edge
    while (trueX < minX || trueX > maxX)
      trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
    double trueY = -1;
    while (trueY < minY || trueY > maxY)
      trueY = rnd.Uniform(centerY-widthY/2,centerY+widthY/2);

    kFitTheta34InDegreesSterile.SetValue(calc4f, trueX);
    kFitTheta24InDegreesSterile.SetValue(calc4f, trueY);
    kFitDelta24InPiUnitsSterile.SetValue(calc4f, rnd.Uniform(0, 2));

    //    SystShifts shifts;
    std::vector<const ISyst*> systlist = getAllNusFDSysts();

    // Update POT to 2017 Analysis POT
    const double pot = kAna2017POT;

    Spectrum obs  = pred.get()->Predict(calc4f);
    obs += cosmic;
    Spectrum mock = obs.MockData(pot);
    TH1 *hMockSpectrum = mock.ToTH1(pot);
    Spectrum mock_spec(hMockSpectrum, pot, intime->Livetime());
    delete hMockSpectrum;

    // Compares (mock)data sample with MC sample
    Double_t meanSinSqTheta23Sterile = 0;
    Double_t sigmaSinSqTheta23Sterile = 0;

    if (kLowerOctant) {
        meanSinSqTheta23Sterile = 0.404;
        sigmaSinSqTheta23Sterile = 0.030;
    } else {    
        meanSinSqTheta23Sterile = 0.603;
        sigmaSinSqTheta23Sterile = 0.030;
    }

    GaussianConstraint s2th23Constraint(&kFitSinSqTheta23Sterile, meanSinSqTheta23Sterile, sigmaSinSqTheta23Sterile);
    GaussianConstraint dmsq32Constraint(&kFitDmSq32Sterile, 2.67e-3, 0.11e-3);

    SingleSampleExperiment nus_expt(pred.get(), mock_spec, cosmic);
    MultiExperiment expt({&s2th23Constraint, &dmsq32Constraint, &nus_expt});

    if (kLowerOctant) {
      kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.404, 0.030));
    } else {
      kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.603, 0.030));
    }
    kFitDmSq32Sterile.SetValue(calc4f, rnd.Gaus(2.67e-3, 0.11e-3));

    MinuitFitter fittrue(&expt, {&kFitTheta23InDegreesSterile, &kFitDmSq32Sterile, &kFitDelta24InPiUnitsSterile}, systlist);
    double trueChi = fittrue.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();

    if (kLowerOctant) {
      kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.404, 0.030));
    } else {
      kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.603, 0.030));
    }
    kFitDmSq32Sterile.SetValue(calc4f, rnd.Gaus(2.67e-3, 0.11e-3));

    MinuitFitter fit(&expt, fitvars, systlist);
    double bestChi = fit.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();

    // The fit set the best fit values in the calculator
    // now we are getting them back
    const double bestX = kFitTheta34InDegreesSterile.GetValue(calc4f);
    const double bestY = kFitTheta24InDegreesSterile.GetValue(calc4f);

    // Now we are setting true oscillation values, chi squared value
    // when fit with true points, best fit oscillation values, and
    // chi squared value for the best fit point
    if (trueChi < bestChi) {
      std::cout << "True Chi less than Best Chi" << std::endl;
      bestChi = trueChi;
    }

    FCPoint pt(trueX, trueY, trueChi,bestX, bestY, bestChi,0); // Don't know seed or rng state
    if (100*iTrial % nTrials == 0){
      std::cout << "Trial: " << iTrial << ", TrueX: " << trueX << ", Y: " << trueY << ", trueChi: " << trueChi << ", bestX: " << bestX << ", Y: " << bestY << ", bestChi: " << bestChi << ", true-best: " << trueChi-bestChi << std::endl;
    }

    fccol.AddPoint(pt);
  }
  
  std::string binname = std::to_string(bin);
  std::string number  = std::to_string(Idx);
  std::string octant  = "_lowerOctant";
  if (!kLowerOctant) {
        octant = "_upperOctant";
  }

  fccol.SaveToFile("FC_nus17_Col_"+binname+"_"+number+octant+".root");
}


void ResetAngles(osc::OscCalcSterile* calc)
{
  ResetSterileCalcToDefault(calc);
  calc->SetDm(4, 0.5);
  calc->SetAngle(3, 4, 0); // 0, icecube
  calc->SetAngle(2, 4, 0); // 0, sk
  calc->SetDm(3, 7.53e-5 + 2.67e-3);
  if (kLowerOctant) {
    calc->SetAngle(2, 3, (asin(sqrt(0.404))));
    calc->SetDelta(1, 3, 1.48*M_PI);
  } else {
    calc->SetAngle(2, 3, (asin(sqrt(0.603))));
    calc->SetDelta(1, 3, 0.74*M_PI);
  }
  return;
}