L1GetHistLimits.cc

Go to the documentation of this file.

// this class header
#include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GetHistLimits.h"

// system include files
#include <iostream>
#include <iomanip>
#include <string>

// user include files
#include "L1Trigger/GlobalTriggerAnalyzer/interface/L1PhiConversion.h"
#include "DataFormats/Math/interface/normalizedPhi.h"

#include "DataFormats/L1GlobalCaloTrigger/interface/L1GctEtMiss.h"
#include "DataFormats/L1GlobalCaloTrigger/interface/L1GctEtTotal.h"
#include "DataFormats/L1GlobalCaloTrigger/interface/L1GctEtHad.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

L1GetHistLimits::Tokens::Tokens(edm::ConsumesCollector iC, bool doEtaOrPhi)
    : m_muPTScaleToken(iC.esConsumes<edm::Transition::BeginRun>()),
      m_etScaleToken(iC.esConsumes<edm::Transition::BeginRun>()),
      m_jetScaleToken(iC.esConsumes<edm::Transition::BeginRun>()),
      m_jetFinderParamsToken(iC.esConsumes<edm::Transition::BeginRun>()),
      m_htMissScaleToken(iC.esConsumes<edm::Transition::BeginRun>()),
      m_hfRingEtScaleToken(iC.esConsumes<edm::Transition::BeginRun>()) {
  if (doEtaOrPhi) {
    m_muScalesToken = iC.esConsumes<edm::Transition::BeginRun>();
    m_caloGeomESHToken = iC.esConsumes<edm::Transition::BeginRun>();
  }
}

// constructor
L1GetHistLimits::L1GetHistLimits(const Tokens& tokens, const edm::EventSetup& evSetup)
    : m_tokens(tokens), m_evSetup(evSetup) {
  //
}

// destructor
L1GetHistLimits::~L1GetHistLimits() {
  // empty
}

