HcalRecAlgos/src/HBHEIsolatedNoiseAlgos.cc

0001 /*
0002 Description: Isolation algorithms used to identify anomalous noise in the HB/HE.
0003              These algorithms will be used to reflag HB/HE rechits as noise.
0004
0005              There are 4 objects implemented here:
0006              1) ObjectValidator
0007          2) PhysicsTowerOrganizer
0008          3) HBHEHitMap
0009          4) HBHEHitMapOrganizer
0010          See comments below for details.
0011
0012 Original Author: John Paul Chou (Brown University)
0013                  Thursday, September 2, 2010
0014 */
0015
0016 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHEIsolatedNoiseAlgos.h"
0017
0018 #include "FWCore/Utilities/interface/EDMException.h"
0019
0020 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0021 #include "DataFormats/METReco/interface/CaloMET.h"
0022 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0023 #include "DataFormats/HcalRecHit/interface/HBHERecHitAuxSetter.h"
0024 #include "DataFormats/HcalRecHit/interface/CaloRecHitAuxSetter.h"
0025 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0026 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0027 #include "DataFormats/TrackReco/interface/Track.h"
0028
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0031 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0032 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0033 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0034 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0035
0036 ////////////////////////////////////////////////////////////
0037 //
0038 // ObjectValidator
0039 //
0040 ////////////////////////////////////////////////////////////
0041
0042 ObjectValidator::ObjectValidator(const edm::ParameterSet& iConfig) {
0043   HBThreshold_ = iConfig.getParameter<double>("HBThreshold");
0044   HESThreshold_ = iConfig.getParameter<double>("HESThreshold");
0045   HEDThreshold_ = iConfig.getParameter<double>("HEDThreshold");
0046   EBThreshold_ = iConfig.getParameter<double>("EBThreshold");
0047   EEThreshold_ = iConfig.getParameter<double>("EEThreshold");
0048
0049   HcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("HcalAcceptSeverityLevel");
0050   EcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("EcalAcceptSeverityLevel");
0051   UseHcalRecoveredHits_ = iConfig.getParameter<bool>("UseHcalRecoveredHits");
0052   UseEcalRecoveredHits_ = iConfig.getParameter<bool>("UseEcalRecoveredHits");
0053   UseAllCombinedRechits_ = iConfig.getParameter<bool>("UseAllCombinedRechits");
0054
0055   MinValidTrackPt_ = iConfig.getParameter<double>("MinValidTrackPt");
0056   MinValidTrackPtBarrel_ = iConfig.getParameter<double>("MinValidTrackPtBarrel");
0057   MinValidTrackNHits_ = iConfig.getParameter<int>("MinValidTrackNHits");
0058
0059   theHcalChStatus_ = nullptr;
0060   theEcalChStatus_ = nullptr;
0061   theHcalSevLvlComputer_ = nullptr;
0062   theEcalSevLvlAlgo_ = nullptr;
0063   theEBRecHitCollection_ = nullptr;
0064   theEERecHitCollection_ = nullptr;
0065
0066   return;
0067 }
0068
0069 ObjectValidator::~ObjectValidator() {}
0070
0071 bool ObjectValidator::validHit(const HBHERecHit& hit) const {
0072   assert(theHcalSevLvlComputer_ != nullptr && theHcalChStatus_ != nullptr);
0073
0074   if (UseAllCombinedRechits_)
0075     if (CaloRecHitAuxSetter::getBit(hit.auxPhase1(), HBHERecHitAuxSetter::OFF_COMBINED))
0076       return true;
0077
0078   // require the hit to pass a certain energy threshold
0079   if (hit.id().subdet() == HcalBarrel && hit.energy() < HBThreshold_)
0080     return false;
0081   else if (hit.id().subdet() == HcalEndcap && hit.id().ietaAbs() <= 20 && hit.energy() < HESThreshold_)
0082     return false;
0083   else if (hit.id().subdet() == HcalEndcap && hit.id().ietaAbs() > 20 && hit.energy() < HEDThreshold_)
0084     return false;
0085
0086   // determine if the hit is good, bad, or recovered
0087   const DetId id = hit.detid();
0088   const uint32_t recHitFlag = hit.flags();
0089   const uint32_t dbStatusFlag = theHcalChStatus_->getValues(id)->getValue();
0090   int severityLevel = theHcalSevLvlComputer_->getSeverityLevel(id, recHitFlag, dbStatusFlag);
0091   bool isRecovered = theHcalSevLvlComputer_->recoveredRecHit(id, recHitFlag);
0092
0093   if (severityLevel == 0)
0094     return true;
0095   if (isRecovered)
0096     return UseHcalRecoveredHits_;
0097   if (severityLevel > static_cast<int>(HcalAcceptSeverityLevel_))
0098     return false;
0099   else
0100     return true;
0101 }
0102
0103 bool ObjectValidator::validHit(const EcalRecHit& hit) const {
0104   assert(theEcalSevLvlAlgo_ != nullptr && theEcalChStatus_ != nullptr);
0105
0106   // require the hit to pass a certain energy threshold
0107   const DetId id = hit.detid();
0108   if (id.subdetId() == EcalBarrel && hit.energy() < EBThreshold_)
0109     return false;
0110   else if (id.subdetId() == EcalEndcap && hit.energy() < EEThreshold_)
0111     return false;
0112
0113   // determine if the hit is good, bad, or recovered
0114   int severityLevel = 999;
0115   if (id.subdetId() == EcalBarrel && theEBRecHitCollection_ != nullptr)
0116     severityLevel = theEcalSevLvlAlgo_->severityLevel(
0117         hit);  //id, *theEBRecHitCollection_, *theEcalChStatus_, 5., EcalSeverityLevelAlgo::kSwissCross, 0.95, 2., 15., 0.999);
0118   else if (id.subdetId() == EcalEndcap && theEERecHitCollection_ != nullptr)
0119     severityLevel = theEcalSevLvlAlgo_->severityLevel(
0120         hit);  //id, *theEERecHitCollection_, *theEcalChStatus_, 5., EcalSeverityLevelAlgo::kSwissCross, 0.95, 2., 15., 0.999);
0121   else
0122     return false;
0123
0124   if (severityLevel == EcalSeverityLevel::kGood)
0125     return true;
0126   if (severityLevel == EcalSeverityLevel::kRecovered)
0127     return UseEcalRecoveredHits_;
0128   if (severityLevel > static_cast<int>(EcalAcceptSeverityLevel_))
0129     return false;
0130   else
0131     return true;
0132 }
0133
0134 bool ObjectValidator::validTrack(const reco::Track& trk) const {
0135   if (trk.pt() < MinValidTrackPt_)
0136     return false;
0137   if (trk.pt() < MinValidTrackPtBarrel_ && std::fabs(trk.momentum().eta()) < 1.479)
0138     return false;
0139   if (trk.numberOfValidHits() < MinValidTrackNHits_)
0140     return false;
0141   return true;
0142 }
0143
0144 ////////////////////////////////////////////////////////////
0145 //
0146 // PhysicsTowerOrganizer
0147 //
0148 ////////////////////////////////////////////////////////////
0149
0150 PhysicsTowerOrganizer::PhysicsTowerOrganizer(
0151     const edm::Handle<HBHERecHitCollection>& hbhehitcoll_h,
0152     const edm::Handle<EcalRecHitCollection>& ebhitcoll_h,
0153     const edm::Handle<EcalRecHitCollection>& eehitcoll_h,
0154     const edm::Handle<std::vector<reco::TrackExtrapolation> >& trackextrapcoll_h,
0155     const ObjectValidatorAbs& objectvalidator,
0156     const CaloTowerConstituentsMap& ctcm,
0157     const CaloGeometry& geo) {
0158   // get some geometries
0159   const CaloSubdetectorGeometry* gEB = geo.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0160   const CaloSubdetectorGeometry* gEE = geo.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0161
0162   // do the HCAL hits
0163   for (HBHERecHitCollection::const_iterator it = hbhehitcoll_h->begin(); it != hbhehitcoll_h->end(); ++it) {
0164     const HBHERecHit* hit = &(*it);
0165
0166     // check that the hit is valid
0167     if (!objectvalidator.validHit(*hit))
0168       continue;
0169
0170     // add the hit to the organizer
0171     CaloTowerDetId tid = ctcm.towerOf(hit->id());
0172     insert_(tid, hit);
0173   }
0174
0175   // do the EB hits
0176   for (EcalRecHitCollection::const_iterator it = ebhitcoll_h->begin(); it != ebhitcoll_h->end(); ++it) {
0177     const EcalRecHit* hit = &(*it);
0178
0179     if (!objectvalidator.validHit(*hit))
0180       continue;
0181     CaloTowerDetId tid = ctcm.towerOf(hit->id());
0182     insert_(tid, hit);
0183   }
0184
0185   // do the EE hits
0186   for (EcalRecHitCollection::const_iterator it = eehitcoll_h->begin(); it != eehitcoll_h->end(); ++it) {
0187     const EcalRecHit* hit = &(*it);
0188
0189     if (!objectvalidator.validHit(*hit))
0190       continue;
0191     CaloTowerDetId tid = ctcm.towerOf(hit->id());
0192     insert_(tid, hit);
0193   }
0194
0195   // do the tracks
0196   for (std::vector<reco::TrackExtrapolation>::const_iterator it = trackextrapcoll_h->begin();
0197        it != trackextrapcoll_h->end();
0198        ++it) {
0199     const reco::TrackExtrapolation* extrap = &(*it);
0200     const reco::Track* track = &(*(extrap->track()));
0201
0202     // validate track
0203     if (!objectvalidator.validTrack(*track))
0204       continue;
0205
0206     // get the point
0207     if (extrap->positions().empty())
0208       continue;
0209     const GlobalPoint point(
0210         extrap->positions().front().x(), extrap->positions().front().y(), extrap->positions().front().z());
0211
0212     if (std::fabs(point.eta()) < 1.479) {
0213       EBDetId cell = gEB->getClosestCell(point);
0214       CaloTowerDetId tid = ctcm.towerOf(cell);
0215       insert_(tid, track);
0216     } else {
0217       EEDetId cell = gEE->getClosestCell(point);
0218       CaloTowerDetId tid = ctcm.towerOf(cell);
0219       insert_(tid, track);
0220     }
0221   }
0222
0223   return;
0224 }
0225
0226 PhysicsTower* PhysicsTowerOrganizer::findTower(const CaloTowerDetId& id) {
0227   // create dummy PhysicsTower
0228   PhysicsTower dummy;
0229
0230   // correct for the merging of the |ieta|=28-29 towers
0231   if (id.ietaAbs() == 29)
0232     dummy.id = CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi());
0233   else
0234     dummy.id = id;
0235
0236   // search on the dummy
0237   std::set<PhysicsTower, towercmp>::iterator it = towers_.find(dummy);
0238
0239   if (it == towers_.end())
0240     return nullptr;
0241
0242   // for whatever reason, I can't get a non-const out of the find method
0243   PhysicsTower& twr = const_cast<PhysicsTower&>(*it);
0244   return &twr;
0245 }
0246
0247 const PhysicsTower* PhysicsTowerOrganizer::findTower(const CaloTowerDetId& id) const {
0248   // create dummy PhysicsTower
0249   PhysicsTower dummy;
0250
0251   // correct for the merging of the |ieta|=28-29 towers
0252   if (id.ietaAbs() == 29)
0253     dummy.id = CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi());
0254   else
0255     dummy.id = id;
0256
0257   // search on the dummy
0258   std::set<PhysicsTower, towercmp>::iterator it = towers_.find(dummy);
0259
0260   if (it == towers_.end())
0261     return nullptr;
0262   return &(*it);
0263 }
0264
0265 const PhysicsTower* PhysicsTowerOrganizer::findTower(int ieta, int iphi) const {
0266   CaloTowerDetId tid(ieta, iphi);
0267   return findTower(tid);
0268 }
0269
0270 PhysicsTower* PhysicsTowerOrganizer::findTower(int ieta, int iphi) {
0271   CaloTowerDetId tid(ieta, iphi);
0272   return findTower(tid);
0273 }
0274
0275 void PhysicsTowerOrganizer::findNeighbors(const CaloTowerDetId& tempid,
0276                                           std::set<const PhysicsTower*>& neighbors) const {
0277   // correct for the merging of the |ieta|=28-29 towers
0278   CaloTowerDetId id(tempid);
0279   if (tempid.ietaAbs() == 29)
0280     id = CaloTowerDetId((tempid.ietaAbs() - 1) * tempid.zside(), tempid.iphi());
0281
0282   std::vector<CaloTowerDetId> ids;
0283   // get the neighbor with higher iphi
0284   if (id.ietaAbs() <= 20) {
0285     if (id.iphi() == 72)
0286       ids.push_back(CaloTowerDetId(id.ieta(), 1));
0287     else
0288       ids.push_back(CaloTowerDetId(id.ieta(), id.iphi() + 1));
0289   } else {
0290     if (id.iphi() == 71)
0291       ids.push_back(CaloTowerDetId(id.ieta(), 1));
0292     else
0293       ids.push_back(CaloTowerDetId(id.ieta(), id.iphi() + 2));
0294   }
0295
0296   // get the neighbor with the lower iphi
0297   if (id.ietaAbs() <= 20) {
0298     if (id.iphi() == 1)
0299       ids.push_back(CaloTowerDetId(id.ieta(), 72));
0300     else
0301       ids.push_back(CaloTowerDetId(id.ieta(), id.iphi() - 1));
0302   } else {
0303     if (id.iphi() == 1)
0304       ids.push_back(CaloTowerDetId(id.ieta(), 71));
0305     else
0306       ids.push_back(CaloTowerDetId(id.ieta(), id.iphi() - 2));
0307   }
0308
0309   // get the neighbor with the higher ietaAbs
0310   if (id.ietaAbs() == 20 && (id.iphi() % 2) == 0)
0311     ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), id.iphi() - 1));
0312   else
0313     ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), id.iphi()));
0314
0315   // get the neighbor(s) with the lower ietaAbs
0316   if (id.ietaAbs() == 21) {
0317     ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi()));
0318     ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi() + 1));
0319   } else if (id.ietaAbs() == 1) {
0320     ids.push_back(CaloTowerDetId(-id.ieta(), id.iphi()));
0321   } else {
0322     ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi()));
0323   }
0324
0325   // get the neighbor with higher ieta and higher iphi
0326   if (id.ietaAbs() <= 19 || (id.ietaAbs() == 20 && (id.iphi() % 2) == 0)) {
0327     if (id.iphi() == 72)
0328       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), 1));
0329     else
0330       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), id.iphi() + 1));
0331   } else if (id.ietaAbs() >= 21) {
0332     if (id.iphi() == 71)
0333       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), 1));
0334     else
0335       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), id.iphi() + 2));
0336   }
0337
0338   // get the neighbor with higher ieta and lower iphi
0339   if (id.ietaAbs() <= 19) {
0340     if (id.iphi() == 1)
0341       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), 72));
0342     else
0343       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), id.iphi() - 1));
0344   } else if (id.ietaAbs() >= 21 || (id.ietaAbs() == 20 && (id.iphi() % 2) == 1)) {
0345     if (id.iphi() == 1)
0346       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), 71));
0347     else
0348       ids.push_back(CaloTowerDetId((id.ietaAbs() + 1) * id.zside(), id.iphi() - 2));
0349   }
0350
0351   // get the neighbor with lower ieta and higher iphi
0352   if (id.ietaAbs() == 1) {
0353     if (id.iphi() == 72)
0354       ids.push_back(CaloTowerDetId(-id.ieta(), 1));
0355     else
0356       ids.push_back(CaloTowerDetId(-id.ieta(), id.iphi() + 1));
0357   } else if (id.ietaAbs() <= 20) {
0358     if (id.iphi() == 72)
0359       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), 1));
0360     else
0361       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi() + 1));
0362   } else if (id.ietaAbs() >= 21) {
0363     if (id.iphi() == 71)
0364       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), 1));
0365     else
0366       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi() + 2));
0367   }
0368
0369   // get the neighbor with lower ieta and lower iphi
0370   if (id.ietaAbs() == 1) {
0371     if (id.iphi() == 1)
0372       ids.push_back(CaloTowerDetId(-id.ieta(), 72));
0373     else
0374       ids.push_back(CaloTowerDetId(-id.ieta(), id.iphi() - 1));
0375   } else if (id.ietaAbs() <= 20) {
0376     if (id.iphi() == 1)
0377       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), 72));
0378     else
0379       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi() - 1));
0380   } else if (id.ietaAbs() >= 22) {
0381     if (id.iphi() == 1)
0382       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), 71));
0383     else
0384       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi() - 2));
0385   } else if (id.ietaAbs() == 21) {
0386     if (id.iphi() == 1)
0387       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), 72));
0388     else
0389       ids.push_back(CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi() - 1));
0390   }
0391
0392   // clear neighbors
0393   neighbors.clear();
0394
0395   // find the neighbors and add them to the eponymous set
0396   for (std::vector<CaloTowerDetId>::const_iterator it = ids.begin(); it != ids.end(); ++it) {
0397     const PhysicsTower* twr = findTower(*it);
0398     if (twr)
0399       neighbors.insert(twr);
0400   }
0401
0402   return;
0403 }
0404
0405 void PhysicsTowerOrganizer::findNeighbors(const PhysicsTower* twr, std::set<const PhysicsTower*>& neighbors) const {
0406   findNeighbors(twr->id, neighbors);
0407   return;
0408 }
0409
0410 void PhysicsTowerOrganizer::findNeighbors(int ieta, int iphi, std::set<const PhysicsTower*>& neighbors) const {
0411   findNeighbors(CaloTowerDetId(ieta, iphi), neighbors);
0412   return;
0413 }
0414
0415 void PhysicsTowerOrganizer::insert_(CaloTowerDetId& id, const HBHERecHit* hit) {
0416   PhysicsTower* twr = findTower(id);
0417   if (twr == nullptr) {
0418     PhysicsTower dummy;
0419     if (id.ietaAbs() == 29)
0420       dummy.id = CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi());
0421     else
0422       dummy.id = id;
0423     dummy.hcalhits.insert(hit);
0424     towers_.insert(dummy);
0425   } else {
0426     twr->hcalhits.insert(hit);
0427   }
0428   return;
0429 }
0430
0431 void PhysicsTowerOrganizer::insert_(CaloTowerDetId& id, const EcalRecHit* hit) {
0432   PhysicsTower* twr = findTower(id);
0433   if (twr == nullptr) {
0434     PhysicsTower dummy;
0435     if (id.ietaAbs() == 29)
0436       dummy.id = CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi());
0437     else
0438       dummy.id = id;
0439     dummy.ecalhits.insert(hit);
0440     towers_.insert(dummy);
0441   } else {
0442     twr->ecalhits.insert(hit);
0443   }
0444   return;
0445 }
0446
0447 void PhysicsTowerOrganizer::insert_(CaloTowerDetId& id, const reco::Track* track) {
0448   PhysicsTower* twr = findTower(id);
0449   if (twr == nullptr) {
0450     PhysicsTower dummy;
0451     if (id.ietaAbs() == 29)
0452       dummy.id = CaloTowerDetId((id.ietaAbs() - 1) * id.zside(), id.iphi());
0453     else
0454       dummy.id = id;
0455     dummy.tracks.insert(track);
0456     towers_.insert(dummy);
0457   } else {
0458     twr->tracks.insert(track);
0459   }
0460   return;
0461 }
0462
0463 ////////////////////////////////////////////////////////////
0464 //
0465 // HBHEHitMap
0466 //
0467 ////////////////////////////////////////////////////////////
0468
0469 HBHEHitMap::HBHEHitMap() {
0470   hitEnergy_ = hitEnergyTrkFid_ = -999.;
0471   nHits_ = -999;
0472   hcalEnergySameTowers_ = ecalEnergySameTowers_ = trackEnergySameTowers_ = -999.;
0473   nHcalHitsSameTowers_ = nEcalHitsSameTowers_ = nTracksSameTowers_ = -999;
0474   hcalEnergyNeighborTowers_ = ecalEnergyNeighborTowers_ = trackEnergyNeighborTowers_ = -999.;
0475   nHcalHitsNeighborTowers_ = nEcalHitsNeighborTowers_ = nTracksNeighborTowers_ = -999;
0476 }
0477
0478 double HBHEHitMap::hitEnergy(void) const {
0479   if (hitEnergy_ < -900)
0480     calcHits_();
0481   return hitEnergy_;
0482 }
0483
0484 int HBHEHitMap::nHits(void) const {
0485   if (nHits_ < -900)
0486     calcHits_();
0487   return nHits_;
0488 }
0489
0490 double HBHEHitMap::hitEnergyTrackFiducial(void) const {
0491   if (hitEnergyTrkFid_ < -900)
0492     calcHits_();
0493   return hitEnergyTrkFid_;
0494 }
0495
0496 double HBHEHitMap::hcalEnergySameTowers(void) const {
0497   if (hcalEnergySameTowers_ < -900)
0498     calcHcalSameTowers_();
0499   return hcalEnergySameTowers_;
0500 }
0501
0502 int HBHEHitMap::nHcalHitsSameTowers(void) const {
0503   if (nHcalHitsSameTowers_ < -900)
0504     calcHcalSameTowers_();
0505   return nHcalHitsSameTowers_;
0506 }
0507
0508 double HBHEHitMap::ecalEnergySameTowers(void) const {
0509   if (ecalEnergySameTowers_ < -900)
0510     calcEcalSameTowers_();
0511   return ecalEnergySameTowers_;
0512 }
0513
0514 int HBHEHitMap::nEcalHitsSameTowers(void) const {
0515   if (nEcalHitsSameTowers_ < -900)
0516     calcEcalSameTowers_();
0517   return nEcalHitsSameTowers_;
0518 }
0519
0520 double HBHEHitMap::trackEnergySameTowers(void) const {
0521   if (trackEnergySameTowers_ < -900)
0522     calcTracksSameTowers_();
0523   return trackEnergySameTowers_;
0524 }
0525
0526 int HBHEHitMap::nTracksSameTowers(void) const {
0527   if (nTracksSameTowers_ < -900)
0528     calcTracksSameTowers_();
0529   return nTracksSameTowers_;
0530 }
0531
0532 void HBHEHitMap::hcalHitsSameTowers(std::set<const HBHERecHit*>& v) const {
0533   v.clear();
0534   for (hitmap_const_iterator it1 = beginHits(); it1 != endHits(); ++it1) {
0535     for (std::set<const HBHERecHit*>::const_iterator it2 = it1->second->hcalhits.begin();
0536          it2 != it1->second->hcalhits.end();
0537          ++it2) {
0538       const HBHERecHit* hit = (*it2);
0539       // if the hit in the tower is already in the hitmap, don't include it
0540       if (findHit(hit) == endHits())
0541         v.insert(hit);
0542     }
0543   }
0544   return;
0545 }
0546
0547 void HBHEHitMap::ecalHitsSameTowers(std::set<const EcalRecHit*>& v) const {
0548   v.clear();
0549   for (hitmap_const_iterator it1 = beginHits(); it1 != endHits(); ++it1) {
0550     v.insert(it1->second->ecalhits.begin(), it1->second->ecalhits.end());
0551   }
0552   return;
0553 }
0554
0555 void HBHEHitMap::tracksSameTowers(std::set<const reco::Track*>& v) const {
0556   v.clear();
0557   for (hitmap_const_iterator it1 = beginHits(); it1 != endHits(); ++it1) {
0558     v.insert(it1->second->tracks.begin(), it1->second->tracks.end());
0559   }
0560   return;
0561 }
0562
0563 double HBHEHitMap::hcalEnergyNeighborTowers(void) const {
0564   if (hcalEnergyNeighborTowers_ < -900)
0565     calcHcalNeighborTowers_();
0566   return hcalEnergyNeighborTowers_;
0567 }
0568
0569 int HBHEHitMap::nHcalHitsNeighborTowers(void) const {
0570   if (nHcalHitsNeighborTowers_ < -900)
0571     calcHcalNeighborTowers_();
0572   return nHcalHitsNeighborTowers_;
0573 }
0574
0575 double HBHEHitMap::ecalEnergyNeighborTowers(void) const {
0576   if (ecalEnergyNeighborTowers_ < -900)
0577     calcEcalNeighborTowers_();
0578   return ecalEnergyNeighborTowers_;
0579 }
0580
0581 int HBHEHitMap::nEcalHitsNeighborTowers(void) const {
0582   if (nEcalHitsNeighborTowers_ < -900)
0583     calcEcalNeighborTowers_();
0584   return nEcalHitsNeighborTowers_;
0585 }
0586
0587 double HBHEHitMap::trackEnergyNeighborTowers(void) const {
0588   if (trackEnergyNeighborTowers_ < -900)
0589     calcTracksNeighborTowers_();
0590   return trackEnergyNeighborTowers_;
0591 }
0592
0593 int HBHEHitMap::nTracksNeighborTowers(void) const {
0594   if (nTracksNeighborTowers_ < -900)
0595     calcTracksNeighborTowers_();
0596   return nTracksNeighborTowers_;
0597 }
0598
0599 void HBHEHitMap::hcalHitsNeighborTowers(std::set<const HBHERecHit*>& v) const {
0600   v.clear();
0601   for (neighbor_const_iterator it1 = beginNeighbors(); it1 != endNeighbors(); ++it1) {
0602     const PhysicsTower* twr = (*it1);
0603     v.insert(twr->hcalhits.begin(), twr->hcalhits.end());
0604   }
0605   return;
0606 }
0607
0608 void HBHEHitMap::ecalHitsNeighborTowers(std::set<const EcalRecHit*>& v) const {
0609   v.clear();
0610   for (neighbor_const_iterator it1 = beginNeighbors(); it1 != endNeighbors(); ++it1) {
0611     const PhysicsTower* twr = (*it1);
0612     v.insert(twr->ecalhits.begin(), twr->ecalhits.end());
0613   }
0614
0615   return;
0616 }
0617
0618 void HBHEHitMap::tracksNeighborTowers(std::set<const reco::Track*>& v) const {
0619   v.clear();
0620   for (neighbor_const_iterator it1 = beginNeighbors(); it1 != endNeighbors(); ++it1) {
0621     const PhysicsTower* twr = (*it1);
0622     v.insert(twr->tracks.begin(), twr->tracks.end());
0623   }
0624   return;
0625 }
0626
0627 void HBHEHitMap::byTowers(std::vector<twrinfo>& v) const { assert(false); }
0628
0629 void HBHEHitMap::insert(const HBHERecHit* hit, const PhysicsTower* twr, std::set<const PhysicsTower*>& neighbors) {
0630   hits_[hit] = twr;
0631   neighbors_.insert(neighbors.begin(), neighbors.end());
0632
0633   // make sure none of the neighbors are also are part of the hitmap
0634   for (hitmap_const_iterator it = beginHits(); it != endHits(); ++it) {
0635     const PhysicsTower* t = it->second;
0636     neighbor_const_iterator find = findNeighbor(t);
0637
0638     // if a hit is also a neighbor, remove the neighbor
0639     if (find != endNeighbors())
0640       neighbors_.erase(find);
0641   }
0642   return;
0643 }
0644
0645 void HBHEHitMap::calcHits_(void) const {
0646   hitEnergy_ = 0;
0647   nHits_ = 0;
0648   hitEnergyTrkFid_ = 0;
0649   for (hitmap_const_iterator it = hits_.begin(); it != hits_.end(); ++it) {
0650     const HBHERecHit* hit = it->first;
0651     if (hit->id().ietaAbs() <= 26)
0652       hitEnergyTrkFid_ += hit->energy();
0653     hitEnergy_ += hit->energy();
0654     ++nHits_;
0655   }
0656   return;
0657 }
0658
0659 void HBHEHitMap::calcHcalSameTowers_(void) const {
0660   hcalEnergySameTowers_ = 0;
0661   nHcalHitsSameTowers_ = 0;
0662   std::set<const HBHERecHit*> v;
0663   hcalHitsSameTowers(v);
0664   for (std::set<const HBHERecHit*>::const_iterator it = v.begin(); it != v.end(); ++it) {
0665     const HBHERecHit* hit = (*it);
0666     hcalEnergySameTowers_ += hit->energy();
0667     ++nHcalHitsSameTowers_;
0668   }
0669   return;
0670 }
0671
0672 void HBHEHitMap::calcEcalSameTowers_(void) const {
0673   ecalEnergySameTowers_ = 0;
0674   nEcalHitsSameTowers_ = 0;
0675   std::set<const EcalRecHit*> v;
0676   ecalHitsSameTowers(v);
0677   for (std::set<const EcalRecHit*>::const_iterator it = v.begin(); it != v.end(); ++it) {
0678     const EcalRecHit* hit = (*it);
0679     ecalEnergySameTowers_ += hit->energy();
0680     ++nEcalHitsSameTowers_;
0681   }
0682   return;
0683 }
0684
0685 void HBHEHitMap::calcTracksSameTowers_(void) const {
0686   trackEnergySameTowers_ = 0;
0687   nTracksSameTowers_ = 0;
0688   std::set<const reco::Track*> v;
0689   tracksSameTowers(v);
0690   for (std::set<const reco::Track*>::const_iterator it = v.begin(); it != v.end(); ++it) {
0691     const reco::Track* trk = (*it);
0692     trackEnergySameTowers_ += trk->p();
0693     ++nTracksSameTowers_;
0694   }
0695   return;
0696 }
0697
0698 void HBHEHitMap::calcHcalNeighborTowers_(void) const {
0699   hcalEnergyNeighborTowers_ = 0;
0700   nHcalHitsNeighborTowers_ = 0;
0701   std::set<const HBHERecHit*> v;
0702   hcalHitsNeighborTowers(v);
0703   for (std::set<const HBHERecHit*>::const_iterator it = v.begin(); it != v.end(); ++it) {
0704     const HBHERecHit* hit = (*it);
0705     hcalEnergyNeighborTowers_ += hit->energy();
0706     ++nHcalHitsNeighborTowers_;
0707   }
0708   return;
0709 }
0710
0711 void HBHEHitMap::calcEcalNeighborTowers_(void) const {
0712   ecalEnergyNeighborTowers_ = 0;
0713   nEcalHitsNeighborTowers_ = 0;
0714   std::set<const EcalRecHit*> v;
0715   ecalHitsNeighborTowers(v);
0716   for (std::set<const EcalRecHit*>::const_iterator it = v.begin(); it != v.end(); ++it) {
0717     const EcalRecHit* hit = (*it);
0718     ecalEnergyNeighborTowers_ += hit->energy();
0719     ++nEcalHitsNeighborTowers_;
0720   }
0721   return;
0722 }
0723
0724 void HBHEHitMap::calcTracksNeighborTowers_(void) const {
0725   trackEnergyNeighborTowers_ = 0;
0726   nTracksNeighborTowers_ = 0;
0727   std::set<const reco::Track*> v;
0728   tracksNeighborTowers(v);
0729   for (std::set<const reco::Track*>::const_iterator it = v.begin(); it != v.end(); ++it) {
0730     const reco::Track* trk = (*it);
0731     trackEnergyNeighborTowers_ += trk->p();
0732     ++nTracksNeighborTowers_;
0733   }
0734   return;
0735 }
0736
0737 ////////////////////////////////////////////////////////////
0738 //
0739 // HBHEHitMapOrganizer
0740 //
0741 ////////////////////////////////////////////////////////////
0742
0743 HBHEHitMapOrganizer::HBHEHitMapOrganizer(const edm::Handle<HBHERecHitCollection>& hbhehitcoll_h,
0744                                          const ObjectValidatorAbs& objvalidator,
0745                                          const PhysicsTowerOrganizer& pto,
0746                                          const HcalFrontEndMap* hfemap)
0747     : hfemap_(hfemap) {
0748   // loop over the hits
0749   for (HBHERecHitCollection::const_iterator it = hbhehitcoll_h->begin(); it != hbhehitcoll_h->end(); ++it) {
0750     const HBHERecHit* hit = &(*it);
0751     if (!objvalidator.validHit(*hit))
0752       continue;
0753
0754     // get the Physics Tower and the neighbors
0755     const PhysicsTower* tower = pto.findTower(hit->id().ieta(), hit->id().iphi());
0756
0757     std::set<const PhysicsTower*> neighbors;
0758     pto.findNeighbors(hit->id().ieta(), hit->id().iphi(), neighbors);
0759
0760     // organize the RBXs
0761     int rbxidnum = hfemap_->lookupRBXIndex(hit->id());
0762     rbxs_[rbxidnum].insert(hit, tower, neighbors);
0763
0764     // organize the HPDs
0765     int hpdidnum = hfemap_->lookupRMIndex(hit->id());
0766     hpds_[hpdidnum].insert(hit, tower, neighbors);
0767
0768     // organize the dihits
0769     std::vector<const HBHERecHit*> hpdneighbors;
0770     getHPDNeighbors(hit, hpdneighbors, pto);
0771
0772     if (hpdneighbors.size() == 1) {
0773       std::vector<const HBHERecHit*> hpdneighborsneighbors;
0774       getHPDNeighbors(hpdneighbors[0], hpdneighborsneighbors, pto);
0775
0776       if (hpdneighborsneighbors.size() == 1 && hpdneighborsneighbors[0] == hit &&
0777           hit->energy() > hpdneighbors[0]->energy()) {
0778         // we've found two hits who are neighbors in the same HPD, but who have no other
0779         // neighbors (in the same HPD) in common.  In order not to double-count, we
0780         // require that the first hit has more energy
0781
0782         const PhysicsTower* tower2 = pto.findTower(hpdneighbors[0]->id().ieta(), hpdneighbors[0]->id().iphi());
0783         std::set<const PhysicsTower*> neighbors2;
0784         pto.findNeighbors(hpdneighbors[0]->id().ieta(), hpdneighbors[0]->id().iphi(), neighbors2);
0785
0786         HBHEHitMap dihit;
0787         dihit.insert(hit, tower, neighbors);
0788         dihit.insert(hpdneighbors[0], tower2, neighbors2);
0789         dihits_.push_back(dihit);
0790       }
0791     } else if (hpdneighbors.empty()) {
0792       // organize the monohits
0793       HBHEHitMap monohit;
0794       monohit.insert(hit, tower, neighbors);
0795       monohits_.push_back(monohit);
0796     }
0797
0798   }  // finished looping over HBHERecHits
0799   return;
0800 }
0801
0802 void HBHEHitMapOrganizer::getRBXs(std::vector<HBHEHitMap>& v, double energy) const {
0803   for (std::map<int, HBHEHitMap>::const_iterator it = rbxs_.begin(); it != rbxs_.end(); ++it) {
0804     const HBHEHitMap& map = it->second;
0805     if (map.hitEnergy() > energy)
0806       v.push_back(map);
0807   }
0808   return;
0809 }
0810
0811 void HBHEHitMapOrganizer::getHPDs(std::vector<HBHEHitMap>& v, double energy) const {
0812   for (std::map<int, HBHEHitMap>::const_iterator it = hpds_.begin(); it != hpds_.end(); ++it) {
0813     const HBHEHitMap& map = it->second;
0814     if (map.hitEnergy() > energy)
0815       v.push_back(map);
0816   }
0817   return;
0818 }
0819
0820 void HBHEHitMapOrganizer::getDiHits(std::vector<HBHEHitMap>& v, double energy) const {
0821   for (std::vector<HBHEHitMap>::const_iterator it = dihits_.begin(); it != dihits_.end(); ++it) {
0822     if (it->hitEnergy() > energy)
0823       v.push_back(*it);
0824   }
0825   return;
0826 }
0827
0828 void HBHEHitMapOrganizer::getMonoHits(std::vector<HBHEHitMap>& v, double energy) const {
0829   for (std::vector<HBHEHitMap>::const_iterator it = monohits_.begin(); it != monohits_.end(); ++it) {
0830     if (it->hitEnergy() > energy)
0831       v.push_back(*it);
0832   }
0833   return;
0834 }
0835
0836 void HBHEHitMapOrganizer::getHPDNeighbors(const HBHERecHit* hit,
0837                                           std::vector<const HBHERecHit*>& neighbors,
0838                                           const PhysicsTowerOrganizer& pto) {
0839   std::set<const PhysicsTower*> temp;
0840   pto.findNeighbors(hit->id().ieta(), hit->id().iphi(), temp);
0841
0842   // make sure to include the same tower that the hit is in
0843   temp.insert(pto.findTower(hit->id().ieta(), hit->id().iphi()));
0844
0845   // loop over the rechits in the temp neighbors
0846   for (std::set<const PhysicsTower*>::const_iterator it1 = temp.begin(); it1 != temp.end(); ++it1) {
0847     for (std::set<const HBHERecHit*>::const_iterator it2 = (*it1)->hcalhits.begin(); it2 != (*it1)->hcalhits.end();
0848          ++it2) {
0849       const HBHERecHit* hit2(*it2);
0850       if (hit != hit2 && hfemap_->lookupRMIndex(hit->id()) == hfemap_->lookupRMIndex(hit2->id())) {
0851         neighbors.push_back(hit2);
0852       }
0853     }
0854   }
0855   return;
0856 }