egammaLayerRecalibTool.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include <iostream>
#include <exception>
#include <cassert>
#include <string>
#include <map>
#include <sstream>
#include <algorithm>
#include <cmath>
#include <limits>
 
#include <TFile.h>
#include <TObjString.h>
 
#include "PathResolver/PathResolver.h"
 
#include "egammaLayerRecalibTool/egammaLayerRecalibTool.h"
#include "xAODEgamma/EgammaxAODHelpers.h"
 
namespace {
const float VALUE_OVERFLOW = std::numeric_limits<float>::max();
 
template <typename TargetPtr, typename SourcePtr>
TargetPtr checked_cast(SourcePtr ptr) {
  // Do we have ptr types
  static_assert(std::is_pointer<TargetPtr>::value,
                "attempt to cast to no ptr object");
  static_assert(std::is_pointer<SourcePtr>::value,
                "attempt to cast from no ptr object");
 
  // nullptr input
  if (!ptr) {
    throw std::runtime_error(
        "Attempt to cast from nullptr in egammaLayerRecalibTool");
  }
 
  // dynamic_cast and check
  TargetPtr obj = dynamic_cast<TargetPtr>(ptr);
  if (not obj) {
    throw std::runtime_error("failed dynamic cast for " +
                             std::string(ptr->GetName()) +
                             " in egammaLayerRecalibTool");
  }
 
  return obj;
}
 
}  // end anonymous namespace
 
float GetAmountHisto1D::operator()(const StdCalibrationInputs& input) const {
  const int bin = m_histo->FindFixBin(input.eta);
  if (m_histo->IsBinUnderflow(bin) or m_histo->IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return m_histo->GetBinContent(bin);
}
 
float GetAmountHisto1DUp::operator()(const StdCalibrationInputs & input) const {
  const int bin = m_histo->FindFixBin(input.eta);
  if (m_histo->IsBinUnderflow(bin) or m_histo->IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return m_histo->GetBinContent(bin) + m_histo->GetBinError(bin);
}
 
float GetAmountHisto1DDown::operator()(const StdCalibrationInputs & input) const {
  const int bin = m_histo->FindFixBin(input.eta);
  if (m_histo->IsBinUnderflow(bin) or m_histo->IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return m_histo->GetBinContent(bin) - m_histo->GetBinError(bin);
}
 
float GetAmountHisto1DErrorUp::operator()(const StdCalibrationInputs & input) const {
  const int bin = m_histo->FindFixBin(input.eta);
  if (m_histo->IsBinUnderflow(bin) or m_histo->IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return m_histo->GetBinError(bin);
}
 
float GetAmountHisto1DErrorDown::operator()(const StdCalibrationInputs & input) const {
  const int bin = m_histo->FindFixBin(input.eta);
  if (m_histo->IsBinUnderflow(bin) or m_histo->IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return -m_histo->GetBinError(bin);
}
 
float GetAmountHisto2D::operator()(const StdCalibrationInputs & input) const {
  const int bin = m_histo.FindFixBin(input.eta, input.phi);
  if (m_histo.IsBinUnderflow(bin) or m_histo.IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return m_histo.GetBinContent(bin);
}
 
float GetAmountHisto2DEtaCaloRunNumber::operator()(const StdCalibrationInputs & input) const {
  const int bin = m_histo.FindFixBin(input.etaCalo, input.RunNumber);
  if (m_histo.IsBinUnderflow(bin) or m_histo.IsBinOverflow(bin)) return VALUE_OVERFLOW;
  return m_histo.GetBinContent(bin);
}
 
float GetAmountFixed::operator()(const StdCalibrationInputs & /*input*/ ) const {
  return m_amount;
}
 
float GetAmountFormula::operator()(const StdCalibrationInputs & input) const {
  return m_formula.Eval(input.eta, input.phi, input.RunNumber);
}
 
float GetAmountHVPSGuillaume::operator()(const StdCalibrationInputs & input) const {
  return m_tool.getCorr(input.RunNumber, input.eta, input.phi);
}
 
float GetAmountHVEMECPS207::operator() (const StdCalibrationInputs & input) const {
  return m_toolEMECPS.getCorr(input.RunNumber, input.eta, input.phi);
}
 
float GetAmountPileupE0::operator()(const StdCalibrationInputs & input) const {
  return m_tool->getCorr(0, input.RunNumber, input.averageInteractionsPerCrossing, input.eta);
}
 
float GetAmountPileupE1::operator()(const StdCalibrationInputs & input) const {
  return m_tool->getCorr(1, input.RunNumber, input.averageInteractionsPerCrossing, input.eta);
}
 
float GetAmountPileupE2::operator()(const StdCalibrationInputs & input) const {
  return m_tool->getCorr(2, input.RunNumber, input.averageInteractionsPerCrossing, input.eta);
}
 
float GetAmountPileupE3::operator()(const StdCalibrationInputs & input) const {
  return m_tool->getCorr(3, input.RunNumber, input.averageInteractionsPerCrossing, input.eta);
}
 
 
CP::CorrectionCode InputModifier::operator()(StdCalibrationInputs& inputs, float amount) const
{
  if (amount == VALUE_OVERFLOW) return CP::CorrectionCode::OutOfValidityRange;
  switch (m_base)
  {
    case SHIFT: shift_inputs(inputs, amount); return CP::CorrectionCode::Ok;
    case SUBTRACT: shift_inputs(inputs, -amount); return CP::CorrectionCode::Ok;
    case SCALE: scale_inputs(inputs, amount); return CP::CorrectionCode::Ok;
    case ZEROBASED: scale_inputs(inputs, 1. + amount); return CP::CorrectionCode::Ok;
    case ONEBASED: scale_inputs(inputs, amount); return CP::CorrectionCode::Ok;
    case ONEBASED_ALPHA: scale_inputs(inputs, 1. / amount); return CP::CorrectionCode::Ok;
    case ZEROBASED_ALPHA: scale_inputs(inputs, 1. / (1. + amount)); return CP::CorrectionCode::Ok;
    default: return CP::CorrectionCode::Error;
  };
}
 
void ScaleE0::scale_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E0raw *= amount; }
void ScaleE1::scale_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E1raw *= amount; }
void ScaleE2::scale_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E2raw *= amount; }
void ScaleE3::scale_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E3raw *= amount; }
 
void ScaleE0::shift_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E0raw += amount; }
void ScaleE1::shift_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E1raw += amount; }
void ScaleE2::shift_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E2raw += amount; }
void ScaleE3::shift_inputs(StdCalibrationInputs & inputs, float amount) const { inputs.E3raw += amount; }
 
 
void ScaleE1overE2::scale_inputs(StdCalibrationInputs & inputs, float amount) const
{
  const double Es1 = inputs.E1raw;
  const double Es2 = inputs.E2raw;
  if (Es1 == 0 and Es2 == 0) {
    inputs.E1raw = -999;
    inputs.E2raw = -999;
    return;
  }
  const double sum = Es1 + Es2;
  const double alpha = amount;
  const double den = (alpha * Es1 + Es2);
  inputs.E1raw = alpha * Es1 * sum / den;
  inputs.E2raw = Es2 * sum / den;
}
 
