dumpNPs.cxx File Reference

#include "CaloGeoHelpers/CaloSampling.h"
#include "AsgTools/AnaToolHandle.h"
#include "EgammaAnalysisInterfaces/IAsgElectronEfficiencyCorrectionTool.h"
#include "ElectronEfficiencyCorrection/AsgElectronEfficiencyCorrectionTool.h"
#include "AsgMessaging/MessageCheck.h"
#include "CreateDummyEl.h"
#include <vector>
#include <map>
#include <tuple>
#include <bitset>
#include <sstream>
#include <type_traits>
#include <cstdlib>
#include <limits>
#include <algorithm>

Go to the source code of this file.

Classes
struct	Domain
struct	NP
struct	Config

Macros
#define	MSGSOURCE   "EgEfficiencyCorr_dumpNPs"

Typedefs
using	map_t = std::map<NP, std::vector<std::size_t>>

Functions
template<typename... Args>
bool	parse_csv_list (std::string val, Args &... args)
bool	scanPhaseSpace (Config &config, map_t &affected_bins)
bool	displayFindings (const Config &config, const map_t &affected_bins)
bool	displayFindings_analysis (const Config &cfg, const map_t &affected_bins)
bool	find_boundaries (const Config &cfg, const map_t::mapped_type &affected_bins, Domain &bounds, bool &abs_eta, bool &holes)
int	get_run_number (const int year)
int	main (int argc, char *argv[])
bool	parse_csv_token (const float x, std::vector< float > &bins)
bool	parse_csv_token (std::string s, std::vector< int8_t > &flags, const std::vector< float > &edges)

Variables
std::vector< float >	eta_edges
std::vector< float >	pt_edges
const std::size_t	subdiv_eta {5}
const std::size_t	subdiv_pt {3}

Macro Definition Documentation

MSGSOURCE

#define MSGSOURCE   "EgEfficiencyCorr_dumpNPs"

Definition at line 23 of file dumpNPs.cxx.

Typedef Documentation

map_t

using map_t = std::map<NP, std::vector<std::size_t>>

Definition at line 58 of file dumpNPs.cxx.

Function Documentation

displayFindings()

bool displayFindings	(	const Config &	config,
		const map_t &	affected_bins )

Definition at line 376 of file dumpNPs.cxx.

{
  using namespace asg::msgUserCode;
  for (const auto& s : cfg.systematics) {
    const NP np{std::get<NP>(s)};
    const std::string name(std::get<CP::SystematicSet>(s).begin()->basename());
    Domain bounds{};
    bool valid{false}, holes{false}, abs_eta{false};
    auto itr = affected_bins.find(np);
    if (itr != affected_bins.end()) {
      valid = find_boundaries(cfg, itr->second, bounds, abs_eta, holes);
    }
    if (!valid) {
      Warning(MSGSOURCE, "%s seems to not affect any bin", name.c_str());
      continue;
    }
    std::string txt(name + " --> ");
    if (holes) txt += "subdomain of ";
    txt += bounds.str(abs_eta) + '.';
    if (!holes) {
      Info(MSGSOURCE, "%s", txt.c_str());
    } else {
      Warning(MSGSOURCE, "%s", txt.c_str());
    }
  }
  return true;
}

displayFindings_analysis()

bool displayFindings_analysis	(	const Config &	cfg,
		const map_t &	affected_bins )

Definition at line 405 of file dumpNPs.cxx.

{
  using namespace asg::msgUserCode;
 
  std::vector<std::pair<std::string, int>> summary;
  NP prev{-888, false};
  for (const auto& kv : affected_bins) {
    const auto& bins{kv.second};
    if (bins.empty()) continue;
    const NP np{kv.first};
    const bool next{prev.uncorr==np.uncorr && np.index==prev.index+1};
    prev = np;
    if (!summary.empty() && next) {
      std::get<int>(summary.back()) = np.index;
      continue;
    }
    const auto& sys{
      std::get<CP::SystematicSet>(
        *std::find_if(
          cfg.systematics.cbegin(), cfg.systematics.cend(), 
          [=](auto& x){ return std::get<NP>(x) == np; }))};
    summary.emplace_back(sys.begin()->basename(), -1);
  }
  for (const auto& x: summary) {
    std::string s("   ++ " + std::get<std::string>(x));
    const int i{std::get<int>(x)};
    if (i >= 0) s += " to " + std::to_string(i);
    Info(MSGSOURCE, "%s", s.c_str());
  }
  return true;
}

