NDDataMC.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

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

#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/MultiVar.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Vars/Vars.h"
#include "NuXAna/Vars/NusVars.h"
#include "NuXAna/Analysis/NusPlotLists.h"
#include "OscLib/IOscCalc.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

using namespace ana;

#include "TFile.h"

#include <iostream>
#include <vector>

void NDDataMC(const std::string predFileName, const bool makeFile=false, 
              const std::string ldrData="", const std::string ldrMC="")
{
  // Variable groups from NusVars
  std::vector<HistDef> nus_defs = getAllNusAna01NDHistDefs();

  // 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_NoCut", kNoCut},
    {"1_NusEventQuality", kNus17EventQuality},
    {"2_NusFiducial", kNus17EventQuality && kNus17NDFiducial},
    {"3_NusPresel", kNus17NDPresel},
    {"4_kNusND", kNus17ND},
    {"5_kNusNDExtra", kNus17NDExtra},
    {"6_kNusNDPtp", kNus17NDExtra && (kPartPtp <= 0.8)}
  };
  
  MultiVar kCutFlowVar(
    [sels](const caf::SRProxy* sr)
    {
      std::vector<double> passes;
      for(unsigned int i = 0; i < sels.size(); ++i){
        if(sels[i].cut(sr)) passes.push_back(i);
      }
      return passes;
    });

  const Binning flowBins = Binning::Simple(6, 0, 6);
  const Var kDummy(
    [](const caf::SRProxy* sr)
    {
      return 1;
    });
  nus_defs.insert(nus_defs.begin(), {"_dummy", {"Dummy", Binning::Simple(3, -0.5, 2.5), kDummy}});
  
  //Save Data and MC spectra. Prediction -> get MC components
  
  if(makeFile){
    SpectrumLoader loaderData(ldrData.c_str());
    SpectrumLoader loaderMC(ldrMC.c_str());
    
    loaderData.SetSpillCut(kStandardSpillCuts);
    loaderMC.SetSpillCut(kStandardSpillCuts);
    
    std::vector<std::vector<Spectrum*> > spects;
    std::vector<std::vector<Spectrum*> > spectsMC;
    std::vector<std::vector<CheatDecomp*> > preds;
    
    for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
      std::vector<Spectrum*>    temp_spec;
      std::vector<Spectrum*>    temp_specMC;
      std::vector<CheatDecomp*> temp_pred;
      
      for(int varIdx = 0; varIdx < (int)nus_defs.size(); ++varIdx){
        const HistAxis& axis = nus_defs[varIdx].axis;
        temp_spec.emplace_back(new Spectrum(loaderData, axis, sels[selIdx].cut, kNoShift));     
        temp_specMC.emplace_back(new Spectrum(loaderMC, axis, sels[selIdx].cut, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt));
        temp_pred.emplace_back(new CheatDecomp(loaderMC,
                                               axis,
                                               sels[selIdx].cut,
                                               kNoShift,
                                               kXSecCVWgt2017*kPPFXFluxCVWgt));
      }
      spects.push_back(temp_spec);
      spectsMC.push_back(temp_specMC);
      preds.push_back(temp_pred);
    }

    auto spfl_data = new Spectrum("Cut Flow", flowBins, loaderData, kCutFlowVar, kNoCut);
    auto spfl_ndmc = new Spectrum("Cut Flow", flowBins, loaderMC, kCutFlowVar, kNoCut);
    
    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();
      for(int varIdx = 0; varIdx < (int)nus_defs.size(); ++varIdx){
        const char* name = nus_defs[varIdx].name.c_str();
        spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
        spectsMC[selIdx][varIdx]->SaveTo(fout, TString::Format("spectMC_%s_%s", name, cut));
        preds[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_%s_%s", name, cut));
        
      }
    }

    spfl_data->SaveTo(fout, "spect_data_cutflow");
    spfl_ndmc->SaveTo(fout, "spect_mc_cutflow");
    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 validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_data_validation.root");
    
    TFile * validFile1 = new TFile(validFileName1,"RECREATE");
    TFile * validFile2 = new TFile(validFileName2,"RECREATE");

    int kNumSels = sels.size();
    int kNumVars = nus_defs.size();
    
    
    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 = nus_defs[varIdx].name.c_str();
        
        auto spec   = Spectrum::LoadFrom(predFile, TString::Format("spect_%s_%s", name, cut));
        auto specMC = Spectrum::LoadFrom(predFile, TString::Format("spectMC_%s_%s", name, cut));
        auto pred   = CheatDecomp::LoadFrom(predFile, TString::Format("pred_%s_%s", name, cut));
        
        
        spec->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}", name,cut));
        specMC->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}", name,cut)); 
        
        // 1D Spectra only
        bool components = true;
        if(components){
          TDirectory* comp = validFile1->GetDirectory("MCComponents");
          if(!comp) comp = validFile1->mkdir("MCComponents");
        
          pred->NCTotalComponent().      SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nc}", name,cut));
          pred->NueComponent().     SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_nue}", name,cut));
          pred->NumuComponent().    SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_numu}", name,cut));
          pred->AntiNueComponent(). SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
          pred->AntiNumuComponent().SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
          
          (*specMC-(pred->TotalMC())).SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=5_other}",name,cut));
        }
      } // end var
      auto spfl = Spectrum::LoadFrom(predFile, "spect_data_cutflow");
      spfl->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}","_cutflow",cut));

      spfl = Spectrum::LoadFrom(predFile, "spect_mc_cutflow");
      spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}","_cutflow",cut));
    
    } // end cut
    std::cout << "Formatted outputs saved to " << validFileName1 << " and " << validFileName2 << std::endl;
      delete validFile1;
      delete validFile2;

  }
}