test_nue2018_fitter.C

Go to the documentation of this file.

#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Fit/MinuitFitter.h"

#include "3FlavorAna/Systs/JointAna2018Systs.h"

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/LoadFromFile.h"

#include "CAFAna/Vars/FitVars.h"

#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"

#include "CAFAna/Experiment/SingleSampleExperiment.h"

#include "Utilities/func/MathUtil.h"

#include "TCanvas.h"
#include "TH1.h"
#include "TFile.h"

#include <iostream>
#include <cmath>

using namespace ana;

void test_nue2018_fitter()
{
  double pot = 9.48e20; // From Nitish

  std::string dirName = "./";

  double dCP = 1.5;
  double ssth23 = 0.402;
  double dmsq32 = 0.0027;

  // random testing calc
  osc::IOscCalcAdjustable* calc = DefaultOscCalc();

  kFitDeltaInPiUnits.SetValue(calc,dCP);
  kFitDmSq32.SetValue(calc,dmsq32);
  kFitSinSqTheta23.SetValue(calc,ssth23);

  std::vector<const IFitVar*> fitvar = {&kFitDmSq32, &kFitDeltaInPiUnits, &kFitSinSqTheta23};

  auto pred = new PredictionSystJoint2018 ( kProportional, calc, getAna2018LargeXsecSysts(kNueAna2018), "fhc",
                                                 false, true, "", // isfakedata, mergePeripheral, 
                                                 false, {}); //mini memo, summingsys

  Spectrum data = pred->Predict(calc).FakeData(pot);
  std::string anaDir = "/nova/ana/nu_e_ana/Ana2018/Predictions/cosmic/cosmic_prediction_real_data.root";
  Spectrum* cosmic = LoadFromFile<Spectrum>(anaDir, "cosmic_spect_fhc" ).release();
  TH1D* hcosmic = cosmic->ToTH1( 1, kLivetime);
  hcosmic->Reset();
  std::cout<<"sanity check:"<<hcosmic->Integral()<<std::endl;
  const SingleSampleExperiment* exp = new SingleSampleExperiment( pred, data, hcosmic);
  
  MinuitFitter fitter(exp, fitvar, {});
  SystShifts bestSysts;

  double counter, chi;

  MinuitFitter::FitOpts prec[4] = {MinuitFitter::FitOpts::kFast, MinuitFitter::FitOpts::kNormal, MinuitFitter::FitOpts::kCareful, MinuitFitter::FitOpts::kGradDesc};
  std::cout<<"Fake data: dcp = "<<calc->GetdCP()<<"; ssth23 = "<<calc->GetTh23()<<"; dmsq32 = "<<calc->GetDmsq32()<<std::endl;
  std::cout<<"Starting the loop..."<<std::endl;
  
  for(auto& p: prec){
    ResetOscCalcToDefault(calc);
    clock_t begin = clock();  
    fitter.SetFitOpts(p);
    chi = fitter.Fit(calc, bestSysts, {}, {})->EvalMetricVal();
    clock_t end = clock();
    counter = (double)(end - begin)/CLOCKS_PER_SEC;
    std::cout<<"Timer: "<<counter<<std::endl;
  }

}