HcalNoiseAlgo.cc

Go to the documentation of this file.

#include "RecoMET/METAlgorithms/interface/HcalNoiseAlgo.h"

CommonHcalNoiseRBXData::CommonHcalNoiseRBXData(const reco::HcalNoiseRBX& rbx,
                                               double minRecHitE,
                                               double minLowHitE,
                                               double minHighHitE,
                                               double TS4TS5EnergyThreshold,
                                               std::vector<std::pair<double, double> > const& TS4TS5UpperCut,
                                               std::vector<std::pair<double, double> > const& TS4TS5LowerCut,
                                               double minRBXRechitR45E)
    : r45Count_(0), r45Fraction_(0), r45EnergyFraction_(0) {
  // energy
  energy_ = rbx.recHitEnergy(minRecHitE);

  // ratio
  e2ts_ = rbx.allChargeHighest2TS();
  e10ts_ = rbx.allChargeTotal();

  // TS4TS5
  TS4TS5Decision_ = true;
  if (energy_ > TS4TS5EnergyThreshold)  // check filter
  {
    std::vector<float> AllCharge = rbx.allCharge();
    double BaseCharge = AllCharge[4] + AllCharge[5];
    if (BaseCharge < 1)
      BaseCharge = 1;
    double TS4TS5 = (AllCharge[4] - AllCharge[5]) / BaseCharge;

    if (CheckPassFilter(BaseCharge, TS4TS5, TS4TS5UpperCut, 1) == false)
      TS4TS5Decision_ = false;
    if (CheckPassFilter(BaseCharge, TS4TS5, TS4TS5LowerCut, -1) == false)
      TS4TS5Decision_ = false;
  } else
    TS4TS5Decision_ = true;

  // Rechit-wide R45
  int rbxHitCount = rbx.numRecHits(minRBXRechitR45E);
  r45Count_ = 0;
  r45Fraction_ = 0;
  r45EnergyFraction_ = 0;
  if (rbxHitCount > 0) {
    r45Count_ = rbx.numRecHitsFailR45(minRBXRechitR45E);
    r45Fraction_ = (double)(r45Count_) / (double)(rbxHitCount);
    r45EnergyFraction_ = rbx.recHitEnergyFailR45(minRBXRechitR45E) / rbx.recHitEnergy(minRBXRechitR45E);
  }

  // # of hits
  numHPDHits_ = 0;
  for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
    int nhpdhits = it1->numRecHits(minRecHitE);
    if (numHPDHits_ < nhpdhits)
      numHPDHits_ = nhpdhits;
  }
  numRBXHits_ = rbx.numRecHits(minRecHitE);
  numHPDNoOtherHits_ = (numHPDHits_ == numRBXHits_) ? numHPDHits_ : 0;

  // # of ADC zeros
  numZeros_ = rbx.totalZeros();

  // timing
  minLowEHitTime_ = minHighEHitTime_ = 99999.;
  maxLowEHitTime_ = maxHighEHitTime_ = -99999.;
  lowEHitTimeSqrd_ = highEHitTimeSqrd_ = 0;
  numLowEHits_ = numHighEHits_ = 0;
  for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
    edm::RefVector<HBHERecHitCollection> rechits = it1->recHits();
    for (edm::RefVector<HBHERecHitCollection>::const_iterator it2 = rechits.begin(); it2 != rechits.end(); ++it2) {
      float energy = (*it2)->energy();
      float time = (*it2)->time();
      if (energy >= minLowHitE) {
        if (minLowEHitTime_ > time)
          minLowEHitTime_ = time;
        if (maxLowEHitTime_ < time)
          maxLowEHitTime_ = time;
        lowEHitTimeSqrd_ += time * time;
        ++numLowEHits_;
      }
      if (energy >= minHighHitE) {
        if (minHighEHitTime_ > time)
          minHighEHitTime_ = time;
        if (maxHighEHitTime_ < time)
          maxHighEHitTime_ = time;
        highEHitTimeSqrd_ += time * time;
        ++numHighEHits_;
      }
    }
  }

  // emf (get the one with minimum value)
  HPDEMF_ = 999.;
  for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
    double eme = it1->caloTowerEmE();
    double hade = it1->recHitEnergy(minRecHitE);
    double emf = (eme + hade) == 0 ? 999 : eme / (eme + hade);
    if (HPDEMF_ > emf)
      HPDEMF_ = emf;
  }
  double eme = rbx.caloTowerEmE();
  RBXEMF_ = (eme + energy_) == 0 ? 999 : eme / (eme + energy_);

  // calotowers
  rbxtowers_.clear();
  JoinCaloTowerRefVectorsWithoutDuplicates join;
  for (std::vector<reco::HcalNoiseHPD>::const_iterator it1 = rbx.HPDsBegin(); it1 != rbx.HPDsEnd(); ++it1) {
    join(rbxtowers_, it1->caloTowers());
  }

  return;
}