void ScaleE1overE2::shift_inputs(StdCalibrationInputs &, float) const
{
  // not very useful, never used
  throw std::runtime_error("not implemented");
}
 
void ScaleEaccordion::scale_inputs(StdCalibrationInputs & inputs, float amount) const
{
  inputs.E1raw *= amount;
  inputs.E2raw *= amount;
  inputs.E3raw *= amount;
}
 
void ScaleEaccordion::shift_inputs(StdCalibrationInputs & inputs, float amount) const
{
  inputs.E1raw += amount;
  inputs.E2raw += amount;
  inputs.E3raw += amount;
}
 
void ScaleEcalorimeter::scale_inputs(StdCalibrationInputs & inputs, float amount) const
{
  inputs.E0raw *= amount;
  inputs.E1raw *= amount;
  inputs.E2raw *= amount;
  inputs.E3raw *= amount;
}
 
void ScaleEcalorimeter::shift_inputs(StdCalibrationInputs & inputs, float amount) const
{
  inputs.E0raw += amount;
  inputs.E1raw += amount;
  inputs.E2raw += amount;
  inputs.E3raw += amount;
}
 
std::string egammaLayerRecalibTool::resolve_alias(const std::string& tune) {
 
  if ("layer1_2012" == tune)         return "layer1_2012_v5";
  if ("layer1_alt_2012" == tune)     return "layer1_alt_2012_v5";
  if ("layer1_2011" == tune)         return "layer1_2011_v5";
  if ("layer1_alt_2011" == tune)     return "layer1_alt_2011_v5";
  if ("layer1_2010" == tune)         return "layer1_2010_v5";
  if ("ps_2012" == tune)             return "ps_2012_v3";
  if ("ps_2011" == tune)             return "ps_2011_v3";
  if ("ps_2010" == tune)             return "ps_2010_v3";
  if ("layer1_2012_up" == tune)      return "layer1_2012_v5_up";
  if ("layer1_2012_down" == tune)    return "layer1_2012_v5_down";
  if ("layer1_2012_errup" == tune)   return "layer1_2012_v5_errup";
  if ("layer1_2012_errdown" == tune) return "layer1_2012_v5_errdown";
  if ("layer1_2011_up" == tune)      return "layer1_2011_v5_up";
  if ("layer1_2011_down" == tune)    return "layer1_2011_v5_down";
  if ("layer1_2011_errup" == tune)   return "layer1_2011_v5_errup";
  if ("layer1_2011_errdown" == tune) return "layer1_2011_v5_errdown";
  if ("layer1_2010_up" == tune)      return "layer1_2010_v5_up";
  if ("layer1_2010_down" == tune)    return "layer1_2010_v5_down";
  if ("layer1_2010_errup" == tune)   return "layer1_2010_v5_errup";
  if ("layer1_2010_errdown" == tune) return "layer1_2010_v5_errdown";
  if ("ps_2012_up" == tune)          return "ps_2012_v3_up";
  if ("ps_2012_down" == tune)        return "ps_2012_v3_down";
  if ("ps_2012_errup" == tune)       return "ps_2012_v3_errup";
  if ("ps_2012_errdown" == tune)     return "ps_2012_v3_errdown";
  if ("ps_2011_up" == tune)          return "ps_2011_v3_up";
  if ("ps_2011_down" == tune)        return "ps_2011_v3_down";
  if ("ps_2011_errup" == tune)       return "ps_2011_v3_errup";
  if ("ps_2011_errdown" == tune)     return "ps_2011_v3_errdown";
  if ("ps_2010_up" == tune)          return "ps_2010_v3_up";
  if ("ps_2010_down" == tune)        return "ps_2010_v3_down";
  if ("ps_2010_errup" == tune)       return "ps_2010_v3_errup";
  if ("ps_2010_errdown" == tune)     return "ps_2010_v3_errdown";
 
  return tune;
}
 
