MuonTriggerScaleFactors.cxx

Go to the documentation of this file.

/*
 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
 */
 
#include "TRandom3.h"
#include "TROOT.h"
#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TKey.h"
 
#include "xAODMuon/MuonContainer.h"
#include "MuonEfficiencyCorrections/MuonTriggerScaleFactors.h"
 
#include "AsgMessaging/StatusCode.h"
#include "PATInterfaces/SystematicRegistry.h"
#include "PATInterfaces/SystematicVariation.h"
#include "PathResolver/PathResolver.h"
#include "AsgDataHandles/ReadHandle.h"
 
#include <iostream>
#include <functional>
#include <sstream>
#include <string>
#include <cmath>
 
namespace ADD_CMODE(CP) {
    constexpr double muon_barrel_endcap_boundary = 1.05;
 
    struct RunYear {
        unsigned int lastRun{};
        unsigned int year{};
    };
 
    constexpr std::array runYears{
        RunYear{284484, 2015},
        RunYear{311481, 2016},
        RunYear{340453, 2017},
        RunYear{364292, 2018},
        RunYear{440613, 2022},
        RunYear{456749, 2023},
        RunYear{486706, 2024},
        RunYear{509849, 2025}
    };
 
    const std::map<unsigned short int, std::vector< std::tuple<unsigned int,unsigned int,std::string_view> > > MuonTriggerScaleFactors::m_year_periods = {
        {2015, {
            {266904,272531,"AC"},
            {276073,276954,"D"},
            {278727,279928,"E"},
            {279932,280422,"F"},
            {280423,281075,"G"},
            {281130,281411,"H"},
            {281662,282482,"I"},
            {282625,284484,"J"}}
        },
        {2016, {
            {296939,300287,"A"},
            {300345,300908,"B"},
            {301912,302393,"C"},
            {302737,302872,"D1D3"},
            {302919,303560,"D4D8"},
            {303638,303892,"E"},
            {303943,304494,"F"},
            {305291,306714,"G"},
            {307124,308084,"I"},
            {309311,309759,"K"},
            {310015,311481,"L"}}
        },
        {2017, {
            {324320,325558,"A"},
            {325713,328393,"B"},
            {329385,330470,"C"},
            {330857,332304,"D"},
            {332720,334779,"E"},
            {334842,335290,"F"},
            {336497,336782,"H"},
            {336832,337833,"I"},
            {338183,340453,"K"}}
        },
        {2018, {
            {348197,348836,"A"},
            {348885,349533,"B"},
            {349534,350220,"C"},
            {350310,352107,"D"},
            {352123,352137,"E"},
            {352274,352514,"F"},
            {354107,354494,"G"},
            {354826,355224,"H"},
            {355261,355273,"I"},
            {355331,355468,"J"},
            {355529,356259,"K"},
            {357050,359171,"L"},
            {359191,360414,"M"},
            {361635,361696,"N"},
            {361738,363400,"O"},
            {363664,364292,"Q"}}
        },
        {2022, {
            {430536,432180,"F"},
            {435816,439927,"H"},
            {440199,440613,"J"}}
        },
        {2023, {
            {451094,455924,"F"},
            {455975,456749,"G"}}
        },
        {2024, {
            {473235,473400,"E"},
            {473617,474271,"F"},
            {474441,474602,"G"},
            {474657,475522,"H"},
            {476060,477048,"I"},
            {479103,480032,"K"},
            {480188,482221,"M"},
            {482374,484799,"N"},
            {484909,486706,"O"}}
        },
        {2025, {
            {497924,498434,"D"},
            {498515,499592,"E"},
            {499618,500348,"F"},
            {502557,505458,"K"},
            {506231,507775,"M"},
            {508160,509452,"O"},
            {509803,509849,"Q"}}
        }
    };
 
    MuonTriggerScaleFactors::MuonTriggerScaleFactors(const std::string& name) :
        asg::AsgTool(name) {}
 
    MuonTriggerScaleFactors::~MuonTriggerScaleFactors() { }
 
