CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoLocalCalo/HcalRecAlgos/src/HBHEIsolatedNoiseAlgos.cc

Go to the documentation of this file.
00001 /*
00002 Description: Isolation algorithms used to identify anomalous noise in the HB/HE.
00003              These algorithms will be used to reflag HB/HE rechits as noise.
00004 
00005              There are 4 objects implemented here:
00006              1) ObjectValidator
00007              2) PhysicsTowerOrganizer
00008              3) HBHEHitMap
00009              4) HBHEHitMapOrganizer
00010              See comments below for details.
00011 
00012 Original Author: John Paul Chou (Brown University)
00013                  Thursday, September 2, 2010
00014 */
00015 
00016 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHEIsolatedNoiseAlgos.h"
00017 
00018 #include "FWCore/Utilities/interface/EDMException.h"
00019 #include "FWCore/Framework/interface/ESHandle.h"
00020 
00021 #include "DataFormats/METReco/interface/CaloMETCollection.h"
00022 #include "DataFormats/METReco/interface/CaloMET.h"
00023 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00024 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00025 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00026 #include "DataFormats/TrackReco/interface/Track.h"
00027 
00028 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00029 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00030 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00031 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
00032 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
00033 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00034 #include "RecoMET/METAlgorithms/interface/HcalHPDRBXMap.h"
00035 
00037 //
00038 // ObjectValidator
00039 //
00041 
00042 ObjectValidator::ObjectValidator(const edm::ParameterSet& iConfig)
00043 {
00044   HBThreshold_  = iConfig.getParameter<double>("HBThreshold");
00045   HESThreshold_ = iConfig.getParameter<double>("HESThreshold");
00046   HEDThreshold_ = iConfig.getParameter<double>("HEDThreshold");
00047   EBThreshold_  = iConfig.getParameter<double>("EBThreshold");
00048   EEThreshold_  = iConfig.getParameter<double>("EEThreshold");
00049 
00050   HcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("HcalAcceptSeverityLevel");
00051   EcalAcceptSeverityLevel_ = iConfig.getParameter<uint32_t>("EcalAcceptSeverityLevel");
00052   UseHcalRecoveredHits_ = iConfig.getParameter<bool>("UseHcalRecoveredHits");
00053   UseEcalRecoveredHits_ = iConfig.getParameter<bool>("UseEcalRecoveredHits");
00054 
00055   MinValidTrackPt_ = iConfig.getParameter<double>("MinValidTrackPt");
00056   MinValidTrackPtBarrel_ = iConfig.getParameter<double>("MinValidTrackPtBarrel");
00057   MinValidTrackNHits_ = iConfig.getParameter<int>("MinValidTrackNHits");
00058 
00059   theHcalChStatus_=0;
00060   theEcalChStatus_=0;
00061   theHcalSevLvlComputer_=0;
00062   theEcalSevLvlAlgo_=0;
00063   theEBRecHitCollection_=0;
00064   theEERecHitCollection_=0;
00065 
00066   return;
00067 }
00068 
00069 ObjectValidator::~ObjectValidator()
00070 {
00071 }
00072 
00073 bool ObjectValidator::validHit(const HBHERecHit& hit) const
00074 {
00075   assert(theHcalSevLvlComputer_!=0 && theHcalChStatus_!=0);
00076 
00077   // require the hit to pass a certain energy threshold
00078   if(hit.id().subdet()==HcalBarrel && hit.energy()<HBThreshold_) return false;
00079   else if(hit.id().subdet()==HcalEndcap && hit.id().ietaAbs()<=20 && hit.energy()<HESThreshold_) return false;
00080   else if(hit.id().subdet()==HcalEndcap && hit.id().ietaAbs()>20 && hit.energy()<HEDThreshold_) return false;
00081 
00082   // determine if the hit is good, bad, or recovered
00083   const DetId id = hit.detid();
00084   const uint32_t recHitFlag = hit.flags();
00085   const uint32_t dbStatusFlag = theHcalChStatus_->getValues(id)->getValue();
00086   int severityLevel = theHcalSevLvlComputer_->getSeverityLevel(id, recHitFlag, dbStatusFlag);
00087   bool isRecovered  = theHcalSevLvlComputer_->recoveredRecHit(id, recHitFlag);
00088 
00089   if(severityLevel==0) return true;
00090   if(isRecovered) return UseHcalRecoveredHits_;
00091   if(severityLevel>static_cast<int>(HcalAcceptSeverityLevel_)) return false;
00092   else return true;
00093 }
00094 
00095 bool ObjectValidator::validHit(const EcalRecHit& hit) const
00096 {
00097   assert(theEcalSevLvlAlgo_!=0 && theEcalChStatus_!=0);
00098 
00099   // require the hit to pass a certain energy threshold
00100   const DetId id = hit.detid();
00101   if(id.subdetId()==EcalBarrel && hit.energy()<EBThreshold_) return false;
00102   else if(id.subdetId()==EcalEndcap && hit.energy()<EEThreshold_) return false;
00103 
00104   // determine if the hit is good, bad, or recovered
00105   int severityLevel = 999;
00106   if     (id.subdetId() == EcalBarrel && theEBRecHitCollection_!=0) severityLevel = theEcalSevLvlAlgo_->severityLevel(id, *theEBRecHitCollection_, *theEcalChStatus_, 5., EcalSeverityLevelAlgo::kSwissCross, 0.95, 2., 15., 0.999);
00107   else if(id.subdetId() == EcalEndcap && theEERecHitCollection_!=0) severityLevel = theEcalSevLvlAlgo_->severityLevel(id, *theEERecHitCollection_, *theEcalChStatus_, 5., EcalSeverityLevelAlgo::kSwissCross, 0.95, 2., 15., 0.999);
00108   else return false;
00109   
00110   if(severityLevel == EcalSeverityLevelAlgo::kGood) return true;
00111   if(severityLevel == EcalSeverityLevelAlgo::kRecovered) return UseEcalRecoveredHits_;
00112   if(severityLevel > static_cast<int>(EcalAcceptSeverityLevel_)) return false;
00113   else return true;
00114 }
00115 
00116 bool ObjectValidator::validTrack(const reco::Track& trk) const
00117 {
00118   if(trk.pt()<MinValidTrackPt_) return false;
00119   if(trk.pt()<MinValidTrackPtBarrel_  && std::fabs(trk.momentum().eta())<1.479) return false;
00120   if(trk.numberOfValidHits()<MinValidTrackNHits_) return false;
00121   return true;
00122 }
00123 
00125 //
00126 // PhysicsTowerOrganizer
00127 //
00129 
00130 PhysicsTowerOrganizer::PhysicsTowerOrganizer(const edm::Event& iEvent,
00131                                              const edm::EventSetup& evSetup,
00132                                              const edm::Handle<HBHERecHitCollection>& hbhehitcoll_h,
00133                                              const edm::Handle<EcalRecHitCollection>& ebhitcoll_h,
00134                                              const edm::Handle<EcalRecHitCollection>& eehitcoll_h,
00135                                              const edm::Handle<std::vector<reco::TrackExtrapolation> >& trackextrapcoll_h,
00136                                              const ObjectValidatorAbs& objectvalidator,
00137                                              const CaloTowerConstituentsMap& ctcm)
00138 {
00139   // get some geometries
00140   edm::ESHandle<CaloGeometry> pG;
00141   evSetup.get<CaloGeometryRecord>().get(pG);
00142   const CaloGeometry* geo = pG.product();
00143   const CaloSubdetectorGeometry* gEB = geo->getSubdetectorGeometry(DetId::Ecal,EcalBarrel);
00144   const CaloSubdetectorGeometry* gEE = geo->getSubdetectorGeometry(DetId::Ecal,EcalEndcap);
00145 
00146   // do the HCAL hits
00147   for(HBHERecHitCollection::const_iterator it=hbhehitcoll_h->begin(); it!=hbhehitcoll_h->end(); ++it) {
00148     const HBHERecHit* hit=&(*it);
00149       
00150     // check that the hit is valid
00151     if(!objectvalidator.validHit(*hit)) continue;
00152       
00153     // add the hit to the organizer
00154     CaloTowerDetId tid = ctcm.towerOf(hit->id());
00155     insert_(tid, hit);
00156   }
00157 
00158   // do the EB hits
00159   for(EcalRecHitCollection::const_iterator it=ebhitcoll_h->begin(); it!=ebhitcoll_h->end(); ++it) {
00160     const EcalRecHit* hit=&(*it);
00161       
00162     if(!objectvalidator.validHit(*hit)) continue;
00163     CaloTowerDetId tid = ctcm.towerOf(hit->id());
00164     insert_(tid, hit);
00165   }
00166 
00167   // do the EE hits
00168   for(EcalRecHitCollection::const_iterator it=eehitcoll_h->begin(); it!=eehitcoll_h->end(); ++it) {
00169     const EcalRecHit* hit=&(*it);
00170     
00171     if(!objectvalidator.validHit(*hit)) continue;
00172     CaloTowerDetId tid = ctcm.towerOf(hit->id());
00173     insert_(tid, hit);
00174   }
00175   
00176   // do the tracks
00177   for(std::vector<reco::TrackExtrapolation>::const_iterator it=trackextrapcoll_h->begin(); it!=trackextrapcoll_h->end(); ++it) {
00178     const reco::TrackExtrapolation* extrap=&(*it);
00179     const reco::Track* track = &(*(extrap->track()));
00180 
00181     // validate track
00182     if(!objectvalidator.validTrack(*track)) continue;
00183     
00184     // need a valid extrapolation point
00185     if(extrap->positions().size()<=0 || !(extrap->isValid().front())) continue;
00186     
00187     // get the point
00188     const GlobalPoint point(extrap->positions().front().x(),
00189                             extrap->positions().front().y(),
00190                             extrap->positions().front().z());
00191     
00192     if(std::fabs(point.eta())<1.479) {
00193       EBDetId cell = gEB->getClosestCell(point);
00194       CaloTowerDetId tid = ctcm.towerOf(cell);
00195       insert_(tid, track);
00196     } else {
00197       EEDetId cell = gEE->getClosestCell(point);
00198       CaloTowerDetId tid = ctcm.towerOf(cell);
00199       insert_(tid, track);
00200     }
00201   }
00202 
00203   return;
00204 }
00205 
00206 PhysicsTower* PhysicsTowerOrganizer::findTower(const CaloTowerDetId& id)
00207 {
00208   // create dummy PhysicsTower
00209   PhysicsTower dummy;
00210 
00211   // correct for the merging of the |ieta|=28-29 towers
00212   if(id.ietaAbs()==29) dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00213   else dummy.id=id;
00214 
00215   // search on the dummy
00216   std::set<PhysicsTower, towercmp>::iterator it=towers_.find(dummy);
00217   
00218   if(it==towers_.end()) return 0;
00219 
00220   // for whatever reason, I can't get a non-const out of the find method
00221   PhysicsTower &twr = const_cast<PhysicsTower&>(*it);
00222   return &twr;
00223 }
00224 
00225 const PhysicsTower* PhysicsTowerOrganizer::findTower(const CaloTowerDetId& id) const
00226 {
00227   // create dummy PhysicsTower
00228   PhysicsTower dummy;
00229 
00230   // correct for the merging of the |ieta|=28-29 towers
00231   if(id.ietaAbs()==29) dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00232   else dummy.id=id;
00233 
00234   // search on the dummy
00235   std::set<PhysicsTower, towercmp>::iterator it=towers_.find(dummy);
00236   
00237   if(it==towers_.end()) return 0;
00238   return &(*it);
00239 }
00240 
00241 const PhysicsTower* PhysicsTowerOrganizer::findTower(int ieta, int iphi) const
00242 {
00243   CaloTowerDetId tid(ieta, iphi);
00244   return findTower(tid);
00245 }
00246 
00247 PhysicsTower* PhysicsTowerOrganizer::findTower(int ieta, int iphi)
00248 {
00249   CaloTowerDetId tid(ieta, iphi);
00250   return findTower(tid);
00251 }
00252 
00253 void PhysicsTowerOrganizer::findNeighbors(const CaloTowerDetId& tempid, std::set<const PhysicsTower*>& neighbors) const
00254 {
00255   // correct for the merging of the |ieta|=28-29 towers
00256   CaloTowerDetId id(tempid);
00257   if(tempid.ietaAbs()==29) id = CaloTowerDetId((tempid.ietaAbs()-1)*tempid.zside(), tempid.iphi());
00258 
00259   std::vector<CaloTowerDetId> ids;
00260   // get the neighbor with higher iphi
00261   if(id.ietaAbs()<=20) {
00262     if(id.iphi()==72) ids.push_back(CaloTowerDetId(id.ieta(), 1));
00263     else              ids.push_back(CaloTowerDetId(id.ieta(), id.iphi()+1));
00264   } else {
00265     if(id.iphi()==71) ids.push_back(CaloTowerDetId(id.ieta(), 1));
00266     else              ids.push_back(CaloTowerDetId(id.ieta(), id.iphi()+2));
00267   }
00268 
00269   // get the neighbor with the lower iphi
00270   if(id.ietaAbs()<=20) {
00271     if(id.iphi()==1) ids.push_back(CaloTowerDetId(id.ieta(), 72));
00272     else             ids.push_back(CaloTowerDetId(id.ieta(), id.iphi()-1));
00273   } else {
00274     if(id.iphi()==1) ids.push_back(CaloTowerDetId(id.ieta(), 71));
00275     else             ids.push_back(CaloTowerDetId(id.ieta(), id.iphi()-2));
00276   }
00277 
00278   // get the neighbor with the higher ietaAbs
00279   if(id.ietaAbs()==20 && (id.iphi()%2)==0)
00280     ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), id.iphi()-1));
00281   else
00282     ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), id.iphi()));
00283 
00284   // get the neighbor(s) with the lower ietaAbs
00285   if(id.ietaAbs()==21) {
00286     ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()));
00287     ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()+1));
00288   } else if(id.ietaAbs()==1) {
00289     ids.push_back(CaloTowerDetId(-id.ieta(), id.iphi()));
00290   } else {
00291     ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()));
00292   }
00293   
00294   // get the neighbor with higher ieta and higher iphi
00295   if(id.ietaAbs()<=19 || (id.ietaAbs()==20 && (id.iphi()%2)==0)) {
00296     if(id.iphi()==72) ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), 1));
00297     else              ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), id.iphi()+1));
00298   } else if(id.ietaAbs()>=21) {
00299     if(id.iphi()==71) ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), 1));
00300     else              ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), id.iphi()+2));
00301   }
00302 
00303   // get the neighbor with higher ieta and lower iphi
00304   if(id.ietaAbs()<=19) {
00305     if(id.iphi()==1) ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), 72));
00306     else             ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), id.iphi()-1));
00307   } else if(id.ietaAbs()>=21 || (id.ietaAbs()==20 && (id.iphi()%2)==1)) {
00308     if(id.iphi()==1) ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), 71));
00309     else             ids.push_back(CaloTowerDetId((id.ietaAbs()+1)*id.zside(), id.iphi()-2));
00310   }
00311 
00312   // get the neighbor with lower ieta and higher iphi
00313   if(id.ietaAbs()==1) {
00314     if(id.iphi()==72) ids.push_back(CaloTowerDetId(-id.ieta(), 1));
00315     else              ids.push_back(CaloTowerDetId(-id.ieta(), id.iphi()+1));
00316   } else if(id.ietaAbs()<=20) {
00317     if(id.iphi()==72) ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), 1));
00318     else              ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()+1));
00319   } else if(id.ietaAbs()>=21) {
00320     if(id.iphi()==71) ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), 1));
00321     else              ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()+2));
00322   }
00323 
00324   // get the neighbor with lower ieta and lower iphi
00325   if(id.ietaAbs()==1) {
00326     if(id.iphi()==1) ids.push_back(CaloTowerDetId(-id.ieta(), 72));
00327     else             ids.push_back(CaloTowerDetId(-id.ieta(), id.iphi()-1));
00328   } else if(id.ietaAbs()<=20) {
00329     if(id.iphi()==1) ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), 72));
00330     else             ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()-1));
00331   } else if(id.ietaAbs()>=22) {
00332     if(id.iphi()==1) ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), 71));
00333     else             ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()-2));
00334   } else if(id.ietaAbs()==21) {
00335     if(id.iphi()==1) ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), 72));
00336     else             ids.push_back(CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi()-1));
00337   }
00338   
00339   // clear neighbors
00340   neighbors.clear();
00341 
00342   // find the neighbors and add them to the eponymous set
00343   for(std::vector<CaloTowerDetId>::const_iterator it=ids.begin(); it!=ids.end(); ++it) {
00344     const PhysicsTower* twr=findTower(*it);
00345     if(twr) neighbors.insert(twr);
00346   }
00347 
00348   return;
00349 }
00350 
00351 void PhysicsTowerOrganizer::findNeighbors(const PhysicsTower* twr, std::set<const PhysicsTower*>& neighbors) const
00352 {
00353   findNeighbors(twr->id, neighbors);
00354   return;
00355 }
00356 
00357 void PhysicsTowerOrganizer::findNeighbors(int ieta, int iphi, std::set<const PhysicsTower*>& neighbors) const
00358 {
00359   findNeighbors(CaloTowerDetId(ieta, iphi), neighbors);
00360   return;
00361 }
00362 
00363 void PhysicsTowerOrganizer::insert_(CaloTowerDetId& id, const HBHERecHit* hit)
00364 {
00365   PhysicsTower* twr=findTower(id);
00366   if(twr==0) {
00367     PhysicsTower dummy;
00368     if(id.ietaAbs()==29)
00369       dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00370     else
00371       dummy.id = id;
00372     dummy.hcalhits.insert(hit);
00373     towers_.insert(dummy);
00374   } else {
00375     twr->hcalhits.insert(hit);
00376   }
00377   return;
00378 }
00379 
00380 void PhysicsTowerOrganizer::insert_(CaloTowerDetId& id, const EcalRecHit* hit)
00381 {
00382   PhysicsTower* twr=findTower(id);
00383   if(twr==0) {
00384     PhysicsTower dummy;
00385     if(id.ietaAbs()==29)
00386       dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00387     else
00388       dummy.id = id;
00389     dummy.ecalhits.insert(hit);
00390     towers_.insert(dummy);
00391   } else {
00392     twr->ecalhits.insert(hit);
00393   }
00394   return;
00395 }
00396 
00397 void PhysicsTowerOrganizer::insert_(CaloTowerDetId& id, const reco::Track* track)
00398 {
00399   PhysicsTower* twr=findTower(id);
00400   if(twr==0) {
00401     PhysicsTower dummy;
00402     if(id.ietaAbs()==29)
00403       dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00404     else
00405       dummy.id = id;
00406     dummy.tracks.insert(track);
00407     towers_.insert(dummy);
00408   } else {
00409     twr->tracks.insert(track);
00410   }
00411   return;
00412 }
00413 
00414 
00416 //
00417 // HBHEHitMap
00418 //
00420 
00421 HBHEHitMap::HBHEHitMap()
00422 {
00423   hitEnergy_=hitEnergyTrkFid_=-999.;
00424   nHits_=-999;
00425   hcalEnergySameTowers_=ecalEnergySameTowers_=trackEnergySameTowers_=-999.;
00426   nHcalHitsSameTowers_=nEcalHitsSameTowers_=nTracksSameTowers_=-999;
00427   hcalEnergyNeighborTowers_=ecalEnergyNeighborTowers_=trackEnergyNeighborTowers_=-999.;
00428   nHcalHitsNeighborTowers_=nEcalHitsNeighborTowers_=nTracksNeighborTowers_=-999;
00429 
00430 }
00431 
00432 double HBHEHitMap::hitEnergy(void) const
00433 {
00434   if(hitEnergy_<-900) calcHits_();
00435   return hitEnergy_;
00436 }
00437 
00438 int HBHEHitMap::nHits(void) const
00439 {
00440   if(nHits_<-900) calcHits_();
00441   return nHits_;
00442 }
00443 
00444 double HBHEHitMap::hitEnergyTrackFiducial(void) const
00445 {
00446   if(hitEnergyTrkFid_<-900) calcHits_();
00447   return hitEnergyTrkFid_;
00448 }
00449 
00450 
00451 double HBHEHitMap::hcalEnergySameTowers(void) const
00452 {
00453   if(hcalEnergySameTowers_<-900) calcHcalSameTowers_();
00454   return hcalEnergySameTowers_;
00455 }
00456 
00457 int HBHEHitMap::nHcalHitsSameTowers(void) const
00458 {
00459   if(nHcalHitsSameTowers_<-900) calcHcalSameTowers_();
00460   return nHcalHitsSameTowers_;
00461 }
00462 
00463 double HBHEHitMap::ecalEnergySameTowers(void) const
00464 {
00465   if(ecalEnergySameTowers_<-900) calcEcalSameTowers_();
00466   return ecalEnergySameTowers_;
00467 }
00468 
00469 int HBHEHitMap::nEcalHitsSameTowers(void) const
00470 {
00471   if(nEcalHitsSameTowers_<-900) calcEcalSameTowers_();
00472   return nEcalHitsSameTowers_;
00473 }
00474 
00475 double HBHEHitMap::trackEnergySameTowers(void) const
00476 {
00477   if(trackEnergySameTowers_<-900) calcTracksSameTowers_();
00478   return trackEnergySameTowers_;
00479 }
00480 
00481 int HBHEHitMap::nTracksSameTowers(void) const
00482 {
00483   if(nTracksSameTowers_<-900) calcTracksSameTowers_();
00484   return nTracksSameTowers_;
00485 }
00486 
00487 void HBHEHitMap::hcalHitsSameTowers(std::set<const HBHERecHit*>& v) const
00488 {
00489   v.clear();
00490   for(hitmap_const_iterator it1=beginHits(); it1!=endHits(); ++it1) {
00491     for(std::set<const HBHERecHit*>::const_iterator it2=it1->second->hcalhits.begin(); it2!=it1->second->hcalhits.end(); ++it2) {
00492       const HBHERecHit* hit=(*it2);
00493       // if the hit in the tower is already in the hitmap, don't include it
00494       if(findHit(hit)==endHits()) v.insert(hit);
00495     }
00496   }
00497   return;
00498 }
00499 
00500 void HBHEHitMap::ecalHitsSameTowers(std::set<const EcalRecHit*>& v) const
00501 {
00502   v.clear();
00503   for(hitmap_const_iterator it1=beginHits(); it1!=endHits(); ++it1) {
00504     v.insert(it1->second->ecalhits.begin(), it1->second->ecalhits.end());
00505   }
00506   return;
00507 }
00508 
00509 void HBHEHitMap::tracksSameTowers(std::set<const reco::Track*>& v) const
00510 {
00511   v.clear();
00512   for(hitmap_const_iterator it1=beginHits(); it1!=endHits(); ++it1) {
00513     v.insert(it1->second->tracks.begin(), it1->second->tracks.end());
00514   }
00515   return;
00516 }
00517 
00518 double HBHEHitMap::hcalEnergyNeighborTowers(void) const
00519 {
00520   if(hcalEnergyNeighborTowers_<-900) calcHcalNeighborTowers_();
00521   return hcalEnergyNeighborTowers_;
00522 }
00523 
00524 int HBHEHitMap::nHcalHitsNeighborTowers(void) const
00525 {
00526   if(nHcalHitsNeighborTowers_<-900) calcHcalNeighborTowers_();
00527   return nHcalHitsNeighborTowers_;
00528 }
00529 
00530 double HBHEHitMap::ecalEnergyNeighborTowers(void) const
00531 {
00532   if(ecalEnergyNeighborTowers_<-900) calcEcalNeighborTowers_();
00533   return ecalEnergyNeighborTowers_;
00534 }
00535 
00536 int HBHEHitMap::nEcalHitsNeighborTowers(void) const
00537 {
00538   if(nEcalHitsNeighborTowers_<-900) calcEcalNeighborTowers_();
00539   return nEcalHitsNeighborTowers_;
00540 }
00541 
00542 double HBHEHitMap::trackEnergyNeighborTowers(void) const
00543 {
00544   if(trackEnergyNeighborTowers_<-900) calcTracksNeighborTowers_();
00545   return trackEnergyNeighborTowers_;
00546 }
00547 
00548 int HBHEHitMap::nTracksNeighborTowers(void) const
00549 {
00550   if(nTracksNeighborTowers_<-900) calcTracksNeighborTowers_();
00551   return nTracksNeighborTowers_;
00552 }
00553 
00554 void HBHEHitMap::hcalHitsNeighborTowers(std::set<const HBHERecHit*>& v) const
00555 {
00556   v.clear();
00557   for(neighbor_const_iterator it1=beginNeighbors(); it1!=endNeighbors(); ++it1) {
00558     const PhysicsTower* twr=(*it1);
00559     v.insert(twr->hcalhits.begin(), twr->hcalhits.end());
00560   }
00561   return;
00562 }
00563 
00564 void HBHEHitMap::ecalHitsNeighborTowers(std::set<const EcalRecHit*>& v) const
00565 {
00566   v.clear();
00567   for(neighbor_const_iterator it1=beginNeighbors(); it1!=endNeighbors(); ++it1) {
00568     const PhysicsTower* twr=(*it1);
00569     v.insert(twr->ecalhits.begin(), twr->ecalhits.end());
00570   }
00571 
00572   return;
00573 }
00574 
00575 void HBHEHitMap::tracksNeighborTowers(std::set<const reco::Track*>& v) const
00576 {
00577   v.clear();
00578   for(neighbor_const_iterator it1=beginNeighbors(); it1!=endNeighbors(); ++it1) {
00579     const PhysicsTower* twr=(*it1);
00580     v.insert(twr->tracks.begin(), twr->tracks.end());
00581   }
00582   return;
00583 }
00584 
00585 void HBHEHitMap::byTowers(std::vector<twrinfo>& v) const
00586 {
00587   assert(0);
00588 }
00589 
00590 void HBHEHitMap::insert(const HBHERecHit* hit, const PhysicsTower* twr, std::set<const PhysicsTower*>& neighbors)
00591 {
00592   hits_[hit]=twr;
00593   neighbors_.insert(neighbors.begin(), neighbors.end());
00594 
00595   // make sure none of the neighbors are also are part of the hitmap
00596   for(hitmap_const_iterator it=beginHits(); it!=endHits(); ++it) {
00597     const PhysicsTower* t=it->second;
00598     neighbor_const_iterator find=findNeighbor(t);
00599 
00600     // if a hit is also a neighbor, remove the neighbor
00601     if(find!=endNeighbors()) neighbors_.erase(find);
00602   }
00603   return;
00604 }
00605 
00606 
00607 void HBHEHitMap::calcHits_(void) const
00608 {
00609   hitEnergy_=0;
00610   nHits_=0;
00611   hitEnergyTrkFid_=0;
00612   for(hitmap_const_iterator it=hits_.begin(); it!=hits_.end(); ++it) {
00613     const HBHERecHit* hit=it->first;
00614     if(hit->id().ietaAbs()<=26) hitEnergyTrkFid_+=hit->energy();
00615     hitEnergy_+=hit->energy();
00616     ++nHits_;
00617   }
00618   return;
00619 }
00620 
00621 void HBHEHitMap::calcHcalSameTowers_(void) const
00622 {
00623   hcalEnergySameTowers_=0;
00624   nHcalHitsSameTowers_=0;
00625   std::set<const HBHERecHit*> v;
00626   hcalHitsSameTowers(v);
00627   for(std::set<const HBHERecHit*>::const_iterator it=v.begin(); it!=v.end(); ++it) {
00628     const HBHERecHit* hit=(*it);
00629     hcalEnergySameTowers_+=hit->energy();
00630     ++nHcalHitsSameTowers_;
00631   }
00632   return;
00633 }
00634 
00635 void HBHEHitMap::calcEcalSameTowers_(void) const
00636 {
00637   ecalEnergySameTowers_=0;
00638   nEcalHitsSameTowers_=0;
00639   std::set<const EcalRecHit*> v;
00640   ecalHitsSameTowers(v);
00641   for(std::set<const EcalRecHit*>::const_iterator it=v.begin(); it!=v.end(); ++it) {
00642     const EcalRecHit* hit=(*it);
00643     ecalEnergySameTowers_+=hit->energy();
00644     ++nEcalHitsSameTowers_;
00645   }
00646   return;
00647 }
00648 
00649 void HBHEHitMap::calcTracksSameTowers_(void) const
00650 {
00651   trackEnergySameTowers_=0;
00652   nTracksSameTowers_=0;
00653   std::set<const reco::Track*> v;
00654   tracksSameTowers(v);
00655   for(std::set<const reco::Track*>::const_iterator it=v.begin(); it!=v.end(); ++it) {
00656     const reco::Track* trk=(*it);
00657     trackEnergySameTowers_+=trk->p();
00658     ++nTracksSameTowers_;
00659   }
00660   return;
00661 }
00662 
00663 void HBHEHitMap::calcHcalNeighborTowers_(void) const
00664 {
00665   hcalEnergyNeighborTowers_=0;
00666   nHcalHitsNeighborTowers_=0;
00667   std::set<const HBHERecHit*> v;
00668   hcalHitsNeighborTowers(v);
00669   for(std::set<const HBHERecHit*>::const_iterator it=v.begin(); it!=v.end(); ++it) {
00670     const HBHERecHit* hit=(*it);
00671     hcalEnergyNeighborTowers_+=hit->energy();
00672     ++nHcalHitsNeighborTowers_;
00673   }
00674   return;
00675 }
00676 
00677 void HBHEHitMap::calcEcalNeighborTowers_(void) const
00678 {
00679   ecalEnergyNeighborTowers_=0;
00680   nEcalHitsNeighborTowers_=0;
00681   std::set<const EcalRecHit*> v;
00682   ecalHitsNeighborTowers(v);
00683   for(std::set<const EcalRecHit*>::const_iterator it=v.begin(); it!=v.end(); ++it) {
00684     const EcalRecHit* hit=(*it);
00685     ecalEnergyNeighborTowers_+=hit->energy();
00686     ++nEcalHitsNeighborTowers_;
00687   }
00688   return;
00689 }
00690 
00691 void HBHEHitMap::calcTracksNeighborTowers_(void) const
00692 {
00693   trackEnergyNeighborTowers_=0;
00694   nTracksNeighborTowers_=0;
00695   std::set<const reco::Track*> v;
00696   tracksNeighborTowers(v);
00697   for(std::set<const reco::Track*>::const_iterator it=v.begin(); it!=v.end(); ++it) {
00698     const reco::Track* trk=(*it);
00699     trackEnergyNeighborTowers_+=trk->p();
00700     ++nTracksNeighborTowers_;
00701   }
00702   return;
00703 }
00704 
00706 //
00707 // HBHEHitMapOrganizer
00708 //
00710 
00711 
00712 HBHEHitMapOrganizer::HBHEHitMapOrganizer(const edm::Handle<HBHERecHitCollection>& hbhehitcoll_h,
00713                                          const ObjectValidatorAbs& objvalidator,
00714                                          const PhysicsTowerOrganizer& pto)
00715 {
00716   // loop over the hits
00717   for(HBHERecHitCollection::const_iterator it=hbhehitcoll_h->begin(); it!=hbhehitcoll_h->end(); ++it) {
00718     const HBHERecHit *hit=&(*it);
00719     if(!objvalidator.validHit(*hit)) continue;
00720     
00721     // get the Physics Tower and the neighbors
00722     const PhysicsTower* tower=pto.findTower(hit->id().ieta(), hit->id().iphi());
00723     
00724     std::set<const PhysicsTower*> neighbors;
00725     pto.findNeighbors(hit->id().ieta(), hit->id().iphi(), neighbors);
00726     
00727     // organize the RBXs
00728     int rbxidnum = HcalHPDRBXMap::indexRBX(hit->id());
00729     rbxs_[rbxidnum].insert(hit, tower, neighbors);
00730     
00731     // organize the HPDs
00732     int hpdidnum = HcalHPDRBXMap::indexHPD(hit->id());
00733     hpds_[hpdidnum].insert(hit, tower, neighbors);
00734     
00735     
00736     // organize the dihits
00737     std::vector<const HBHERecHit*> hpdneighbors;
00738     getHPDNeighbors(hit, hpdneighbors, pto);
00739     
00740     if(hpdneighbors.size()==1) {
00741       std::vector<const HBHERecHit*> hpdneighborsneighbors;
00742       getHPDNeighbors(hpdneighbors[0], hpdneighborsneighbors, pto);
00743       
00744       if(hpdneighborsneighbors.size()==1 && hpdneighborsneighbors[0]==hit && hit->energy()>hpdneighbors[0]->energy()) {
00745         // we've found two hits who are neighbors in the same HPD, but who have no other
00746         // neighbors (in the same HPD) in common.  In order not to double-count, we
00747         // require that the first hit has more energy
00748         
00749         const PhysicsTower* tower2=pto.findTower(hpdneighbors[0]->id().ieta(), hpdneighbors[0]->id().iphi());
00750         std::set<const PhysicsTower*> neighbors2;
00751         pto.findNeighbors(hpdneighbors[0]->id().ieta(), hpdneighbors[0]->id().iphi(), neighbors2);
00752         
00753         HBHEHitMap dihit;
00754         dihit.insert(hit, tower, neighbors);
00755         dihit.insert(hpdneighbors[0], tower2, neighbors2);
00756         dihits_.push_back(dihit);
00757       }
00758     } else if(hpdneighbors.size()==0) {
00759       
00760       // organize the monohits
00761       HBHEHitMap monohit;
00762       monohit.insert(hit, tower, neighbors);
00763       monohits_.push_back(monohit);
00764     }
00765     
00766   } // finished looping over HBHERecHits
00767   return;
00768 }
00769 
00770 void HBHEHitMapOrganizer::getRBXs(std::vector<HBHEHitMap>& v, double energy) const
00771 {
00772   for(std::map<int, HBHEHitMap>::const_iterator it=rbxs_.begin(); it!=rbxs_.end(); ++it) {
00773     const HBHEHitMap &map=it->second;
00774     if(map.hitEnergy()>energy) v.push_back(map);
00775   }
00776   return;
00777 }
00778 
00779 void HBHEHitMapOrganizer::getHPDs(std::vector<HBHEHitMap>& v, double energy) const
00780 {
00781   for(std::map<int, HBHEHitMap>::const_iterator it=hpds_.begin(); it!=hpds_.end(); ++it) {
00782     const HBHEHitMap &map=it->second;
00783     if(map.hitEnergy()>energy) v.push_back(map);
00784   }
00785   return;
00786 }
00787 
00788 void HBHEHitMapOrganizer::getDiHits(std::vector<HBHEHitMap>& v, double energy) const
00789 {
00790   for(std::vector<HBHEHitMap>::const_iterator it=dihits_.begin(); it!=dihits_.end(); ++it) {
00791     if(it->hitEnergy()>energy) v.push_back(*it);
00792   }
00793   return;
00794 }
00795 
00796 void HBHEHitMapOrganizer::getMonoHits(std::vector<HBHEHitMap>& v, double energy) const
00797 {
00798   for(std::vector<HBHEHitMap>::const_iterator it=monohits_.begin(); it!=monohits_.end(); ++it) {
00799     if(it->hitEnergy()>energy) v.push_back(*it);
00800   }
00801   return;
00802 }
00803 
00804 
00805 
00806 void HBHEHitMapOrganizer::getHPDNeighbors(const HBHERecHit* hit, std::vector<const HBHERecHit*>& neighbors, const PhysicsTowerOrganizer& pto)
00807 {
00808   std::set<const PhysicsTower*> temp;
00809   pto.findNeighbors(hit->id().ieta(), hit->id().iphi(), temp);
00810 
00811   // make sure to include the same tower that the hit is in
00812   temp.insert(pto.findTower(hit->id().ieta(), hit->id().iphi()));
00813 
00814   // loop over the rechits in the temp neighbors
00815   for(std::set<const PhysicsTower*>::const_iterator it1=temp.begin(); it1!=temp.end(); ++it1) {
00816     for(std::set<const HBHERecHit*>::const_iterator it2=(*it1)->hcalhits.begin(); it2!=(*it1)->hcalhits.end(); ++it2) {
00817       const HBHERecHit* hit2(*it2);
00818       if(hit!=hit2 && HcalHPDRBXMap::indexHPD(hit->id())==HcalHPDRBXMap::indexHPD(hit2->id())) {
00819         neighbors.push_back(hit2);
00820       }
00821     }
00822   }
00823   return;
00824 }
00825 
00826 
00827