void L1GetHistLimits::getHistLimits(const L1GtObject& l1GtObject, const std::string& quantity) {
  m_l1HistLimits.nrBins = 0;
  m_l1HistLimits.lowerBinValue = 0;
  m_l1HistLimits.upperBinValue = 0;
  m_l1HistLimits.binThresholds.clear();

  // number of objects is independent of the object type
  // and hardcoded in the actual version
  if (quantity == "NrObjects") {
    m_l1HistLimits.nrBins = 15;
    m_l1HistLimits.lowerBinValue = -0.5;
    m_l1HistLimits.upperBinValue = 14.5;

    m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

    for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
      m_l1HistLimits.binThresholds[iBin] =
          m_l1HistLimits.lowerBinValue +
          iBin * (m_l1HistLimits.upperBinValue - m_l1HistLimits.lowerBinValue) / m_l1HistLimits.nrBins;
    }

    // set last bin upper edge
    m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

    return;
  }

  switch (l1GtObject) {
    case Mu: {
      if (quantity == "PT") {
        const L1MuTriggerPtScale& muPtScale = m_evSetup.getData(m_tokens.m_muPTScaleToken);

        m_l1HistLimits.nrBins = muPtScale.getPtScale()->getNBins();
        m_l1HistLimits.lowerBinValue = muPtScale.getPtScale()->getScaleMin();
        m_l1HistLimits.upperBinValue = muPtScale.getPtScale()->getScaleMax();

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = muPtScale.getPtScale()->getValue(iBin);
        }

        // last limit for muon is set too high (10^6) - resize the last bin

        float lastBinSize = m_l1HistLimits.upperBinValue - m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins - 1];

        // limit the size of the last bin to maximum 200 GeV
        float maxLaxtBinsize = 200;
        if (lastBinSize >= maxLaxtBinsize) {
          m_l1HistLimits.upperBinValue = m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins - 1] +
                                         2. * (m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins - 1] -
                                               m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins - 2]);
          LogDebug("L1GetHistLimits") << "\n L1ExtraMuon: PT histogram"
                                      << "\nm_l1HistLimits.upperBinValue truncated to = "
                                      << m_l1HistLimits.upperBinValue << std::endl;
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        const L1MuTriggerScales& muScales = m_evSetup.getData(m_tokens.m_muScalesToken);

        if (quantity == "eta") {
          // eta scale defined for positive values - need to be symmetrized
          int histNrBinsHalf = muScales.getGMTEtaScale()->getNBins();
          m_l1HistLimits.lowerBinValue = muScales.getGMTEtaScale()->getScaleMin();
          m_l1HistLimits.upperBinValue = muScales.getGMTEtaScale()->getScaleMax();

          m_l1HistLimits.nrBins = 2 * histNrBinsHalf;
          m_l1HistLimits.lowerBinValue = -m_l1HistLimits.upperBinValue;

          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          int iBin = 0;
          for (int j = histNrBinsHalf; j > 0; j--, iBin++) {
            m_l1HistLimits.binThresholds[iBin] = (-1) * muScales.getGMTEtaScale()->getValue(j);
          }
          for (int j = 0; j <= histNrBinsHalf; j++, iBin++) {
            m_l1HistLimits.binThresholds[iBin] = muScales.getGMTEtaScale()->getValue(j);
          }

          // set last bin upper edge
          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

        } else {
          m_l1HistLimits.nrBins = muScales.getPhiScale()->getNBins();
          m_l1HistLimits.lowerBinValue = rad2deg(muScales.getPhiScale()->getScaleMin());
          m_l1HistLimits.upperBinValue = rad2deg(muScales.getPhiScale()->getScaleMax());

          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin <= m_l1HistLimits.nrBins; iBin++) {
            m_l1HistLimits.binThresholds[iBin] = rad2deg(muScales.getPhiScale()->getValue(iBin));
          }

          // set last bin upper edge
          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;
        }
      }

    } break;
    case NoIsoEG:
    case IsoEG: {
      // common scales for NoIsoEG and IsoEG
      if (quantity == "ET") {
        const L1CaloEtScale& emScale = m_evSetup.getData(m_tokens.m_etScaleToken);

        std::vector<double> emThresholds = emScale.getThresholds();
        m_l1HistLimits.nrBins = emThresholds.size();
        m_l1HistLimits.lowerBinValue = emThresholds.at(0);

        // FIXME high edge retrieval in the scale definition
        // now, last bin has the same width like the last but one
        m_l1HistLimits.upperBinValue =
            emThresholds[m_l1HistLimits.nrBins - 1] +
            (emThresholds[m_l1HistLimits.nrBins - 1] - emThresholds[m_l1HistLimits.nrBins - 2]);

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = static_cast<float>(emThresholds[iBin]);
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        const L1CaloGeometry* caloGeomScales = &m_evSetup.getData(m_tokens.m_caloGeomESHToken);

        if (quantity == "eta") {
          m_l1HistLimits.nrBins =
              2 * (caloGeomScales->numberGctCentralEtaBinsPerHalf() + caloGeomScales->numberGctForwardEtaBinsPerHalf());
          m_l1HistLimits.lowerBinValue = caloGeomScales->globalEtaBinLowEdge(0);
          m_l1HistLimits.upperBinValue = caloGeomScales->globalEtaBinHighEdge(m_l1HistLimits.nrBins - 1);

          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
            m_l1HistLimits.binThresholds[iBin] = caloGeomScales->globalEtaBinLowEdge(iBin);
          }

          // set last bin upper edge
          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

        } else {
          m_l1HistLimits.nrBins = caloGeomScales->numberGctEmJetPhiBins();
          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
            m_l1HistLimits.binThresholds[iBin] = rad2deg(caloGeomScales->emJetPhiBinLowEdge(iBin));

            // treat correctly the 10 deg anti-clockwise rotation
            if (rad2deg(caloGeomScales->emJetPhiBinHighEdge(iBin)) < m_l1HistLimits.binThresholds[iBin]) {
              m_l1HistLimits.binThresholds[iBin] = m_l1HistLimits.binThresholds[iBin] - 360.;
            }
          }

          m_l1HistLimits.lowerBinValue = m_l1HistLimits.binThresholds[0];

          // last bin upper limit
          m_l1HistLimits.upperBinValue = rad2deg(caloGeomScales->emJetPhiBinHighEdge(m_l1HistLimits.nrBins - 1));

          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;
        }
      }

    } break;
    case CenJet:
    case ForJet:
    case TauJet: {
      // common scales for all jets
      if (quantity == "ET") {
        const L1CaloEtScale& jetScale = m_evSetup.getData(m_tokens.m_jetScaleToken);

        std::vector<double> jetThresholds = jetScale.getThresholds();
        m_l1HistLimits.nrBins = jetThresholds.size();
        m_l1HistLimits.lowerBinValue = jetThresholds.at(0);

        // FIXME high edge retrieval in the scale definition
        // now, last bin has the same width like the last but one
        m_l1HistLimits.upperBinValue =
            jetThresholds[m_l1HistLimits.nrBins - 1] +
            (jetThresholds[m_l1HistLimits.nrBins - 1] - jetThresholds[m_l1HistLimits.nrBins - 2]);

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = static_cast<float>(jetThresholds[iBin]);
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        const L1CaloGeometry* caloGeomScales = &m_evSetup.getData(m_tokens.m_caloGeomESHToken);

        if (quantity == "eta") {
          m_l1HistLimits.nrBins =
              2 * (caloGeomScales->numberGctCentralEtaBinsPerHalf() + caloGeomScales->numberGctForwardEtaBinsPerHalf());
          m_l1HistLimits.lowerBinValue = caloGeomScales->globalEtaBinLowEdge(0);
          m_l1HistLimits.upperBinValue = caloGeomScales->globalEtaBinHighEdge(m_l1HistLimits.nrBins - 1);

          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
            m_l1HistLimits.binThresholds[iBin] = caloGeomScales->globalEtaBinLowEdge(iBin);
          }

          // set last bin upper edge
          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

        } else {
          m_l1HistLimits.nrBins = caloGeomScales->numberGctEmJetPhiBins();
          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
            m_l1HistLimits.binThresholds[iBin] = rad2deg(caloGeomScales->emJetPhiBinLowEdge(iBin));

            // treat correctly the 10 deg anti-clockwise rotation
            if (rad2deg(caloGeomScales->emJetPhiBinHighEdge(iBin)) < m_l1HistLimits.binThresholds[iBin]) {
              m_l1HistLimits.binThresholds[iBin] = m_l1HistLimits.binThresholds[iBin] - 360.;
            }
          }

          m_l1HistLimits.lowerBinValue = m_l1HistLimits.binThresholds[0];

          // last bin upper limit
          m_l1HistLimits.upperBinValue = rad2deg(caloGeomScales->emJetPhiBinHighEdge(m_l1HistLimits.nrBins - 1));

          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;
        }
      }
    } break;
    case ETM: {
      if (quantity == "ET") {
        const L1CaloEtScale& etMissScale = m_evSetup.getData(m_tokens.m_jetScaleToken);

        const double etSumLSB = etMissScale.linearLsb();

        m_l1HistLimits.nrBins = L1GctEtMiss::kEtMissMaxValue;

        m_l1HistLimits.lowerBinValue = 0;
        m_l1HistLimits.upperBinValue = (m_l1HistLimits.nrBins + 1) * etSumLSB;

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = iBin * etSumLSB;
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        const L1CaloGeometry* caloGeomScales = &m_evSetup.getData(m_tokens.m_caloGeomESHToken);

        if (quantity == "eta") {
          // do nothing, eta is not defined for ETM

        } else {
          m_l1HistLimits.nrBins = caloGeomScales->numberGctEtSumPhiBins();
          m_l1HistLimits.lowerBinValue = rad2deg(caloGeomScales->etSumPhiBinLowEdge(0));
          m_l1HistLimits.upperBinValue = rad2deg(caloGeomScales->etSumPhiBinHighEdge(m_l1HistLimits.nrBins - 1));

          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
            m_l1HistLimits.binThresholds[iBin] = rad2deg(caloGeomScales->etSumPhiBinLowEdge(iBin));
          }

          // last bin upper limit
          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;
        }
      }

    } break;
    case ETT: {
      if (quantity == "ET") {
        const L1CaloEtScale& etMissScale = m_evSetup.getData(m_tokens.m_jetScaleToken);

        const double etSumLSB = etMissScale.linearLsb();

        m_l1HistLimits.nrBins = L1GctEtTotal::kEtTotalMaxValue;

        m_l1HistLimits.lowerBinValue = 0;
        m_l1HistLimits.upperBinValue = (m_l1HistLimits.nrBins + 1) * etSumLSB;

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = iBin * etSumLSB;
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        // do nothing, eta and phi are not defined for ETT
      }

    } break;
    case HTT: {
      if (quantity == "ET") {
        const L1GctJetFinderParams& jetFinderParams = m_evSetup.getData(m_tokens.m_jetFinderParamsToken);
        double htSumLSB = jetFinderParams.getHtLsbGeV();

        m_l1HistLimits.nrBins = L1GctEtHad::kEtHadMaxValue;

        m_l1HistLimits.lowerBinValue = 0;
        m_l1HistLimits.upperBinValue = (m_l1HistLimits.nrBins + 1) * htSumLSB;

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = iBin * htSumLSB;
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        // do nothing, eta and phi are not defined for HTT
      }
    } break;
    case HTM: {
      if (quantity == "ET") {
        const L1CaloEtScale& htMissScale = m_evSetup.getData(m_tokens.m_htMissScaleToken);

        const std::vector<double>& htThresholds = htMissScale.getThresholds();
        m_l1HistLimits.nrBins = htThresholds.size();
        m_l1HistLimits.lowerBinValue = htThresholds[0];

        // FIXME high edge retrieval in the scale definition
        // now, last bin has the same width like the last but one
        m_l1HistLimits.upperBinValue =
            htThresholds[m_l1HistLimits.nrBins - 1] +
            (htThresholds[m_l1HistLimits.nrBins - 1] - htThresholds[m_l1HistLimits.nrBins - 2]);

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = static_cast<float>(htThresholds[iBin]);
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        const L1CaloGeometry* caloGeomScales = &m_evSetup.getData(m_tokens.m_caloGeomESHToken);

        if (quantity == "eta") {
          // do nothing, eta is not defined for HTM

        } else {
          m_l1HistLimits.nrBins = caloGeomScales->numberGctHtSumPhiBins();
          m_l1HistLimits.lowerBinValue = rad2deg(caloGeomScales->htSumPhiBinLowEdge(0));
          m_l1HistLimits.upperBinValue = rad2deg(caloGeomScales->htSumPhiBinHighEdge(m_l1HistLimits.nrBins - 1));

          m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

          for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
            m_l1HistLimits.binThresholds[iBin] = rad2deg(caloGeomScales->htSumPhiBinLowEdge(iBin));
          }

          // last bin upper limit
          m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;
        }
      }
    } break;
    case JetCounts: {
    } break;
    case HfBitCounts: {
      // there are no scales for HfBitCounts, so one implements a fixed scale
      // use same values as GCT for 3 bits
      const unsigned int R3BINS = 8;
      const float R3MIN = -0.5;
      const float R3MAX = 7.5;

      m_l1HistLimits.nrBins = R3BINS;
      m_l1HistLimits.lowerBinValue = R3MIN;
      m_l1HistLimits.upperBinValue = R3MAX;

      m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

      for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
        m_l1HistLimits.binThresholds[iBin] = R3MIN + iBin * (R3MAX - R3MIN) / R3BINS;
      }

      // last bin upper limit
      m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

    } break;
    case HfRingEtSums: {
      if (quantity == "ET") {
        const L1CaloEtScale& hfRingEtScale = m_evSetup.getData(m_tokens.m_hfRingEtScaleToken);

        const std::vector<double>& hfRingEtThresholds = hfRingEtScale.getThresholds();
        m_l1HistLimits.nrBins = hfRingEtThresholds.size();
        m_l1HistLimits.lowerBinValue = hfRingEtThresholds[0];

        // FIXME high edge retrieval in the scale definition
        // now, last bin has the same width like the last but one
        m_l1HistLimits.upperBinValue =
            hfRingEtThresholds[m_l1HistLimits.nrBins - 1] +
            (hfRingEtThresholds[m_l1HistLimits.nrBins - 1] - hfRingEtThresholds[m_l1HistLimits.nrBins - 2]);

        m_l1HistLimits.binThresholds.resize(m_l1HistLimits.nrBins + 1);

        for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
          m_l1HistLimits.binThresholds[iBin] = static_cast<float>(hfRingEtThresholds[iBin]);
        }

        // set last bin upper edge
        m_l1HistLimits.binThresholds[m_l1HistLimits.nrBins] = m_l1HistLimits.upperBinValue;

      } else if (quantity == "eta" || quantity == "phi") {
        // do nothing, eta and phi are not defined for HfRingEtSums
      }
    } break;
    case TechTrig:
    case Castor:
    case BPTX:
    default: {
      // do nothing, for these cases ET/PT, eta and phi are not defined

    } break;
  }
}

