test_genie_systs.C

Go to the documentation of this file.

#include "OscLib/OscCalcPMNSOpt.h"

#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Vars/Vars.h"

#include "CAFAna/Systs/Systs.h"

#include "TCanvas.h"
#include "TGraph.h"
#include "TH1.h"

#include <cmath>
#include <iostream>

#include "Utilities/func/MathUtil.h"

using namespace ana;

void test_genie_systs()
{
  // NB: this is only about 1% of the available stats
  SpectrumLoader loader("$NOVA_PROD/mc/S14-05-12/fardet/genie/prod_pid_caf_loosened_presel_neutrino2014nightmare_besttoavoid/fardet_genie_fhc_nonswap_none_3000_r0100000?_s01*.caf.root");
  SpectrumLoader loaderSwap("$NOVA_PROD/mc/S14-05-12/fardet/genie/prod_pid_caf_loosened_presel_neutrino2014nightmare_besttoavoid/fardet_genie_fhc_swap_none_3000_r0100000?_s01*.caf.root");

  // Make a calculator. Any of the variants is fine
  osc::OscCalcPMNSOpt calc;

  // Standard oscillation parameters
  calc.SetL(810);
  calc.SetRho(2.75);
  calc.SetDmsq21(7.6e-5);
  calc.SetDmsq32(2.35e-3);
  calc.SetTh12(asin(sqrt(.87))/2);
  calc.SetTh13(asin(sqrt(.10))/2);
  calc.SetTh23(TMath::Pi()/4);
  calc.SetdCP(0);

  HistAxis axis("Calorimetric Energy (GeV)", 20, 0, 5, kCaloE);

  // Nominal prediction
  PredictionNoExtrap predNom(loader, loaderSwap, axis.label, axis.bins, axis.var, kLEM > 0.8, kNoShift);

  // M_A QE +1 sigma
  PredictionNoExtrap predUp(loader, loaderSwap, axis.label, axis.bins, axis.var, kLEM > 0.8, SystShifts(GetGenieKnobSyst(rwgt::fReweightMaCCQE), +1));

  // M_A QE -1 sigma
  PredictionNoExtrap predDn(loader, loaderSwap, axis.label, axis.bins, axis.var, kLEM > 0.8, SystShifts(GetGenieKnobSyst(rwgt::fReweightMaCCQE), -1));


  // Load the files
  loader.Go();
  loaderSwap.Go();

  // Draw the plot

  TH1* hNom = predNom.Predict(&calc).ToTH1(18e20);
  hNom->Draw("hist");

  TH1* hUp = predUp.Predict(&calc).ToTH1(18e20);
  hUp->SetLineColor(kRed);
  hUp->Draw("hist same");

  TH1* hDn = predDn.Predict(&calc).ToTH1(18e20);
  hDn->SetLineColor(kBlue);
  hDn->Draw("hist same");

  new TCanvas;
  TH1* rUp = (TH1*)hUp->Clone();
  rUp->Divide(hNom);
  rUp->Draw("hist ][");
  rUp->GetYaxis()->SetTitle("Ratio to nominal");
  rUp->GetYaxis()->SetRangeUser(.8, 1.2);

  TH1* rDn = (TH1*)hDn->Clone();
  rDn->Divide(hNom);
  rDn->Draw("hist same ][");

  TGraph* one = new TGraph;
  one->SetPoint(0, 0, 1);
  one->SetPoint(1, 100, 1);
  one->SetLineWidth(2);
  one->SetLineStyle(7);
  one->Draw("same");
}