find_boundaries()

bool find_boundaries	(	const Config &	cfg,
		const map_t::mapped_type &	affected_bins,
		Domain &	bounds,
		bool &	abs_eta,
		bool &	holes )

Definition at line 438 of file dumpNPs.cxx.

{
  if (affected_bins.empty()) return false;
  constexpr float inf{std::numeric_limits<float>::max()};
  constexpr float inf_{std::numeric_limits<float>::lowest()};
  Domain bAC[2] = {{inf, inf_, inf, inf_}, {inf, inf_, inf, inf_}};
  auto update = [&](const int side, const auto& dom,
                    const float etamin, const float etamax) {
    auto& b = bAC[side];
    b.etamin = std::min(b.etamin, etamin);
    b.etamax = std::max(b.etamax, etamax);
    b.ptmin = std::min(b.ptmin, dom.ptmin);
    b.ptmax = std::max(b.ptmax, dom.ptmax);
  };
 
  for (const int bin : affected_bins) {
    const Domain& dom{cfg.domains.at(bin)};
    if (dom.etamax > 0.f) update(1, dom, std::max(0.f, dom.etamin), dom.etamax);
    if (dom.etamin < 0.f) update(0, dom, dom.etamin, std::min(0.f, dom.etamax));
  }
 
  holes = false;
  const int ac{2*(bAC[1].ptmin!=inf) + (bAC[0].ptmin!=inf)};
  if (ac == 3) {
    abs_eta = (bAC[0].ptmin == bAC[1].ptmin) && (bAC[0].ptmax == bAC[1].ptmax) &&
      (std::fabs(bAC[0].etamax + bAC[1].etamin) < 1e-4f) && 
      (std::fabs(bAC[0].etamin + bAC[1].etamax) < 1e-4f);
    if (abs_eta) {
      bounds = bAC[1];
    } else {
      holes = true;
      bounds.etamin = bAC[0].etamin;
      bounds.etamax = bAC[1].etamax;
      bounds.ptmin = std::min(bAC[0].ptmin, bAC[1].ptmin);
      bounds.ptmax = std::min(bAC[0].ptmax, bAC[1].ptmax);
    }
  } else if (ac) bounds = bAC[ac - 1];
  else return false;
 
  if (!holes) {
    // Search for overlapping domains not affected by this NP.
    const std::size_t imax{cfg.domains.size()};
    auto f = [=](const float x) { return abs_eta? std::fabs(x) : x; };
    for (std::size_t i{0}; i<imax && !holes; ++i) {
      if (std::count(affected_bins.cbegin(), affected_bins.cend(), i)) continue;
      const auto& dom{cfg.domains[i]};
      holes = (dom.ptmax > bounds.ptmin) && 
        (dom.ptmin < bounds.ptmax) && 
        (f(dom.etamax) > bounds.etamin) &&
        (f(dom.etamin) < bounds.etamax);
    }
  }
  return true;
}

get_run_number()

int get_run_number

(

const int

year

)

Definition at line 576 of file dumpNPs.cxx.

{
  using namespace asg::msgUserCode;
  constexpr int random_runs[] = {280423, 302737, 332720, 354826};
  switch(year) {
  case 2015: case 2016: case 2017: case 2018: 
    return random_runs[year - 2015];
  default:
    Error(MSGSOURCE,
          "Invalid year specified, it should be between 2015 and 2018.");
    std::abort();
  }
}

main()

