3FlavorAna2020Systs.cxx

Go to the documentation of this file.

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

#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Systs/XSecSysts.h"
#include "CAFAna/Systs/MECSysts.h"
#include "CAFAna/Systs/RPASysts.h"
#include "CAFAna/Systs/FSISysts.h"
#include "CAFAna/Systs/RESSysts.h"
#include "CAFAna/Systs/DISSysts.h"
#include "CAFAna/Systs/BeamSysts.h"

#include "CAFAna/Experiment/SingleSampleExperiment.h" // kCosmicBkgScale

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

#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/EnergySysts2018.h" // kNeutronVisEScalePrimariesSyst2018
#include "3FlavorAna/Systs/EnergySysts2020.h"

#include "3FlavorAna/Systs/AngleSysts.h"
#include "3FlavorAna/Systs/GeniePCASyst.h"
#include "3FlavorAna/Systs/MichelTaggingSyst.h"
#include "3FlavorAna/Systs/WrongSignSysts.h"
#include "3FlavorAna/Systs/NueAcceptSysts.h"


namespace ana
{
  std::string AnaTypeToString(const EAnaType2020 ana)
  {
    std::string ret = "";
    if(ana == EAnaType2020::kNueAna)   ret = "Nue2020";
    if(ana == EAnaType2020::kNumuAna)  ret = "Numu2020";
    if(ana == EAnaType2020::k3FlavorAna) ret = "3Flavor2020";
    return ret;
  }

  std::vector<const ISyst *> get3FlavorAna2020LargeXsecSysts(const EAnaType2020 ana)
  {
    std::vector<const ISyst *> ret =
    {
      &kZExpEV1Syst2020,
      &kZExpEV2Syst2020,
      &kZExpEV3Syst2020,
      &kZExpEV4Syst2020,

      GetGenieKnobSyst(rwgt::fReweightMaCCRES),
      GetGenieKnobSyst(rwgt::fReweightMvCCRES),
      GetGenieKnobSyst(rwgt::fReweightMaNCRES),
      GetGenieKnobSyst(rwgt::fReweightMvNCRES),

      GetGenieKnobSyst(rwgt::fReweightZNormCCQE),

      &kRPACCQEEnhSyst2020,
      &kRPACCQESuppSyst2020,

      &kRESLowQ2SuppressionSyst2020,

      &kDISvnCC1pi_2020,

      &khNFSISyst2020_MFP,
      &khNFSISyst2020_EV1,

      &kMECEnuShapeSyst2020Nu,
      &kMECEnuShapeSyst2020AntiNu,
      &kMECShapeSyst2020Nu,
      &kMECShapeSyst2020AntiNu,
      &kMECInitStateNPFracSyst2020Nu,
      &kMECInitStateNPFracSyst2020AntiNu,
    };

    //nue only or joint
    if(ana != EAnaType2020::kNumuAna)
    {
      ret.push_back(&kRadCorrNue);
      ret.push_back(&kRadCorrNuebar);
      ret.push_back(&k2ndClassCurrs);
    }
    //numu only or joint
    if(ana != EAnaType2020::kNueAna)
    {
    }
    return ret;
  }


  void Add3FlavorAna2020XSecSysts(std::vector<const ISyst*> &systs,
                                const EAnaType2020 ana,
                                bool smallgenie)
  {
    auto temp = get3FlavorAna2020LargeXsecSysts(ana);
    for (auto const& s:temp)
    {
      systs.push_back(s);
    }
    if(smallgenie)
    {
      // 12 Genie PCs for 2020
      for (int i =0;i<12; ++i)
      {
        systs.push_back(GetGeniePrincipals2020Small(i));
      }
    }
    systs.push_back(&kTauScaleSyst);
  }

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

  void Add3FlavorAna2020CalibSysts( std::vector<const ISyst*>& systs )
  {
    systs.push_back( &kAnaCalibrationSyst );
    systs.push_back( &kAnaRelativeCalibSyst );
    systs.push_back( &kAnaCalibShapeSyst );
    systs.push_back( &kAnaCalibDriftSyst );
  }

  void Add3FlavorAna2020LightSysts( std::vector<const ISyst*>& systs )
  {
    systs.push_back( &kAnaLightlevelFDSyst );
    systs.push_back( &kAnaLightlevelNDSyst );
    systs.push_back( &kAnaCherenkovSyst );
  }

  void Add3FlavorAna2020FileSysts( std::vector<const ISyst*> &systs )
  {
    Add3FlavorAna2020CalibSysts( systs );
    Add3FlavorAna2020LightSysts( systs );
  }

  void Add3FlavorAna2020MuEnergySysts( std::vector<const ISyst*> &systs,
                                       const EAnaType2020 ana )
  {
    systs.push_back( &kUnCorrNDMuEScaleSyst2020 );
    systs.push_back( &kUnCorrMuCatMuESyst2020 );
    systs.push_back( &kPileupMuESyst2020 );
    systs.push_back( &kCorrMuEScaleSyst2020 );
    if ( ana != EAnaType2020::kNueAna )
    {
      systs.push_back( &kUnCorrFDMuEScaleSyst2020 );
    }
  }

  void Add3FlavorAna2020LeptonAngleSysts( std::vector<const ISyst*> &systs,
                                          const bool ptExtrap )
  {
    if ( ptExtrap )
    {
      systs.push_back( &kLeptonAngleSystNDXZ2020 );
      systs.push_back( &kLeptonAngleSystNDYZ2020 );
      systs.push_back( &kLeptonAngleSystFDXZ2020 );
      systs.push_back( &kLeptonAngleSystFDYZ2020 );
    }
  }

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

