make_nus_fc_surfs.C

Go to the documentation of this file.

// Ran at Caltech

#include "OscLib/OscCalcSterile.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/Fit/MinuitFitter.h"
#include "CAFAna/Vars/FitVars.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "NuXAna/Systs/NusSysts.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Analysis/Calcs.h"

#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSPlotFunctions.h"

using namespace ana;

#include "TCanvas.h"
#include "TFile.h"
#include "TH2.h"
#include "TLegend.h"
#include "TMarker.h"
#include "TRandom3.h"

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

void make_nus_fc_surfs(int nTrials = 100, int bin = 0, int Idx = 0)
{
  //std::cout << "Bin: " << bin << ", nTrials: " << nTrials << std::endl;
  // Load in the NuMI data spectrum
  Spectrum* intime = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/fout_nus_ana01_box_opening_restricted.root","spec_nus_cale_numi").release();
  //Spectrum* intime = LoadFromFile<Spectrum>("fout_nus_ana01_box_opening_restricted.root","spec_nus_cale_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;
  

  // Load in the out-of-NuMI-time weighted spectrum (cosmics)
  Spectrum* cosmic_scale = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/AnaResults.root","cosmic__CalE").release();
  //Spectrum* cosmic_scale = LoadFromFile<Spectrum>("AnaResults.root","cosmic__CalE").release();
  TH1D* hcosmic = (TH1D*)cosmic_scale->ToTH1(cosmic_scale->Livetime(), kLivetime);
  Spectrum cosmics(hcosmic, intime->POT(), intime->Livetime());
  double totcosmics = cosmics.ToTH1(intime->POT())->GetSumOfWeights();
  std::cout << "Total cosmic count: " << totcosmics << std::endl;

  // Load in MC prediction
  TFile input("$NOVA_ANA/steriles/Ana01/nus_ana01_prediction.root", "READ");
  //TFile input("nus_ana01_prediction.root", "READ");
  
  std::unique_ptr<PredictionInterp> pred_p1p2e3b =
    PredictionInterp::LoadFrom((TDirectory*)input.Get("nus_pred_p1p2e3b"));
  
  std::unique_ptr<PredictionInterp> pred_e3c     =
    PredictionInterp::LoadFrom((TDirectory*)input.Get("nus_pred_e3c"));
  
  const double pot_p1p2e3b =
    kSecondAnaPeriod1POT
    + kSecondAnaPeriod2POT
    + kSecondAnaEpoch3bPOT;
  
  const double pot_e3ce3d  = kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT;
  
  // Combine period1/period2/epoch3b + epoch3c, weighted appropriately
  PredictionCombinePeriods pred({
      {&(*pred_p1p2e3b),pot_p1p2e3b},
        {&(*pred_e3c),pot_e3ce3d}
    });
  
  
  osc::OscCalcSterile* calc4f = DefaultSterileCalc(4);
  calc4f->SetNFlavors(4);
  calc4f->SetDm(4,0.5);       // #Delta m^{2}_{41} = 0.5 eV^2
  calc4f->SetAngle(2, 4, 0.); // #theta_{24} = 0 degrees
  calc4f->SetAngle(3, 4, 0.); // #theta_{34} = 0 degrees

  // 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}
  std::vector<const IFitVar*> fitvars;
  fitvars.push_back(&kFitTheta34InDegreesSterile);
  fitvars.push_back(&kFitTheta24InDegreesSterile);
  fitvars.push_back(&kFitTheta23InDegreesSterile);

  FCCollection fccol;
  //Zero sets the seed to be a unique value in space and time
  TRandom3 rnd(0); // zero seeds randomly

  //Progress prog("Performing mock experiments");

  for(int iTrial = 0; iTrial < nTrials; ++iTrial){
    ResetSterileCalcToDefault(calc4f);
    ResetAngles(calc4f);
    
    std::cout << "On trial: " << iTrial << std::endl;

    
    //For each trial, getting random location in true X, Y space.
    // Limits are the edges of 2D plot you want contours drawn on.
    // 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);

    //const double trueX = rnd.Uniform(0.0, 50.);
    //const double trueY = rnd.Uniform(0.0, 45.);

    kFitTheta34InDegreesSterile.SetValue(calc4f, trueX);
    kFitTheta24InDegreesSterile.SetValue(calc4f, trueY);

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

    // for(const auto& syst:systlist) {
    //   shifts.SetShift(syst, rnd.Gaus(0, 1));
    // }

    const double pot = kSecondAnaPOT;
    
    Spectrum obs  = pred.Predict(calc4f);
    obs += cosmics;
    Spectrum mock = obs.MockData(pot);

    // Compares (mock)data sample with MC sample
    GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile,
                                      45.8,
                                      3.2);
    CountingExperiment nus_expt(&pred, mock, cosmics);
    MultiExperiment    expt({&th23Constraint, &nus_expt});

    double trueChi;
    kFitTheta23InDegreesSterile.SetValue(calc4f, rnd.Gaus(45.8,3.2));
    
    MinuitFitter fittrue(&expt, {&kFitTheta23InDegreesSterile}, systlist);
    trueChi = fittrue.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();

    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);
  
  fccol.SaveToFile("FCCol_"+binname+"_"+number+".root");
}


void ResetAngles(osc::OscCalcSterile* calc)
{
  calc->SetAngle(2, 3, M_PI/4); // 45 degress, maximal mixing
  calc->SetAngle(3, 4, 0.);     // 0 degrees
  calc->SetAngle(2, 4, 0.);     // 0 degrees

  return;
}