const L1GetHistLimits::L1HistLimits& L1GetHistLimits::l1HistLimits(const L1GtObject& l1GtObject,
                                                                   const std::string& quantity) {
  getHistLimits(l1GtObject, quantity);

  if (edm::isDebugEnabled()) {
    LogDebug("L1GetHistLimits") << "\n Histogram limits for L1GtObject " << l1GtObject << " and quantity " << quantity
                                << "\n  Number of bins:           " << m_l1HistLimits.nrBins
                                << "\n  Lower limit of first bin: " << m_l1HistLimits.lowerBinValue
                                << "\n  Upper limit of last bin:  " << m_l1HistLimits.upperBinValue << std::endl;

    int binThreshSize = static_cast<int>(m_l1HistLimits.binThresholds.size());

    if (binThreshSize != (m_l1HistLimits.nrBins + 1)) {
      LogTrace("L1GetHistLimits") << "\n Warning: inconsistent nrBins and binThresholds size"
                                  << "\n   Number of bins nrBins = " << m_l1HistLimits.nrBins
                                  << "\n   binThresholds size =    " << binThreshSize
                                  << "\n Please fix the L1GetLimits class.\n\n"
                                  << std::endl;
    }

    for (int iBin = 0; iBin < binThreshSize; ++iBin) {
      LogTrace("L1GetHistLimits") << " Bin " << std::right << std::setw(5) << iBin << ":  "
                                  << m_l1HistLimits.binThresholds[iBin] << std::endl;
    }
  }

  return m_l1HistLimits;
}

