make_nue_xsec_pred.C

Go to the documentation of this file.

#include "CAFAna/Analysis/PredictionWriter.h"
//#include "3FlavorAna/Ana2017/nue/SystVariations.h"
#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Systs/XSecSystLists.h"  // getAllXsecSysts_2017
#include "CAFAna/Vars/GenieWeights.h"  // kXSecCVWgt2017
#include "CAFAna/Vars/PPFXWeights.h" // kPPFXFluxCVWgt
#include "CAFAna/Analysis/Prod3Loaders.h"  // kNueConcat, kNueDecaf
#include "CAFAna/Vars/HistAxes.h"  // kNumuNonQEAxisFirstAna
//#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"  // kShwCalE, kHadCalE
#include "3FlavorAna/Vars/NueVars.h"  // kNueEnergy2017
#include "CAFAna/Vars/Vars.h"  // kCVNe
#include "3FlavorAna/Cuts/NueCuts2017.h"  // kNue2017Binning, kNue2017AnaBin
#include "CAFAna/Prediction/IPrediction.h"
#include "3FlavorAna/Prediction/PredictionGenerator3Flavor.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Extrap/ModularExtrap.h"
//#include "CAFAna/Core/Registry.h"
//#include "CAFAna/Systs/XSecSysts.h"
//#include "CAFAna/Analysis/Calcs.h"
//#include "CAFAna/Cuts/Cuts.h"
//#include "CAFAna/Cuts/NumuCuts.h"
//#include "CAFAna/Cuts/SpillCuts.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "TString.h"
#include "TFile.h"

#include <vector>
#include <iostream>
#include <string>

using namespace ana;

// Structure to hold weights --> only way to combine weights easily
struct Weight
{
  const TString name;
  Var weight;
};


