test_nue2018_surface.C

Go to the documentation of this file.

#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Analysis/Calcs.h"
#include "OscLib/IOscCalc.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"

using namespace ana;

#include <iostream>

void test_nue2018_surface(bool create){

  if(create){
    double pot = 6E20;
    osc::IOscCalcAdjustable* calc = DefaultOscCalc();
    
    auto pred = new PredictionSystJoint2018 ( kProportional, calc,
                                        {},
                                        "fhc",
                                        false, true, "", // isfakedata, mergePeripheral, 
                                        true, {}); //mini-memo, summingsys

    auto data = pred->Predict(calc).FakeData(pot);
    SingleSampleExperiment expt (pred, data);
    
    FrequentistSurface surf(&expt, calc,
                 &kFitSinSqTheta23, 30,0.3,0.75,
                 &kFitDmSq32, 30, 2e-3, 3.5e-3,
                 {&kFitDeltaInPiUnits}, {}, {}, {},false);
    TFile * file = new TFile ("test_surf.root","recreate");
    surf.SaveTo(file, "surf");
    file->Close();
  }
  else{
    TFile * file = new TFile ("test_surf.root","read");
    
    std::unique_ptr<FrequentistSurface> surf = ISurface::LoadFrom(file, "surf");

    surf->Draw();

    surf->DrawBestFit(kRed);
    std::cout << "MinChi " << surf->BestLikelihood() << std::endl;
    new TCanvas;
    surf->DrawContour(Gaussian90Percent2D(*surf),kSolid, kRed);
    new TCanvas;
    auto h2 = surf->ToTH2();
    h2->Draw("colz");

    auto graphs = surf->GetGraphs(Gaussian2Sigma2D(*surf));
    for (auto gr:graphs){
      gr->SetLineColor(kBlue);
      gr->Draw("C same");
    }

    for (auto h:(surf->GetMarginalizedHists())){
      new TCanvas;
      h->Draw();
    }
    
    
  }


}