HcalNoiseAlgo::HcalNoiseAlgo(const edm::ParameterSet& iConfig) {
  pMinERatio_ = iConfig.getParameter<double>("pMinERatio");
  pMinEZeros_ = iConfig.getParameter<double>("pMinEZeros");
  pMinEEMF_ = iConfig.getParameter<double>("pMinEEMF");

  minERatio_ = iConfig.getParameter<double>("minERatio");
  minEZeros_ = iConfig.getParameter<double>("minEZeros");
  minEEMF_ = iConfig.getParameter<double>("minEEMF");

  pMinE_ = iConfig.getParameter<double>("pMinE");
  pMinRatio_ = iConfig.getParameter<double>("pMinRatio");
  pMaxRatio_ = iConfig.getParameter<double>("pMaxRatio");
  pMinHPDHits_ = iConfig.getParameter<int>("pMinHPDHits");
  pMinRBXHits_ = iConfig.getParameter<int>("pMinRBXHits");
  pMinHPDNoOtherHits_ = iConfig.getParameter<int>("pMinHPDNoOtherHits");
  pMinZeros_ = iConfig.getParameter<int>("pMinZeros");
  pMinLowEHitTime_ = iConfig.getParameter<double>("pMinLowEHitTime");
  pMaxLowEHitTime_ = iConfig.getParameter<double>("pMaxLowEHitTime");
  pMinHighEHitTime_ = iConfig.getParameter<double>("pMinHighEHitTime");
  pMaxHighEHitTime_ = iConfig.getParameter<double>("pMaxHighEHitTime");
  pMaxHPDEMF_ = iConfig.getParameter<double>("pMaxHPDEMF");
  pMaxRBXEMF_ = iConfig.getParameter<double>("pMaxRBXEMF");

  if (iConfig.existsAs<int>("pMinRBXRechitR45Count"))
    pMinRBXRechitR45Count_ = iConfig.getParameter<int>("pMinRBXRechitR45Count");
  else
    pMinRBXRechitR45Count_ = 0;
  if (iConfig.existsAs<double>("pMinRBXRechitR45Fraction"))
    pMinRBXRechitR45Fraction_ = iConfig.getParameter<double>("pMinRBXRechitR45Fraction");
  else
    pMinRBXRechitR45Fraction_ = 0;
  if (iConfig.existsAs<double>("pMinRBXRechitR45EnergyFraction"))
    pMinRBXRechitR45EnergyFraction_ = iConfig.getParameter<double>("pMinRBXRechitR45EnergyFraction");
  else
    pMinRBXRechitR45EnergyFraction_ = 0;

  lMinRatio_ = iConfig.getParameter<double>("lMinRatio");
  lMaxRatio_ = iConfig.getParameter<double>("lMaxRatio");
  lMinHPDHits_ = iConfig.getParameter<int>("lMinHPDHits");
  lMinRBXHits_ = iConfig.getParameter<int>("lMinRBXHits");
  lMinHPDNoOtherHits_ = iConfig.getParameter<int>("lMinHPDNoOtherHits");
  lMinZeros_ = iConfig.getParameter<int>("lMinZeros");
  lMinLowEHitTime_ = iConfig.getParameter<double>("lMinLowEHitTime");
  lMaxLowEHitTime_ = iConfig.getParameter<double>("lMaxLowEHitTime");
  lMinHighEHitTime_ = iConfig.getParameter<double>("lMinHighEHitTime");
  lMaxHighEHitTime_ = iConfig.getParameter<double>("lMaxHighEHitTime");

  if (iConfig.existsAs<std::vector<double> >("lRBXRecHitR45Cuts"))
    lMinRBXRechitR45Cuts_ = iConfig.getParameter<std::vector<double> >("lRBXRecHitR45Cuts");
  else {
    double defaultCut[4] = {-999, -999, -999, -999};
    lMinRBXRechitR45Cuts_ = std::vector<double>(defaultCut, defaultCut + 4);
  }

  tMinRatio_ = iConfig.getParameter<double>("tMinRatio");
  tMaxRatio_ = iConfig.getParameter<double>("tMaxRatio");
  tMinHPDHits_ = iConfig.getParameter<int>("tMinHPDHits");
  tMinRBXHits_ = iConfig.getParameter<int>("tMinRBXHits");
  tMinHPDNoOtherHits_ = iConfig.getParameter<int>("tMinHPDNoOtherHits");
  tMinZeros_ = iConfig.getParameter<int>("tMinZeros");
  tMinLowEHitTime_ = iConfig.getParameter<double>("tMinLowEHitTime");
  tMaxLowEHitTime_ = iConfig.getParameter<double>("tMaxLowEHitTime");
  tMinHighEHitTime_ = iConfig.getParameter<double>("tMinHighEHitTime");
  tMaxHighEHitTime_ = iConfig.getParameter<double>("tMaxHighEHitTime");

  if (iConfig.existsAs<std::vector<double> >("tRBXRecHitR45Cuts"))
    tMinRBXRechitR45Cuts_ = iConfig.getParameter<std::vector<double> >("tRBXRecHitR45Cuts");
  else {
    double defaultCut[4] = {-999, -999, -999, -999};
    tMinRBXRechitR45Cuts_ = std::vector<double>(defaultCut, defaultCut + 4);
  }

  hlMaxHPDEMF_ = iConfig.getParameter<double>("hlMaxHPDEMF");
  hlMaxRBXEMF_ = iConfig.getParameter<double>("hlMaxRBXEMF");
}

