NusSysts.cxx

Go to the documentation of this file.

#include "NuXAna/Systs/NusSysts.h"

#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/Vars.h"

#include "StandardRecord/Proxy/SRProxy.h"
#include "Utilities/func/MathUtil.h"

#include "TFile.h"
#include "TH1.h"

#include <cassert>
#include <cmath>


namespace ana
{
  //const NusSystFromHist kNusOscParamSysts  (kNusAna01SystFAlt, "EX", "OscParam",  "Oscillation Parameters");
  const NusSystFromHist kNusBeamSysts      (kNusAna01SystFile, "EX", "Beam",      "All Beam");
  const NusSystFromHist kNusBirksSyst      (kNusAna01SystFile, "EX", "Birks",     "Birks C");
  const NusSystFromHist kNusCalibFlatSyst  (kNusAna01SystFile, "EX", "CalFlat",   "Flat Miscalibration");
  const NusSystFromHist kNusCalibRelSyst   (kNusAna01SystFile, "EX", "CalRel",    "Relative Detector Miscalibration");
  const NusSystFromHist kNusCalibSlopeXSyst(kNusAna01SystFile, "EX", "CalSlopeX", "Sloped Miscalibration, X");
  const NusSystFromHist kNusCalibSlopeYSyst(kNusAna01SystFile, "EX", "CalSlopeY", "Sloped Miscalibration, Y");
  const NusSystFromHist kNusGENIESmallSysts(kNusAna01SystFile, "EX", "GENIESm",   "Summed small GENIE Systs");
  const NusSystFromHist kNusNDRockSyst     (kNusAna01SystFile, "EX", "NDNDRock",    "ND Rock");
  const NusSystFromHist kNusNueCCSyst      (kNusAna01SystFile, "EX", "NueCC",     "#nu_{e} CC Background");
  const NusSystFromHist kNusNumuCCSyst     (kNusAna01SystFile, "EX", "NumuCC",    "#nu_{#mu} CC Background");

  const NusSystFromHist kNusNDBeamSysts      (kNusAna01SystFile, "ND", "NDBeam",      "All Beam");
  const NusSystFromHist kNusNDBirksSyst      (kNusAna01SystFile, "ND", "NDBirks",     "Birks C");
  const NusSystFromHist kNusNDCalibFlatSyst  (kNusAna01SystFile, "ND", "NDCalFlat",   "Flat Miscalibration");
  const NusSystFromHist kNusNDCalibSlopeXSyst(kNusAna01SystFile, "ND", "NDCalSlopeX", "Sloped Miscalibration, X");
  const NusSystFromHist kNusNDCalibSlopeYSyst(kNusAna01SystFile, "ND", "NDCalSlopeY", "Sloped Miscalibration, Y");
  const NusSystFromHist kNusNDGENIESmallSysts(kNusAna01SystFile, "ND", "NDGENIESm",   "Summed small GENIE Systs");
  const NusSystFromHist kNusNDNDRockSyst     (kNusAna01SystFile, "ND", "NDRock",    "ND Rock");
  const NusSystFromHist kNusNDNueCCSyst      (kNusAna01SystFile, "ND", "NDNueCC",     "#nu_{e} CC Background");
  const NusSystFromHist kNusNDNumuCCSyst     (kNusAna01SystFile, "ND", "NDNumuCC",    "#nu_{#mu} CC Background");

  // Systematics applied as flat shifts
  const NusFlatSyst kNusMCStatsSyst("mcstat", "MC Stats", 2.0, 4.8);
  const NusFlatSyst kNusNDContSyst ("ndcont", "ND Containment", 1.0, 0.6);
  const NusFlatSyst kNusNormSyst   ("normNus","Normalization", 4.9, 4.9);

  //----------------------------------------------------------------------
  NusSystFromHist::NusSystFromHist(const std::string &fname,
                                   const std::string &det,
                                   const std::string &shortname,
                                   const std::string &latexname)
    : ISyst(shortname, latexname),
      fFileName(fname),
      fDet(det)
  {
  }

  //----------------------------------------------------------------------
  NusSystFromHist::~NusSystFromHist()
  {
  }

  //----------------------------------------------------------------------
  void NusSystFromHist::LoadHists() const
  {
    // Someone already called us
    if(!fHists.empty()) return;

    TFile fin(fFileName.c_str(), "read");
    if(fin.IsZombie()) {
      std::cout << "Warning: couldn't open " << fFileName
                << ". Crashing" << std::endl;
      abort();
    }

    std::vector<std::string> channels = {"NC", "BG"};

    const std::vector<int> sigmas = {-1, 0, +1};
    const std::vector<std::string> sigstr = {"-1", "0", "+1"};

    for(int i_chan = 0; i_chan < (int)channels.size(); ++i_chan) {
      std::vector< std::pair<int,TH1D*> > curHists;

      for(int i_sig = 0; i_sig < (int)sigmas.size(); ++i_sig) {
        // Find histogram for given NusChannel at given sigma
        std::string hName = TString::Format("h%s_%s_%s_%s",
                                            channels[i_chan].c_str(),
                                            fDet.c_str(),
                                            ShortName().c_str(),
                                            sigstr[i_sig].c_str()).Data();
        TH1D* h = (TH1D*)fin.Get(hName.c_str());
        if(!h) {
          // Surely it's a bad thing if we can't find nom, +-1 sigma hists
          std::cout << "Error: can't find necessary " << hName
                    << " histogram in file " << fFileName
                    << ". Crashing" << std::endl;
          abort();
        }

        h->SetDirectory(0);
        curHists.emplace_back(sigmas[i_sig], h);
      }

      fHists.push_back(curHists);
    }
  }

