EgammaCoreTools/src/EcalTools.cc

0001
0002 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0003 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0004 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 float EcalTools::swissCross(const DetId& id,
0010                             const EcalRecHitCollection& recHits,
0011                             float recHitThreshold,
0012                             bool avoidIeta85) {
0013   // compute swissCross
0014   if (id.subdetId() == EcalBarrel) {
0015     EBDetId ebId(id);
0016     // avoid recHits at |eta|=85 where one side of the neighbours is missing
0017     // (may improve considering also eta module borders, but no
0018     // evidence for the time being that there the performance is
0019     // different)
0020     if (abs(ebId.ieta()) == 85 && avoidIeta85)
0021       return 0;
0022     // select recHits with Et above recHitThreshold
0023     if (recHitApproxEt(id, recHits) < recHitThreshold)
0024       return 0;
0025     float s4 = 0;
0026     float e1 = recHitE(id, recHits);
0027     // protect against nan (if 0 threshold is given above)
0028     if (e1 == 0)
0029       return 0;
0030     s4 += recHitE(id, recHits, 1, 0);
0031     s4 += recHitE(id, recHits, -1, 0);
0032     s4 += recHitE(id, recHits, 0, 1);
0033     s4 += recHitE(id, recHits, 0, -1);
0034     return 1 - s4 / e1;
0035   } else if (id.subdetId() == EcalEndcap) {
0036     EEDetId eeId(id);
0037     // select recHits with E above recHitThreshold
0038     float e1 = recHitE(id, recHits);
0039     if (e1 < recHitThreshold)
0040       return 0;
0041     float s4 = 0;
0042     // protect against nan (if 0 threshold is given above)
0043     if (e1 == 0)
0044       return 0;
0045     s4 += recHitE(id, recHits, 1, 0);
0046     s4 += recHitE(id, recHits, -1, 0);
0047     s4 += recHitE(id, recHits, 0, 1);
0048     s4 += recHitE(id, recHits, 0, -1);
0049     return 1 - s4 / e1;
0050   }
0051   return 0;
0052 }
0053
0054 bool EcalTools::isNextToDead(const DetId& id, const EcalNextToDeadChannel& dch) {
0055   EcalNextToDeadChannel::const_iterator chIt = dch.find(id);
0056
0057   if (chIt != dch.end()) {
0058     return *chIt;
0059   } else {
0060     edm::LogError("EcalDBError") << "No NextToDead  status found for xtal " << id.rawId();
0061   }
0062
0063   return false;
0064 }
0065
0066 bool EcalTools::isNextToDeadFromNeighbours(const DetId& id, const EcalChannelStatus& chs, int chStatusThreshold) {
0067   if (deadNeighbour(id, chs, chStatusThreshold, 1, 0))
0068     return true;
0069   if (deadNeighbour(id, chs, chStatusThreshold, -1, 0))
0070     return true;
0071   if (deadNeighbour(id, chs, chStatusThreshold, 0, 1))
0072     return true;
0073   if (deadNeighbour(id, chs, chStatusThreshold, 0, -1))
0074     return true;
0075   if (deadNeighbour(id, chs, chStatusThreshold, 1, -1))
0076     return true;
0077   if (deadNeighbour(id, chs, chStatusThreshold, 1, 1))
0078     return true;
0079   if (deadNeighbour(id, chs, chStatusThreshold, -1, 1))
0080     return true;
0081   if (deadNeighbour(id, chs, chStatusThreshold, -1, -1))
0082     return true;
0083
0084   return false;
0085 }
0086
0087 bool EcalTools::deadNeighbour(const DetId& id, const EcalChannelStatus& chs, int chStatusThreshold, int dx, int dy) {
0088   DetId nid;
0089   if (id.subdetId() == EcalBarrel)
0090     nid = EBDetId::offsetBy(id, dx, dy);
0091   else if (id.subdetId() == EcalEndcap)
0092     nid = EEDetId::offsetBy(id, dx, dy);
0093
0094   if (!nid)
0095     return false;
0096
0097   EcalChannelStatus::const_iterator chIt = chs.find(nid);
0098   uint16_t dbStatus = 0;
0099   if (chIt != chs.end()) {
0100     // note that
0101     dbStatus = chIt->getStatusCode();
0102   } else {
0103     edm::LogError("EcalDBError") << "No channel status found for xtal " << id.rawId()
0104                                  << "! something wrong with EcalChannelStatus in your DB? ";
0105   }
0106
0107   return (dbStatus >= chStatusThreshold);
0108 }
0109
0110 float EcalTools::recHitE(const DetId id, const EcalRecHitCollection& recHits, int di, int dj) {
0111   // in the barrel:   di = dEta   dj = dPhi
0112   // in the endcap:   di = dX     dj = dY
0113
0114   DetId nid;
0115   if (id.subdetId() == EcalBarrel)
0116     nid = EBDetId::offsetBy(id, di, dj);
0117   else if (id.subdetId() == EcalEndcap)
0118     nid = EEDetId::offsetBy(id, di, dj);
0119
0120   return (nid == DetId(0) ? 0 : recHitE(nid, recHits));
0121 }
0122
0123 float EcalTools::recHitE(const DetId id, const EcalRecHitCollection& recHits) {
0124   if (id == DetId(0)) {
0125     return 0;
0126   } else {
0127     EcalRecHitCollection::const_iterator it = recHits.find(id);
0128     if (it != recHits.end())
0129       return (*it).energy();
0130   }
0131   return 0;
0132 }
0133
0134 float EcalTools::recHitApproxEt(const DetId id, const EcalRecHitCollection& recHits) {
0135   // for the time being works only for the barrel
0136   if (id.subdetId() == EcalBarrel) {
0137     return recHitE(id, recHits) / cosh(EBDetId::approxEta(id));
0138   }
0139   return 0;
0140 }
0141
0142 bool isNextToBoundary(const DetId& id) {
0143   if (id.subdetId() == EcalBarrel)
0144     return EBDetId::isNextToBoundary(EBDetId(id));
0145   else if (id.subdetId() == EcalEndcap)
0146     return EEDetId::isNextToBoundary(EEDetId(id));
0147
0148   return false;
0149 }