int main	(	int	argc,
		char *	argv[] )

Definition at line 85 of file dumpNPs.cxx.

{
  using namespace asg::msgUserCode;
 
  #ifdef XAOD_STANDALONE
    StatusCode::enableFailure();
    ANA_CHECK(xAOD::Init());
    xAOD::TEvent event{xAOD::TEvent::kClassAccess};
  #endif
 
  CP::CorrectionCode::enableFailure();
  //copy negated elements onto end of vector without duplicating the first element (which is zero)
  std::transform(eta_edges.begin()+1,eta_edges.end(), std::back_inserter(eta_edges),std::negate());
  //
  std::sort(eta_edges.begin(), eta_edges.end());
  std::sort(pt_edges.begin(), pt_edges.end());
  const bool duplicates_eta{
    std::adjacent_find(eta_edges.begin(), eta_edges.end()) != eta_edges.end()};
  const bool duplicates_pt{
    std::adjacent_find(pt_edges.begin(), pt_edges.end()) != pt_edges.end()};
  if (duplicates_eta || duplicates_pt) {
    Error(MSGSOURCE, "Duplicated values in the specified bin edges.");
    std::abort();
  }
 
  // Setup tool using command-line arguments
  Config cfg;
  ANA_CHECK( cfg.initialize(argc, argv) );
  ANA_CHECK( cfg.tool.initialize() );
 
  // List nuisance parameters for systematics
  const CP::SystematicSet allsysts(cfg.tool->recommendedSystematics());
  try {
    std::set<std::string> recorded;
    for (const auto& var : allsysts) {
      std::string name(var.basename());
      if (!recorded.emplace(name).second) continue;
      NP np;
      if (name.find("UncorrUncertainty") != std::string::npos) {
        np.uncorr = true;
      } else if (name.find("CorrUncertainty") != std::string::npos) {
        np.uncorr = false;
      } else {
        throw std::invalid_argument("");
      }
      np.index = std::stoul(name.substr(name.find("NP") + 2));
      cfg.systematics.emplace_back(CP::SystematicSet{var}, np);
    }
  } catch (std::invalid_argument&) {
    Error(MSGSOURCE, "The pattern of the SystematicVariation names seems to have changed.");
    std::abort();
  }
  std::sort(cfg.systematics.begin(), cfg.systematics.end(),
    [](auto& x, auto& y) { return std::get<NP>(x) < std::get<NP>(y); });
 
  if (cfg.analysis_mode) {
    for (const int year : {2015, 2016, 2017, 2018}) {
       map_t affected_bins;
       cfg.run_number = get_run_number(year);
      ANA_CHECK( scanPhaseSpace(cfg, affected_bins) );
      Info(MSGSOURCE, "%s", "");
      Info(MSGSOURCE, "Relevant nuisance parameters for the year %i:", year);
      ANA_CHECK( displayFindings_analysis(cfg, affected_bins) );
    }
  } else {
    map_t affected_bins;
    ANA_CHECK( scanPhaseSpace(cfg, affected_bins) );
    ANA_CHECK( displayFindings(cfg, affected_bins) );
  }
 
  return 0;
}

parse_csv_list()

template<typename... Args>

bool parse_csv_list	(	std::string	val,
		Args &...	args )

Definition at line 562 of file dumpNPs.cxx.

{
  std::replace(val.begin(), val.end(), ',', ' ');
  std::stringstream ss{val};
  while (true) {
    std::conditional_t<(sizeof...(Args)<2), float, std::string> x;
    ss >> x;
    if (ss.fail()) break;
    if (!parse_csv_token(x, args...)) return false;
  }
  return ss.eof();
}

parse_csv_token() [1/2]

bool parse_csv_token	(	const float	x,
		std::vector< float > &	bins )

Definition at line 513 of file dumpNPs.cxx.

{
  bins.push_back(x);
  return true;
}

parse_csv_token() [2/2]

bool parse_csv_token	(	std::string	s,
		std::vector< int8_t > &	flags,
		const std::vector< float > &	edges )