const L1GetHistLimits::L1HistLimits& L1GetHistLimits::l1HistLimits(const L1GtObject& l1GtObject,
                                                                   const std::string& quantity,
                                                                   const double histMinValue,
                                                                   const double histMaxValue) {
  getHistLimits(l1GtObject, quantity);

  bool foundLowerBinValue = false;
  bool foundUpperBinValue = false;

  for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
    if (m_l1HistLimits.binThresholds[iBin] <= histMinValue) {
      m_l1HistLimits.lowerBinValue = m_l1HistLimits.binThresholds[iBin];
      foundLowerBinValue = true;
      break;
    }
  }

  for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
    if (m_l1HistLimits.binThresholds[iBin] > histMaxValue) {
      m_l1HistLimits.upperBinValue = m_l1HistLimits.binThresholds[iBin];
      foundUpperBinValue = true;
      break;
    }
  }

  if (foundLowerBinValue && foundUpperBinValue) {
    int countBins = -1;
    std::vector<float> binThresh;
    binThresh.reserve(m_l1HistLimits.binThresholds.size());

    for (int iBin = 0; iBin < m_l1HistLimits.nrBins; ++iBin) {
      if ((m_l1HistLimits.binThresholds[iBin] >= histMinValue) && m_l1HistLimits.binThresholds[iBin] < histMaxValue) {
        m_l1HistLimits.upperBinValue = m_l1HistLimits.binThresholds[iBin];

        countBins++;
        binThresh.push_back(m_l1HistLimits.binThresholds[iBin]);
      }
    }

    m_l1HistLimits.nrBins = countBins;
    m_l1HistLimits.binThresholds.clear();
    m_l1HistLimits.binThresholds = binThresh;

    // FIXME last bin untested.

  } else {
    m_l1HistLimits.nrBins = 0;
    m_l1HistLimits.lowerBinValue = 0;
    m_l1HistLimits.upperBinValue = 0;
    m_l1HistLimits.binThresholds.clear();

    LogDebug("L1GetHistLimits") << "\n Histogram limits for L1GtObject" << l1GtObject << " and quantity " << quantity
                                << " within the required range [" << histMinValue << ", " << histMaxValue
                                << "] not found."
                                << "\n The range is not included in the original histogram range." << std::endl;

    return m_l1HistLimits;
  }

  if (edm::isDebugEnabled()) {
    LogDebug("L1GetHistLimits") << "\n Histogram limits for L1GtObject" << l1GtObject << " and quantity " << quantity
                                << "\n  Number of bins:           " << m_l1HistLimits.nrBins
                                << "\n  Lower limit of first bin: " << m_l1HistLimits.lowerBinValue
                                << "\n  Upper limit of last bin:  " << m_l1HistLimits.upperBinValue << std::endl;

    int binThreshSize = static_cast<int>(m_l1HistLimits.binThresholds.size());

    if (binThreshSize != (m_l1HistLimits.nrBins + 1)) {
      LogTrace("L1GetHistLimits") << "\n Warning: inconsistent nrBins and binThresholds size"
                                  << "\n   Number of bins nrBins = " << m_l1HistLimits.nrBins
                                  << "\n   binThresholds size =    " << binThreshSize
                                  << "\n Please fix the L1GetLimits class.\n\n"
                                  << std::endl;
    }

    for (int iBin = 0; iBin < binThreshSize; ++iBin) {
      LogTrace("L1GetHistLimits") << " Bin " << std::right << std::setw(5) << iBin << ":  "
                                  << m_l1HistLimits.binThresholds[iBin] << std::endl;
    }
  }

  return m_l1HistLimits;
}

const int L1GetHistLimits::l1HistNrBins(const L1GtObject& l1GtObject, const std::string& quantity) {
  getHistLimits(l1GtObject, quantity);
  return m_l1HistLimits.nrBins;
}

const double L1GetHistLimits::l1HistLowerBinValue(const L1GtObject& l1GtObject, const std::string& quantity) {
  getHistLimits(l1GtObject, quantity);
  return m_l1HistLimits.lowerBinValue;
}

const double L1GetHistLimits::l1HistUpperBinValue(const L1GtObject& l1GtObject, const std::string& quantity) {
  getHistLimits(l1GtObject, quantity);
  return m_l1HistLimits.upperBinValue;
}

const std::vector<float>& L1GetHistLimits::l1HistBinThresholds(const L1GtObject& l1GtObject,
                                                               const std::string& quantity) {
  getHistLimits(l1GtObject, quantity);
  return m_l1HistLimits.binThresholds;
}