Nus17Systs.cxx

Go to the documentation of this file.

#include "NuXAna/Systs/Nus17Systs.h"

#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/Vars.h"
#include "NuXAna/Vars/NusVars.h"
#include "StandardRecord/Proxy/SRProxy.h"
#include "Utilities/func/MathUtil.h"

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

#include <cassert>
#include <cmath>


namespace ana
{
  //const Nus17SystFromHist kNus17OscParamSysts          (kNusAna01SystFAlt, "EX", "OscParam",  "Oscillation Parameters");
/*   const Nus17SystFromHist kNus17BeamSysts             (kNusAna17SystFile, "EX", "Beam17",      "All Beam");
   const Nus17SystFromHist kNus17CherenkovSyst    (kNusAna17SystFile, "EX", "Cherenkov17",     "Cherenkov");
   const Nus17SystFromHist kNus17Cherenkov_mSyst  (kNusAna17SystFile, "EX", "Cherenkov_m17",     "Light level method");
   const Nus17SystFromHist kNus17CalibFlatSyst           (kNusAna17SystFile, "EX", "CalFlat17",   "Flat Miscalibration");
   const Nus17SystFromHist kNus17CalibRelFDSyst   (kNusAna17SystFile, "EX", "CalRelFD17",    "Relative Detector FD Miscalibration");
   const Nus17SystFromHist kNus17CalibShapeSyst  (kNusAna17SystFile, "EX", "CalShape17", "Calibration shape");
   const Nus17SystFromHist kNus17GENIESmallSysts         (kNusAna17SystFile, "EX", "GENIESm17",   "Summed small GENIE Systs");
   const Nus17SystFromHist kNus17LightLSyst              (kNusAna17SystFile, "EX", "CalLight17",    "Lightlevel syst");
   const Nus17SystFromHist kNus17NueCCSyst       (kNusAna17SystFile, "EX", "NueCC17",     "#nu_{e} CC Background");
   const Nus17SystFromHist kNus17NumuCCSyst      (kNusAna17SystFile, "EX", "NumuCC17",    "#nu_{#mu} CC Background");
   const Nus17SystFromHist kNus17NCSyst          (kNusAna17SystFile, "EX", "NC17",    "NC Background");
   const Nus17SystFromHist kNus17AccSyst          (kNusAna17SystFile, "EX", "NDAcc17",    "NC Background"); //Caveat: Latex name is "acceptance syst". To match with tat of ISyst object created, the wrong Latex name is retained. But it is not troublesome. Will update soon.
   const Nus17SystFromHist kNus17CalibRelNDSyst   (kNusAna17SystFile, "EX", "CalRelND17",    "Relative Detector NDMiscalibration");
*/

/*  const Nus17SystFromHist kNus17NDBeamSysts        (kNusAna17SystFileND, "ND", "NDBeam17",      "All Beam");
  const Nus17SystFromHist kNus17NDCherenkovSyst    (kNusAna17SystFileND, "ND", "NDCherenkov17",     "Cherenkov");
  const Nus17SystFromHist kNus17NDCalibFlatSyst    (kNusAna17SystFileND, "ND", "NDCalFlat17",   "Flat Miscalibration");
  const Nus17SystFromHist kNus17NDCalibShapeSyst   (kNusAna17SystFileND, "ND", "NDCalShape17", "Calibration shape");
  const Nus17SystFromHist kNus17NDGENIESmallSysts  (kNusAna17SystFileND, "ND", "NDGENIESm17",   "Summed small GENIE Systs");
  const Nus17SystFromHist kNus17NDCalibRelSyst     (kNusAna17SystFileND, "ND", "NDCalRel17",    "Relative systematics");
  const Nus17SystFromHist kNus17NDNueCCSyst        (kNusAna17SystFileND, "ND", "NDNueCC17",     "#nu_{e} CC Background");
  const Nus17SystFromHist kNus17NDNumuCCSyst       (kNusAna17SystFileND, "ND", "NDNumuCC17",    "#nu_{#mu} CC Background");
  const Nus17SystFromHist kNus17NDNCSyst           (kNusAna17SystFileND, "ND", "NDNC17",    "NC Background");  
  const Nus17SystFromHist kNus17NDCherenkov_mSyst  (kNusAna17SystFileND, "ND", "NDCherenkov_m17",     "Light level method");
  const Nus17SystFromHist kNus17NDLightLSyst       (kNusAna17SystFileND, "ND", "NDCalLight17",    "Lightlevel syst");
*/

