EcalClusterSeverityLevelAlgo.cc

Go to the documentation of this file.

#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterSeverityLevelAlgo.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"

float EcalClusterSeverityLevelAlgo::goodFraction(const reco::CaloCluster& cluster,
                                                 const EcalRecHitCollection& recHits,
                                                 const EcalSeverityLevelAlgo& sevlv) {
  float fraction = 0.;
  std::vector<std::pair<DetId, float> > hitsAndFracs = cluster.hitsAndFractions();
  std::vector<std::pair<DetId, float> >::const_iterator it;
  for (it = hitsAndFracs.begin(); it != hitsAndFracs.end(); ++it) {
    DetId id = (*it).first;
    EcalRecHitCollection::const_iterator jrh = recHits.find(id);
    if (jrh == recHits.end()) {
      edm::LogError("EcalClusterSeverityLevelAlgo")
          << "The cluster DetId " << id.rawId() << " is not in the recHit collection!!";
      return -1;
    }

    uint32_t sev = sevlv.severityLevel(id, recHits);
    //                if ( sev == EcalSeverityLevelAlgo::kBad ) ++recoveryFailed;
    if (sev == EcalSeverityLevel::kProblematic || sev == EcalSeverityLevel::kRecovered ||
        sev == EcalSeverityLevel::kBad) {
      //            std::cout << "[goodFraction] Found a problematic channel " << EBDetId(id) << " " << flag << " energy: " <<  (*jrh).energy() << std::endl;
      fraction += (*jrh).energy() * (*it).second / cluster.energy();
    }
  }
  return 1. - fraction;
}

float EcalClusterSeverityLevelAlgo::fractionAroundClosestProblematic(const reco::CaloCluster& cluster,
                                                                     const EcalRecHitCollection& recHits,
                                                                     const CaloTopology* topology,
                                                                     const EcalSeverityLevelAlgo& sevlv) {
  DetId closestProb = closestProblematic(cluster, recHits, topology, sevlv);
  //  std::cout << "%%%%%%%%%%% Closest prob is " << EBDetId(closestProb) << std::endl;
  if (closestProb.null())
    return 0.;

  std::vector<DetId> neighbours = topology->getWindow(closestProb, 3, 3);
  std::vector<DetId>::const_iterator itn;

  std::vector<std::pair<DetId, float> > hitsAndFracs = cluster.hitsAndFractions();
  std::vector<std::pair<DetId, float> >::const_iterator it;

  float fraction = 0.;

  for (itn = neighbours.begin(); itn != neighbours.end(); ++itn) {
    //      std::cout << "Checking detId " << EBDetId((*itn)) << std::endl;
    for (it = hitsAndFracs.begin(); it != hitsAndFracs.end(); ++it) {
      DetId id = (*it).first;
      if (id != (*itn))
        continue;
      //      std::cout << "Is in cluster detId " << EBDetId(id) << std::endl;
      EcalRecHitCollection::const_iterator jrh = recHits.find(id);
      if (jrh == recHits.end()) {
        edm::LogError("EcalClusterSeverityLevelAlgo")
            << "The cluster DetId " << id.rawId() << " is not in the recHit collection!!";
        return -1;
      }

      fraction += (*jrh).energy() * (*it).second / cluster.energy();
    }
  }
  //  std::cout << "%%%%%%%%%%% Fraction is " << fraction << std::endl;
  return fraction;
}

DetId EcalClusterSeverityLevelAlgo::closestProblematic(const reco::CaloCluster& cluster,
                                                       const EcalRecHitCollection& recHits,
                                                       const CaloTopology* topology,
                                                       const EcalSeverityLevelAlgo& sevlv) {
  DetId seed = EcalClusterTools::getMaximum(cluster, &recHits).first;
  if ((seed.det() != DetId::Ecal) || (EcalSubdetector(seed.subdetId()) != EcalBarrel)) {
    //method not supported if not in Barrel
    edm::LogError("EcalClusterSeverityLevelAlgo") << "The cluster seed is not in the BARREL";
    return DetId(0);
  }

  int minDist = 9999;
  DetId closestProb(0);
  //Get a window of DetId around the seed crystal
  std::vector<DetId> neighbours = topology->getWindow(seed, 51, 11);

  for (std::vector<DetId>::const_iterator it = neighbours.begin(); it != neighbours.end(); ++it) {
    EcalRecHitCollection::const_iterator jrh = recHits.find(*it);
    if (jrh == recHits.end())
      continue;
    //Now checking rh flag
    uint32_t sev = sevlv.severityLevel(*it, recHits);
    if (sev == EcalSeverityLevel::kGood)
      continue;
    //      std::cout << "[closestProblematic] Found a problematic channel " << EBDetId(*it) << " " << flag << std::endl;
    //Find the closest DetId in eta,phi space (distance defined by deta^2 + dphi^2)
    int deta = EBDetId::distanceEta(EBDetId(seed), EBDetId(*it));
    int dphi = EBDetId::distancePhi(EBDetId(seed), EBDetId(*it));
    double r = sqrt(deta * deta + dphi * dphi);
    if (r < minDist) {
      closestProb = *it;
      minDist = r;
    }
  }

  return closestProb;
}

std::pair<int, int> EcalClusterSeverityLevelAlgo::etaphiDistanceClosestProblematic(const reco::CaloCluster& cluster,
                                                                                   const EcalRecHitCollection& recHits,
                                                                                   const CaloTopology* topology,
                                                                                   const EcalSeverityLevelAlgo& sevlv) {
  DetId seed = EcalClusterTools::getMaximum(cluster, &recHits).first;
  if ((seed.det() != DetId::Ecal) || (EcalSubdetector(seed.subdetId()) != EcalBarrel)) {
    edm::LogError("EcalClusterSeverityLevelAlgo") << "The cluster seed is not in the BARREL";
    //method not supported if not in Barrel
    return std::pair<int, int>(-1, -1);
  }

  DetId closestProb = closestProblematic(cluster, recHits, topology, sevlv);

  if (!closestProb.null())
    return std::pair<int, int>(EBDetId::distanceEta(EBDetId(seed), EBDetId(closestProb)),
                               EBDetId::distancePhi(EBDetId(seed), EBDetId(closestProb)));
  else
    return std::pair<int, int>(-1, -1);
}