JointAna2018Systs.cxx

Go to the documentation of this file.

#include "3FlavorAna/Systs/JointAna2018Systs.h"

#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/MECSysts.h"
#include "CAFAna/Systs/RPASysts.h"
#include "CAFAna/Systs/BeamSysts.h"

#include "3FlavorAna/Systs/DummyRockScaleSyst.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/EnergySysts2017.h"
#include "3FlavorAna/Systs/EnergySysts2018.h"
#include "3FlavorAna/Systs/GeniePCASyst.h"
#include "3FlavorAna/Systs/MichelTaggingSyst.h"
#include "3FlavorAna/Systs/WrongSignSysts.h"

#include "3FlavorAna/Systs/DummySystStorage.h"
#include "3FlavorAna/Systs/NueAcceptSysts.h"


namespace ana
{
  std::string AnaTypeToString(const EAnaType2018 ana)
  {
    std::string ret = "";
    if(ana == kNueAna2018)   ret = "Nue2018";
    if(ana == kNumuAna2018)  ret = "Numu2018";
    if(ana == kJointAna2018) ret = "Joint2018";
    return ret;
  }

  std::vector<const ISyst *> getAna2018LargeXsecSysts(const EAnaType2018 ana)
  {
    std::vector<const ISyst *> ret =
      {
        &kMECEnuShapeSyst2018Nu,
        &kMECEnuShapeSyst2018AntiNu,
        &kMECShapeSyst2018Nu,
        &kMECShapeSyst2018AntiNu,
        &kMECInitStateNPFracSyst2018Nu,
        &kMECInitStateNPFracSyst2018AntiNu,
        &kMAQEGenieReducedSyst2018,
        &kRPARESSyst2018,
        &kRPACCQEEnhSyst2018,
        &kCOHCCScaleSyst2018,
        &kCOHNCScaleSyst2018,
        GetGenieKnobSyst(rwgt::fReweightMaCCRES),
        GetGenieKnobSyst(rwgt::fReweightMaNCRES),
        GetGenieKnobSyst(rwgt::fReweightMvCCRES),
        GetGenieKnobSyst(rwgt::fReweightCCQEPauliSupViaKF),
        //GetGenieKnobSyst(rwgt::fReweightFrElas_N),
        GetGenieKnobSyst(rwgt::fReweightFrCEx_N),
        GetGenieKnobSyst(rwgt::fReweightFrAbs_N),
        GetGenieKnobSyst(rwgt::fReweightFrInel_pi),
        GetGenieKnobSyst(rwgt::fReweightFormZone),
      };
    //nue only or joint
    if(ana != kNumuAna2018){
      ret.push_back(&kRadCorrNue);
      ret.push_back(&kRadCorrNuebar);
      ret.push_back(&k2ndClassCurrs);
    }
    //numu only or joint
    if(ana != kNueAna2018){
    }
    return ret;
  }

  
  void AddJointAna2018XSecSysts(std::vector<const ISyst*> &systs,
                                const EAnaType2018 ana,
                                bool smallgenie)
  {
    auto temp = getAna2018LargeXsecSysts(ana);
    for (auto const& s:temp){
      systs.push_back(s);
    }
    if(smallgenie){
      for (int i =0;i<5; ++i){
        systs.push_back(GetGeniePrincipals2018Small(i));
      }
      systs.push_back(&kTauScaleSyst);
    }
  }

  void AddJointAna2018BeamSysts (std::vector<const ISyst*> &systs,
                                 const EAnaType2018 ana)
  {
    for (int i =0;i<5; ++i){
      systs.push_back(GetFluxPrincipals2018(i));
    }
  }

  void AddJointAna2018FileSysts (std::vector<const ISyst*> &systs,
                                 const EAnaType2018 ana, const BeamType2018 beam){
    systs.push_back(&kAnaCalibrationSyst);
    systs.push_back(&kAnaRelativeCalibSyst);
    systs.push_back(&kAnaCalibShapeSyst);
    systs.push_back(&kAnaLightlevelSyst);
    systs.push_back(&kAnaCherenkovSyst); 
  }

  void AddJointAna2018NotCorrelSysts(std::vector<const ISyst*> &systs,
                                     const EAnaType2018 ana, const BeamType2018 beam)
  {

    if(ana != kNumuAna2018){
      if( beam != kRHC) {
        systs.push_back(&kNueAcceptSystSignalKin2018FHC);
        systs.push_back(&kNueAcceptSystBkg2018FHC);
        systs.push_back(&kMichelTaggingSyst2018);
      }
      if( beam != kFHC) {
        systs.push_back(&kNueAcceptSystSignalKin2018RHC);
        systs.push_back(&kNueAcceptSystBkg2018RHC);
      }
    }
    if(ana != kNueAna2018){
      systs.push_back(&kRelMuEScaleSyst2017);                   
    }

  }