bool HcalNoiseAlgo::isProblematic(const CommonHcalNoiseRBXData& data) const {
  if (data.energy() > pMinE_)
    return true;
  if (data.validRatio() && data.energy() > pMinERatio_ && data.ratio() < pMinRatio_)
    return true;
  if (data.validRatio() && data.energy() > pMinERatio_ && data.ratio() > pMaxRatio_)
    return true;
  if (data.numHPDHits() >= pMinHPDHits_)
    return true;
  if (data.numRBXHits() >= pMinRBXHits_)
    return true;
  if (data.numHPDNoOtherHits() >= pMinHPDNoOtherHits_)
    return true;
  if (data.numZeros() >= pMinZeros_ && data.energy() > pMinEZeros_)
    return true;
  if (data.minLowEHitTime() < pMinLowEHitTime_)
    return true;
  if (data.maxLowEHitTime() > pMaxLowEHitTime_)
    return true;
  if (data.minHighEHitTime() < pMinHighEHitTime_)
    return true;
  if (data.maxHighEHitTime() > pMaxHighEHitTime_)
    return true;
  if (data.HPDEMF() < pMaxHPDEMF_ && data.energy() > pMinEEMF_)
    return true;
  if (data.RBXEMF() < pMaxRBXEMF_ && data.energy() > pMinEEMF_)
    return true;
  if (data.r45Count() >= pMinRBXRechitR45Count_)
    return true;
  if (data.r45Fraction() >= pMinRBXRechitR45Fraction_)
    return true;
  if (data.r45EnergyFraction() >= pMinRBXRechitR45EnergyFraction_)
    return true;

  return false;
}

bool HcalNoiseAlgo::passLooseNoiseFilter(const CommonHcalNoiseRBXData& data) const {
  return (passLooseRatio(data) && passLooseHits(data) && passLooseZeros(data) && passLooseTiming(data) &&
          passLooseRBXRechitR45(data));
}

