nue_data_mc_validation.C

Go to the documentation of this file.

// Make Data and MC Spectra for validation
// makeFile = true saves the spectra for ldrData, ldrMC
// makeFile = false reruns formatting for generic validation

// SA nue cuts and variables
// Same macro should work for different datasets,
// as long as all vars and cuts are proper

// TODO: add logic for FD
// TODO: add logic to split data and MC

#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "3FlavorAna/Plotting/NuePlotLists.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2018Alt.h"
#include "OscLib/IOscCalc.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include "CAFAna/Vars/CVNFinalStates.h"
#include "3FlavorAna/Vars/DecompWeights.h"

using namespace ana;

#include "TFile.h"

#include <iostream>
#include <vector>

void nue_data_mc_validation(const std::string beamMode,
                            const std::string predFileName,
                            const bool makeFile=false) 
{
  // Variable groups from NueVarsExtra
  // Variables from NueVars.h are used here or in the cuts
  assert(beamMode == "rhc" || beamMode == "fhc");

  std::string ldrData;
  std::string ldrMC;
  if(beamMode == "rhc") {
    ldrData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns";
    ldrMC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
  }
  else {
    ldrData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
    ldrMC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
  }

  std::vector <HistDef> defs;
  
  AddNueHistDefBasic(defs); 
  AddNueHistDefVertexND(defs); 
  AddNueHistDefShowerND(defs);
  AddNueHistDefShower(defs);
  AddNueHistDefPIDs(defs);     
  AddNueHistDefForCvn(defs);   
  AddNueHistDefSelectionExtras(defs);
  AddNueHistDefProngCVN(defs); 

  // Cut flow. Change the names/cuts here only if new spectra are going to be saved
  struct CutDef
  {
    std::string name;
    Cut cut;
  };

  std::vector <CutDef> sels = {
      {"0_Presel",kNue2018NDPresel},
      {"1_CVNeLow",kNue2018NDCVNSsb},
  };
  if(beamMode ==  "rhc")
    sels.push_back({"2_CVNeHigh",kNue2018NDPresel && kCVNSSe >= .98});
  else
    sels.push_back({"2_CVNeHigh",kNue2018NDPresel && kCVNSSe >= .96});

  std::string decompFName = "nue2018_decomp_ratios_rhc_prop.root";
  std::string decompShortName = "prop";
  if(beamMode == "fhc") {
    decompFName = "nue2018_decomp_ratios_fhc.root";
    decompShortName = "combo";
  }

  const Var kDecompWeight(DecompWeightFunc(FindCAFAnaDir()+
                                           "/data/decomp/" + decompFName,
                                           decompShortName,
                                           kNue2018AnaBin));

  if(makeFile){
    SpectrumLoader loaderData(ldrData.c_str());
    SpectrumLoader loaderMC(ldrMC.c_str());
    
    loaderData.SetSpillCut(kStandardSpillCuts);
    loaderMC.SetSpillCut(kStandardDQCuts);
    
    std::vector < std::vector<Spectrum*> > spects;
    std::vector < std::vector<IPrediction*> > preds_mc;
    std::vector < std::vector<IPrediction*> > preds_mc_decomp;
    
    for(int selIdx = 0; selIdx<(int)sels.size(); ++selIdx){
      std::vector<Spectrum *> temp_spec;
      std::vector<IPrediction *> temp_pred_mc;
      std::vector<IPrediction *> temp_pred_mc_decomp;
      
      for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
            const HistAxis& axis = defs[varIdx].axis;
            temp_spec.emplace_back( new Spectrum(loaderData, axis, sels[selIdx].cut, kNoShift));
            temp_pred_mc.emplace_back(new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
                                      axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
                                                             sels[selIdx].cut, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt)); //TODO Broken w new cafs
            temp_pred_mc_decomp.emplace_back(new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
                                                                  axis.GetLabels()[0], axis.GetBinnings()[0],
                                                                  axis.GetVars()[0], sels[selIdx].cut,
                                                                  kNoShift,
                                                                  kXSecCVWgt2018*kPPFXFluxCVWgt*kDecompWeight));
      }
      spects.push_back(temp_spec);
      preds_mc.push_back(temp_pred_mc);
      preds_mc_decomp.push_back(temp_pred_mc_decomp);
    }
    
    loaderMC.Go();
    loaderData.Go();
    
    TFile* fout = new TFile(predFileName.c_str(), "RECREATE");
    for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
      const char* cut = sels[selIdx].name.c_str();
      std::cout << cut << std::endl;
      for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
        const char* name = defs[varIdx].name.c_str();
        std::string tempname = defs[varIdx].name;
        spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
        preds_mc[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_mc_%s_%s", name, cut));
        preds_mc_decomp[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_mc_decomp_%s_%s", name, cut));
        
      }
    }
    delete fout;
    return;
  } //end makeFile
  else{
  // Play with the format for the generic validation scripts
  // Here I'm taking the one file and splitting data and MC
  // Data and MC will be treated as "datasets" = separate files
  // Alternatively, you could treat them as folders, or categories in a group
  // Your mileage may vary
  
  TFile * predFile = new TFile (predFileName.c_str(), "READ");
  TString validFileName1 (predFileName); validFileName1.ReplaceAll(".root","_mc_validation.root");
  TString validFileName3 (predFileName); validFileName3.ReplaceAll(".root","_data_validation.root");
  TString validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_mc_decomp_validation.root");

  TFile * validFile1 = new TFile(validFileName1,"RECREATE");
  TFile * validFile3 = new TFile(validFileName3,"RECREATE");
  TFile * validFile2 = new TFile(validFileName2,"RECREATE");

  int kNumSels = sels.size();
  int kNumVars = defs.size();

  // Don't make all the canvases and .json right now
  // kNumSels = 4;
  //  kNumVars = 2;
  
  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
    const char* cut = sels[selIdx].name.c_str();
    for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
      const char* name = defs[varIdx].name.c_str();
      
      auto spec = Spectrum::LoadFrom(predFile, TString::Format("spect_%s_%s", name, cut));
      auto pred_mc = PredictionNoExtrap::LoadFrom (predFile,TString::Format("pred_mc_%s_%s", name, cut).Data());
      auto pred_mc_decomp = PredictionNoExtrap::LoadFrom(predFile, TString::Format("pred_mc_decomp_%s_%s", name, cut).Data());

      
      // Subfolders -> new pages (not shown here)
      // Data and MC for variable "name" will be plotted together 
      // Plots will have selector for cut level {group=Cut,cat=%s}
      // Plots will have selector for MC component component; total MC = all_sum will be created as well
      // Add functions to operate on variable, cut or axis labels 

      spec->SaveTo(validFile3, TString::Format("%s{group=Cut,cat=%s}", name,cut));

      pred_mc->ComponentCC(12,12).  Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nue}", name,cut));
      pred_mc->ComponentCC(14,14).  Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_numu}", name,cut));
      pred_mc->ComponentCC(-12,-12).Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
      pred_mc->ComponentCC(-14,-14).Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
      pred_mc->ComponentNCTotal().                      SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_nc}", name,cut));
      

      pred_mc_decomp->ComponentCC(12,12).  Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nue}", name,cut));
      pred_mc_decomp->ComponentCC(14,14).  Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_numu}", name,cut));
      pred_mc_decomp->ComponentCC(-12,-12).Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
      pred_mc_decomp->ComponentCC(-14,-14).Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
      pred_mc_decomp->ComponentNCTotal().                      SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_nc}", name,cut));
      
    }
  }
  std::cout << "Formatted outputs saved to " << validFileName1 << ", " << validFileName2 << " and " << validFileName3 << std::endl;
  delete validFile1;
  delete validFile3;
  
  }
}