  StatusCode MuonTriggerScaleFactors::LoadTriggerMap(unsigned int year) {
        std::string fileName;
        if (m_customInputFilePerYear.value().contains(year)) {
            fileName = m_customInputFilePerYear.value().at(year);
            if (fileName.empty()) {
                ANA_MSG_INFO("Skipping loading trigger SF for year " << year << " since the provided file name is empty");
                return StatusCode::SUCCESS;
            } else {
                ATH_MSG_INFO("Loading trigger SF from user specified file " << fileName << " for year " << year);
            }
        }
 
        if (fileName.empty()) {
          if (year == 2015) fileName = "muontrigger_sf_2015_mc20a_v3.root";
          else if (year == 2016) fileName = "muontrigger_sf_2016_mc20a_v3.root";
          else if (year == 2017) fileName = "muontrigger_sf_2017_mc20d_v3.root";
          else if (year == 2018) fileName = "muontrigger_sf_2018_mc20e_v3.root";
          else if (year == 2022) fileName = "muontrigger_sf_2022_mc23a_v5.root";
          else if (year == 2023) fileName = "muontrigger_sf_2023_mc23d_v5.root";
          else if (year == 2024) fileName = "muontrigger_sf_2024_mc23e_v5.root";
          else if (year == 2025) fileName = "muontrigger_sf_2025_mc23g_v2.root";
          else {
            ATH_MSG_WARNING("There is no SF file for year " << year << " yet");
            return StatusCode::SUCCESS;
          }
        }
 
        TDirectory* origDir = gDirectory;
 
        std::string filePath;
 
        if (m_customInputFolder.empty()) {
            filePath = PathResolverFindCalibFile(Form("MuonEfficiencyCorrections/%s/%s", m_calibrationVersion.value().c_str(), fileName.c_str()));
            if (filePath.empty()) {
                ATH_MSG_ERROR("Unable to resolve the input file " << fileName << " via PathResolver.");
            }
        }
        else {
            ATH_MSG_INFO("Note: setting up with user specified input file location " << m_customInputFolder.value() << " - this is not encouraged!");
            filePath = PathResolverFindCalibFile(Form("%s/%s", m_customInputFolder.value().c_str(), fileName.c_str()));
        }
 
        TFile* file = TFile::Open(filePath.c_str());
 
        if (file == nullptr || !file->IsOpen()) {
            ATH_MSG_FATAL("MuonTriggerScaleFactors::initialize Couldn't open file " << filePath);
            return StatusCode::FAILURE;
        }
        TDirectory* tempDir = getTemporaryDirectory();
        tempDir->cd();
 
        static const std::vector<std::string> type { "data", "mc" };
        static const std::vector<std::string> region { "barrel", "endcap" };
        static const std::vector<std::string> systematic { "nominal", "stat_up", "stat_down", "syst_up", "syst_down" };
        if (m_muonQuality.value().compare("LowPt") == 0) {
          m_muonQuality = "Medium";
        }
        TDirectory* qualityDirectory = file->GetDirectory(m_muonQuality.value().c_str());
        if (qualityDirectory == nullptr) {
            ATH_MSG_FATAL("MuonTriggerScaleFactors::initialize cannot find directory with selected quality");
            return StatusCode::FAILURE;
        }
        TKey* periodKey;
        TIter nextPeriod(qualityDirectory->GetListOfKeys());
        while ((periodKey = (TKey*) nextPeriod())) {
            if (not periodKey->IsFolder()) continue;
            TDirectory* periodDirectory = qualityDirectory->GetDirectory(periodKey->GetName());
            std::string periodName = std::string(periodKey->GetName());
            std::string periodName_InMap = periodName.substr(std::string("Period").size(), periodName.size());
            YearPeriod period = YearPeriod(year, periodName_InMap);
            TKey* triggerKey;
            TIter nextTrigger(periodDirectory->GetListOfKeys());
            while ((triggerKey = (TKey*) nextTrigger())) {
                if (not triggerKey->IsFolder()) continue;
                TDirectory* triggerDirectory = periodDirectory->GetDirectory(triggerKey->GetName());
                std::string triggerName = std::string(triggerKey->GetName());
                if(!std::set<std::string>{"HLT_mu26_ivarmedium", "HLT_mu50", "HLT_mu26_ivarmedium_OR_HLT_mu50"}.count(triggerName) && m_binning == "coarse"){
                  ATH_MSG_DEBUG("Coarse binning not supported for di-muon trigger legs at the moment");
                  continue;
                }
                for (const auto& iregion : region) {
                    bool isBarrel = iregion.find("barrel") != std::string::npos;
                    for (const auto& itype : type) {
                        bool isData = itype.find("data") != std::string::npos;
                        std::string histname = ("_MuonTrigEff_" + periodName + "_" + triggerName + "_" + m_muonQuality.value() + "_" + "_EtaPhi_" + m_binning + "_" + iregion + "_" + itype);
                        for (const auto& isys : systematic) {
                            if (itype.find("data") != std::string::npos && isys.find("syst") != std::string::npos) continue;
                            std::string path = "eff_etaphi_" + m_binning + "_" + iregion + "_" + itype + "_" + isys;
                            TH2* hist = dynamic_cast<TH2*>(triggerDirectory->Get(path.c_str()));
                            if (not hist) {
 
                                ATH_MSG_FATAL("MuonTriggerScaleFactors::initialize " << path << " not found under trigger " << triggerName << " and period " << periodName << " for year: " << year);
                                continue;
                            }
                            hist->SetDirectory(0);
 
                            EffiHistoIdent HistoId = EffiHistoIdent(period, encodeHistoName(periodName_InMap, triggerName, isData, isys, isBarrel));
                            if (m_efficiencyMap.find(HistoId) != m_efficiencyMap.end()) {
                                ATH_MSG_FATAL("MuonTriggerScaleFactors::initialize(): histogram " << path << " is duplicated for year" << year << " in period " << periodName);
                                return StatusCode::FAILURE;
                            }
                            m_efficiencyMap.insert(std::pair<EffiHistoIdent, TH1_Ptr>(HistoId, std::shared_ptr < TH1 > (hist)));
 
                        }
                        //If the trigger is chosen for toy evaluation, generate all the replicas from
                        // NOMINAL with STAT variations stored in the data hist, load them in corresponding vector
                        if (m_replicaSet.find(triggerName) != m_replicaSet.end() && itype.find("data") != std::string::npos) {
 
                            TH1_Ptr Nominal_H = getEfficiencyHistogram(year, periodName, triggerName, isData, "nominal", isBarrel);
                            TH1_Ptr StatUp_H = getEfficiencyHistogram(year, periodName, triggerName, isData, "stat_up", isBarrel);
 
                            TH1_Ptr tmp_h2 = TH1_Ptr(dynamic_cast<TH2F*>(Nominal_H->Clone(Form("tmp_h2_%s", Nominal_H->GetName()))));
                            const int xbins = tmp_h2->GetNbinsX(), ybins = tmp_h2->GetNbinsY();
                            for (int x_i = 0; x_i <= xbins; ++x_i) {
                                for (int y_i = 0; y_i <= ybins; ++y_i) {
                                    double statErr = std::abs(tmp_h2->GetBinContent(x_i, y_i) - StatUp_H->GetBinContent(x_i, y_i));
                                    tmp_h2->SetBinError(x_i, y_i, statErr);
                                }
                            }
                            m_efficiencyMapReplicaArray[EffiHistoIdent(period, encodeHistoName(periodName, triggerName, isData, "repl", isBarrel))] = generateReplicas(tmp_h2, m_nReplicas, m_ReplicaRandomSeed);
                        }
                    }
                }
            }
        }
        file->Close();
        delete file;
        origDir->cd();
        return StatusCode::SUCCESS;
    }
    // ==================================================================================
    // == MuonTriggerScaleFactors::initialize()
    // ==================================================================================
    StatusCode MuonTriggerScaleFactors::initialize() {
 
        ATH_MSG_INFO("MuonQuality = '" << m_muonQuality.value() << "'");
        ATH_MSG_INFO("Binning = '" << m_binning << "'");
        ATH_MSG_INFO("Campaign = '" << m_campaign << "'");
        ATH_MSG_INFO("CalibrationRelease = '" << m_calibrationVersion.value() << "'");
        ATH_MSG_INFO("CustomInputFolder = '" << m_customInputFolder.value() << "'");
        ATH_MSG_INFO("AllowZeroSF = " << m_allowZeroSF);
        ATH_MSG_INFO("experimental = " << m_experimental);
 
        ATH_CHECK(m_eventInfo.initialize());
 
        if (registerSystematics() != StatusCode::SUCCESS) {
            return StatusCode::FAILURE;
        }
 
        if (applySystematicVariation(CP::SystematicSet()) != StatusCode::SUCCESS) {
            ATH_MSG_ERROR("Could not configure for nominal settings");
            return StatusCode::FAILURE;
        }
        // Initialize indexes of replicas for trigges which are asked
        for (auto trigToy : m_replicaTriggerList)
            m_replicaSet.insert(trigToy);
 
 
        if (m_campaign.empty()) {
            constexpr auto years_to_run = std::to_array<int>({2015, 2016, 2017, 2018, 2022, 2023, 2024, 2025});
            for (const int &year: years_to_run) {
                ATH_CHECK(LoadTriggerMap(year));
            }
        } else if (m_campaign.value() == "mc20a") {
            ATH_CHECK(LoadTriggerMap(2015));
            ATH_CHECK(LoadTriggerMap(2016));
        } else if (m_campaign.value() == "mc20d") {
            ATH_CHECK(LoadTriggerMap(2017));
        } else if (m_campaign.value() == "mc20e") {
            ATH_CHECK(LoadTriggerMap(2018));
        } else if (m_campaign.value() == "mc23a") {
            ATH_CHECK(LoadTriggerMap(2022));
        } else if (m_campaign.value() == "mc23d") {
            ATH_CHECK(LoadTriggerMap(2023));
        } else if (m_campaign.value() == "mc23e") {
            ATH_CHECK(LoadTriggerMap(2024));
        } else if (m_campaign.value() == "mc23g") {
            ATH_CHECK(LoadTriggerMap(2025));
        } else {
            ATH_MSG_ERROR("Campaign " << m_campaign.value() << " is not supported. Please choose a valid campaign or leave empty to load all years.");
        }
        return StatusCode::SUCCESS;
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::getTriggerScaleFactor(const xAOD::Muon& muon, Double_t& triggersf, const std::string& trigger) const {
      if(!m_experimental){
    ATH_MSG_ERROR("MuonTriggerScaleFactors::getTriggerScaleFactor This is an experimental function. If you really know what you are doing set UseExperimental property.");
      return CP::CorrectionCode::Error;
      }
 
      if (trigger.empty()) {
    ATH_MSG_ERROR("MuonTriggerScaleFactors::getTriggerScaleFactor Trigger must have value.");
    return CP::CorrectionCode::Error;
      }
 
      TrigMuonEff::Configuration configuration;
 
      if (trigger == "HLT_mu8noL1")
    ATH_MSG_WARNING("What you are trying to do is not correct. For di-muon triggers you should get the efficiency with getTriggerEfficiency and compute the SF by yourself.");
      else if (trigger.find("HLT_2mu10") != std::string::npos || trigger.find("HLT_2mu14") != std::string::npos)
    ATH_MSG_WARNING("Di-muon trigger scale factors for single reco muons are not supported!");
      else
    return GetTriggerSF(triggersf, configuration, muon, trigger);
      return CP::CorrectionCode::Ok;
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::getTriggerScaleFactor(const xAOD::MuonContainer& mucont, Double_t& triggersf, const std::string& trigger) const{
        if (trigger.empty()) {
            ATH_MSG_ERROR("MuonTriggerScaleFactors::getTriggerScaleFactor Trigger must have value.");
            return CP::CorrectionCode::Error;
        }
 
        TrigMuonEff::Configuration configuration;
 
        if (trigger == "HLT_mu8noL1") {
            ATH_MSG_WARNING("What you are trying to do is not correct. For di-muon triggers you should get the efficiency with getTriggerEfficiency and compute the SF by yourself.");
        } else if (trigger.find("HLT_2mu10") != std::string::npos || trigger.find("HLT_2mu14") != std::string::npos) {
            CP::CorrectionCode cc = GetTriggerSF_dimu(triggersf, configuration, mucont, trigger);
            return cc;
        } else {
            CP::CorrectionCode cc = GetTriggerSF(triggersf, configuration, mucont, trigger);
            return cc;
        }
        return CP::CorrectionCode::Ok;
    }
 
    // ==================================================================================
    // == MuonTriggerScaleFactors::getReplica_index
    // ==================================================================================
    // Gets  replica index correponding to the toy.
    // Also checks if the sys_name contains "MCTOY" and if the trigger has replicas generated.
    // Returns -1 if conditions are note satisfied
    int MuonTriggerScaleFactors::getReplica_index(const std::string& sysBaseName, const std::string& trigStr) const{
        if (m_replicaSet.find(trigStr) == m_replicaSet.end()) return -1; //No toys for this trigger
        std::size_t pos = sysBaseName.find("MCTOY");
        if (pos == std::string::npos) return -1; //sys variation not affected by TOYS
        return atoi(sysBaseName.substr(pos + 5, pos + 8).c_str()); //toys for this trigger are around get the 3-digit number
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::getTriggerEfficiency(const xAOD::Muon& mu, Double_t& efficiency, const std::string& trigger, Bool_t dataType) const{
        if (trigger.empty()) {
            ATH_MSG_ERROR("MuonTriggerScaleFactors::getTriggerEfficiency Trigger must have value.");
            return CP::CorrectionCode::Error;
        }
        TrigMuonEff::Configuration configuration;
        configuration.isData = dataType;
        configuration.replicaIndex = -1;
        Int_t threshold;
        CP::CorrectionCode result = getThreshold(threshold, trigger);
        if (result != CP::CorrectionCode::Ok) return result;
        if (mu.pt() < threshold) {
            efficiency = 0;
            return CP::CorrectionCode::Ok;
        }
 
        // Pre-define uncertainty variations
        static const CP::SystematicVariation stat_up("MUON_EFF_TrigStatUncertainty", 1);
        static const CP::SystematicVariation stat_down("MUON_EFF_TrigStatUncertainty", -1);
        static const CP::SystematicVariation syst_up("MUON_EFF_TrigSystUncertainty", 1);
        static const CP::SystematicVariation syst_down("MUON_EFF_TrigSystUncertainty", -1);
 
        std::string systype = "";
        if (appliedSystematics().matchSystematic(syst_down) && !dataType) {
            systype = "syst_down";
        } else if (appliedSystematics().matchSystematic(syst_up) && !dataType) {
            systype = "syst_up";
        } else if (appliedSystematics().matchSystematic(stat_down)) {
            systype = "stat_down";
        } else if (appliedSystematics().matchSystematic(stat_up)) {
            systype = "stat_up";
        } else {
            systype = "nominal";
        }
 
        // Toys, if found, will overwrite the data hists stat with the generated toy
        //+++++++++++++
        // The best way is the use of filterByName with the 000MCTOY at the end. See:
        // if( !(appliedSystematics().filterByBaseName("MUON_EFF_Trig_MCTOY000")).empty()){//The following is a hack!!!
        //++++++++++THE FOLLOWING IS A PARTIAL HACK!!!
        if (!appliedSystematics().empty() && configuration.isData == true) {
            configuration.replicaIndex = getReplica_index(appliedSystematics().begin()->basename(), trigger);
            if (configuration.replicaIndex != -1) systype = "replicas";
        }
        CP::CorrectionCode cc = getMuonEfficiency(efficiency, configuration, mu, trigger, systype);
        return cc;
    }

    // Private functions //
 
    // ==================================================================================
    // == MuonTriggerScaleFactors::generateReplicas
    // ==================================================================================
    // Generate replicas of h for Toys with each bin of h varied with Gaussian distribution
    // with mean from bin content and sigma from bin error
    std::vector<TH1_Ptr> MuonTriggerScaleFactors::generateReplicas(TH1_Ptr h, int nrep, int seed) const {
        TRandom3 Rndm(seed);
        std::vector<TH1_Ptr> replica_v;
        const int xbins = h->GetNbinsX(), ybins = h->GetNbinsY();
 
        for (int t = 0; t < nrep; ++t) {
            TH2* replica = dynamic_cast<TH2*>(h->Clone(Form("rep%d_%s", t, h->GetName())));
 
            for (int x_i = 0; x_i <= xbins; ++x_i) {
                for (int y_i = 0; y_i <= ybins; ++y_i) {
                    replica->SetBinContent(x_i, y_i, Rndm.Gaus(h->GetBinContent(x_i, y_i), h->GetBinError(x_i, y_i)));
                }
            }
            replica_v.push_back(TH1_Ptr(replica));
        }
        return replica_v;
    }
 
   bool MuonTriggerScaleFactors::isTriggerSupported(const std::string& trigger) const{
       TH1_Ptr H1 = getEfficiencyHistogram(trigger, true, "nominal");
       return H1.get() != nullptr;
    }
 
  int MuonTriggerScaleFactors::getBinNumber(const xAOD::Muon& muon, const std::string& trigger) const{
    if(!m_experimental){
      ATH_MSG_ERROR("MuonTriggerScaleFactors::getTriggerScaleFactor This is an experimental function. If you really know what you are doing set UseExperimental property.");
      return CP::CorrectionCode::Error;
    }
 
    const double mu_eta = muon.eta();
    const double mu_phi = muon.phi();
    bool isBarrel = std::abs(mu_eta) < muon_barrel_endcap_boundary;
    TH1_Ptr cit = getEfficiencyHistogram(trigger, true, "nominal", isBarrel);
    if(!cit.get()){
      if(!m_allowZeroSF)
    ATH_MSG_ERROR("Could not find efficiency map for muon with eta: " << mu_eta << " and phi: " << mu_phi << ". Something is inconsistent. Please check your settings for year, mc and trigger." );
      return -1;
    }
    auto eff_h2 = cit;
    double mu_phi_corr = mu_phi;
    if (mu_phi_corr < eff_h2->GetYaxis()->GetXmin()) mu_phi_corr += 2.0 * M_PI;
    if (mu_phi_corr > eff_h2->GetYaxis()->GetXmax()) mu_phi_corr -= 2.0 * M_PI;
    return eff_h2->FindFixBin(mu_eta, mu_phi_corr);
  }
 
    unsigned int MuonTriggerScaleFactors::encodeHistoName(const std::string_view period, const std::string& Trigger, bool isData, const std::string& Systematic, bool isBarrel) const {
        //keep the string as short as possible
        const std::string histName = std::string{period} + "_" + Trigger + "_" + (isBarrel ? "b" : "e") + "_" + (isData ? "data" : "mc") + Systematic;
        return std::hash<std::string>()(histName);
    }
 
    unsigned int MuonTriggerScaleFactors::encodeHistoName(const std::string& Trigger, const TrigMuonEff::Configuration& configuration, const std::string& Systematic, bool isBarrel) const {
        //keep the string as short as possible
        return encodeHistoName(getDataPeriod(), Trigger, configuration.isData, Systematic, isBarrel);
 
    }
    TH1_Ptr MuonTriggerScaleFactors::getEfficiencyHistogram(unsigned int year, const std::string_view period, const std::string& trigger, bool isData, const std::string& Systematic, bool isBarrel) const {
        EffiHistoIdent Ident = EffiHistoIdent(YearPeriod(year, period), encodeHistoName(period, trigger, isData, Systematic, isBarrel));
        EfficiencyMap::const_iterator Itr = m_efficiencyMap.find(Ident);
 
        if (Itr == m_efficiencyMap.end()) {
            return TH1_Ptr();
        }
        return Itr->second;
    }
    TH1_Ptr MuonTriggerScaleFactors::getEfficiencyHistogram(const std::string& trigger, bool isData, const std::string& Systematic, bool isBarrel) const {
        unsigned int run = getRunNumber();
        return getEfficiencyHistogram(getYear(run), getDataPeriod(run), trigger, isData, Systematic, isBarrel);
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::getMuonEfficiency(Double_t& eff, const TrigMuonEff::Configuration& configuration, const xAOD::Muon& muon, const std::string& trigger, const std::string& systematic) const{
        const double mu_eta = muon.eta();
        const double mu_phi = muon.phi();
        bool isBarrel = std::abs(mu_eta) < muon_barrel_endcap_boundary;
 
        TH1_Ptr eff_h2 = nullptr;
        if (configuration.replicaIndex >= 0) { //Only look into the replicas if asking for them
 
            unsigned int run = getRunNumber();
            EffiHistoIdent Ident = EffiHistoIdent(YearPeriod(getYear(run), getDataPeriod(run)), encodeHistoName(getDataPeriod(run), trigger, configuration.isData, "repl", isBarrel));
            std::map<EffiHistoIdent, std::vector<TH1_Ptr> >::const_iterator cit = m_efficiencyMapReplicaArray.find(Ident);
            if (cit == m_efficiencyMapReplicaArray.end()) {
                if (m_allowZeroSF) {
                    ATH_MSG_WARNING("Could not find what you are looking for in the efficiency map. The trigger you are looking for, year and mc are not consistent, or the trigger is unavailable in this data period. Returning efficiency = 0.");
                    eff = 0.;
                    return CP::CorrectionCode::Ok;
                }
 
                else {
                    ATH_MSG_ERROR("Could not find what you are looking for in the efficiency map. The trigger you are looking for, year and mc are not consistent, or the trigger is unavailable in this data period. Please check how you set up the tool.");
                    return CP::CorrectionCode::OutOfValidityRange;
                }
            }
 
            if (configuration.replicaIndex >= (int) cit->second.size()) {
                ATH_MSG_ERROR("MuonTriggerScaleFactors::getMuonEfficiency ; index for replicated histograms is out of range.");
                return CP::CorrectionCode::OutOfValidityRange;
            }
 
            eff_h2 = cit->second[configuration.replicaIndex];
        } else { //Standard case, look into the usual eff map
            TH1_Ptr cit = getEfficiencyHistogram(trigger, configuration.isData, systematic, isBarrel);
            if (cit.get() == nullptr) {
                if (m_allowZeroSF) {
                    ATH_MSG_WARNING("Could not find what you are looking for in the efficiency map. The trigger you are looking for, year and mc are not consistent, or the trigger is unavailable in this data period. Returning efficiency = 0.");
                    eff = 0.;
                    return CP::CorrectionCode::Ok;
                } else {
                    ATH_MSG_ERROR("Could not find what you are looking for in the efficiency map. The trigger you are looking for, year and mc are not consistent, or the trigger is unavailable in this data period. Please check how you set up the tool.");
                    return CP::CorrectionCode::OutOfValidityRange;
                }
            }
            eff_h2 = cit;
        }
 
        double mu_phi_corr = mu_phi;
        if (mu_phi_corr < eff_h2->GetYaxis()->GetXmin()) mu_phi_corr += 2.0 * M_PI;
        if (mu_phi_corr > eff_h2->GetYaxis()->GetXmax()) mu_phi_corr -= 2.0 * M_PI;
 
        const int bin = eff_h2->FindFixBin(mu_eta, mu_phi_corr);
        const double efficiency = eff_h2->GetBinContent(bin);
 
        eff = efficiency;
 
        ATH_MSG_DEBUG("getMuonEfficiency [eta,phi,phi_corr]=[" << mu_eta << "," << mu_phi << "," << mu_phi_corr << "], ibin=" << bin << " eff=" << eff);
 
        return CP::CorrectionCode::Ok;
 
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::GetTriggerSF_dimu(Double_t& TriggerSF, TrigMuonEff::Configuration& configuration, const xAOD::MuonContainer& mucont, const std::string& trigger) const{
 
        if (mucont.size() != 2) {
            ATH_MSG_FATAL("MuonTriggerScaleFactors::GetTriggerSF;Currently dimuon trigger chains only implemented for events with exactly 2 muons.");
        }
        ATH_MSG_DEBUG("The trigger that you choose : " << trigger);
 
        Double_t eff_data = 0;
        Double_t eff_mc = 0;
 
        std::string data_err = "";
        std::string mc_err = "";
 
        if (appliedSystematics().matchSystematic(CP::SystematicVariation("MUON_EFF_TrigSystUncertainty", -1))) {
            data_err = "nominal";
            mc_err = "syst_up";
        } else if (appliedSystematics().matchSystematic(CP::SystematicVariation("MUON_EFF_TrigSystUncertainty", 1))) {
            data_err = "nominal";
            mc_err = "syst_down";
        } else if (appliedSystematics().matchSystematic(CP::SystematicVariation("MUON_EFF_TrigStatUncertainty", -1))) {
            data_err = "stat_down";
            mc_err = "nominal";
        } else if (appliedSystematics().matchSystematic(CP::SystematicVariation("MUON_EFF_TrigStatUncertainty", 1))) {
            data_err = "stat_up";
            mc_err = "nominal";
        } else {
            data_err = "nominal";
            mc_err = "nominal";
        }
 
        //Toys, if found, will overwrite the data hists with the sys generated with one toy
        //+++++++++++++
        //The best way is the use of filterByName with the 000MCTOY at the end. See:
        // if( !(appliedSystematics().filterByBaseName("MUON_EFF_Trig_MCTOY000")).empty()){//The following is a hack!!!
        //++++++++++THE FOLLOWING IS A PARTIAL HACK!!!
        if (!appliedSystematics().empty()) {
            configuration.replicaIndex = getReplica_index(appliedSystematics().begin()->basename(), trigger);
            if (configuration.replicaIndex != -1) data_err = "replicas";
        }
 
        configuration.isData = true;
        CP::CorrectionCode result = getDimuonEfficiency(eff_data, configuration, mucont, trigger, data_err);
        if (result != CP::CorrectionCode::Ok) return result;
 
        configuration.isData = false;
        configuration.replicaIndex = -1;
        result = getDimuonEfficiency(eff_mc, configuration, mucont, trigger, mc_err);
        if (result != CP::CorrectionCode::Ok) return result;
 
        double event_SF = 1.;
 
        if (std::abs(1. - eff_mc) > 0.0001) {
            event_SF = eff_data / eff_mc;
        }
 
        TriggerSF = event_SF;
        return CP::CorrectionCode::Ok;
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::GetTriggerSF(Double_t& TriggerSF, TrigMuonEff::Configuration& configuration, const xAOD::MuonContainer& mucont, const std::string& trigger) const{
        Int_t threshold;
        CP::CorrectionCode result = getThreshold(threshold, trigger);
        if (result != CP::CorrectionCode::Ok) return result;
 
        double rate_not_fired_data = 1.;
        double rate_not_fired_mc = 1.;
 
        for (const auto mu : mucont) {
 
            double eff_data = 0., eff_mc = 0.;
 
            if (mu->pt() < threshold) {
 
                eff_data = 0.;
                eff_mc = 0.;
 
            } else {
 
                std::string muon_trigger_name = trigger;
                std::string data_err = "";
                std::string mc_err = "";
 
                // Pre-define uncertainty variations
                static const CP::SystematicVariation stat_up("MUON_EFF_TrigStatUncertainty", 1);
                static const CP::SystematicVariation stat_down("MUON_EFF_TrigStatUncertainty", -1);
                static const CP::SystematicVariation syst_up("MUON_EFF_TrigSystUncertainty", 1);
                static const CP::SystematicVariation syst_down("MUON_EFF_TrigSystUncertainty", -1);
 
                if (appliedSystematics().matchSystematic(syst_down)) {
                    data_err = "nominal";
                    mc_err = "syst_up";
                } else if (appliedSystematics().matchSystematic(syst_up)) {
                    data_err = "nominal";
                    mc_err = "syst_down";
                } else if (appliedSystematics().matchSystematic(stat_down)) {
                    data_err = "stat_down";
                    mc_err = "nominal";
                } else if (appliedSystematics().matchSystematic(stat_up)) {
                    data_err = "stat_up";
                    mc_err = "nominal";
                } else {
                    data_err = "nominal";
                    mc_err = "nominal";
                }
 
                //Toys, if found, will overwrite the data hists, on which toys for stat uncertainty have been generated
                //+++++++++++++
                //The best way is the use of filterByName with the 000MCTOY at the end. See:
                // if( !(appliedSystematics().filterByBaseName("MUON_EFF_Trig_MCTOY000")).empty()){//The following is a hack!!!
                //++++++++++The following is a hack!!!
                if (!appliedSystematics().empty()) {
                    configuration.replicaIndex = getReplica_index(appliedSystematics().begin()->basename(), trigger);
                    if (configuration.replicaIndex != -1) data_err = "replicas";
                }
 
                configuration.isData = true;
                CP::CorrectionCode result_data = getMuonEfficiency(eff_data, configuration, *mu, muon_trigger_name, data_err);
                if (result_data != CP::CorrectionCode::Ok) return result_data;
                configuration.isData = false;
                configuration.replicaIndex = -1;
                CP::CorrectionCode result_mc = getMuonEfficiency(eff_mc, configuration, *mu, muon_trigger_name, mc_err);
                if (result_mc != CP::CorrectionCode::Ok) return result_mc;
            }
            rate_not_fired_data *= (1. - eff_data);
            rate_not_fired_mc *= (1. - eff_mc);
        }
 
        double event_SF = 1.;
        if (1 - rate_not_fired_data == 0) event_SF = 0;
        if ((mucont.size()) and (std::abs(1. - rate_not_fired_mc) > 0.0001)) {
 
            event_SF = (1. - rate_not_fired_data) / (1. - rate_not_fired_mc);
        }
        TriggerSF = event_SF;
 
        return CP::CorrectionCode::Ok;
    }
 
  CP::CorrectionCode MuonTriggerScaleFactors::GetTriggerSF(Double_t& TriggerSF, TrigMuonEff::Configuration& configuration, const xAOD::Muon& mu, const std::string& trigger) const{
        Int_t threshold;
        CP::CorrectionCode result = getThreshold(threshold, trigger);
        if (result != CP::CorrectionCode::Ok)
      return result;
 
    double eff_data = 0., eff_mc = 0.;
 
    if (mu.pt() < threshold) {
      eff_data = 0.;
      eff_mc = 0.;
      TriggerSF = 1.;
      return CP::CorrectionCode::Ok;
    }
 
    std::string muon_trigger_name = trigger;
    std::string data_err = "";
    std::string mc_err = "";
 
    // Pre-define uncertainty variations
    static const CP::SystematicVariation stat_up("MUON_EFF_TrigStatUncertainty", 1);
    static const CP::SystematicVariation stat_down("MUON_EFF_TrigStatUncertainty", -1);
    static const CP::SystematicVariation syst_up("MUON_EFF_TrigSystUncertainty", 1);
    static const CP::SystematicVariation syst_down("MUON_EFF_TrigSystUncertainty", -1);
 
    if (appliedSystematics().matchSystematic(syst_down)) {
      data_err = "nominal";
      mc_err = "syst_up";
    } else if (appliedSystematics().matchSystematic(syst_up)) {
      data_err = "nominal";
      mc_err = "syst_down";
    } else if (appliedSystematics().matchSystematic(stat_down)) {
      data_err = "stat_down";
      mc_err = "nominal";
    } else if (appliedSystematics().matchSystematic(stat_up)) {
      data_err = "stat_up";
      mc_err = "nominal";
    } else {
      data_err = "nominal";
      mc_err = "nominal";
    }
 
    if (!appliedSystematics().empty()) {
      configuration.replicaIndex = getReplica_index(appliedSystematics().begin()->basename(), trigger);
      if (configuration.replicaIndex != -1) data_err = "replicas";
    }
 
    configuration.isData = true;
    CP::CorrectionCode result_data = getMuonEfficiency(eff_data, configuration, mu, muon_trigger_name, data_err);
    if (result_data != CP::CorrectionCode::Ok)
      return result_data;
    configuration.isData = false;
    configuration.replicaIndex = -1;
    CP::CorrectionCode result_mc = getMuonEfficiency(eff_mc, configuration, mu, muon_trigger_name, mc_err);
    if (result_mc != CP::CorrectionCode::Ok)
      return result_mc;
    if (eff_data == 0)
      TriggerSF =  0;
        if (std::abs(eff_mc) > 0.0001)
      TriggerSF = eff_data / eff_mc;
        return CP::CorrectionCode::Ok;
  }
 
    CP::CorrectionCode MuonTriggerScaleFactors::getDimuonEfficiency(Double_t& eff, const TrigMuonEff::Configuration& configuration, const xAOD::MuonContainer& mucont, const std::string& chain, const std::string& systematic) const{
 
        std::string trigger = getTriggerCorrespondingToDimuonTrigger(chain);
        Int_t threshold;
        CP::CorrectionCode result = getThreshold(threshold, trigger);
        if (result != CP::CorrectionCode::Ok) return result;
 
        xAOD::MuonContainer::const_iterator mu1 = mucont.begin();
        xAOD::MuonContainer::const_iterator mu2 = mucont.begin() + 1;
 
        Double_t eff1 = 0;
        if ((**mu1).pt() > threshold) {
            CP::CorrectionCode result1 = getMuonEfficiency(eff1, configuration, (**mu1), trigger, systematic);
            if (result1 != CP::CorrectionCode::Ok) return result1;
        }
        Double_t eff2 = 0;
        if ((**mu2).pt() > threshold) {
            CP::CorrectionCode result2 = getMuonEfficiency(eff2, configuration, (**mu2), trigger, systematic);
            if (result2 != CP::CorrectionCode::Ok) return result2;
        }
 
        eff = eff1 * eff2;
        return CP::CorrectionCode::Ok;
    }
 
    std::string MuonTriggerScaleFactors::getTriggerCorrespondingToDimuonTrigger(const std::string& trigger) const {
        if (trigger.find("2mu10") != std::string::npos) return "HLT_mu10";
        if (trigger.find("2mu14") != std::string::npos) return "HLT_mu14";
        throw std::runtime_error("Unknown dimuon trigger");
    }
 
    CP::CorrectionCode MuonTriggerScaleFactors::getThreshold(Int_t& threshold, const std::string& trigger) const{
        std::size_t index = trigger.find("HLT_mu");
        if (index != std::string::npos) {
            std::string rawNumber = trigger.substr(index + 6);
            if (!rawNumber.empty() && isdigit(rawNumber[0])) {
                std::stringstream(rawNumber) >> threshold;
                if (threshold < 10) threshold = 10000;
                else threshold = (threshold + 1) * 1000;
                return CP::CorrectionCode::Ok;
            }
        }
        ATH_MSG_ERROR("MuonTriggerScaleFactors::getThreshold Could not extract threshold for trigger " << trigger);
        return CP::CorrectionCode::Error;
    }
 
    unsigned int MuonTriggerScaleFactors::getYear(unsigned int run) const {
        if(m_forceYear != -1) return m_forceYear;
 
        auto runCheck = RunYear{run,2050};
        auto lower = std::ranges::lower_bound(runYears, runCheck, [](const auto &s1, const auto &s2) { return s1.lastRun < s2.lastRun; });
        return (lower != runYears.end()) ? lower->year : runCheck.year;
    }
 
    const std::string_view MuonTriggerScaleFactors::getDataPeriod() const {
        return getDataPeriod(getRunNumber());
    }
 
    const std::string_view MuonTriggerScaleFactors::getDataPeriod(unsigned int run) const {
        return getDataPeriod(run, getYear(run));
    }
 
    const std::string_view MuonTriggerScaleFactors::getDataPeriod(unsigned int runNumber, unsigned year) const {
        if(!m_forcePeriod.empty()) return std::string_view{m_forcePeriod};
 
        auto yit = m_year_periods.find(year);
        if (yit != m_year_periods.end()) {
            auto runit = std::find_if(yit->second.begin(), yit->second.end(), [&runNumber](const auto &item){return runNumber >= std::get<0>(item) && runNumber <= std::get<1>(item);});
            if (runit != yit->second.end()) return std::get<2>(*runit);
            else {
                ATH_MSG_FATAL("RunNumber " << runNumber << " does not belong to any data period in " << year);
                throw std::invalid_argument{""};
            }
        }
        else {
            ATH_MSG_FATAL("No data were collected during " << year);
            throw std::invalid_argument{""};
        }
    }
 
    unsigned int MuonTriggerScaleFactors::getRunNumber() const {
        static const SG::ConstAccessor<unsigned int> acc_rnd("RandomRunNumber");
        SG::ReadHandle<xAOD::EventInfo> info(m_eventInfo);
        if (info.operator->()==nullptr) {
            ATH_MSG_FATAL("Could not retrieve the xAOD::EventInfo with name: " << m_eventInfo.key() << ". Exiting the code.");
            throw std::invalid_argument{""};
        }
        if (!info->eventType(xAOD::EventInfo::IS_SIMULATION)) {
            ATH_MSG_DEBUG("The current event is a data event. Return runNumber instead.");
            return info->runNumber();
        }
        if (!acc_rnd.isAvailable(*info)) {
            if(m_forceYear == -1 && m_forcePeriod == "")
            ATH_MSG_FATAL("Failed to find the RandomRunNumber decoration. Please call the apply() method from the PileupReweightingTool beforehand in order to get period dependent SFs");
            throw std::invalid_argument{""};
        } else if (acc_rnd(*info) == 0) {
            ATH_MSG_FATAL("Pile up tool has given runNumber 0. Exiting the code.");
            throw std::invalid_argument{""};
        }
 
        // standard behaviour for MC, get the random RunNumber
        return acc_rnd(*info);
    }
 
    TDirectory* MuonTriggerScaleFactors::getTemporaryDirectory(void) const {
        gROOT->cd();
        TDirectory* tempDir = 0;
        int counter = 0;
        while (not tempDir) {
            std::stringstream dirname;
            dirname << "MuonTriggerScaleFactorsTempDir_%i" << counter;
            if (gROOT->GetDirectory((dirname.str()).c_str())) {
                ++counter;
                continue;
            }
            tempDir = gROOT->mkdir((dirname.str()).c_str());
            if (not tempDir) {
                ATH_MSG_ERROR("getTemporaryDirectory::Temporary directory could not be created");
            }
        }
        return tempDir;
    }
 
    //=======================================================================
    //   Systematics Interface
    //=======================================================================
    bool MuonTriggerScaleFactors::isAffectedBySystematic(const CP::SystematicVariation& systematic) const {
        if (!systematic.empty()) {
            CP::SystematicSet sys = affectingSystematics();
            return sys.find(systematic) != sys.end();
        }
        return true;
    }

    CP::SystematicSet MuonTriggerScaleFactors::affectingSystematics() const {
        CP::SystematicSet mySysSet;
 
        mySysSet.insert(CP::SystematicVariation("MUON_EFF_TrigSystUncertainty", 1));
        mySysSet.insert(CP::SystematicVariation("MUON_EFF_TrigSystUncertainty", -1));
 
        //Consider full statUncertainty if  TOY replicas are not used
        if (m_replicaTriggerList.size() == 0) {
            mySysSet.insert(CP::SystematicVariation("MUON_EFF_TrigStatUncertainty", 1));
            mySysSet.insert(CP::SystematicVariation("MUON_EFF_TrigStatUncertainty", -1));
        } else {
            for (int i = 0; i < m_nReplicas; ++i) { //TOFIX Hack with just up variations! needs ASG reserved words for a clean handling//+++++++
                mySysSet.insert(CP::SystematicVariation(Form("MUON_EFF_Trig_MCTOY%03d", i), 1));
            }
        }
 
        return mySysSet;
    }
 
    // Register the systematics with the registry and add them to the recommended list
    StatusCode MuonTriggerScaleFactors::registerSystematics() {
        CP::SystematicRegistry& registry = CP::SystematicRegistry::getInstance();
        if (registry.registerSystematics(*this) != StatusCode::SUCCESS) {
            ATH_MSG_ERROR("Failed to add systematic to list of recommended systematics.");
            return StatusCode::FAILURE;
        }
        return StatusCode::SUCCESS;
    }

    CP::SystematicSet MuonTriggerScaleFactors::recommendedSystematics() const {
        return affectingSystematics();
    }
 
    StatusCode MuonTriggerScaleFactors::applySystematicVariation(const CP::SystematicSet& systConfig){
        // First, check if we already know this systematic configuration
        auto itr = m_systFilter.find(systConfig);
 
        // If it's a new input set, we need to filter it
        if (itr == m_systFilter.end()) {
 
            // New systematic. We need to parse it.
            static const CP::SystematicSet affectingSys = affectingSystematics();
            CP::SystematicSet filteredSys;
            if (!CP::SystematicSet::filterForAffectingSystematics(systConfig, affectingSys, filteredSys)) {
                ATH_MSG_ERROR("Unsupported combination of systematics passed to the tool!");
                return StatusCode::FAILURE;
            }
 
            // Insert filtered set into the map
            itr = m_systFilter.insert(std::make_pair(systConfig, filteredSys)).first;
        }
 
        CP::SystematicSet& mySysConf = itr->second;
 
        // Check to see if the set of variations tries to add in the uncertainty up and down. Since the errors
        // are symetric this would result in 0 and so should not be done.
        static const CP::SystematicVariation stat_up("MUON_EFF_TrigStatUncertainty", 1);
        static const CP::SystematicVariation stat_down("MUON_EFF_TrigStatUncertainty", -1);
        static const CP::SystematicVariation syst_up("MUON_EFF_TrigSystUncertainty", 1);
        static const CP::SystematicVariation syst_down("MUON_EFF_TrigSystUncertainty", -1);
 
        if ((mySysConf.matchSystematic(stat_up) && mySysConf.matchSystematic(stat_down)) || (mySysConf.matchSystematic(syst_up) && mySysConf.matchSystematic(syst_down))) {
            return StatusCode::FAILURE;
        }
 
        m_appliedSystematics = &mySysConf;
        return StatusCode::SUCCESS;
    }
 
} /* namespace ADD_CMODE(CP) */