  //----------------------------------------------------------------------
  NusChannel NusSystFromHist::GetNusChannel(caf::SRProxy* sr) const
  {
    if(sr->mc.nnu == 0)    { return NusChannel::kNC; }
    if(!sr->mc.nu[0].iscc) { return NusChannel::kNC; }
    if( sr->mc.nu[0].iscc) { return NusChannel::kBG; }
    assert(0 && "Unknown Oscillation Channel");
  }

  //----------------------------------------------------------------------
  double NusSystFromHist::WeightFor(NusChannel chan,
                                    double sigma,
                                    double calE) const
  {
    LoadHists();

    const int bin = fHists[chan][0].second->FindBin(calE);

    int LowIdx = 0;
    if (sigma < fHists[chan].front().first)
      LowIdx = 0;
    else if (sigma >= fHists[chan].back().first)
      LowIdx = fHists[chan].size()-2;
    else{
      for (int i = 0; i < (int)fHists[chan].size()-1; i++){
        if (sigma >= fHists[chan][i].first){
          LowIdx = i;
          break;
        }
      }
    }

    // Why would we have templates differing by more than 1 sigma?
    // fracpart below assumes this
    assert(fHists[chan][LowIdx+1].first - fHists[chan][LowIdx].first == 1);

    const double fracpart = sigma - fHists[chan][LowIdx].first;
    const double ret = fracpart*   fHists[chan][LowIdx+1].second->GetBinContent(bin) +
                      (1-fracpart)*fHists[chan][LowIdx]  .second->GetBinContent(bin);

    return std::max(0., ret); // Keep the LL from blowing up
  }

  //----------------------------------------------------------------------
  void NusSystFromHist::Shift(double sigma,
                              caf::SRProxy* sr, double& weight) const
  {
    if(sr->hdr.det == caf::kFARDET &&
       fDet.compare("FD") != 0 &&
       fDet.compare("EX") != 0) { return; }
    if(sr->hdr.det == caf::kNEARDET &&
       fDet.compare("ND") != 0) { return; }

    // Find the energy estimate for this slice
    double calE = kCaloE(sr);

    // Find the signal or background oscillation channel
    NusChannel chan = GetNusChannel(sr);

    weight *= WeightFor(chan, sigma, calE);
  }

  //----------------------------------------------------------------------
  NusFlatSyst::NusFlatSyst(const std::string &shortname,
                           const std::string &latexname,
                           const double& ncWei,
                           const double& bgWei)
    : ISyst(shortname, latexname),
      fNCWei(ncWei/100.),
      fBGWei(bgWei/100.)
  {
  }

  //----------------------------------------------------------------------
  NusFlatSyst::~NusFlatSyst()
  {
  }

  //----------------------------------------------------------------------
  void NusFlatSyst::Shift(double sigma,
                          caf::SRProxy* sr, double& weight) const
  {
    // Takes in output from extrap - only shift FD spectra
    if(sr->hdr.det != caf::kFARDET) return;
    if(sr->mc.nnu == 0) return;
    if(sigma == 0) return;

    weight *= ((sr->mc.nu[0].iscc) ?
               1. + fBGWei*sigma :
               1. + fNCWei*sigma);

    return;
  }

  //----------------------------------------------------------------------
  /// Get a vector of all the nus group systs
  std::vector<const ISyst*> getAllNusSysts()
  {
    std::vector<const ISyst*> systs;

   // systs.push_back(&kNusOscParamSysts);
    systs.push_back(&kNusBeamSysts);
    systs.push_back(&kNusBirksSyst);
    systs.push_back(&kNusCalibFlatSyst);
    systs.push_back(&kNusCalibRelSyst);
    systs.push_back(&kNusCalibSlopeXSyst);
    systs.push_back(&kNusCalibSlopeYSyst);
    systs.push_back(&kNusGENIESmallSysts);
    systs.push_back(&kNusMCStatsSyst);
    systs.push_back(&kNusNDContSyst);
    systs.push_back(&kNusNDRockSyst);
    systs.push_back(&kNusNormSyst);
    systs.push_back(&kNusNueCCSyst);
    systs.push_back(&kNusNumuCCSyst);

    return systs;
  }
}