EgammaCoreTools/src/EcalClusterSeverityLevelAlgo.cc

0001 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0002 #include "FWCore/ServiceRegistry/interface/Service.h"
0003 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterSeverityLevelAlgo.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0007 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0008 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0009
0010 float EcalClusterSeverityLevelAlgo::goodFraction(const reco::CaloCluster& cluster,
0011                                                  const EcalRecHitCollection& recHits,
0012                                                  const EcalSeverityLevelAlgo& sevlv) {
0013   float fraction = 0.;
0014   std::vector<std::pair<DetId, float> > hitsAndFracs = cluster.hitsAndFractions();
0015   std::vector<std::pair<DetId, float> >::const_iterator it;
0016   for (it = hitsAndFracs.begin(); it != hitsAndFracs.end(); ++it) {
0017     DetId id = (*it).first;
0018     EcalRecHitCollection::const_iterator jrh = recHits.find(id);
0019     if (jrh == recHits.end()) {
0020       edm::LogError("EcalClusterSeverityLevelAlgo")
0021           << "The cluster DetId " << id.rawId() << " is not in the recHit collection!!";
0022       return -1;
0023     }
0024
0025     uint32_t sev = sevlv.severityLevel(id, recHits);
0026     //                if ( sev == EcalSeverityLevelAlgo::kBad ) ++recoveryFailed;
0027     if (sev == EcalSeverityLevel::kProblematic || sev == EcalSeverityLevel::kRecovered ||
0028         sev == EcalSeverityLevel::kBad) {
0029       //            std::cout << "[goodFraction] Found a problematic channel " << EBDetId(id) << " " << flag << " energy: " <<  (*jrh).energy() << std::endl;
0030       fraction += (*jrh).energy() * (*it).second / cluster.energy();
0031     }
0032   }
0033   return 1. - fraction;
0034 }
0035
0036 float EcalClusterSeverityLevelAlgo::fractionAroundClosestProblematic(const reco::CaloCluster& cluster,
0037                                                                      const EcalRecHitCollection& recHits,
0038                                                                      const CaloTopology* topology,
0039                                                                      const EcalSeverityLevelAlgo& sevlv) {
0040   DetId closestProb = closestProblematic(cluster, recHits, topology, sevlv);
0041   //  std::cout << "%%%%%%%%%%% Closest prob is " << EBDetId(closestProb) << std::endl;
0042   if (closestProb.null())
0043     return 0.;
0044
0045   std::vector<DetId> neighbours = topology->getWindow(closestProb, 3, 3);
0046   std::vector<DetId>::const_iterator itn;
0047
0048   std::vector<std::pair<DetId, float> > hitsAndFracs = cluster.hitsAndFractions();
0049   std::vector<std::pair<DetId, float> >::const_iterator it;
0050
0051   float fraction = 0.;
0052
0053   for (itn = neighbours.begin(); itn != neighbours.end(); ++itn) {
0054     //      std::cout << "Checking detId " << EBDetId((*itn)) << std::endl;
0055     for (it = hitsAndFracs.begin(); it != hitsAndFracs.end(); ++it) {
0056       DetId id = (*it).first;
0057       if (id != (*itn))
0058         continue;
0059       //      std::cout << "Is in cluster detId " << EBDetId(id) << std::endl;
0060       EcalRecHitCollection::const_iterator jrh = recHits.find(id);
0061       if (jrh == recHits.end()) {
0062         edm::LogError("EcalClusterSeverityLevelAlgo")
0063             << "The cluster DetId " << id.rawId() << " is not in the recHit collection!!";
0064         return -1;
0065       }
0066
0067       fraction += (*jrh).energy() * (*it).second / cluster.energy();
0068     }
0069   }
0070   //  std::cout << "%%%%%%%%%%% Fraction is " << fraction << std::endl;
0071   return fraction;
0072 }
0073
0074 DetId EcalClusterSeverityLevelAlgo::closestProblematic(const reco::CaloCluster& cluster,
0075                                                        const EcalRecHitCollection& recHits,
0076                                                        const CaloTopology* topology,
0077                                                        const EcalSeverityLevelAlgo& sevlv) {
0078   DetId seed = EcalClusterTools::getMaximum(cluster, &recHits).first;
0079   if ((seed.det() != DetId::Ecal) || (EcalSubdetector(seed.subdetId()) != EcalBarrel)) {
0080     //method not supported if not in Barrel
0081     edm::LogError("EcalClusterSeverityLevelAlgo") << "The cluster seed is not in the BARREL";
0082     return DetId(0);
0083   }
0084
0085   int minDist = 9999;
0086   DetId closestProb(0);
0087   //Get a window of DetId around the seed crystal
0088   std::vector<DetId> neighbours = topology->getWindow(seed, 51, 11);
0089
0090   for (std::vector<DetId>::const_iterator it = neighbours.begin(); it != neighbours.end(); ++it) {
0091     EcalRecHitCollection::const_iterator jrh = recHits.find(*it);
0092     if (jrh == recHits.end())
0093       continue;
0094     //Now checking rh flag
0095     uint32_t sev = sevlv.severityLevel(*it, recHits);
0096     if (sev == EcalSeverityLevel::kGood)
0097       continue;
0098     //      std::cout << "[closestProblematic] Found a problematic channel " << EBDetId(*it) << " " << flag << std::endl;
0099     //Find the closest DetId in eta,phi space (distance defined by deta^2 + dphi^2)
0100     int deta = EBDetId::distanceEta(EBDetId(seed), EBDetId(*it));
0101     int dphi = EBDetId::distancePhi(EBDetId(seed), EBDetId(*it));
0102     double r = sqrt(deta * deta + dphi * dphi);
0103     if (r < minDist) {
0104       closestProb = *it;
0105       minDist = r;
0106     }
0107   }
0108
0109   return closestProb;
0110 }
0111
0112 std::pair<int, int> EcalClusterSeverityLevelAlgo::etaphiDistanceClosestProblematic(const reco::CaloCluster& cluster,
0113                                                                                    const EcalRecHitCollection& recHits,
0114                                                                                    const CaloTopology* topology,
0115                                                                                    const EcalSeverityLevelAlgo& sevlv) {
0116   DetId seed = EcalClusterTools::getMaximum(cluster, &recHits).first;
0117   if ((seed.det() != DetId::Ecal) || (EcalSubdetector(seed.subdetId()) != EcalBarrel)) {
0118     edm::LogError("EcalClusterSeverityLevelAlgo") << "The cluster seed is not in the BARREL";
0119     //method not supported if not in Barrel
0120     return std::pair<int, int>(-1, -1);
0121   }
0122
0123   DetId closestProb = closestProblematic(cluster, recHits, topology, sevlv);
0124
0125   if (!closestProb.null())
0126     return std::pair<int, int>(EBDetId::distanceEta(EBDetId(seed), EBDetId(closestProb)),
0127                                EBDetId::distancePhi(EBDetId(seed), EBDetId(closestProb)));
0128   else
0129     return std::pair<int, int>(-1, -1);
0130 }