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

#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "NuXAna/Cuts/NusCuts17.h"
#include "CAFAna/Vars/Vars.h"
#include "NuXAna/Vars/NusVars.h"
#include "NuXAna/Vars/NusVarsTemp.h"
#include "OscLib/IOscCalc.h"
#include "NuXAna/Analysis/NusPlotLists.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"

using namespace ana;

#include "TFile.h"

#include <iostream>
#include <vector>

void FDDataMC(const std::string predFileName,
              const bool makeFile=false, 
              const std::string ldrFDnon="",
              const std::string ldrFDswp="",
              const std::string ldrCosmic="")
{
  
  std::vector<HistDef> defs = getAllNusAna01FDHistDefs();
  
  // 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_Veto", kApplySecondAnalysisMask && kNus17Veto},
    {"1_NusEventQuality", kApplySecondAnalysisMask && kNus17Veto && kNus17EventQuality},
    {"2_NusPresel", kNus17FDPresel},
    {"3_kNusBackCut", kNus17FDPresel && kNus17BackwardCut},
    {"4_kNusCosRej", kNus17FDPresel && kNus17BackwardCut && kNus17NCCosRejAlt},
    {"5_kNusCVN", kNus17FDSansTimePtpAlt},
    {"6_kNusTime", kNus17FDSansTimePtpAlt && !(kNus17SlcTimeGap && kNus17SlcDist)},
    {"7_kNusFD", kNus17FDAlt},
    {"8_kNusFDNhit", kNus17FDExtraAlt}
  };
  
  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;});
  defs.insert(defs.begin(),{"_dummy", {"Dummy", Binning::Simple(3, -0.5,2.5), kDummy}});
  if(makeFile){
    SpectrumLoader loaderMCnon(ldrFDnon.c_str());
    SpectrumLoader loaderMCswp(ldrFDswp.c_str());
    SpectrumLoader loaderCosmic(ldrCosmic.c_str(), DataSource::kCosmic);
    loaderMCnon.SetSpillCut(kStandardSpillCuts);
    loaderMCswp.SetSpillCut(kStandardSpillCuts);
    loaderCosmic.SetSpillCut(kStandardSpillCuts);
    
    std::vector<std::vector<Spectrum*> > spects;
    std::vector<std::vector<Spectrum*> > spectsMC;
    std::vector<std::vector<IPrediction*> > preds;
    
    for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
      std::vector<Spectrum *>    temp_spec_non;
      std::vector<Spectrum *>    temp_specMC;
      std::vector<IPrediction *> temp_pred;
      
      for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
        const HistAxis& axis = defs[varIdx].axis;
        temp_spec_non.emplace_back(new Spectrum(loaderCosmic,
                                                axis,
                                                sels[selIdx].cut&&kInTimingSideband,
                                                kNoShift));

        temp_specMC.emplace_back(new Spectrum(loaderMCnon,
                                              axis,
                                              sels[selIdx].cut && kInBeamSpill,
                                              kNoShift,
                                              kXSecCVWgt2017*kPPFXFluxCVWgt));

        temp_pred.emplace_back(new PredictionNoExtrap(loaderMCnon,
                                                      loaderMCswp,
                                                      axis.GetLabels()[0],
                                                      axis.GetBinnings()[0],
                                                      axis.GetVars()[0],
                                                      sels[selIdx].cut&&kInBeamSpill,
                                                      kNoShift,
                                                      kXSecCVWgt2017*kPPFXFluxCVWgt));
      }
      spects.push_back(temp_spec_non);
      preds.push_back(temp_pred);
    }
    
    loaderCosmic.Go();
    loaderMCnon.Go(); 
    loaderMCswp.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)defs.size(); ++varIdx){
        const char* name = defs[varIdx].name.c_str();
        spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
        preds[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_%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","_fhc_mc_validation.root");
    TString validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_cosmic_validation.root");
    
    TFile * validFile1 = new TFile(validFileName1,"RECREATE");
    TFile * validFile2 = new TFile(validFileName2,"RECREATE");
    
    int kNumSels = sels.size();
    int kNumVars = defs.size();
    
    osc::IOscCalcAdjustable* calc = DefaultOscCalc();
    
    
    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 = PredictionNoExtrap::LoadFrom(predFile, TString::Format("pred_%s_%s", name, cut));
        
        spec->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}", name,cut));
        

        pred->PredictComponent(calc, ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth).       SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=6_mc_total}", name,cut));
        pred->PredictComponent(calc, ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth).         SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nc}", name,cut));
        pred->PredictComponent(calc, ana::Flavors::kNuEToNuE, ana::Current::kCC, ana::Sign::kNu).      SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_nue_surv}", name,cut));
        pred->PredictComponent(calc, ana::Flavors::kNuEToNuE, ana::Current::kCC, ana::Sign::kAntiNu).  SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_antinue_surv}", name,cut));
        pred->PredictComponent(calc, ana::Flavors::kNuMuToNuMu, ana::Current::kCC, ana::Sign::kNu).    SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_numu_surv}", name,cut));
        pred->PredictComponent(calc, ana::Flavors::kNuMuToNuMu, ana::Current::kCC, ana::Sign::kAntiNu).SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu_surv}", name,cut));
        
        
      }
    }
    std::cout << "Formatted outputs saved to "
              << validFileName1
              << " and "
              << validFileName2
              << std::endl;
    delete validFile1;
    delete validFile2;
    
  }
}