bool HcalNoiseAlgo::passTightNoiseFilter(const CommonHcalNoiseRBXData& data) const {
  return (passTightRatio(data) && passTightHits(data) && passTightZeros(data) && passTightTiming(data) &&
          passTightRBXRechitR45(data));
}

bool HcalNoiseAlgo::passHighLevelNoiseFilter(const CommonHcalNoiseRBXData& data) const {
  if (passEMFThreshold(data)) {
    if (data.HPDEMF() < hlMaxHPDEMF_)
      return false;
    if (data.RBXEMF() < hlMaxRBXEMF_)
      return false;
  }
  return true;
}

bool HcalNoiseAlgo::passLooseRatio(const CommonHcalNoiseRBXData& data) const {
  if (passRatioThreshold(data)) {
    if (data.validRatio() && data.ratio() < lMinRatio_)
      return false;
    if (data.validRatio() && data.ratio() > lMaxRatio_)
      return false;
  }
  return true;
}

bool HcalNoiseAlgo::passLooseHits(const CommonHcalNoiseRBXData& data) const {
  if (data.numHPDHits() >= lMinHPDHits_)
    return false;
  if (data.numRBXHits() >= lMinRBXHits_)
    return false;
  if (data.numHPDNoOtherHits() >= lMinHPDNoOtherHits_)
    return false;
  return true;
}

bool HcalNoiseAlgo::passLooseZeros(const CommonHcalNoiseRBXData& data) const {
  if (passZerosThreshold(data)) {
    if (data.numZeros() >= lMinZeros_)
      return false;
  }
  return true;
}

bool HcalNoiseAlgo::passLooseTiming(const CommonHcalNoiseRBXData& data) const {
  if (data.minLowEHitTime() < lMinLowEHitTime_)
    return false;
  if (data.maxLowEHitTime() > lMaxLowEHitTime_)
    return false;
  if (data.minHighEHitTime() < lMinHighEHitTime_)
    return false;
  if (data.maxHighEHitTime() > lMaxHighEHitTime_)
    return false;
  return true;
}

bool HcalNoiseAlgo::passLooseRBXRechitR45(const CommonHcalNoiseRBXData& data) const {
  int Count = data.r45Count();
  double Fraction = data.r45Fraction();
  double EnergyFraction = data.r45EnergyFraction();

  for (int i = 0; i + 3 < (int)lMinRBXRechitR45Cuts_.size(); i = i + 4) {
    double Value = Count * lMinRBXRechitR45Cuts_[i] + Fraction * lMinRBXRechitR45Cuts_[i + 1] +
                   EnergyFraction * lMinRBXRechitR45Cuts_[i + 2] + lMinRBXRechitR45Cuts_[i + 3];
    if (Value > 0)
      return false;
  }

  return true;
}

bool HcalNoiseAlgo::passTightRatio(const CommonHcalNoiseRBXData& data) const {
  if (passRatioThreshold(data)) {
    if (data.validRatio() && data.ratio() < tMinRatio_)
      return false;
    if (data.validRatio() && data.ratio() > tMaxRatio_)
      return false;
  }
  return true;
}

bool HcalNoiseAlgo::passTightHits(const CommonHcalNoiseRBXData& data) const {
  if (data.numHPDHits() >= tMinHPDHits_)
    return false;
  if (data.numRBXHits() >= tMinRBXHits_)
    return false;
  if (data.numHPDNoOtherHits() >= tMinHPDNoOtherHits_)
    return false;
  return true;
}

bool HcalNoiseAlgo::passTightZeros(const CommonHcalNoiseRBXData& data) const {
  if (passZerosThreshold(data)) {
    if (data.numZeros() >= tMinZeros_)
      return false;
  }
  return true;
}

bool HcalNoiseAlgo::passTightTiming(const CommonHcalNoiseRBXData& data) const {
  if (data.minLowEHitTime() < tMinLowEHitTime_)
    return false;
  if (data.maxLowEHitTime() > tMaxLowEHitTime_)
    return false;
  if (data.minHighEHitTime() < tMinHighEHitTime_)
    return false;
  if (data.maxHighEHitTime() > tMaxHighEHitTime_)
    return false;
  return true;
}