  // const Nus17SystFromHist kNus17FDBeamSysts        (kNusAna17SystFileND, "FD", "NDBeam17",      "All Beam");
  // const Nus17SystFromHist kNus17FDCherenkovSyst    (kNusAna17SystFileND, "FD", "NDCherenkov17",     "Cherenkov");
  // const Nus17SystFromHist kNus17FDCalibFlatSyst    (kNusAna17SystFileND, "FD", "NDCalFlat17",   "Flat Miscalibration");
  // const Nus17SystFromHist kNus17FDCalibShapeSyst   (kNusAna17SystFileND, "FD", "NDCalShape17", "Calibration shape");
  // const Nus17SystFromHist kNus17FDGENIESmallSysts  (kNusAna17SystFileND, "FD", "NDGENIESm17",   "Summed small GENIE Systs");
  // const Nus17SystFromHist kNus17FDCalibRelSyst     (kNusAna17SystFileND, "FD", "NDCalRel17",    "Relative systematics");
  // const Nus17SystFromHist kNus17FDCherenkov_mSyst  (kNusAna17SystFileND, "FD", "NDCherenkov_m17",     "Light level method");
  // const Nus17SystFromHist kNus17FDLightLSyst       (kNusAna17SystFileND, "FD", "NDCalLight17",    "Lightlevel syst");




  // 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 Nus17FlatSyst kNus17NormSyst   ("normNus17","Normalization", 4.9, 4.9);
   const Nus17FlatSyst kNus17OscParamSysts  ("oscparm", "Oscillation_Parameter",0.7,10.7);
   const Nus17FlatSyst kNusFDPPFXSyst  ("ppfx_fd","Fluxweight_fd",0.96,2.8);

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

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

  //----------------------------------------------------------------------
  void Nus17SystFromHist::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 Nus17SystFromHist::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 Nus17SystFromHist::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 Nus17SystFromHist::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 = kNus17Energy(sr);


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

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

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

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

  //----------------------------------------------------------------------
  void Nus17FlatSyst::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*> getAllNusFDSysts()
   {
     std::vector<const ISyst*> systs;

     systs.push_back(&kNus17OscParamSysts);
     systs.push_back(&kNusFDPPFXSyst);
     systs.push_back(&kNus17BeamSysts);
     systs.push_back(&kNus17CherenkovSyst);
     systs.push_back(&kNus17CalibFlatSyst);
     systs.push_back(&kNus17CalibRelNDSyst);
     systs.push_back(&kNus17CalibRelFDSyst);
     systs.push_back(&kNus17CalibShapeSyst);
     systs.push_back(&kNus17Cherenkov_mSyst);
     systs.push_back(&kNus17GENIESmallSysts);
     systs.push_back(&kNus17LightLSyst);
     systs.push_back(&kNus17NumuCCSyst);
     systs.push_back(&kNus17NCSyst);
     systs.push_back(&kNus17AccSyst);
     return systs;
   }*/
//..............................................................................................
/* std::vector<const ISyst*> getAllNusNDSysts()
  {
    std::vector<const ISyst*> systs;

    systs.push_back(&kNus17OscParamSysts);
    systs.push_back(&kNus17NDBeamSysts);
    systs.push_back(&kNus17NDCherenkovSyst);
    systs.push_back(&kNus17NDCalibFlatSyst);
    systs.push_back(&kNus17NDCalibRelSyst);
    systs.push_back(&kNus17NDCalibShapeSyst);
    systs.push_back(&kNus17NDCherenkov_mSyst);
    systs.push_back(&kNus17NDGENIESmallSysts);
    systs.push_back(&kNus17NDLightLSyst);
    systs.push_back(&kNus17NDNueCCSyst);
    systs.push_back(&kNus17NDNumuCCSyst);
    systs.push_back(&kNus17NDNCSyst);
    systs.push_back(&kNus17NormSyst);
    return systs;
  }*/
 //................................................................................................................
//    std::vector<const ISyst*> getAllNusFDRSysts()
//   {
//     std::vector<const ISyst*> systs;

//     systs.push_back(&kNus17OscParamSysts);
//     systs.push_back(&kNus17FDBeamSysts);
//     systs.push_back(&kNus17FDCherenkovSyst);
//     systs.push_back(&kNus17FDCalibFlatSyst);
//     systs.push_back(&kNus17FDCalibRelSyst);
//     systs.push_back(&kNus17FDCalibShapeSyst);
//     systs.push_back(&kNus17FDCherenkov_mSyst);
//     systs.push_back(&kNus17FDGENIESmallSysts);
//     systs.push_back(&kNus17FDLightLSyst);
//     systs.push_back(&kNus17NormSyst);
//     return systs;
// } 


}