  void Add3FlavorAna2020NueAcceptSysts( std::vector<const ISyst*> &systs,
                                        const EAnaType2020 ana,
                                        const BeamType2020 beam,
                                        const bool ptExtrap )
  {
    if ( ana != EAnaType2020::kNumuAna )
    {
      if ( beam != BeamType2020::kRHC )
      {
        if ( ptExtrap ) systs.push_back( &kNueAcceptSystSignalKinPtExtrap2020FHC );
        else            systs.push_back( &kNueAcceptSystSignalKin2020FHC );
      }
      if ( beam != BeamType2020::kFHC )
      {
        if ( ptExtrap ) systs.push_back( &kNueAcceptSystSignalKinPtExtrap2020RHC );
        else            systs.push_back( &kNueAcceptSystSignalKin2020RHC );
      }
    }
  }

  void Add3FlavorAna2020MichelTagSysts( std::vector<const ISyst*> &systs,
                                        const EAnaType2020 ana,
                                        const BeamType2020 beam )
  {
    if ( ana != EAnaType2020::kNumuAna && beam != BeamType2020::kRHC )
    {
      systs.push_back( &kMichelTaggingSyst2020 );
    }
  }

  void Add3FlavorAna2020NotCorrelSysts( std::vector<const ISyst*> &systs,
                                        const EAnaType2020 ana,
                                        const BeamType2020 beam,
                                        const bool ptExtrap )
  {
    Add3FlavorAna2020NueAcceptSysts( systs, ana, beam, ptExtrap );
    Add3FlavorAna2020MichelTagSysts( systs, ana, beam );
    if ( ana != EAnaType2020::kNueAna )
    {
      systs.push_back( &kUnCorrFDMuEScaleSyst2020 );
    }
  }

  void Add3FlavorAna2020NormSysts( std::vector<const ISyst*> &systs,
                                   const BeamType2020 beam )
  {
    systs.push_back(&kAna2020NormHornCorr);
    if( beam != BeamType2020::kRHC ) systs.push_back(&kAna2020NormFHC);
    if( beam != BeamType2020::kFHC ) systs.push_back(&kAna2020NormRHC);
  }

  // not for 2020?
  void AddNonLoadable2020Systs(std::vector<const ISyst*> &systs,
                                 const EAnaType2020 ana)
  {
    if(ana != EAnaType2020::kNumuAna)
    {
      systs.push_back(&kRockScaleSyst);
    }
    systs.push_back(&kCosmicBkgScaleSyst); 
  }

  void Add3FlavorAna2020OtherSysts( std::vector<const ISyst*> &systs,
                                    const EAnaType2020 ana,
                                    const BeamType2020 beam,
                                    const bool ptExtrap )
  {
    Add3FlavorAna2020MuEnergySysts( systs, ana );
    Add3FlavorAna2020LeptonAngleSysts( systs, ptExtrap );
    Add3FlavorAna2020NeutronSysts( systs );
    Add3FlavorAna2020NormSysts( systs, beam );
    Add3FlavorAna2020NueAcceptSysts( systs, ana, beam, ptExtrap );
    Add3FlavorAna2020MichelTagSysts( systs, ana, beam );
  }

//////////////////////////////////////////////////////////////////////
  std::vector<const ISyst*> get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap )
  {
    std::vector<const ISyst*> ret;
    
    Add3FlavorAna2020XSecSysts (ret, ana, smallgenie);
    Add3FlavorAna2020BeamSysts (ret, ana);
    Add3FlavorAna2020FileSysts (ret);
    Add3FlavorAna2020OtherSysts(ret, ana, beam, ptExtrap);
    if(isFit) AddNonLoadable2020Systs(ret, ana);
    
    return ret;
  }


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

    Add3FlavorAna2020XSecSysts (ret, ana, smallgenie);
    Add3FlavorAna2020BeamSysts (ret, ana);
    Add3FlavorAna2020FileSysts (ret);
    Add3FlavorAna2020OtherSysts(ret, ana, beam, ptExtrap);
    // deprecated AddWrongSignSysts(ret, ana);
    if(isFit) AddNonLoadable2020Systs(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 *> get3FlavorAna2020SmallXsecSysts(const EAnaType2020 ana)
  {
    auto allgenie = getAllXsecSysts_2020();
    auto largegenie = get3FlavorAna2020LargeXsecSysts(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 * get3FlavorAna2020SummedSmallXsecSysts(const EAnaType2020 ana)
  {
    auto systs =  get3FlavorAna2020SmallXsecSysts(ana);
    auto anastr = AnaTypeToString(ana);
    return new SummedSyst(systs, "SumSmallXSec"+anastr,
                          "Summed small XSec systs "+anastr);
   }

  //////////////////////////////////////////////////////////////////////
  std::vector<const ISyst *> getAna2020NueOnly( const BeamType2020 beam, const bool ptExtrap )
  {
    std::vector<const ISyst*> ret;
    ret.push_back(&kRadCorrNue);
    ret.push_back(&kRadCorrNuebar);
    ret.push_back(&k2ndClassCurrs);
    Add3FlavorAna2020NueAcceptSysts( ret, EAnaType2020::kNueAna, beam, ptExtrap );
    Add3FlavorAna2020MichelTagSysts( ret, EAnaType2020::kNueAna, beam );
    AddNonLoadable2020Systs( ret, EAnaType2020::kNueAna );

    return ret;
  }

  std::vector<const ISyst *> getAna2020NumuOnly( const BeamType2020 beam )
  {
    std::vector<const ISyst*> ret;
    ret.push_back( &kUnCorrFDMuEScaleSyst2020 );
    AddNonLoadable2020Systs( ret, EAnaType2020::kNumuAna );
    return ret;
  }
}