 void AddJointAna2018OtherSysts(std::vector<const ISyst*> &systs,
                                 const EAnaType2018 ana, const BeamType2018 beam){

    systs.push_back(&kNeutronVisEScalePrimariesSyst2018);
    systs.push_back(&kMuEScaleSyst2017);
    if( beam != kRHC ) systs.push_back(&kAna2018NormFHC);                           // temp location
    if( beam != kFHC ) systs.push_back(&kAna2018NormRHC);                           // temp location

  }
 void AddNonLoadable2018Systs(std::vector<const ISyst*> &systs,
                                 const EAnaType2018 ana)
  {
    if(ana != kNumuAna2018){
       systs.push_back(&kRockScaleSyst);
     }
  }


  void AddWrongSignSysts(std::vector<const ISyst*> &systs,
                         const EAnaType2018 ana)
  {
      systs.push_back(&kWrongSignScale);
  }

//////////////////////////////////////////////////////////////////////
  std::vector<const ISyst*> getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
  {
    std::vector<const ISyst*> ret;
    
    AddJointAna2018XSecSysts (ret, ana, smallgenie);
    AddJointAna2018BeamSysts (ret, ana);
    AddJointAna2018FileSysts (ret, ana, beam);
    AddJointAna2018OtherSysts(ret, ana, beam);
    AddJointAna2018NotCorrelSysts(ret, ana, beam);
    if(isFit) AddNonLoadable2018Systs(ret, ana);
    
    return ret;
  }


  //////////////////////////////////////////////////////////////////////
  //test for CPT analysis
  //should ideally go into the std analyses
  std::vector<const ISyst*> getAllCPTSysts(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
  {
    std::vector<const ISyst*> ret;

    AddJointAna2018XSecSysts (ret, ana, smallgenie);
    AddJointAna2018BeamSysts (ret, ana);
    AddJointAna2018FileSysts (ret, ana, beam);
    AddJointAna2018OtherSysts(ret, ana, beam);
    AddJointAna2018NotCorrelSysts(ret, ana, beam);
    AddWrongSignSysts(ret, ana);
    if(isFit) AddNonLoadable2018Systs(ret, ana);

    return ret;
  }


  //////////////////////////////////////////////////////////////////////
  //don't really need these functions for analysis, but there are some parts in ana that call them
  
  std::vector<const ISyst *> getAna2018SmallXsecSysts(const EAnaType2018 ana)
  {
    auto allgenie = getAllXsecSysts_2018();
    auto largegenie = getAna2018LargeXsecSysts(ana);

    for (const auto & s:largegenie){
      for (const auto & a:allgenie){
        if( a == s ){
          auto it = std::remove(allgenie.begin(),allgenie.end(),a);
          allgenie.erase(it);
          continue;
        }
      }
    }
    return allgenie;
  }

  const SummedSyst * getAna2018SummedSmallXsecSysts(const EAnaType2018 ana)
  {
    auto systs =  getAna2018SmallXsecSysts(ana);
    auto anastr = AnaTypeToString(ana);
    return new SummedSyst(systs, "SumSmallXSec"+anastr,
                          "Summed small XSec systs "+anastr);
   }

  //////////////////////////////////////////////////////////////////////
  std::vector<const ISyst *> getAna2018NueOnly( const BeamType2018 beam)
  {
    std::vector<const ISyst*> ret;
    ret.push_back(&kRadCorrNue);
    ret.push_back(&kRadCorrNuebar);
    ret.push_back(&k2ndClassCurrs);
    AddJointAna2018NotCorrelSysts(ret, kNueAna2018, beam);
    AddNonLoadable2018Systs(ret, kNueAna2018);

    if(beam != kFHC) ret.push_back(&kAna2018NormRHC);
    if(beam != kRHC) ret.push_back(&kAna2018NormFHC);

    return ret;
  }

  std::vector<const ISyst *> getAna2018NumuOnly( const BeamType2018 beam)
  {
    std::vector<const ISyst*> ret;
    AddJointAna2018NotCorrelSysts(ret, kNumuAna2018, kBoth);
    AddNonLoadable2018Systs(ret, kNumuAna2018);
    return ret;
  }



}