Definition at line 520 of file dumpNPs.cxx.

{
  const bool accept_aliases{&edges == &eta_edges};
  if (flags.size() != edges.size() - 1) return false;
  int alias{0};
  if (s == "nocrack") alias = 1;
  else if (s == "barrel") alias = 2;
  else if (s == "endcap") alias = 3;
  else if (s == "sym") alias = 4;
  if (!accept_aliases && alias>0) return false;
  if (alias>=1 && alias<=3) {
    for (std::size_t i{0}; i<edges.size()-1; ++i) {
      const float eta{std::fabs(edges[i])};
      if (alias==1 && eta>=1.37f && eta<1.52f) flags[i] = -1;
      if (alias==2 && eta<1.37f) flags[i] = 1;
      if (alias==3 && eta>1.52f) flags[i] = 1;
    }
    return true;
  } else if (alias == 4) {
    for (std::size_t i=0;i<edges.size();++i) {
      if (edges[i] < 0) flags[i] = 2;
    }
    return true;
  }
  auto pos = s.find("..");
  if (pos == std::string::npos) return false;
  std::stringstream ss{s.replace(pos, 2, 1, ' ')};
  float xmin, xmax;
  ss >> xmin >> xmax;
  if (ss.fail() || !(ss >> std::ws).eof()) return false;
  for (std::size_t i{0}; i<edges.size()-1; ++i) {
    if (xmin<edges[i+1] && xmax>edges[i] && flags[i]>=0) {
      flags[i] = 1;
    }
  }
  return true;
}

scanPhaseSpace()

bool scanPhaseSpace	(	Config &	config,
		map_t &	affected_bins )

Definition at line 279 of file dumpNPs.cxx.

