nue_fd_mc_validation.C

Go to the documentation of this file.

// 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

#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/MultiVar.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "OscLib/IOscCalc.h"
#include "OscLib/OscCalcDumb.h"
#include "CAFAna/Vars/GenieWeights.h"

using namespace ana;

#include "TFile.h"

#include <iostream>
#include <vector>

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

  std::vector <HistDef> defs;
  
  AddNueHistDefSecondAnaBasic(defs); 
  AddNueHistDefVertexFD(defs); 
  AddNueHistDefShowerFD(defs); 
  AddNueHistDefShower(defs);   
  AddNueHistDefPIDs(defs);     
  AddNueHistDefForLid(defs);   
  AddNueHistDefForLem(defs);   
  AddNueHistDefForCvn(defs);   
  AddNueHistDefConfusion(defs);
  AddNueHistDefSelectionExtras(defs);
  AddNueHistDefTruth(defs);
  // Variables that shouldn't be reweighted
  // Set to -1 for no reweights anywhere
  // const int kIdUnweight = -1;/
  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_NueFDSecondAnaDecaf",kNueFDSecondAnaDecafCut},
    {"3_NueSecondAnaFullPresel", kNueSecondAnaFullPresel},
    {"4_NueSAFDCVNSsb" ,kNueSAFDCVNSsb},
    {"5_NueSAFDCVNSb", kNueSecondAnaFullPresel && 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);
                          }
                          if(sr->mc.nnu==0)passes.push_back(sels.size()+1);
                          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}});

  if(makeFile){

    SpectrumLoader loaderNonSwap(ldrNonSwap);
    SpectrumLoader loaderFlxSwap(ldrFluxSwap);

    NullLoader loaderTauSwap;
    std::cerr << "WARNING: not trying to load tau swap" << std::endl;

    loaderNonSwap.SetSpillCut(kStandardSpillCuts);
    loaderFlxSwap.SetSpillCut(kStandardSpillCuts);
    loaderTauSwap.SetSpillCut(kStandardSpillCuts);
    
    std::vector < std::vector<Spectrum *> > other_nonswap;
    std::vector < std::vector<Spectrum *> > other_flxswap;
    std::vector < std::vector<IPrediction*> > preds;
    
    for(int selIdx = 0; selIdx<(int)sels.size(); ++selIdx){
      std::vector<Spectrum *> temp_nonswap;
      std::vector<Spectrum *> temp_flxswap;
      std::vector<IPrediction *> temp_pred;
      
      for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
        const HistAxis& axis = defs[varIdx].axis;
        Var weight = (varIdx < kIdUnweight) ? kkWeight : kUnweighted;

        auto sp_nonswap = new Spectrum(loaderNonSwap, axis, sels[selIdx].cut && !kHasNeutrino, kNoShift, weight);
        auto sp_flxswap = new Spectrum(loaderFlxSwap, axis, sels[selIdx].cut && !kHasNeutrino, kNoShift, weight);
        temp_nonswap.emplace_back(sp_nonswap);
        temp_flxswap.emplace_back(sp_flxswap);

        temp_pred.emplace_back
          (new PredictionNoExtrap(loaderNonSwap,loaderFlxSwap,loaderTauSwap,
                                  axis.GetLabels()[0],axis.GetBinnings()[0],axis.GetVars()[0], 
                                  sels[selIdx].cut, kNoShift, weight));
      }
      other_nonswap.push_back(temp_nonswap);
      other_flxswap.push_back(temp_flxswap);
      preds.push_back(temp_pred);
    }

    auto spfl_nonswap = new Spectrum ("Cut Flow",flowBins,loaderNonSwap, kCutFlowVar, kNoCut);
    auto spfl_flxswap = new Spectrum ("Cut Flow",flowBins,loaderFlxSwap, kCutFlowVar, kNoCut);
   
    loaderNonSwap.Go();
    loaderFlxSwap.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();

        other_nonswap[selIdx][varIdx]->SaveTo(fout, TString::Format("other_nonswap_%s_%s", name, cut));
        other_flxswap[selIdx][varIdx]->SaveTo(fout, TString::Format("other_flxswap_%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_nonswap->SaveTo(fout, "spect_nonswap_cutflow");
    spfl_flxswap->SaveTo(fout, "spect_flxswap_cutflow");
    prog.Done();
    delete fout;
    std::cout << "Saved prediction to " << predFileName << "\n";
    return;
  } //end makeFile
  else{
    
    TFile * predFile = new TFile (predFileName.c_str(), "READ");
    TString validFileName1 (predFileName); 
    validFileName1.ReplaceAll(".root","_mc_validation.root");
    TFile * validFile1 = new TFile(validFileName1,"RECREATE");
    
    int kNumSels = sels.size();
    int kNumVars = defs.size();
    // Don't make all the canvases and .json right now
//     kNumSels = 3;
    //   kNumVars = 5;

    bool components = 1;
    TDirectory* comp = 0;
    if(components) comp = validFile1->mkdir("MCComponents");
            
    Progress prog("Arranging spectra...");
    int nsteps =0;
  
    osc::IOscCalc * noosc = new osc::NoOscillations();
    osc::IOscCalc * dumb = new osc::OscCalcDumb();
    osc::IOscCalc * calc;
    const char* oscName;

    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 pred = PredictionNoExtrap::LoadFrom(
          predFile->GetDirectory(TString::Format("pred_%s_%s", name, cut)));
        auto other_nonswap= Spectrum::LoadFrom(
          predFile->GetDirectory(TString::Format("other_nonswap_%s_%s", name, cut)));
        auto other_flxswap = Spectrum::LoadFrom(
          predFile->GetDirectory(TString::Format("other_flxswap_%s_%s", name, cut)));
        auto other = *other_nonswap + *other_flxswap;

        for(int i=0;i<2;i++){
          if(i==0){
            calc = noosc; 
            oscName = "NoOscillations";
          }
          if(i==1) {
            calc = dumb;
            oscName="DumbOscillations";
          }
          
          (pred->Predict(calc)+other).SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Osc,cat=%s}", name,cut,oscName));
          
          if (components){

            pred->PredictComponent(calc, Flavors::kNuMuToNuE,Current::kCC,Sign::kBoth).
              SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=0_sig}", name,cut,oscName));
            pred->PredictComponent(calc, Flavors::kNuEToNuE,Current::kCC,Sign::kBoth).
              SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=1_beam}", name,cut,oscName));
            pred->PredictComponent(calc, Flavors::kAll,Current::kNC,Sign::kBoth).
              SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=2_nc}", name,cut,oscName));
            pred->PredictComponent(calc, Flavors::kAllNuMu,Current::kCC,Sign::kBoth).
              SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=3_numu}", name,cut,oscName));
            pred->PredictComponent(calc, Flavors::kAllNuTau,Current::kCC,Sign::kBoth).
              SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=4_tau}", name,cut,oscName));
            other.
              SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=5_other}", name,cut,oscName));
          }
          nsteps++;
          prog.SetProgress((double)nsteps/(kNumSels*kNumVars*2.));
        }
      } //end var
      
      for(auto oscName:{"NoOscillations","DumbOscillations"}){
        auto spfl = Spectrum::LoadFrom(predFile, "spect_nonswap_cutflow");
        spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Osc,cat=%s}", "_cutflow_nonswap",cut,oscName));
        spfl = Spectrum::LoadFrom(predFile, "spect_flxswap_cutflow");
        spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Osc,cat=%s}", "_cutflow_flxswap",cut,oscName));
      }
    }// end cut
    prog.Done();
    delete validFile1;
    std::cout << "Formatted outputs saved to " << validFileName1  << std::endl;

    
  }
}