00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
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
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
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
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
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
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
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
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
00147 for(HBHERecHitCollection::const_iterator it=hbhehitcoll_h->begin(); it!=hbhehitcoll_h->end(); ++it) {
00148 const HBHERecHit* hit=&(*it);
00149
00150
00151 if(!objectvalidator.validHit(*hit)) continue;
00152
00153
00154 CaloTowerDetId tid = ctcm.towerOf(hit->id());
00155 insert_(tid, hit);
00156 }
00157
00158
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
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
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
00182 if(!objectvalidator.validTrack(*track)) continue;
00183
00184
00185 if(extrap->positions().size()<=0 || !(extrap->isValid().front())) continue;
00186
00187
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
00209 PhysicsTower dummy;
00210
00211
00212 if(id.ietaAbs()==29) dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00213 else dummy.id=id;
00214
00215
00216 std::set<PhysicsTower, towercmp>::iterator it=towers_.find(dummy);
00217
00218 if(it==towers_.end()) return 0;
00219
00220
00221 PhysicsTower &twr = const_cast<PhysicsTower&>(*it);
00222 return &twr;
00223 }
00224
00225 const PhysicsTower* PhysicsTowerOrganizer::findTower(const CaloTowerDetId& id) const
00226 {
00227
00228 PhysicsTower dummy;
00229
00230
00231 if(id.ietaAbs()==29) dummy.id = CaloTowerDetId((id.ietaAbs()-1)*id.zside(), id.iphi());
00232 else dummy.id=id;
00233
00234
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
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
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
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
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
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
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
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
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
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
00340 neighbors.clear();
00341
00342
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
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
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
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
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
00708
00710
00711
00712 HBHEHitMapOrganizer::HBHEHitMapOrganizer(const edm::Handle<HBHERecHitCollection>& hbhehitcoll_h,
00713 const ObjectValidatorAbs& objvalidator,
00714 const PhysicsTowerOrganizer& pto)
00715 {
00716
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
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
00728 int rbxidnum = HcalHPDRBXMap::indexRBX(hit->id());
00729 rbxs_[rbxidnum].insert(hit, tower, neighbors);
00730
00731
00732 int hpdidnum = HcalHPDRBXMap::indexHPD(hit->id());
00733 hpds_[hpdidnum].insert(hit, tower, neighbors);
00734
00735
00736
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
00746
00747
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
00761 HBHEHitMap monohit;
00762 monohit.insert(hit, tower, neighbors);
00763 monohits_.push_back(monohit);
00764 }
00765
00766 }
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
00812 temp.insert(pto.findTower(hit->id().ieta(), hit->id().iphi()));
00813
00814
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