Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:51 

0001 /*
0002 Description: "Reflags" HB/HE hits based on their ECAL, HCAL, and tracking isolation.
0003 
0004 Original Author: John Paul Chou (Brown University)
0005                  Thursday, September 2, 2010
0006 */
0007 
0008 #include "HBHEIsolatedNoiseReflagger.h"
0009 
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0012 #include "DataFormats/JetReco/interface/TrackExtrapolation.h"
0013 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0014 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0015 #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h"
0016 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0017 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0018 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0019 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0020 #include "CondFormats/DataRecord/interface/HcalFrontEndMapRcd.h"
0021 #include "CondFormats/HcalObjects/interface/HcalFrontEndMap.h"
0022 
0023 HBHEIsolatedNoiseReflagger::HBHEIsolatedNoiseReflagger(const edm::ParameterSet& iConfig)
0024     :
0025 
0026       LooseHcalIsol_(iConfig.getParameter<double>("LooseHcalIsol")),
0027       LooseEcalIsol_(iConfig.getParameter<double>("LooseEcalIsol")),
0028       LooseTrackIsol_(iConfig.getParameter<double>("LooseTrackIsol")),
0029       TightHcalIsol_(iConfig.getParameter<double>("TightHcalIsol")),
0030       TightEcalIsol_(iConfig.getParameter<double>("TightEcalIsol")),
0031       TightTrackIsol_(iConfig.getParameter<double>("TightTrackIsol")),
0032 
0033       LooseRBXEne1_(iConfig.getParameter<double>("LooseRBXEne1")),
0034       LooseRBXEne2_(iConfig.getParameter<double>("LooseRBXEne2")),
0035       LooseRBXHits1_(iConfig.getParameter<int>("LooseRBXHits1")),
0036       LooseRBXHits2_(iConfig.getParameter<int>("LooseRBXHits2")),
0037       TightRBXEne1_(iConfig.getParameter<double>("TightRBXEne1")),
0038       TightRBXEne2_(iConfig.getParameter<double>("TightRBXEne2")),
0039       TightRBXHits1_(iConfig.getParameter<int>("TightRBXHits1")),
0040       TightRBXHits2_(iConfig.getParameter<int>("TightRBXHits2")),
0041 
0042       LooseHPDEne1_(iConfig.getParameter<double>("LooseHPDEne1")),
0043       LooseHPDEne2_(iConfig.getParameter<double>("LooseHPDEne2")),
0044       LooseHPDHits1_(iConfig.getParameter<int>("LooseHPDHits1")),
0045       LooseHPDHits2_(iConfig.getParameter<int>("LooseHPDHits2")),
0046       TightHPDEne1_(iConfig.getParameter<double>("TightHPDEne1")),
0047       TightHPDEne2_(iConfig.getParameter<double>("TightHPDEne2")),
0048       TightHPDHits1_(iConfig.getParameter<int>("TightHPDHits1")),
0049       TightHPDHits2_(iConfig.getParameter<int>("TightHPDHits2")),
0050 
0051       LooseDiHitEne_(iConfig.getParameter<double>("LooseDiHitEne")),
0052       TightDiHitEne_(iConfig.getParameter<double>("TightDiHitEne")),
0053       LooseMonoHitEne_(iConfig.getParameter<double>("LooseMonoHitEne")),
0054       TightMonoHitEne_(iConfig.getParameter<double>("TightMonoHitEne")),
0055 
0056       RBXEneThreshold_(iConfig.getParameter<double>("RBXEneThreshold")),
0057 
0058       debug_(iConfig.getUntrackedParameter<bool>("debug", true)),
0059       objvalidator_(iConfig) {
0060   tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"));
0061   tok_EB_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ebInput"));
0062   tok_EE_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("eeInput"));
0063   tok_trackExt_ =
0064       consumes<std::vector<reco::TrackExtrapolation> >(iConfig.getParameter<edm::InputTag>("trackExtrapolationInput"));
0065 
0066   // ES tokens
0067   ecalChStatusToken_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0068   hcalChStatusToken_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
0069   hcalSevToken_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
0070   ecalSevToken_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0071   ctcmToken_ = esConsumes<CaloTowerConstituentsMap, CaloGeometryRecord>();
0072   hfemapToken_ = esConsumes<HcalFrontEndMap, HcalFrontEndMapRcd>();
0073   geoToken_ = esConsumes();
0074 
0075   produces<HBHERecHitCollection>();
0076 }
0077 
0078 HBHEIsolatedNoiseReflagger::~HBHEIsolatedNoiseReflagger() {}
0079 
0080 void HBHEIsolatedNoiseReflagger::produce(edm::Event& iEvent, const edm::EventSetup& evSetup) {
0081   // get the ECAL channel status map
0082   const EcalChannelStatus* dbEcalChStatus = &evSetup.getData(ecalChStatusToken_);
0083 
0084   // get the HCAL channel status map
0085   const HcalChannelQuality* dbHcalChStatus = &evSetup.getData(hcalChStatusToken_);
0086 
0087   // get the severity level computers
0088   const HcalSeverityLevelComputer* hcalSevLvlComputer = &evSetup.getData(hcalSevToken_);
0089   const EcalSeverityLevelAlgo* ecalSevLvlAlgo = &evSetup.getData(ecalSevToken_);
0090 
0091   // get the calotower mappings
0092   const CaloTowerConstituentsMap& ctcm = evSetup.getData(ctcmToken_);
0093 
0094   // get hcal frontend map
0095   const HcalFrontEndMap* hfemap = &evSetup.getData(hfemapToken_);
0096 
0097   // get the HB/HE hits
0098   edm::Handle<HBHERecHitCollection> hbhehits_h;
0099   iEvent.getByToken(tok_hbhe_, hbhehits_h);
0100 
0101   // get the ECAL hits
0102   edm::Handle<EcalRecHitCollection> ebhits_h;
0103   iEvent.getByToken(tok_EB_, ebhits_h);
0104   edm::Handle<EcalRecHitCollection> eehits_h;
0105   iEvent.getByToken(tok_EE_, eehits_h);
0106 
0107   // get the tracks
0108   edm::Handle<std::vector<reco::TrackExtrapolation> > trackextraps_h;
0109   iEvent.getByToken(tok_trackExt_, trackextraps_h);
0110 
0111   // set the status maps and severity level computers for the hit validator
0112   objvalidator_.setHcalChannelQuality(dbHcalChStatus);
0113   objvalidator_.setEcalChannelStatus(dbEcalChStatus);
0114   objvalidator_.setHcalSeverityLevelComputer(hcalSevLvlComputer);
0115   objvalidator_.setEcalSeverityLevelAlgo(ecalSevLvlAlgo);
0116   objvalidator_.setEBRecHitCollection(&(*ebhits_h));
0117   objvalidator_.setEERecHitCollection(&(*eehits_h));
0118 
0119   // organizer the hits
0120   PhysicsTowerOrganizer pto(
0121       hbhehits_h, ebhits_h, eehits_h, trackextraps_h, objvalidator_, ctcm, evSetup.getData(geoToken_));
0122   HBHEHitMapOrganizer organizer(hbhehits_h, objvalidator_, pto, hfemap);
0123 
0124   // get the rbxs, hpds, dihits, and monohits
0125   std::vector<HBHEHitMap> rbxs;
0126   std::vector<HBHEHitMap> hpds;
0127   std::vector<HBHEHitMap> dihits;
0128   std::vector<HBHEHitMap> monohits;
0129   organizer.getRBXs(rbxs, LooseRBXEne1_ < TightRBXEne1_ ? LooseRBXEne1_ : TightRBXEne1_);
0130   organizer.getHPDs(hpds, LooseHPDEne1_ < TightHPDEne1_ ? LooseHPDEne1_ : TightHPDEne1_);
0131   organizer.getDiHits(dihits, LooseDiHitEne_ < TightDiHitEne_ ? LooseDiHitEne_ : TightDiHitEne_);
0132   organizer.getMonoHits(monohits, LooseMonoHitEne_ < TightMonoHitEne_ ? LooseMonoHitEne_ : TightMonoHitEne_);
0133 
0134   if (debug_ && (!rbxs.empty() || !hpds.empty() || !dihits.empty() || !monohits.empty())) {
0135     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "RBXs:" << std::endl;
0136     DumpHBHEHitMap(rbxs);
0137     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "\nHPDs:" << std::endl;
0138     DumpHBHEHitMap(hpds);
0139     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "\nDiHits:" << std::endl;
0140     DumpHBHEHitMap(dihits);
0141     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "\nMonoHits:" << std::endl;
0142     DumpHBHEHitMap(monohits);
0143   }
0144 
0145   //  bool result=true;
0146 
0147   // determine which hits are noisy
0148   std::set<const HBHERecHit*> noisehits;
0149   for (int i = 0; i < static_cast<int>(rbxs.size()); i++) {
0150     int nhits = rbxs[i].nHits();
0151     double ene = rbxs[i].hitEnergy();
0152     double trkfide = rbxs[i].hitEnergyTrackFiducial();
0153     double isolhcale = rbxs[i].hcalEnergySameTowers() + rbxs[i].hcalEnergyNeighborTowers();
0154     double isolecale = rbxs[i].ecalEnergySameTowers();
0155     double isoltrke = rbxs[i].trackEnergySameTowers() + rbxs[i].trackEnergyNeighborTowers();
0156     //
0157     // RBX mistag reduction
0158     bool isLooseIso = false;
0159     bool isTightIso = false;
0160     if (ene > RBXEneThreshold_ && ene > 0) {  // New absolute iso-cut for high energy RBX clusters
0161       if (isolhcale < LooseHcalIsol_ * RBXEneThreshold_ && isolecale < LooseEcalIsol_ * RBXEneThreshold_ &&
0162           isoltrke < LooseTrackIsol_ * RBXEneThreshold_)
0163         isLooseIso = true;
0164       if (isolhcale < TightHcalIsol_ * RBXEneThreshold_ && isolecale < TightEcalIsol_ * RBXEneThreshold_ &&
0165           isoltrke < TightTrackIsol_ * RBXEneThreshold_)
0166         isTightIso = true;
0167     }
0168     if (ene <= RBXEneThreshold_ && ene > 0) {  // Old relative iso-cut for low energy RBX clusters
0169       if (isolhcale / ene < LooseHcalIsol_ && isolecale / ene < LooseEcalIsol_ && isoltrke / ene < LooseTrackIsol_)
0170         isLooseIso = true;
0171       if (isolhcale / ene < TightHcalIsol_ && isolecale / ene < TightEcalIsol_ && isoltrke / ene < TightTrackIsol_)
0172         isTightIso = true;
0173     }
0174     //
0175     if ((isLooseIso && ((trkfide > LooseRBXEne1_ && nhits >= LooseRBXHits1_) ||
0176                         (trkfide > LooseRBXEne2_ && nhits >= LooseRBXHits2_))) ||
0177         (isTightIso && ((trkfide > TightRBXEne1_ && nhits >= TightRBXHits1_) ||
0178                         (trkfide > TightRBXEne2_ && nhits >= TightRBXHits2_)))) {
0179       for (HBHEHitMap::hitmap_const_iterator it = rbxs[i].beginHits(); it != rbxs[i].endHits(); ++it)
0180         noisehits.insert(it->first);
0181       //      result=false;
0182     }
0183   }
0184 
0185   for (int i = 0; i < static_cast<int>(hpds.size()); i++) {
0186     int nhits = hpds[i].nHits();
0187     double ene = hpds[i].hitEnergy();
0188     double trkfide = hpds[i].hitEnergyTrackFiducial();
0189     double isolhcale = hpds[i].hcalEnergySameTowers() + hpds[i].hcalEnergyNeighborTowers();
0190     double isolecale = hpds[i].ecalEnergySameTowers();
0191     double isoltrke = hpds[i].trackEnergySameTowers() + hpds[i].trackEnergyNeighborTowers();
0192     if ((ene > 0 && isolhcale / ene < LooseHcalIsol_ && isolecale / ene < LooseEcalIsol_ &&
0193          isoltrke / ene < LooseTrackIsol_ &&
0194          ((trkfide > LooseHPDEne1_ && nhits >= LooseHPDHits1_) ||
0195           (trkfide > LooseHPDEne2_ && nhits >= LooseHPDHits2_))) ||
0196         (ene > 0 && isolhcale / ene < TightHcalIsol_ && isolecale / ene < TightEcalIsol_ &&
0197          isoltrke / ene < TightTrackIsol_ &&
0198          ((trkfide > TightHPDEne1_ && nhits >= TightHPDHits1_) ||
0199           (trkfide > TightHPDEne2_ && nhits >= TightHPDHits2_)))) {
0200       for (HBHEHitMap::hitmap_const_iterator it = hpds[i].beginHits(); it != hpds[i].endHits(); ++it)
0201         noisehits.insert(it->first);
0202       //      result=false;
0203     }
0204   }
0205 
0206   for (int i = 0; i < static_cast<int>(dihits.size()); i++) {
0207     double ene = dihits[i].hitEnergy();
0208     double trkfide = dihits[i].hitEnergyTrackFiducial();
0209     double isolhcale = dihits[i].hcalEnergySameTowers() + dihits[i].hcalEnergyNeighborTowers();
0210     double isolecale = dihits[i].ecalEnergySameTowers();
0211     double isoltrke = dihits[i].trackEnergySameTowers() + dihits[i].trackEnergyNeighborTowers();
0212     if ((ene > 0 && isolhcale / ene < LooseHcalIsol_ && isolecale / ene < LooseEcalIsol_ &&
0213          isoltrke / ene < LooseTrackIsol_ && trkfide > 0.99 * ene && trkfide > LooseDiHitEne_) ||
0214         (ene > 0 && isolhcale / ene < TightHcalIsol_ && isolecale / ene < TightEcalIsol_ &&
0215          isoltrke / ene < TightTrackIsol_ && ene > TightDiHitEne_)) {
0216       for (HBHEHitMap::hitmap_const_iterator it = dihits[i].beginHits(); it != dihits[i].endHits(); ++it)
0217         noisehits.insert(it->first);
0218       //      result=false;
0219     }
0220   }
0221 
0222   for (int i = 0; i < static_cast<int>(monohits.size()); i++) {
0223     double ene = monohits[i].hitEnergy();
0224     double trkfide = monohits[i].hitEnergyTrackFiducial();
0225     double isolhcale = monohits[i].hcalEnergySameTowers() + monohits[i].hcalEnergyNeighborTowers();
0226     double isolecale = monohits[i].ecalEnergySameTowers();
0227     double isoltrke = monohits[i].trackEnergySameTowers() + monohits[i].trackEnergyNeighborTowers();
0228     if ((ene > 0 && isolhcale / ene < LooseHcalIsol_ && isolecale / ene < LooseEcalIsol_ &&
0229          isoltrke / ene < LooseTrackIsol_ && trkfide > 0.99 * ene && trkfide > LooseMonoHitEne_) ||
0230         (ene > 0 && isolhcale / ene < TightHcalIsol_ && isolecale / ene < TightEcalIsol_ &&
0231          isoltrke / ene < TightTrackIsol_ && ene > TightMonoHitEne_)) {
0232       for (HBHEHitMap::hitmap_const_iterator it = monohits[i].beginHits(); it != monohits[i].endHits(); ++it)
0233         noisehits.insert(it->first);
0234       //      result=false;
0235     }
0236   }
0237 
0238   // prepare the output HBHE RecHit collection
0239   auto pOut = std::make_unique<HBHERecHitCollection>();
0240   // loop over rechits, and set the new bit you wish to use
0241   for (HBHERecHitCollection::const_iterator it = hbhehits_h->begin(); it != hbhehits_h->end(); ++it) {
0242     const HBHERecHit* hit = &(*it);
0243     HBHERecHit newhit(*hit);
0244     if (noisehits.end() != noisehits.find(hit)) {
0245       newhit.setFlagField(1, HcalCaloFlagLabels::HBHEIsolatedNoise);
0246     }
0247     pOut->push_back(newhit);
0248   }
0249 
0250   iEvent.put(std::move(pOut));
0251 
0252   return;
0253 }
0254 
0255 void HBHEIsolatedNoiseReflagger::DumpHBHEHitMap(std::vector<HBHEHitMap>& i) const {
0256   for (std::vector<HBHEHitMap>::const_iterator it = i.begin(); it != i.end(); ++it) {
0257     edm::LogInfo("HBHEIsolatedNoiseReflagger")
0258         << "hit energy=" << it->hitEnergy() << "; # hits=" << it->nHits()
0259         << "; hcal energy same=" << it->hcalEnergySameTowers() << "; ecal energy same=" << it->ecalEnergySameTowers()
0260         << "; track energy same=" << it->trackEnergySameTowers()
0261         << "; neighbor hcal energy=" << it->hcalEnergyNeighborTowers() << std::endl;
0262     edm::LogInfo("HBHEIsolatedNoiseReflagger") << "hits:" << std::endl;
0263     for (HBHEHitMap::hitmap_const_iterator it2 = it->beginHits(); it2 != it->endHits(); ++it2) {
0264       const HBHERecHit* hit = it2->first;
0265       edm::LogInfo("HBHEIsolatedNoiseReflagger")
0266           << "RBX #=" << hfemap->lookupRBX(hit->id()) << "; HPD #=" << hfemap->lookupRMIndex(hit->id()) << "; "
0267           << (*hit) << std::endl;
0268     }
0269   }
0270   return;
0271 }
0272 
0273 //define this as a plug-in
0274 DEFINE_FWK_MODULE(HBHEIsolatedNoiseReflagger);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!