{
  using namespace asg::msgUserCode;
  
  static int display_count{-1};
  if (!++display_count) {
    Info(MSGSOURCE,
         "Scanning the eta/pt plane, this will take a bit of time...");
  }
 
  xAOD::TStore store;
 
  using RealTool = const AsgElectronEfficiencyCorrectionTool;
  const auto *asgtool = dynamic_cast<RealTool*>(&*cfg.tool);
  const std::size_t jmax{subdiv_eta + 1}, kmax{subdiv_pt + 1};
  const std::size_t dmax{cfg.domains.size()};
  for (std::size_t d{0}; d<dmax; ++d) {
    const auto& dom = cfg.domains[d];
    std::bitset<1024> zero_uncorr{}, nonzero_uncorr{}, expected_uncorr{};
    std::bitset<256> zero_corr{}, nonzero_corr{};
    bool multi_np{false};
    for (std::size_t j{1}; j<jmax; ++j) {
      const float eta{((jmax-j)*dom.etamin + j*dom.etamax) / jmax};
      for (std::size_t k{1}; k<kmax; ++k) {
        const float pt{((kmax-k)*dom.ptmin + k*dom.ptmax) / kmax};
        ANA_CHECK( getElectrons({{pt, eta}}, cfg.run_number, store) );
        const xAOD::ElectronContainer* electrons{nullptr};
        ANA_CHECK( store.retrieve(electrons, "MyElectrons") );
        const xAOD::Electron& electron{*electrons->at(0)};
        double sf_nom{}, sf{};
        ANA_CHECK( cfg.tool->applySystematicVariation({}) );
        ANA_CHECK( cfg.tool->getEfficiencyScaleFactor(electron, sf_nom) );
        int vi{asgtool->systUncorrVariationIndex(electron)};
        expected_uncorr.set(vi);
        unsigned n_nonzero_uncorr{0};
        for (const auto& s : cfg.systematics) {
          const auto& sys = std::get<CP::SystematicSet>(s);
          ANA_CHECK( cfg.tool->applySystematicVariation(sys) );
          ANA_CHECK( cfg.tool->getEfficiencyScaleFactor(electron, sf) );
          const double delta{sf - sf_nom};
          const NP np{std::get<NP>(s)};
          if (delta != 0.) {
            if (np.uncorr) {
              nonzero_uncorr.set(np.index);
              ++n_nonzero_uncorr;
            }
            else nonzero_corr.set(np.index);
          } else {
            if (np.uncorr) zero_uncorr.set(np.index);
            else zero_corr.set(np.index);
          }
        }
        multi_np = multi_np || (n_nonzero_uncorr > 1);
        store.clear();
      }
    }
    if (multi_np) {
      std::string w("several uncorrelated NPs seem to affect the bin at "
        + dom.str() + ".");
      Warning(MSGSOURCE, "%s", w.c_str());
    } else if (nonzero_uncorr.count() > 1) {
      std::string w{"the binning used for the scan is unadapted at " + 
        dom.str() + " (too coarse?)."};
      Warning(MSGSOURCE, "%s", w.c_str());
    }
    if ((zero_uncorr & nonzero_uncorr).any() ||  
        (zero_corr & nonzero_corr).any()) {
      std::string w("the binning used for the scan is unadapted at " + 
        dom.str() + " (wrong boundaries?).");
      Warning(MSGSOURCE, "%s", w.c_str());
    }
    if ((nonzero_uncorr!=expected_uncorr) && nonzero_uncorr.any()) {
      std::string snz, se;
      for (std::size_t i{0}; i<nonzero_uncorr.size(); ++i) {
        std::string x{std::to_string(i) + ','};
        if (nonzero_uncorr[i] && !expected_uncorr[i]) snz += x;
        if (!nonzero_uncorr[i] && expected_uncorr[i]) se += x;
      }
      if (snz.length()) snz.pop_back();
      if (se.length()) se.pop_back();
      std::string w{"systUncorrVariationIndex() at " + dom.str() + 
        " indicates different NP(s) than those found in the scan, "
        "(%s) vs (%s)."};
      Warning(MSGSOURCE, w.c_str(), se.c_str(), snz.c_str());
    }
    for (const auto& s : cfg.systematics) {
      const NP np{std::get<NP>(s)};
      const bool impact{
        np.uncorr? nonzero_uncorr.test(np.index)
        : nonzero_corr.test(np.index)};
      if (impact) affected_bins[np].emplace_back(d);
    }
  }
  return true;
}

Variable Documentation

eta_edges

std::vector<float> eta_edges

Initial value:

= { 
  0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 
  1.15f, 1.2f, 1.3f, 1.37f, 1.52f, 1.6f, 1.65f, 1.7f, 1.8f, 1.9f, 2.0f, 
  2.1f, 2.2f, 2.3f, 2.37f, 2.4f, 2.47f}

Definition at line 71 of file dumpNPs.cxx.

                             { // eta<0 bins are auto-added in main()
  0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 
  1.15f, 1.2f, 1.3f, 1.37f, 1.52f, 1.6f, 1.65f, 1.7f, 1.8f, 1.9f, 2.0f, 
  2.1f, 2.2f, 2.3f, 2.37f, 2.4f, 2.47f};

pt_edges

std::vector<float> pt_edges

Initial value:

= { 
  45e2f, 7e3f, 1e4f, 15e3f, 2e4f, 25e3f, 3e4f, 35e3f, 4e4f, 45e3f, 5e4f,
  6e4f, 8e4f, 15e4f, 25e4f, 5e5f, 1e7f}

Definition at line 76 of file dumpNPs.cxx.

                            { // in MeV
  45e2f, 7e3f, 1e4f, 15e3f, 2e4f, 25e3f, 3e4f, 35e3f, 4e4f, 45e3f, 5e4f,
  6e4f, 8e4f, 15e4f, 25e4f, 5e5f, 1e7f};

subdiv_eta

const std::size_t subdiv_eta {5}

Definition at line 81 of file dumpNPs.cxx.

{5};

subdiv_pt

const std::size_t subdiv_pt {3}

Definition at line 82 of file dumpNPs.cxx.

{3};