void egammaLayerRecalibTool::add_scale(const std::string& tuneIn)
{
  ATH_MSG_INFO("using scale " << tuneIn);
  std::string tune = resolve_alias(tuneIn);
 
  if (tune.empty()) { }
  else if ("es2025_run3_extrapolate_gnn_v0" == tune) {
    add_scale("run3_partial_ofc_extrapolate_gnn_v0");
  }
  else if ("es2024_run3_extrapolate_v0" == tune) {
    add_scale("run3_partial_ofc_extrapolate_v0");
  }
  // R22 layer tune with fixed E1E2 and repeated acc
  else if ("es2022_22.0_Precision_v1" == tune) {
    add_scale("run2_alt_with_layer2_r22_Precision_v1");
  }
  else if ("es2022_22.0_Precision" == tune) {
    add_scale("run2_alt_with_layer2_r22_Precision");
  }
  else if ("es2018_21.0_v0" == tune) {
    add_scale("run2_alt_with_layer2_r21_v1");
  }
  else if ("es2017_21.0_v0" == tune) {
    add_scale("run2_alt_with_layer2_r21_v0");
  }
  else if ("es2017_20.7_final" == tune) {
    add_scale("pileup_20.7");
    add_scale("run2_alt_with_layer2_modif");
  }
  else if ("es2017_20.7_improved" == tune) {
    add_scale("pileup_20.7");                // new pileup correction Guillaume for 20.7
    //TEMPORARY HACK REMOVED (two intermediate tags with this ES model use different layer corrections)
    add_scale("2012_alt_with_layer2_modif"); // temporary old corrections from run1 + EMECPS HV
  }
  else if ("pileup_20.7" == tune) {
    m_pileup_tool = new corr_pileupShift();
    add_scale(new ScaleE0(InputModifier::SUBTRACT), new GetAmountPileupE0(m_pileup_tool));
    add_scale(new ScaleE1(InputModifier::SUBTRACT), new GetAmountPileupE1(m_pileup_tool));
    add_scale(new ScaleE2(InputModifier::SUBTRACT), new GetAmountPileupE2(m_pileup_tool));
    add_scale(new ScaleE3(InputModifier::SUBTRACT), new GetAmountPileupE3(m_pileup_tool));
  }
  // Run3 2022+2023
  else if ("run3_partial_ofc_extrapolate_gnn_v0" == tune) {
    add_scale("layer2_run3_ofc_extrapolate_v0");
    add_scale("ps_run3_ofc_extrapolate_v0");
    if(m_doSaccCorrections) add_scale("acc_zee_run3_gnn_v0");
  }
  else if ("run3_partial_ofc_extrapolate_v0" == tune) {
    add_scale("layer2_run3_ofc_extrapolate_v0");
    add_scale("ps_run3_ofc_extrapolate_v0");
    if(m_doSaccCorrections) add_scale("acc_zee_run3_v0");
  }
  //Run 2 release 22 with fixed E1E2 and repeated acc 
  else if ("run2_alt_with_layer2_r22_Precision_v1"==tune) {
    add_scale("layer2_alt_el_mu_comb_r21_v0_fix");
    add_scale("ps_mu_r21_v0");
    if(m_doSaccCorrections) add_scale("acc_zee_r22_v1");
  }
  //Run 2 release 22
  else if ("run2_alt_with_layer2_r22_Precision"==tune) {
    add_scale("layer2_alt_el_mu_comb_r21_v0");
    add_scale("ps_mu_r21_v0");
    if(m_doSaccCorrections) add_scale("acc_zee_r22_v0");
  }
  else if ("run2_alt_with_layer2_r21_v1"==tune) {
    add_scale("layer2_alt_run2_r21_v1");
    add_scale("ps_2016_r21_v0");
  }
  else if ("run2_alt_with_layer2_r21_v0"==tune) {
    add_scale("layer2_alt_run2_r21_v0");
    add_scale("ps_2016_r21_v0");
  }
  else if("run2_alt_with_layer2_modif" == tune) {
    add_scale("ps_EMECHV1");
    add_scale("layer2_alt_run2_v1");
    add_scale("ps_2016");
  }
  // 2012
  else if ("2012" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012");
    add_scale("ps_2012");
  }
  else if("2012_with_layer2" == tune) {
    add_scale("ps_HV1");
    add_scale("layer2_2012_v5");
    add_scale("ps_2012");
  }
  else if ("2012_alt" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_alt_2012");
    add_scale("ps_2012");
  }
  else if("2012_alt_with_layer2" == tune) {
    add_scale("ps_HV1");
    add_scale("layer2_alt_2012_v5");
    add_scale("ps_2012");
  }
  else if("2012_alt_with_layer2_modif" == tune) {
    add_scale("ps_HV1");
    add_scale("ps_EMECHV1");
    add_scale("layer2_alt_2012_v5");
    add_scale("ps_2012");
  }
  else if("2010_with_layer2" == tune) {
    add_scale("layer2_2010_v5");
    add_scale("ps_2010");
  }
  else if ("2012_layer1_up" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012_up");
    add_scale("ps_2012");
  }
  else if ("2012_layer1_down" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012_down");
    add_scale("ps_2012");
  }
  else if ("2012_layer1_errup" == tune) {
    add_scale("layer1_2012_errup");
  }
  else if ("2012_layer1_errdown" == tune) {
    add_scale("layer1_2012_errdown");
  }
  else if ("2012_ps_down" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012");
    add_scale("ps_2012_down");
  }
  else if ("2012_ps_up" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012");
    add_scale("ps_2012_up");
  }
  else if ("2012_ps_errdown" == tune) {
    add_scale("ps_2012_errdown");
  }
  else if ("2012_ps_errup" == tune) {
    add_scale("ps_2012_errup");
  }
  else if ("2012_up" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012_up");
    add_scale("ps_2012_up");
  }
  else if ("2012_down" == tune) {
    add_scale("ps_HV1");
    add_scale("layer1_2012_down");
    add_scale("ps_2012_down");
  }
  else if ("2012_errup" == tune) {
    add_scale("layer1_2012_errup");
    add_scale("ps_2012_errup");
  }
  else if ("2012_errdown" == tune) {
    add_scale("layer1_2012_errdown");
    add_scale("ps_2012_errdown");
  }
  // 2011
  else if ("2011" == tune) {
    add_scale("layer1_2011");
    add_scale("ps_2011");
  }
  else if("2011_with_layer2" == tune) {
    add_scale("layer2_2011_v5");
    add_scale("ps_2011");
  }
  else if ("2011_alt" == tune) {
    add_scale("layer1_alt_2011");
    add_scale("ps_2011");
  }
  else if("2011_alt_with_layer2" == tune) {
    add_scale("layer2_alt_2011_v5");
    add_scale("ps_2011");
  }
  else if ("2011_layer1_up" == tune) {
    add_scale("layer1_2011_up");
    add_scale("ps_2011");
  }
  else if ("2011_layer1_down" == tune) {
    add_scale("layer1_2011_down");
    add_scale("ps_2011");
  }
  else if ("2011_layer1_errup" == tune) {
    add_scale("layer1_2011_errup");
  }
  else if ("2011_layer1_errdown" == tune) {
    add_scale("layer1_2011_errdown");
  }
  else if ("2011_ps_down" == tune) {
    add_scale("layer1_2011");
    add_scale("ps_2011_down");
  }
  else if ("2011_ps_up" == tune) {
    add_scale("layer1_2011");
    add_scale("ps_2011_up");
  }
  else if ("2011_ps_errdown" == tune) {
    add_scale("ps_2011_errdown");
  }
  else if ("2011_ps_errup" == tune) {
    add_scale("ps_2011_errup");
  }
  else if ("2011_up" == tune) {
    add_scale("layer1_2011_up");
    add_scale("ps_2011_up");
  }
  else if ("2011_down" == tune) {
    add_scale("layer1_2011_down");
    add_scale("ps_2011_down");
  }
  else if ("2011_errup" == tune) {
    add_scale("layer1_2011_errup");
    add_scale("ps_2011_errup");
  }
  else if ("2011_errdown" == tune) {
    add_scale("layer1_2011_errdown");
    add_scale("ps_2011_errdown");
  }
  // 2010
  else if ("2010" == tune) {
    add_scale("layer1_2010");
    add_scale("ps_2010");
  }
  else if ("2010_layer1_up" == tune) {
    add_scale("layer1_2010_up");
    add_scale("ps_2010");
  }
  else if ("2010_layer1_down" == tune) {
    add_scale("layer1_2010_down");
    add_scale("ps_2010");
  }
  else if ("2010_layer1_errup" == tune) {
    add_scale("layer1_2010_errup");
  }
  else if ("2010_layer1_errdown" == tune) {
    add_scale("layer1_2010_errdown");
  }
  else if ("2010_ps_down" == tune) {
    add_scale("layer1_2010");
    add_scale("ps_2010_down");
  }
  else if ("2010_ps_up" == tune) {
    add_scale("layer1_2010");
    add_scale("ps_2010_up");
  }
  else if ("2010_ps_errdown" == tune) {
    add_scale("ps_2010_errdown");
  }
  else if ("2010_ps_errup" == tune) {
    add_scale("ps_2010_errup");
  }
  else if ("2010_up" == tune) {
    add_scale("layer1_2010_up");
    add_scale("ps_2010_up");
  }
  else if ("2010_down" == tune) {
    add_scale("layer1_2010_down");
    add_scale("ps_2010_down");
  }
  else if ("2010_errup" == tune) {
    add_scale("layer1_2010_errup");
    add_scale("ps_2010_errup");
  }
  else if ("2010_errdown" == tune) {
    add_scale("layer1_2010_errdown");
    add_scale("ps_2010_errdown");
  }
  else if ("ps_HV1" == tune) {
    add_scale(new ScaleE0(InputModifier::ONEBASED), new GetAmountHVPSGuillaume());
  }
  else if ("ps_EMECHV1" == tune) {
    add_scale(new ScaleE0(InputModifier::ONEBASED), new GetAmountHVEMECPS207());
  }
  else if ("test1" == tune) {
    TH1F h_presampler("h_presampler", "h_presampler", 10, -2.5, 2.5);
    // just as an example, correct E0 by 0.1 * sign(eta)
    // and E1 by 1%
    for (int ibin = 1; ibin <= 5; ++ibin) {
      h_presampler.SetBinContent(ibin, -0.1);
      h_presampler.SetBinContent(ibin + 5, 0.1);
    }
    add_scale(new ScaleE0(InputModifier::ZEROBASED), new GetAmountHisto1D(h_presampler));
    add_scale(new ScaleE1(InputModifier::ZEROBASED), new GetAmountFixed(0.01));
  }
  else if ("acc_zee_run3_gnn_v0" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v14/egammaLayerRecalibTunes_transformerTune.root");
    TFile f(file.c_str());
    TH2F* histo_acc = static_cast<TH2F*>(f.Get("hACC_Zee_rel23_gnn"));
    assert(histo_acc);
    add_scale(new ScaleEaccordion(InputModifier::ZEROBASED_ALPHA),
              new GetAmountHisto2DEtaCaloRunNumber(*histo_acc));    
  }
  else if ("acc_zee_run3_v0" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v13/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH2F* histo_acc = static_cast<TH2F*>(f.Get("hACC_Zee_rel23"));
    assert(histo_acc);
    add_scale(new ScaleEaccordion(InputModifier::ZEROBASED_ALPHA),
              new GetAmountHisto2DEtaCaloRunNumber(*histo_acc));    
  }
  // repeated acc scale based on layer2_alt_el_mu_comb_r21_v0_fix 
  else if ("acc_zee_r22_v1" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v12/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH2F* histo_acc = static_cast<TH2F*>(f.Get("hACC_Zee_rel22"));
    assert(histo_acc);
    add_scale(new ScaleEaccordion(InputModifier::ZEROBASED_ALPHA),
              new GetAmountHisto2DEtaCaloRunNumber(*histo_acc));
  }
  else if ("acc_zee_r22_v0" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v11/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH2F* histo_acc = static_cast<TH2F*>(f.Get("hACC_Zee_rel22"));
    assert(histo_acc);
    add_scale(new ScaleEaccordion(InputModifier::ZEROBASED_ALPHA),
              new GetAmountHisto2DEtaCaloRunNumber(*histo_acc));
  }
  else if ("layer1_1" == tune) {
    TFormula f("formula_layer1_1", "(abs(x)<1.425) ? 0.97 : 1");
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA), new GetAmountFormula(f));
  }
  else if ("layer1_2" == tune) {
    TFormula f("formula_layer1_2", "(abs(x)<1.425) ? 0.97 : 1.05");
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA), new GetAmountFormula(f));
  }
  else if ("layer1_alt_2012_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_alt_2012"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo));
  }
  else if ("layer1_2012_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo));
  }
  else if ("layer1_2012_v5_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DUp(*histo));
  }
  else if ("layer1_2012_v5_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DDown(*histo));
  }
  else if ("layer1_2012_v5_errdown" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE1(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorUp(*histo));
  }
  else if ("layer1_2012_v5_errup" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE1(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorDown(*histo));
  }
  else if ("layer1_alt_2011_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_alt_2011"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo));
  }
  else if ("layer1_2011_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo));
  }
  else if ("layer1_2011_v5_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DUp(*histo));
  }
  else if ("layer1_2011_v5_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DDown(*histo));
  }
  else if ("layer1_2011_v5_errdown" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE1(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorUp(*histo));
  }
  else if ("layer1_2011_v5_errup" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE1(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorDown(*histo));
  }
  else if ("layer1_2010_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo));
  }
  else if ("layer1_2010_v5_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DUp(*histo));
  }
  else if ("layer1_2010_v5_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE1(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DDown(*histo));
  }
  else if ("layer1_2010_v5_errdown" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE1(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorUp(*histo));
  }
  else if ("layer1_2010_v5_errup" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE1(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorDown(*histo));
  }
  else if ("layer2_run3_ofc_extrapolate_v0"==tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v12/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1D* histo = static_cast<TH1D*>(f.Get("hE1E2_emu_run2_rel21_v1_run3ofc"));
    assert(histo);
    add_scale(new ScaleE2(InputModifier::ONEBASED),
         new GetAmountHisto1D(*histo));    
  }
  // fix E1E2 scale from R21 precision model
  else if("layer2_alt_el_mu_comb_r21_v0_fix"==tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v12/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1D* histo = static_cast<TH1D*>(f.Get("hE1E2_emu_run2_rel21_v0_fix"));
    assert(histo);
    add_scale(new ScaleE2(InputModifier::ONEBASED),
         new GetAmountHisto1D(*histo));
  }
  else if("layer2_alt_el_mu_comb_r21_v0"==tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v11/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1D* histo = static_cast<TH1D*>(f.Get("hE1E2_emu_run2_rel21_v0"));
    assert(histo);
    add_scale(new ScaleE2(InputModifier::ONEBASED),
         new GetAmountHisto1D(*histo));
  }
  else if("layer2_alt_run2_r21_v1"==tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v6/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2_mu_run2_rel21_v1"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
              new GetAmountHisto1D(*histo));
  }
  else if("layer2_alt_run2_r21_v0"==tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v5/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2mu_2016_rel21_v1"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
              new GetAmountHisto1D(*histo));
  }
  else if("layer2_alt_run2_v1" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v3/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2mu_2016_v1"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
              new GetAmountHisto1D(*histo));
  }
  else if("layer2_alt_2012_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_alt_2012"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1D(*histo));
  }
  else if("layer2_2012_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1D(*histo));
  }
  else if("layer2_2012_v5_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DDown(*histo));
  }
  else if("layer2_2012_v5_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DUp(*histo));
  }
  else if ("layer2_2012_v5_errdown" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DErrorDown(*histo));
  }
  else if ("layer2_2012_v5_errup" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2012"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DErrorUp(*histo));
  }
  else if("layer2_alt_2011_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_alt_2011"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1D(*histo));
  }
  else if("layer2_2011_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1D(*histo));
  }
  else if("layer2_2011_v5_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DDown(*histo));
  }
  else if("layer2_2011_v5_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DUp(*histo));
  }
  else if ("layer2_2011_v5_errdown" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DErrorDown(*histo));
  }
  else if ("layer2_2011_v5_errup" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2011"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DErrorUp(*histo));
  }
  else if("layer2_2010_v5" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1D(*histo));
  }
  else if("layer2_2010_v5_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DDown(*histo));
  }
  else if("layer2_2010_v5_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DUp(*histo));
  }
  else if ("layer2_2010_v5_errdown" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DErrorDown(*histo));
  }
  else if ("layer2_2010_v5_errup" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo = checked_cast<TH1*>(f.Get("hE1E2ave_2010"));
    add_scale(new ScaleE2(InputModifier::ONEBASED),
          new GetAmountHisto1DErrorUp(*histo));
  }
  else if ("ps_2016_r21_v0" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v5/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2016_rel21"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
              new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_run3_ofc_extrapolate_v0" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v12/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1F* histo_ps_tot_error = static_cast<TH1F*>(f.Get("hPS_MuonLowMu_rel21_run3ofc"));
    assert(histo_ps_tot_error);
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
              new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_mu_r21_v0" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v11/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1F* histo_ps_tot_error = static_cast<TH1F*>(f.Get("hPS_MuonLowMu_rel21"));
    assert(histo_ps_tot_error);
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
              new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_2016_v1" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v4/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2016"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_2012_v3" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2012"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_2012_v3_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2012"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DUp(*histo_ps_tot_error));
  }
  else if ("ps_2012_v3_up" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2012"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DDown(*histo_ps_tot_error));
  }
  else if ("ps_2012_v3_errdown" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2012"));
    add_scale(new ScaleE0(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorUp(*histo_ps_tot_error));
  }
  else if ("ps_2012_v3_errup" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2012"));
    add_scale(new ScaleE0(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorDown(*histo_ps_tot_error));
  }
  else if ("ps_2011_v3" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2011"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_2011_v3_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2011"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DUp(*histo_ps_tot_error));
  }
  else if ("ps_2011_v3_up" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2011"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DDown(*histo_ps_tot_error));
  }
  else if ("ps_2011_v3_errdown" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2011"));
    add_scale(new ScaleE0(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorUp(*histo_ps_tot_error));
  }
  else if ("ps_2011_v3_errup" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2011"));
    add_scale(new ScaleE0(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorDown(*histo_ps_tot_error));
  }
  // 2010
  else if ("ps_2010_v3" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2010"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1D(*histo_ps_tot_error));
  }
  else if ("ps_2010_v3_down" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2010"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DUp(*histo_ps_tot_error));
  }
  else if ("ps_2010_v3_up" == tune) {
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2010"));
    add_scale(new ScaleE0(InputModifier::ONEBASED_ALPHA),
          new GetAmountHisto1DDown(*histo_ps_tot_error));
  }
  else if ("ps_2010_v3_errdown" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2010"));
    add_scale(new ScaleE0(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorUp(*histo_ps_tot_error));
  }
  else if ("ps_2010_v3_errup" == tune){
    const std::string file = PathResolverFindCalibFile("egammaLayerRecalibTool/v1/egammaLayerRecalibTunes.root");
    TFile f(file.c_str());
    TH1* histo_ps_tot_error = checked_cast<TH1*>(f.Get("hPS_2010"));
    add_scale(new ScaleE0(InputModifier::ZEROBASED_ALPHA),
          new GetAmountHisto1DErrorDown(*histo_ps_tot_error));
  }
  else {
    throw std::runtime_error(tune+" is not a valid tune");
  }
}
 
egammaLayerRecalibTool::egammaLayerRecalibTool(const std::string& name, const std::string& tune, int SaccEnable)
  : asg::AsgMessaging(name), m_tune(tune), m_doSaccCorrections(SaccEnable)
{
  add_scale(tune);
}
 
 
egammaLayerRecalibTool::egammaLayerRecalibTool(const std::string& tune, int SaccEnable)
  : egammaLayerRecalibTool("egammaLayerRecalibTool", tune, SaccEnable) { }
 
 
void egammaLayerRecalibTool::add_scale(InputModifier* modifier, GetAmountBase* amount)
{
  m_modifiers.emplace_back(modifier, amount);
}
 
CP::CorrectionCode egammaLayerRecalibTool::scale_inputs(StdCalibrationInputs & inputs) const
{
  CP::CorrectionCode status = CP::CorrectionCode::Ok;
  for (const auto& modifier : m_modifiers)  {
    const float amount = (*modifier.second)(inputs);
    const auto s = (*modifier.first)(inputs, amount);
    ATH_MSG_DEBUG("  after E0|E1|E2|E3 = " << inputs.E0raw << "|" << inputs.E1raw << "|" << inputs.E2raw << "|" << inputs.E3raw);
    if (s != CP::CorrectionCode::Ok) {
      if (status != CP::CorrectionCode::Error) { status = s; }
    }
  }
  return status;
}
 
CP::CorrectionCode egammaLayerRecalibTool::read_and_scale_inputs( const xAOD::Egamma& particle,
                                                                  const xAOD::EventInfo& event_info,
                                                                  StdCalibrationInputs& inputs,
                                                                  bool& isData,
                                                                  std::string& fixT,
                                                                  double& addE2,
                                                                  double& addE3 ) const
{
  const xAOD::CaloCluster* cluster = particle.caloCluster();
  if (!cluster) {
    ATH_MSG_ERROR("egamma particle without CaloCluster");
    return CP::CorrectionCode::Error;
  }
 
  if (m_aodFixMissingCells &&
      event_info.runNumber() > m_Run2Run3runNumberTransition) {
    fixT = "_egFixForTopoTimingCut";
    unsigned short stat =
      xAOD::EgammaHelpers::energyInMissingCells(particle,addE2,addE3);
    if (stat) {
      ATH_MSG_WARNING("Fix for missing cells required"
              " but some layer info is not available,"
              " from L2 : " << stat%2 << " from L3 : " << stat/2);
    }
  }
 
  double eta_calo;
  static const SG::AuxElement::Accessor<float> accEtaCalo("etaCalo");
  if(particle.author() == xAOD::EgammaParameters::AuthorFwdElectron){
    eta_calo = cluster->eta();
  }
  else if (cluster->retrieveMoment(xAOD::CaloCluster::ETACALOFRAME, eta_calo)){
 
  }
  else if (accEtaCalo.isAvailable(*cluster)) {
    eta_calo = accEtaCalo(*cluster);
  }
  else{
    ATH_MSG_ERROR("etaCalo not available as auxilliary variable,"
          " using cluster eta as eta calo!");
    eta_calo=cluster->eta();
  }
 
  inputs = StdCalibrationInputs{
    event_info.averageInteractionsPerCrossing(),
    event_info.runNumber(),
    cluster->eta(),
    cluster->phi(),
    cluster->energyBE(0),
    cluster->energyBE(1),
    cluster->energyBE(2) + addE2,
    cluster->energyBE(3) + addE3,
    eta_calo};
 
  isData = !event_info.eventType(xAOD::EventInfo::IS_SIMULATION);
  CP::CorrectionCode status = CP::CorrectionCode::Ok;
  if (isData || m_scaleMC)
    status = scale_inputs(inputs);
 
  return status;
}
 
CP::CorrectionCode egammaLayerRecalibTool::applyCorrection(xAOD::Egamma& particle, const xAOD::EventInfo& event_info) const
{
 
  StdCalibrationInputs inputs{};
  bool isData = true;
  std::string fixT;
  double addE2 = 0.0, addE3 = 0.0;
 
  CP::CorrectionCode status = read_and_scale_inputs(particle, event_info, inputs, isData, fixT, addE2, addE3);
 
  if ( !m_doLayerclEdecoration )
    return status;
 
  const xAOD::CaloCluster* cluster = particle.caloCluster();
 
  static const SG::AuxElement::Decorator<double> deco_E0("correctedcl_Es0");
  static const SG::AuxElement::Decorator<double> deco_E1("correctedcl_Es1");
  static const SG::AuxElement::Decorator<double> deco_E2("correctedcl_Es2");
  static const SG::AuxElement::Decorator<double> deco_E3("correctedcl_Es3");
  static const SG::AuxElement::Decorator<std::string>
    deco_layer_correction("layer_correction");
 
  if (status == CP::CorrectionCode::Ok) {
    ATH_MSG_DEBUG("decorating cluster with corrected layer energies");
    deco_E0(*cluster) = m_doPSCorrections  ?
      inputs.E0raw : cluster->energyBE(0);
    deco_E1(*cluster) = m_doS12Corrections or m_doSaccCorrections ?
      inputs.E1raw : cluster->energyBE(1);
    deco_E2(*cluster) = m_doS12Corrections or m_doSaccCorrections ?
      inputs.E2raw : cluster->energyBE(2) + addE2;
    deco_E3(*cluster) = m_doSaccCorrections ?
      inputs.E3raw : cluster->energyBE(3) + addE3;
    deco_layer_correction(*cluster) = isData ? m_tune+fixT : fixT;
    return status;
  }
 
  ATH_MSG_DEBUG("cannot correct layer energies:"
        " decorating particle with non-corrected layer energies");
  // this is done for safety, since when a particle is decorated
  // all the particle in the container are decorated
  // it is not possible to distinguish between decorated / non-decorated
  // since all are decorated
  deco_E0(*cluster) = cluster->energyBE(0);
  deco_E1(*cluster) = cluster->energyBE(1);
  deco_E2(*cluster) = cluster->energyBE(2) + addE2;
  deco_E3(*cluster) = cluster->energyBE(3) + addE3;
  deco_layer_correction(*cluster) = isData ? m_tune + "_Err" + fixT : fixT;
  return status;
 
}
 
std::array<double,4> egammaLayerRecalibTool::getLayerCorrections(const xAOD::Egamma& particle, const xAOD::EventInfo& event_info) const
{
  StdCalibrationInputs inputs{};
  bool isData = true;
  std::string fixT;
  double addE2 = 0.0, addE3 = 0.0;
 
  CP::CorrectionCode status = read_and_scale_inputs(particle, event_info, inputs, isData, fixT, addE2, addE3);
 
  if (status != CP::CorrectionCode::Ok) {
    if (status == CP::CorrectionCode::Error) {
      ATH_MSG_ERROR("Failed to read inputs, returning 1 for all layers, please check if this is expected!");
    }
    else if (status == CP::CorrectionCode::OutOfValidityRange) {
      ATH_MSG_WARNING("Some inputs are out of validity range, returning 1 for all layers; object eta: " << inputs.eta);
    }
    return {1.,1.,1.,1.};
  }
  
  auto safe_divide = [](float a, float b) { return b != 0 ? a/b : 1.f; };
  return {
    safe_divide(inputs.E0raw, particle.caloCluster()->energyBE(0)),
    safe_divide(inputs.E1raw, particle.caloCluster()->energyBE(1)),
    safe_divide(inputs.E2raw, particle.caloCluster()->energyBE(2) + addE2),
    safe_divide(inputs.E3raw, particle.caloCluster()->energyBE(3) + addE3) 
  };
 
}
 
void egammaLayerRecalibTool::clear_corrections()
{
  for (auto modifier :  m_modifiers) {
    delete modifier.first;
    delete modifier.second;
  }
  m_modifiers.clear();
}
 
 
// helper
std::map<std::string, std::string> parse(const std::string& list)
{
  std::cout << "list: '" << list << "'" << std::endl;
  std::map<std::string, std::string> result;
  TIter next(TString(list).Tokenize(","));
  while (TObjString* sObj = (TObjString*) next())
    {
      const TString& item(sObj->GetString());
      std::cout << "item: '" << item << "'" << std::endl;
      TObjArray* item_list = TString(item).Tokenize(":");
      std::string key;
      std::string value;
      if (item_list->GetEntries() == 1) {
    key = "amount";
    value = static_cast<TObjString*>(item_list->At(0))->GetString().Data();
      }
      else if (item_list->GetEntries() == 2) {
    key = static_cast<TObjString*>(item_list->At(0))->GetString().Data();
    value = static_cast<TObjString*>(item_list->At(1))->GetString().Data();
      }
      else {
    std::cerr << "invalid string " << item << std::endl;
      }
      if (result.find(key) != result.end()) {
    std::cerr << "trying to insert two times key " << key << std::endl;
    assert(false);
      }
      result.insert(std::make_pair(key, value));
    }
  return result;
}
 
 
std::pair<std::string, egammaLayerRecalibTool*>
egammaLayerRecalibTool::create(const std::string& type, const std::string& args)
{
  std::map<std::string, std::string> args_map = parse(args);
  egammaLayerRecalibTool* tool = new egammaLayerRecalibTool("egammaLayerRecalibTool", "");
  std::string name = "";
  std::string amount_name = "";
  std::string type_name = "";
 
  GetAmountBase* amount_getter = nullptr;
  InputModifier* modifier = nullptr;
 
  if (args_map.find("amount") != args_map.end()) {
    std::string amount_str = args_map["amount"];
    bool perc = false;
    if (amount_str.back()=='%') {
      perc = true;
      amount_str.pop_back();
    }
    const float amount = TString(amount_str).Atof() * (perc ? 0.01 : 1);
 
    amount_getter = new GetAmountFixed(amount);
    std::stringstream amount_stream;
    amount_stream << amount;
    amount_name = amount_stream.str();
    std::replace(amount_name.begin(), amount_name.end(), '-', 'n');
    std::replace(amount_name.begin(), amount_name.end(), '.', 'p');
  }
  else if (args_map.find("name") != args_map.end()) {
    name = args_map["name"];
  }
  else if (args_map.find("histo") != args_map.end()) {
    int dim = 0;
    if (args_map.find("file") == args_map.end()) {
      std::cerr << "with histo you must specify file" << std::endl;
      assert(false);
    }
    if (args_map.find("formulax") != args_map.end()) dim = 1;
 
    if (dim == 0)
      {
    std::cerr << "with histo you must specify formulax" << std::endl;
    assert(false);
      }
    if (dim == 1) {
      TFile f(args_map["file"].c_str());
      std::cout << "opening histo " << args_map["histo"] << " from file " << args_map["file"] << std::endl;
      TH1F* histo = dynamic_cast<TH1F*>(f.Get(args_map["histo"].c_str()));
 
      if(histo){
    histo->SetDirectory(nullptr);
    amount_getter = new GetAmountHisto1D(*histo);
      }
      else{assert(false); }
    }
    else { assert(false); }
  }
  else {
    std::cerr << "cannot understand argument " << args << std::endl;
    assert(false);
  }
 
  if ("bias-E0" == type)                { modifier = new ScaleE0(InputModifier::ZEROBASED); type_name = "E0"; }
  else if ("bias-E1" == type)           { modifier = new ScaleE1(InputModifier::ZEROBASED); type_name = "E1"; }
  else if ("bias-E2" == type)           { modifier = new ScaleE2(InputModifier::ZEROBASED); type_name = "E2"; }
  else if ("bias-E3" == type)           { modifier = new ScaleE3(InputModifier::ZEROBASED); type_name = "E3"; }
  else if ("bias-E1overE2" == type)     { modifier = new ScaleE1overE2(InputModifier::ZEROBASED); type_name = "E1overE2"; }
  else if ("bias-Eaccordion" == type)   { modifier = new ScaleEaccordion(InputModifier::ZEROBASED); type_name = "Eaccordion"; }
  else if ("bias-Ecalorimeter" == type) { modifier = new ScaleEcalorimeter(InputModifier::ZEROBASED); type_name = "Ecalorimeter"; }
 
  if (not type_name.empty() and not amount_name.empty()) {
    name = type_name + "_" + amount_name;
  }
 
  if (name.empty()) {
    std::cerr << "you need to specify a name for the bias with type " << type << std::endl;
  }
 
  if (modifier and amount_getter) {
    tool->add_scale(modifier, amount_getter);
  }
  else{
    tool->add_scale(type);
    //release resources, if modifier false need to release amount_getter and vice versa
    //since they are not passed to the tool
    if(modifier) delete modifier;
    if(amount_getter) delete  amount_getter;
  }
 
  return {name, tool};
}