NusAna2020Systs.cxx

Go to the documentation of this file.

#include "CAFAna/Core/Registry.h"
#include "CAFAna/Core/Sample.h"
#include "CAFAna/Vars/TruthVars.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/BeamSysts.h"
#include "3FlavorAna/Systs/CovMxSysts.h"
#include "3FlavorAna/Systs/EnergySysts2018.h"
#include "3FlavorAna/Systs/EnergySysts2020.h"
#include "NuXAna/Systs/CovMxSysts.h"
#include "NuXAna/Systs/NusAna2020Systs.h"

#include "TFile.h"
#include "TKey.h"

#include <cmath>
#include <cassert>
#include <iostream>

using std::cout;
using std::endl;
using std::string;
using std::vector;

using ana::covmx::Sample;

namespace ana
{

  const covmx::NormSyst kNusAna2020NormPOT(           "POT",           "POT uncertainty",           0.0055 );
  const covmx::NormSyst kNusAna2020NormNDMass(        "NDMass",        "Near detector mass",        0.003  );
  const covmx::NormSyst kNusAna2020NormFDMass(        "FDMass",        "Far detector mass",         0.003  );
  const covmx::NormSyst kNusAna2020NormMissingLepton( "MissingLepton", "Missing lepton bug",        0.008  );
  const covmx::NormSyst kNusAna2020NormNDNumuOverlay( "NDNumuOverlay", "ND numu selection overlay", 0.0021 );
  const covmx::NormSyst kNusAna2020NormNDNCOverlay(   "NDNCOverlay",   "ND NC selection overlay",   0.01   );

  const NusAna2020KaonSyst kNusAna2020KaonSyst;
  const NusAna2020TauSyst kNusAna2020TauSyst;

  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DNorm_1(    kGauss2DNorm_1Nux   , "Norm_1Nux"     );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ0_1(  kGauss2DMeanQ0_1Nux , "MeanQ0_1Nux"   );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ3_1(  kGauss2DMeanQ3_1Nux , "MeanQ3_1Nux"   );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ0_1( kGauss2DSigmaQ0_1Nux, "SigmaQ0_1Nux"  );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ3_1( kGauss2DSigmaQ3_1Nux, "SigmaQ3_1Nux"  );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DCorr_1(    kGauss2DCorr_1Nux   , "Corr_1Nux"     );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DNorm_2(    kGauss2DNorm_2Nux   , "Norm_2Nux"     );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ0_2(  kGauss2DMeanQ0_2Nux , "MeanQ0_2Nux"   );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ3_2(  kGauss2DMeanQ3_2Nux , "MeanQ3_2Nux"   );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ0_2( kGauss2DSigmaQ0_2Nux, "SigmaQ0_2Nux"  );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ3_2( kGauss2DSigmaQ3_2Nux, "SigmaQ3_2Nux"  );
  const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DCorr_2(    kGauss2DCorr_2Nux   , "Corr_2Nux"     );
  const MECDoubleGaussEnhSystNux kNusAna2020MECBaseline(         kBaselineNux        , "Baseline_2Nux" );

  string InputPath() {
    string ret;
    if (kUseCVMFS) {
      const char* cpath = getenv("NUSDATA_DIR");
      if (!cpath) 
        throw std::runtime_error("nusdata directory not found! Exiting.");
      ret = string(cpath) + "/nus20";
    } else {
      ret = kNus20Path;
    }
    if (kUseXRootD) ret = pnfs2xrootd(ret);
    return ret;
  }

  // Check whether a key syst is one-sided
  bool IsOneSided(string systKey) {
    for (string key : kOneSidedSysts) {
      if (systKey.find(key) != string::npos) return true;
    }
    return false;
  } // function IsOneSided

  // Check whether a syst should be uncorrelated between detectors
  bool IsDetUncorrelated(string systKey) {
    for (string key : kUncorrelatedSysts) {
      if (systKey.find(key) != string::npos) return true;
    }
    return false;
  } // function IsDetUncorrelated

  vector<string> GetKeySystNames() {
    vector<string> ret;
    for (string name : kDetectorSystNames) {
      if (IsDetUncorrelated(name)) {
        ret.push_back("key_neardet_"+name);
        ret.push_back("key_fardet_"+name);
      } else {
        ret.push_back("key_"+name);
      }
    }
    return ret;
  } // function GetKeySystNames

  vector<const KeySyst*> GetKeySysts() {
    vector<const KeySyst*> ret;
    vector<string> keys = GetKeySystNames();
    for (string key : keys) {
      const KeySyst* keySyst = nullptr;
      const ISyst* syst = Registry<ISyst>::ShortNameToPtr(key, true);
      if (!syst) {
        keySyst = new KeySyst(key, key);
      } else {
        keySyst = dynamic_cast<const KeySyst*>(syst);
      }
      ret.push_back(keySyst);
    }
    return ret;
  } // function GetKeySysts