void make_nue_xsec_pred(TString suffix="test", bool useConcat=true,
                        int startIdx=0, int stopIdx=-1)
{
  // Filename for output prediction file
  TString filename = "pred_xsec_fhc_" + suffix + ".root";

  // Make PredictionWriter instance
  PredictionWriter* writer = new PredictionWriter(filename);

  // Set analysis name
  TString analysisName = "pred_propDecomp";

  // Get vector of xsec systematics
  std::vector<const ISyst*> systs = getAllXsecSysts_2017();
  std::cout << "[make_nue_xsec_pred] Found " << systs.size()
    << " xsec systematics!" << std::endl;
  if (stopIdx < 0)  // assume all systematics, set as length of vector
    stopIdx = systs.size();
  std::cout << "[make_nue_xsec_pred] Processing systematics:" << std::endl;
  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
       (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
    std::cout << std::distance(systs.begin(), syst) << ": "
      << (*syst)->ShortName() << std::endl;

  // Set sigmas - only want 0 (no shift) once, so handle separately
  std::vector<int> sigmas = {-2, -1, 1, 2};

  // Create vector/map of weights
  std::vector<Weight> weights;
  weights.push_back({"ppfx_xsec", kPPFXFluxCVWgt * kXSecCVWgt2017});
  weights.push_back({"ppfx_only", kPPFXFluxCVWgt});

  // Setup Loaders
  ECAFType type;
  if (useConcat)
    type = kNueConcat;
  else
    type = kDecaf;
  Prod3NomLoaders loaders(type);

  // Define axes and energy variables for predictions
  const HistAxis axisNumu = kNumuNonQEAxisFirstAna;

  struct HistDef
  {
    const TString name;
    HistAxis axis;
  };
  std::vector<HistDef> hists;
//  hists.push_back({"shwCalE",
//                   {"Shower calE (GeV)",
//                    Binning::Simple(30.0, 0.0, 3.0), kShwCalE}});
//  hists.push_back({"hadCalE",
//                   {"calE_{had} (GeV)",
//                    Binning::Simple(30.0, 0.0,  3.0), kHadCalE}});
  //hists.push_back({"recoNuE",
//                   {"reco E_{#nu} (GeV)",
//                    Binning::Simple(50.0, 0.0,  5.0), kNueEnergy2017}});
  //hists.push_back({"id_cvne",
  //                 {"CVN e", Binning::Simple(40.0, 0.70, 1.1), kCVNe}});
  hists.push_back({"EnergyCVN2D",
                   {"Nue Energy / Selection Bin",
                    kNue2017Binning, kNue2017AnaBin}});  // peripheral binning

  // Create vector/map to store predictions
  std::vector<std::pair<IPrediction*, TString> > predictions;

  //TODO: add vector of pointers to Pred Generators
  //* No extrap
  //* proportional
  //* combo
  // try pred no extrap with numu cut for proof of principle

  // First do Nominal
  // Loop over axes and energy Vars
  for (auto hist:hists)
  {
    for (auto weight:weights)
    {
      TString path;  // In output prediction root file
      // Extrap first
      path = writer->ReadPrediction(
        analysisName, "Nominal", 0, true, hist.name, weight.name);

      auto extrapGen = NuePropExtrapGenerator(
        hist.axis, axisNumu, kNue2017FDAllSamples, kNue2017NDCVNSsb,
        kNumuND, kNoShift, weight.weight);
      auto predProp = extrapGen.Generate(loaders);  // no shift
      predictions.push_back({predProp.release(), path});

      // No Extrap
      path = writer->ReadPrediction(
        analysisName, "Nominal", 0, false, hist.name, weight.name);

      auto noExtrapGen = NoExtrapGenerator(
        hist.axis, kNue2017FDAllSamples, weight.weight);
      auto predNoExtrap = noExtrapGen.Generate(loaders);  // no shift
      predictions.push_back({predNoExtrap.release(), path});

    }// Loop over weights
  }// Loop over axes and energy Vars

  // Loop over systematics
  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
       (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
  {
    // Loop over sigmas
    for (auto sigma:sigmas)
    {
      SystShifts shift = SystShifts((*syst), static_cast<double>(sigma));
      // Loop over axes and energy Vars
      for (auto hist:hists)
      {
        // Loop over weights
        for (auto weight:weights)
        {
          TString path;  // In output prediction root file
          // Extrap first
          path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
                                        sigma, true, hist.name, weight.name);

          auto extrapGen = NuePropExtrapGenerator(
            hist.axis, axisNumu, kNue2017FDAllSamples, kNue2017NDCVNSsb,
            kNumuND,  kNoShift, weight.weight);
          auto predProp = extrapGen.Generate(loaders, shift);
          predictions.push_back({predProp.release(), path});

          // No extrap
          path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
                                        sigma, false, hist.name, weight.name);

          auto noExtrapGen = NoExtrapGenerator(
            hist.axis, kNue2017FDAllSamples, weight.weight);
          auto predNoExtrap = noExtrapGen.Generate(loaders, shift);
          predictions.push_back({predNoExtrap.release(), path});

        }// Loop over weights
      }// Loop over axes and energy Vars
    }// Loop over sigmas
  }// Loop over systNames

  loaders.Go();

  // Loop over again to write predictions
  std::vector<std::pair<IPrediction*, TString> >::iterator iPrediction =
    predictions.begin();

  // Do Nominal first
  // Loop over axes and energy Vars
  for (auto hist:hists)
  {
    // Loop over weights
    for (auto weight:weights)
    {
      // Write extrap prediction
      TString path;
      path = writer->ReadPrediction(
        analysisName, "Nominal", 0, true, hist.name, weight.name);

      // Check indexing
      if (!path.EqualTo(iPrediction->second))  // incorrect save location
        std::cout << "[make_nue_xsec_pred] WARNING: prediction path mismatch!"
          << std::endl;

      writer->WritePrediction(iPrediction->first, analysisName, "Nominal",
                              0, true, hist.name, weight.name);
      ++iPrediction;  // increment

      // Write no extrap prediction
      path = writer->ReadPrediction(
        analysisName, "Nominal", 0, false, hist.name, weight.name);

      // Check indexing
      if (!path.EqualTo(iPrediction->second))  // incorrect save location
        std::cout << "[make_nue_xsec_pred] WARNING: prediction path mismatch!"
          << std::endl;

      writer->WritePrediction(iPrediction->first, analysisName, "Nominal",
                              0, false, hist.name, weight.name);
      ++iPrediction;  // increment
    }// Loop over weights
  }// Loop over axes and energy Vars

  // Loop over systematics
  for (std::vector<const ISyst*>::iterator syst = systs.begin()+startIdx;
       (syst != systs.begin()+stopIdx) && (syst != systs.end()); ++syst)
  {
    // Loop over sigmas
    for (auto sigma:sigmas)
    {
      // Loop over axes and energy Vars
      for (auto hist:hists)
      {
        // Loop over weights
        for (auto weight:weights)
        {
          // Write extrap prediction
          TString path;
          path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
                                        sigma, true, hist.name, weight.name);

          // Check indexing
          if (!path.EqualTo(iPrediction->second))  // incorrect save locatio
            std::cout << "[make_nue_xsec_pred] WARNING: "
              << "prediction path mismatch!" << std::endl;

          writer->WritePrediction(iPrediction->first, analysisName,
                                  (*syst)->ShortName(), sigma,
                                  true, hist.name, weight.name);
          ++iPrediction;  // increment

          // Write no extrap prediction
          path = writer->ReadPrediction(analysisName, (*syst)->ShortName(),
                                        sigma, false, hist.name, weight.name);

          // Check indexing
          if (!path.EqualTo(iPrediction->second))  // incorrect save locatio
            std::cout << "[make_nue_xsec_pred] WARNING: "
              << "prediction path mismatch!" << std::endl;

          writer->WritePrediction(iPrediction->first, analysisName,
                                  (*syst)->ShortName(), sigma,
                                  false, hist.name, weight.name);
          ++iPrediction;
        }// Loop over weights
      }// Loop over axes and energy Vars
    }// Loop over sigmas
  }// Loop over systNames
}// make_nue_xsec_pred