bool HcalNoiseAlgo::passTightRBXRechitR45(const CommonHcalNoiseRBXData& data) const {
  int Count = data.r45Count();
  double Fraction = data.r45Fraction();
  double EnergyFraction = data.r45EnergyFraction();

  for (int i = 0; i + 3 < (int)tMinRBXRechitR45Cuts_.size(); i = i + 4) {
    double Value = Count * tMinRBXRechitR45Cuts_[i] + Fraction * tMinRBXRechitR45Cuts_[i + 1] +
                   EnergyFraction * tMinRBXRechitR45Cuts_[i + 2] + tMinRBXRechitR45Cuts_[i + 3];
    if (Value > 0)
      return false;
  }

  return true;
}

bool HcalNoiseAlgo::passRatioThreshold(const CommonHcalNoiseRBXData& data) const {
  return (data.energy() > minERatio_);
}

bool HcalNoiseAlgo::passZerosThreshold(const CommonHcalNoiseRBXData& data) const {
  return (data.energy() > minEZeros_);
}

bool HcalNoiseAlgo::passEMFThreshold(const CommonHcalNoiseRBXData& data) const { return (data.energy() > minEEMF_); }

void JoinCaloTowerRefVectorsWithoutDuplicates::operator()(edm::RefVector<CaloTowerCollection>& v1,
                                                          const edm::RefVector<CaloTowerCollection>& v2) const {
  // combines them first into a set to get rid of duplicates and then puts them into the first vector

  // sorts them first to get rid of duplicates, then puts them into another RefVector
  twrrefset_t twrrefset;
  for (edm::RefVector<CaloTowerCollection>::const_iterator it = v1.begin(); it != v1.end(); ++it)
    twrrefset.insert(*it);
  for (edm::RefVector<CaloTowerCollection>::const_iterator it = v2.begin(); it != v2.end(); ++it)
    twrrefset.insert(*it);

  // clear the original refvector and put them back in
  v1.clear();
  for (twrrefset_t::const_iterator it = twrrefset.begin(); it != twrrefset.end(); ++it) {
    v1.push_back(*it);
  }
  return;
}

bool CommonHcalNoiseRBXData::CheckPassFilter(double Charge,
                                             double Discriminant,
                                             std::vector<std::pair<double, double> > const& Cuts,
                                             int Side) {
  //
  // Checks whether Discriminant value passes Cuts for the specified Charge.  True if pulse is good.
  //
  // The "Cuts" pairs are assumed to be sorted in terms of size from small to large,
  //    where each "pair" = (Charge, Discriminant)
  // "Side" is either positive or negative, which determines whether to discard the pulse if discriminant
  //    is greater or smaller than the cut value
  //

  if (Cuts.empty())  // safety check that there are some cuts defined
    return true;

  if (Charge <= Cuts[0].first)  // too small to cut on
    return true;

  int IndexLargerThanCharge = -1;  // find the range it is falling in
  for (int i = 1; i < (int)Cuts.size(); i++) {
    if (Cuts[i].first > Charge) {
      IndexLargerThanCharge = i;
      break;
    }
  }

  double limit = 1000000;

  if (IndexLargerThanCharge == -1)  // if charge is greater than the last entry, assume flat line
    limit = Cuts[Cuts.size() - 1].second;
  else  // otherwise, do a linear interpolation to find the cut position
  {
    double C1 = Cuts[IndexLargerThanCharge].first;
    double C2 = Cuts[IndexLargerThanCharge - 1].first;
    double L1 = Cuts[IndexLargerThanCharge].second;
    double L2 = Cuts[IndexLargerThanCharge - 1].second;

    limit = (Charge - C1) / (C2 - C1) * (L2 - L1) + L1;
  }

  if (Side > 0 && Discriminant > limit)
    return false;
  if (Side < 0 && Discriminant < limit)
    return false;

  return true;
}