  // This function should be in CovMxSysts after the CommonAna migration
  vector<const ISyst*> LoadSystsFromFile(string filePath, string dirName)
  {

    DontAddDirectory guard;

    // Open up input directory
    TFile* inFile = TFile::Open(filePath.c_str(), "read");
    TDirectory* dir = inFile->GetDirectory(dirName.c_str());

    // Loop over everything in this directory and load it as a systematic
    std::vector<const ISyst*> systs;
    TIter next(dir->GetListOfKeys());
    TKey* key;
    while ((key = (TKey*)next())) {
      // Check if this systematic has already been loaded
      std::string shortName = ((TObjString*)((TDirectory*)key->ReadObj())->Get("name"))
        ->GetString().Data();

      const ISyst* syst = Registry<ISyst>::ShortNameToPtr(shortName, true);
      // If it's already been loaded, just return the existing version
      if (syst) {
        systs.push_back(syst);
      } else {
        // If it hasn't already been loaded, load it
        cout << "Syst loading code works, but requires an overhaul." << endl;
        TObjString* str = (TObjString*)((TDirectory*)key->ReadObj())->Get("type");
        if (str->GetString() == "NueSyst")
          systs.push_back(NueSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
        else if (str->GetString() == "NumuSyst")
          systs.push_back(NumuSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
        else if (str->GetString() == "NCSyst")
          systs.push_back(NCSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
        else
          throw std::runtime_error("Syst " + str->GetString() + " not recognised!");
        //systs.push_back(NuISyst::LoadFrom((TDirectory*)key->ReadObj()).release());
      }
    }

    inFile->Close();
    return systs;

  } // function LoadSystsFromFile

  void AddNusAna2020NormSysts(vector<const ISyst*> &systs, Sample& s)
  {
    systs.push_back(&kNusAna2020NormPOT);
    systs.push_back(&kNusAna2020NormMissingLepton);
    if (s.detector == covmx::kNearDet) {
      systs.push_back(&kNusAna2020NormNDMass);
      if (s.selection == covmx::kCCNumu)
        systs.push_back(&kNusAna2020NormNDNumuOverlay);
      if (s.selection == covmx::kNC)
        systs.push_back(&kNusAna2020NormNDNCOverlay);
    }
    if (s.detector == covmx::kFarDet)
      systs.push_back(&kNusAna2020NormFDMass);
  }

  void AddNusAna2020NormSysts(vector<const ISyst*> &systs)
  {
    systs.push_back(&kNusAna2020NormPOT);
    systs.push_back(&kNusAna2020NormMissingLepton);
    systs.push_back(&kNusAna2020NormNDMass);
    systs.push_back(&kNusAna2020NormFDMass);
    systs.push_back(&kNusAna2020NormNDNumuOverlay);
    systs.push_back(&kNusAna2020NormNDNCOverlay);
  }

  void AddNusAna2020FileSysts(vector<const ISyst*> &systs, Sample& s)
  {
    // Open systfile and get directory
    string filePath = InputPath() + "/systs/isysts.root";
    string dirName = "isysts_" + s.GetTag();

    for (const ISyst* syst : LoadSystsFromFile(filePath, dirName)) {
      systs.push_back(syst);
    }
  }

  void AddNusAna2020XSecSysts(vector<const ISyst*> &systs)
  {
    for (const ISyst* syst : getAllXsecNuTruthSysts_2020()) {
      // Skip MEC systs that are tuned to ND data
      if (syst->ShortName() == "MECShape2020Nu" ||
        syst->ShortName() == "MECShape2020AntiNu") continue;
      systs.push_back(syst);
    }
  }

  void AddNusAna2020MECSysts ( vector<const ISyst*> &systs)
  {
    systs.push_back(&kNusAna2020MECGauss2DNorm_1);
    systs.push_back(&kNusAna2020MECGauss2DMeanQ0_1);
    systs.push_back(&kNusAna2020MECGauss2DMeanQ3_1);
    systs.push_back(&kNusAna2020MECGauss2DSigmaQ0_1);
    systs.push_back(&kNusAna2020MECGauss2DSigmaQ3_1);
    systs.push_back(&kNusAna2020MECGauss2DCorr_1);
    systs.push_back(&kNusAna2020MECGauss2DNorm_2);
    systs.push_back(&kNusAna2020MECGauss2DMeanQ0_2);
    systs.push_back(&kNusAna2020MECGauss2DMeanQ3_2);
    systs.push_back(&kNusAna2020MECGauss2DSigmaQ0_2);
    systs.push_back(&kNusAna2020MECGauss2DSigmaQ3_2);
    systs.push_back(&kNusAna2020MECGauss2DCorr_2);
    systs.push_back(&kNusAna2020MECBaseline);
  }

  void AddNusAna2020BeamSysts (vector<const ISyst*> &systs)
  {
    // keeping the first 5 flux PCs for 2020
    for (int i = 0; i < 5; ++i)
    {
      systs.push_back(GetFluxPrincipals2020(i));
    }
  }

  void AddNusAna2020KaonSyst(vector<const ISyst*> &systs)
  {
    systs.push_back(&kNusAna2020KaonSyst);
  }

  void AddNusAna2020NeutronSyst(vector<const ISyst*> &systs)
  {
    systs.push_back(&kNeutronVisEScalePrimariesSyst2018);
  }

  void AddNusAna2020TauSyst(vector<const ISyst*> &systs)
  {
    systs.push_back(&kNusAna2020TauSyst);
  }

  void AddNusAna2020NumuSysts(vector<const ISyst*> &systs)
  {
    systs.push_back(&kCorrMuEScaleSyst2020);
    systs.push_back(&kUnCorrMuCatMuESyst2020);
    systs.push_back(&kUnCorrNDMuEScaleSyst2020);
    systs.push_back(&kUnCorrFDMuEScaleSyst2020);
    systs.push_back(&kPileupMuESyst2020);
  }

  //////////////////////////////////////////////////////////////////////
  vector<const ISyst*> getNusAna2020AllSysts(covmx::Sample& s)
  {
      
    vector<const ISyst*> ret;
    AddNusAna2020NormSysts(ret, s);
    AddNusAna2020FileSysts(ret, s);
    AddNusAna2020XSecSysts(ret);
    AddNusAna2020BeamSysts(ret);
    AddNusAna2020MECSysts(ret);
    AddNusAna2020KaonSyst(ret);
    AddNusAna2020NeutronSyst(ret);
    AddNusAna2020TauSyst(ret);
    if (s.selection == covmx::kCCNumu) AddNusAna2020NumuSysts(ret);
       
    return ret;
  }

  //////////////////////////////////////////////////////////////////////
  vector<const ISyst*> getNusAna2020AllSysts()
  {
      
    vector<const ISyst*> ret;

    AddNusAna2020NormSysts(ret);
    for (const KeySyst* syst : GetKeySysts())
      ret.push_back(dynamic_cast<const ISyst*>(syst));
    AddNusAna2020XSecSysts(ret);
    AddNusAna2020BeamSysts(ret);
    AddNusAna2020MECSysts(ret);
    AddNusAna2020KaonSyst(ret);
    AddNusAna2020NeutronSyst(ret);
    AddNusAna2020TauSyst(ret);
    AddNusAna2020NumuSysts(ret);
       
    return ret;
  }

  // Set systematic aliases for a given sample
  void SetSystAlias(Sample& sample) {

    // Mask off normalisation systs for samples that don't apply
    if (sample.detector == covmx::kNearDet)
      sample.SetSystAlias(&kNusAna2020NormFDMass, nullptr);
    if (sample.detector == covmx::kFarDet)
      sample.SetSystAlias(&kNusAna2020NormNDMass, nullptr);
    if (sample.detector == covmx::kFarDet || sample.selection == covmx::kCCNumu)
      sample.SetSystAlias(&kNusAna2020NormNDNCOverlay, nullptr);
    if (sample.detector == covmx::kFarDet || sample.selection == covmx::kNC)
      sample.SetSystAlias(&kNusAna2020NormNDNumuOverlay, nullptr);

    // Loop over systematic names
    for (string systName : kDetectorSystNames) {
      // Get systematic for this sample

      // If this is detector-uncorrelated, then we need the key for
      // THIS detector to work, and the key for the OTHER detector
      // to be a null pointer
      if (IsDetUncorrelated(systName)) {
        // ND syst
        string keyND = "key_neardet_"+systName; // ND key
        const ISyst* keySyst = Registry<ISyst>::ShortNameToPtr(keyND, true);
        if (!keySyst) keySyst = new KeySyst(keyND, keyND);
        const ISyst* childSyst;
        if (sample.detector == covmx::kNearDet) {
          childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
        } else {
          childSyst = nullptr;
        }
        sample.SetSystAlias(keySyst, childSyst);

        std::string keyFD = "key_fardet_"+systName; // FD key
        keySyst = Registry<ISyst>::ShortNameToPtr(keyFD, true);
        if (!keySyst) keySyst = new KeySyst(keyFD, keyFD);
        if (sample.detector == covmx::kFarDet) {
          childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
        } else {
          childSyst = nullptr;
        }
        sample.SetSystAlias(keySyst, childSyst);
      } else {
        // Get key systematic
        const ISyst* childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
        std::string key = "key_"+systName;
        const ISyst* keySyst = Registry<ISyst>::ShortNameToPtr(key, true);
        if (!keySyst) keySyst = new KeySyst(key, key);
        // Set alias from key to child
        sample.SetSystAlias(keySyst, childSyst); // This should automatically set nullptr if the key is true
      }
    } // for syst name

    // Turn off numu systs for non-numu samples
    if (sample.selection != covmx::kCCNumu) {
      vector<const ISyst*> numuSysts;
      AddNusAna2020NumuSysts(numuSysts);
      for (const ISyst* keySyst : numuSysts) sample.SetSystAlias(keySyst, nullptr);
    }

  } // function SetSystAlias

}