nue_data_mc_validation.C

Go to the documentation of this file.

// Make Data and MC Spectra for validation

#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 "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "OscLib/IOscCalc.h"
#include "CAFAna/Vars/GenieWeights.h"

using namespace ana;

#include "TFile.h"

#include <iostream>
#include <vector>

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

  std::vector <HistDef> defs;
  AddNueHistDefSecondAnaBasic(defs); 
  AddNueHistDefVertexND(defs); 
  AddNueHistDefShowerND(defs); 
  AddNueHistDefShower(defs);   
  AddNueHistDefPIDs(defs);     
  AddNueHistDefForLid(defs);   
  AddNueHistDefForLem(defs);   
  AddNueHistDefForCvn(defs);   
  AddNueHistDefConfusion(defs);
  AddNueHistDefSelectionExtras(defs);
  AddNueHistDefTruth(defs);

  //Don't apply weights to next / all variables
  //const int kIdUnweight = -1;//defs.size();
  const int kIdUnweight = defs.size();
  AddNueHistDefWeight(defs);
  const Var kkWeight = kPPFXFluxCVWgt *  kXSecCVWgt2017;

  // 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_NueSecondAnaDQ",kNueSecondAnaDQ},
    {"2_NueNDSecondAnaDecaf",kNueNDSecondAnaDecafCut},
    {"3_NueNDSecondAnaPresel", kNueNDSecondAnaPresel},
    {"4_NueNDCVNSsb" ,kNueNDCVNSsb},
    {"5_NueNDCVNSb",kNueNDSecondAnaPresel && kNueSecondAnaCVNSb},
  };
  
  //Vars for a quick visual comparison of cutflow
  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(10,0,10);
  const Var kDummy ( [] (const caf::SRProxy *sr){return 1;});
  defs.insert(defs.begin(),{"_dummy", {"Dummy", Binning::Simple(3, -0.5,2.5), kDummy}});
 
  //Save Data and MC spectra

  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<Spectrum*> > specMC; //useful when truth info is missing
    std::vector < std::vector<IDecomp*> > preds;
    
    for(int selIdx = 0; selIdx<(int)sels.size(); ++selIdx){
      std::vector<Spectrum *> temp_spec;
      std::vector<Spectrum *> temp_spMC;
      std::vector<IDecomp *> temp_pred;
      for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
        const HistAxis& axis = defs[varIdx].axis;
        Var weight = (varIdx < kIdUnweight) ? kkWeight:kUnweighted;

        temp_spec.emplace_back(new Spectrum (loaderData, axis, sels[selIdx].cut, kNoShift));
        temp_spMC.emplace_back(new Spectrum   (loaderMC, axis, sels[selIdx].cut, kNoShift, weight)); 
        temp_pred.emplace_back(new CheatDecomp(loaderMC, axis, sels[selIdx].cut, kNoShift, weight)); 
      }
      spects.push_back(temp_spec);
      specMC.push_back(temp_spMC);
      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");
    Progress prog("Saving...");
    int nsteps=0;
    for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
      const char* cut = sels[selIdx].name.c_str();
      for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
        const char* name = defs[varIdx].name.c_str();
        spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
        specMC[selIdx][varIdx]->SaveTo(fout, TString::Format("speMC_%s_%s", name, cut));
        preds[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_%s_%s", name, cut));
        nsteps ++;
        prog.SetProgress((double)nsteps/(sels.size()*defs.size()+1));
      }
    }
    spfl_data->SaveTo(fout, "spect_data_cutflow");
    spfl_ndmc->SaveTo(fout, "spect_mc_cutflow");
    prog.Done();

    delete fout;
    std::cout << "Saved prediction to " << predFileName << "\n";
    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 = defs.size();

  // Don't make all the canvases right now
//  kNumSels = 2;
//  kNumVars = 2;
  bool components = 0;
  
  Progress prog ("Arranging spectra...");
  int nsteps = 0;

  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).Data());
      auto spMC = Spectrum::LoadFrom(predFile, TString::Format("speMC_%s_%s", name, cut).Data());
      auto pred = CheatDecomp::LoadFrom (predFile, TString::Format("pred_%s_%s", name, cut).Data());
      
      // Data and MC for variable "name" will be plotted together 
      // Plots will have selector for cut level {group=Cut,cat=%s}
      // Optional MCComponents folder will have extra selector; total MC = all_sum will be created as well

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

      if (components){
        TDirectory* comp = validFile1->GetDirectory("MCComponents");
        if(!comp) comp = validFile1->mkdir("MCComponents");
        pred->NueComponent().     SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=0_nue}", name,cut));
        pred->NumuComponent().    SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=1_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));
        pred->NCTotalComponent().      SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=2_nc}", name,cut));
        (*spMC-(pred->TotalMC())).   SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=5_other}", name,cut));
      } 
      prog.SetProgress((double)nsteps/(kNumSels*kNumVars));
      nsteps++;
    } //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
  prog.Done();
  std::cout << "Closing files...\n";
  delete validFile1;
  delete validFile2;
  std::cout << "Formatted outputs saved to " << validFileName1 << " and " << validFileName2 << std::endl;
  
  }
}