datamc_ND_numu_kinematics.C

Go to the documentation of this file.

#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "StandardRecord/Proxy/SRProxy.h" 

// Core
#include "CAFAna/Core/EventList.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/Cut.h"

// Cuts
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/NueCutsFirstAna.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

// Vars
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/PPFXWeights.h"

// Systs
#include "CAFAna/Systs/Systs.h" 
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Decomp/IDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"

#include "OscLib/OscCalcPMNSOpt.h"


#include "CAFAna/Core/SystShifts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "TFile.h"
#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Core/Spectrum.h"
#include "TH1.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "TCanvas.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"

#include "TStopwatch.h"
#include "CAFAna/Core/Utilities.h"



#include <fstream>
#include <iostream>
#include <cmath>
#include <string>

using namespace ana;

#include "OscLib/IOscCalc.h"

// ROOT
#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TLegend.h"
#include "TPad.h"
#include "TMath.h"
#include "TH2.h"
#include "TLine.h"
#include "TText.h"
#include "TString.h"

#include "Utilities/func/MathUtil.h"
#include "3FlavorAna/Vars/NumuVars.h"


void datamc_ND_numu_kinematics()
{
    
  std::string nameData = "defname: prod_decaf_R17-03-01-prod3reco.d_nd_numi_fhc_full_nue_or_numu_or_nus_contain_v1_goodruns";
  std::string nameMC = "defname: prod_decaf_R17-03-01-prod3reco.d_nd_genie_nonswap_fhc_nova_v08_full_nue_or_numu_or_nus_contain_v1";
  
  //choose your samweb datasets or concat files: 
  SpectrumLoader loaderData(nameData);
  SpectrumLoader loaderMC(nameMC);

  loaderData.SetSpillCut(kStandardSpillCuts);
  loaderMC.SetSpillCut(kStandardDQCuts);

  //which variables do you want to look at?
  const int kNumVars = 7;
  
  const HistDef defs[kNumVars] = {
      {"recoE", {"E_{reco} (GeV)", Binning::Simple(100, 0, 5), kCCE}},
      {"trueQ2", {"true Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kTrueQ2}},
      {"recoQ2", {"reco Q^{2} (GeV^{2})", Binning::Simple(500, 0, 5), kRecoQ2}},
      {"trueW2", {"true W^{2} (GeV^{2})", Binning::Simple(100, 0, 2), kTrueW}},
      {"trueE", {"true E (GeV)", Binning::Simple(100, 0, 5), kTrueE}},
      {"PtP", {"p_{t}/p", Binning::Simple(500,0,1), kPtP}},
      {"CosNumi", {"CosNumi", Binning::Simple(500,0,1), kCosNumi}}
  };
  
  const Cut kIsRock(
       [](const caf::SRProxy* sr)
       {
         if (sr->mc.nnu == 0)
            return false;
         if (fabs(sr->mc.nu[0].vtx.X())>180 ||
             fabs(sr->mc.nu[0].vtx.Y())>180 ||
             sr->mc.nu[0].vtx.Z()<0 ||
             sr->mc.nu[0].vtx.Z()>1250) 
            return true;
         return false;
       });
   
  const int kNumSels = 1;
  const Cut sels[kNumSels] = {
      kNumuNDCvn 
      };
  const std::string selNames[kNumSels] = {
      "numuND"
      };

  Spectrum* spects[kNumSels][kNumVars];
  IPrediction* preds[kNumSels][kNumVars];

  //loop over the selectors and variables to create a full selection
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const HistAxis& axis = defs[varIdx].axis;
      spects[selIdx][varIdx] = new Spectrum(loaderData, axis, sels[selIdx], kNoShift);
      preds[selIdx][varIdx] = new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
                                                     axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
                                                     sels[selIdx], kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt);           
    }
  }

  loaderMC.Go();
  loaderData.Go();
  
  TString fname = "datamc_ND_numu_kinematics.root";
  TFile* fout = new TFile(fname, "RECREATE");
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    TDirectory* d = fout->mkdir(selNames[selIdx].c_str());
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* name = defs[varIdx].name.c_str();
      spects[selIdx][varIdx]->SaveTo(d, TString::Format("spect_%s", name));
      preds[selIdx][varIdx]->SaveTo(d, TString::Format("pred_%s", name));
    }
  }
}