diff --git a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py index e627f103846db..cafddd10560d7 100644 --- a/HLTrigger/Configuration/python/customizeHLTforCMSSW.py +++ b/HLTrigger/Configuration/python/customizeHLTforCMSSW.py @@ -164,10 +164,41 @@ def customiseFor2017DtUnpacking(process): return process +def customiseFor29049(process) : + + listHltPFRecHitHBHE=['hltParticleFlowRecHitHBHE', + 'hltParticleFlowRecHitHBHEForEgamma', + 'hltParticleFlowRecHitHBHEForEgammaUnseeded', + 'hltParticleFlowRecHitHBHEForMuons', + 'hltParticleFlowRecHitHBHERegForMuons'] + for att in listHltPFRecHitHBHE: + if hasattr(process,att): + prod = getattr(process, att) + pset_navi = prod.navigator + if hasattr(pset_navi, "sigmaCut"): delattr(pset_navi,'sigmaCut') + if hasattr(pset_navi, "timeResolutionCalc"): delattr(pset_navi,'timeResolutionCalc') + pset_navi.name = cms.string("PFRecHitHCALDenseIdNavigator") + pset_navi.hcalEnums = cms.vint32(1,2) + + listHltPFRecHitHF=['hltParticleFlowRecHitHF', + 'hltParticleFlowRecHitHFForEgammaUnseeded'] + for att in listHltPFRecHitHF: + if hasattr(process,att): + prod = getattr(process, att) + pset_navi = prod.navigator + if hasattr(pset_navi, "barrel"): delattr(pset_navi,'barrel') + if hasattr(pset_navi, "endcap"): delattr(pset_navi,'endcap') + pset_navi.name = cms.string("PFRecHitHCALDenseIdNavigator") + pset_navi.hcalEnums = cms.vint32(4) + + return process + # CMSSW version specific customizations def customizeHLTforCMSSW(process, menuType="GRun"): # add call to action function in proper order: newest last! # process = customiseFor12718(process) + process = customiseFor29049(process) + return process diff --git a/RecoParticleFlow/PFClusterProducer/interface/HGCRecHitNavigator.h b/RecoParticleFlow/PFClusterProducer/interface/HGCRecHitNavigator.h index 583a2a933df5a..08895e7a6623f 100644 --- a/RecoParticleFlow/PFClusterProducer/interface/HGCRecHitNavigator.h +++ b/RecoParticleFlow/PFClusterProducer/interface/HGCRecHitNavigator.h @@ -59,13 +59,13 @@ class HGCRecHitNavigator : public PFRecHitNavigatorBase { } } - void beginEvent(const edm::EventSetup& iSetup) override { + void init(const edm::EventSetup& iSetup) override { if (nullptr != eeNav_) - eeNav_->beginEvent(iSetup); + eeNav_->init(iSetup); if (nullptr != hefNav_) - hefNav_->beginEvent(iSetup); + hefNav_->init(iSetup); if (nullptr != hebNav_) - hebNav_->beginEvent(iSetup); + hebNav_->init(iSetup); } void associateNeighbours(reco::PFRecHit& hit, diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFECALHashNavigator.h b/RecoParticleFlow/PFClusterProducer/interface/PFECALHashNavigator.h index 7c208762827ce..5e3dce958e0cb 100644 --- a/RecoParticleFlow/PFClusterProducer/interface/PFECALHashNavigator.h +++ b/RecoParticleFlow/PFClusterProducer/interface/PFECALHashNavigator.h @@ -26,7 +26,7 @@ class PFECALHashNavigator : public PFRecHitNavigatorBase { neighbourmapcalculated_ = false; } - void beginEvent(const edm::EventSetup& iSetup) override { + void init(const edm::EventSetup& iSetup) override { edm::ESHandle geoHandle; iSetup.get().get(geoHandle); diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFHCALDenseIdNavigator.h b/RecoParticleFlow/PFClusterProducer/interface/PFHCALDenseIdNavigator.h new file mode 100644 index 0000000000000..a81442514be8a --- /dev/null +++ b/RecoParticleFlow/PFClusterProducer/interface/PFHCALDenseIdNavigator.h @@ -0,0 +1,192 @@ +#ifndef RecoParticleFlow_PFClusterProducer_PFHCALDenseIdNavigator_h +#define RecoParticleFlow_PFClusterProducer_PFHCALDenseIdNavigator_h + +#include "FWCore/Framework/interface/ESWatcher.h" + +#include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" + +#include "RecoCaloTools/Navigation/interface/CaloNavigator.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" +#include "DataFormats/EcalDetId/interface/ESDetId.h" +#include "DataFormats/HcalDetId/interface/HcalDetId.h" + +#include "Geometry/CaloTopology/interface/EcalEndcapTopology.h" +#include "Geometry/CaloTopology/interface/EcalBarrelTopology.h" +#include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h" +#include "Geometry/CaloTopology/interface/HcalTopology.h" + +#include "Geometry/CaloTopology/interface/CaloTowerTopology.h" +#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h" + +template +class PFHCALDenseIdNavigator : public PFRecHitNavigatorBase { +public: + ~PFHCALDenseIdNavigator() override { + if (!ownsTopo) { + topology_.release(); + } + } + + PFHCALDenseIdNavigator(const edm::ParameterSet& iConfig) { + vhcalEnum_ = iConfig.getParameter>("hcalEnums"); + } + + void init(const edm::EventSetup& iSetup) override { + bool check = theRecNumberWatcher_.check(iSetup); + if (!check) + return; + + edm::ESHandle hcalTopology; + iSetup.get().get(hcalTopology); + topology_.release(); + topology_.reset(hcalTopology.product()); + + // Fill a vector of valid denseid's + edm::ESHandle hGeom; + iSetup.get().get(hGeom); + const CaloGeometry& caloGeom = *hGeom; + + std::vector vecHcal; + std::vector vDenseIdHcal; + neighboursHcal_.clear(); + for (auto hcalSubdet : vhcalEnum_) { + std::vector vecDetIds(caloGeom.getValidDetIds(DetId::Hcal, hcalSubdet)); + vecHcal.insert(vecHcal.end(), vecDetIds.begin(), vecDetIds.end()); + } + for (auto hDetId : vecHcal) { + vDenseIdHcal.push_back(topology_.get()->detId2denseId(hDetId)); + } + std::sort(vDenseIdHcal.begin(), vDenseIdHcal.end()); + + // Fill a vector of cell neighbours + denseIdHcalMax_ = *max_element(vDenseIdHcal.begin(), vDenseIdHcal.end()); + denseIdHcalMin_ = *min_element(vDenseIdHcal.begin(), vDenseIdHcal.end()); + neighboursHcal_.resize(denseIdHcalMax_ - denseIdHcalMin_ + 1); + + for (auto denseid : vDenseIdHcal) { + DetId N(0); + DetId E(0); + DetId S(0); + DetId W(0); + DetId NW(0); + DetId NE(0); + DetId SW(0); + DetId SE(0); + std::vector neighbours(9, DetId(0)); + + // the centre + unsigned denseid_c = denseid; + DetId detid_c = topology_.get()->denseId2detId(denseid_c); + CaloNavigator navigator(detid_c, topology_.get()); + + // Using enum in Geometry/CaloTopology/interface/CaloDirection.h + // Order: CENTER(NONE),SOUTH,SOUTHEAST,SOUTHWEST,EAST,WEST,NORTHEAST,NORTHWEST,NORTH + neighbours.at(NONE) = detid_c; + + navigator.home(); + N = navigator.north(); + neighbours.at(NORTH) = N; + if (N != DetId(0)) { + NE = navigator.east(); + } else { + navigator.home(); + E = navigator.east(); + NE = navigator.north(); + } + neighbours.at(NORTHEAST) = NE; + + navigator.home(); + S = navigator.south(); + neighbours.at(SOUTH) = S; + if (S != DetId(0)) { + SW = navigator.west(); + } else { + navigator.home(); + W = navigator.west(); + SW = navigator.south(); + } + neighbours.at(SOUTHWEST) = SW; + + navigator.home(); + E = navigator.east(); + neighbours.at(EAST) = E; + if (E != DetId(0)) { + SE = navigator.south(); + } else { + navigator.home(); + S = navigator.south(); + SE = navigator.east(); + } + neighbours.at(SOUTHEAST) = SE; + + navigator.home(); + W = navigator.west(); + neighbours.at(WEST) = W; + if (W != DetId(0)) { + NW = navigator.north(); + } else { + navigator.home(); + N = navigator.north(); + NW = navigator.west(); + } + neighbours.at(NORTHWEST) = NW; + + unsigned index = getIdx(denseid_c); + neighboursHcal_[index] = neighbours; + } + } + + void associateNeighbours(reco::PFRecHit& hit, + std::unique_ptr& hits, + edm::RefProd& refProd) override { + DetId detid(hit.detId()); + unsigned denseid = topology_.get()->detId2denseId(detid); + + std::vector neighbours(9, DetId(0)); + + if (denseid < denseIdHcalMin_ || denseid > denseIdHcalMax_) { + edm::LogWarning("PFRecHitHCALCachedNavigator") << " DenseId for this cell is out of the range." << std::endl; + } else if (!validNeighbours(denseid)) { + edm::LogWarning("PFRecHitHCALCachedNavigator") + << " DenseId for this cell does not have the neighbour information." << std::endl; + } else { + unsigned index = getIdx(denseid); + neighbours = neighboursHcal_.at(index); + } + + associateNeighbour(neighbours.at(NORTH), hit, hits, refProd, 0, 1, 0); // N + associateNeighbour(neighbours.at(NORTHEAST), hit, hits, refProd, 1, 1, 0); // NE + associateNeighbour(neighbours.at(SOUTH), hit, hits, refProd, 0, -1, 0); // S + associateNeighbour(neighbours.at(SOUTHWEST), hit, hits, refProd, -1, -1, 0); // SW + associateNeighbour(neighbours.at(EAST), hit, hits, refProd, 1, 0, 0); // E + associateNeighbour(neighbours.at(SOUTHEAST), hit, hits, refProd, 1, -1, 0); // SE + associateNeighbour(neighbours.at(WEST), hit, hits, refProd, -1, 0, 0); // W + associateNeighbour(neighbours.at(NORTHWEST), hit, hits, refProd, -1, 1, 0); // NW + } + + bool validNeighbours(const unsigned int denseid) const { + bool ok = true; + unsigned index = getIdx(denseid); + if (neighboursHcal_.at(index).size() != 9) + ok = false; // the neighbour vector size should be 3x3 + return ok; + } + + unsigned int getIdx(const unsigned int denseid) const { + unsigned index = denseid - denseIdHcalMin_; + return index; + } + +protected: + edm::ESWatcher theRecNumberWatcher_; + std::unique_ptr topology_; + std::vector vhcalEnum_; + std::vector> neighboursHcal_; + unsigned int denseIdHcalMax_; + unsigned int denseIdHcalMin_; +}; + +#endif diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitDualNavigator.h b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitDualNavigator.h index 8e0f84384e401..2915ca3240a1a 100644 --- a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitDualNavigator.h +++ b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitDualNavigator.h @@ -14,9 +14,9 @@ class PFRecHitDualNavigator : public PFRecHitNavigatorBase { endcapNav_ = new endcap(iConfig.getParameter("endcap")); } - void beginEvent(const edm::EventSetup& iSetup) override { - barrelNav_->beginEvent(iSetup); - endcapNav_->beginEvent(iSetup); + void init(const edm::EventSetup& iSetup) override { + barrelNav_->init(iSetup); + endcapNav_->init(iSetup); } void associateNeighbours(reco::PFRecHit& hit, diff --git a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h index 47ca3b21823e6..40960ee85a147 100644 --- a/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h +++ b/RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h @@ -32,7 +32,7 @@ class PFRecHitNavigatorBase { virtual ~PFRecHitNavigatorBase() = default; - virtual void beginEvent(const edm::EventSetup&) = 0; + virtual void init(const edm::EventSetup&) = 0; virtual void associateNeighbours(reco::PFRecHit&, std::unique_ptr&, edm::RefProd&) = 0; diff --git a/RecoParticleFlow/PFClusterProducer/plugins/Navigators.cc b/RecoParticleFlow/PFClusterProducer/plugins/Navigators.cc index e68bc0c3130bc..1dbf5c78b820b 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/Navigators.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/Navigators.cc @@ -6,6 +6,7 @@ #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitDualNavigator.h" #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigator.h" #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigatorWithTime.h" +#include "RecoParticleFlow/PFClusterProducer/interface/PFHCALDenseIdNavigator.h" #include "RecoParticleFlow/PFClusterProducer/interface/PFECALHashNavigator.h" #include "RecoParticleFlow/PFClusterProducer/interface/HGCRecHitNavigator.h" @@ -13,7 +14,7 @@ class PFRecHitEcalBarrelNavigatorWithTime : public PFRecHitCaloNavigatorWithTime public: PFRecHitEcalBarrelNavigatorWithTime(const edm::ParameterSet& iConfig) : PFRecHitCaloNavigatorWithTime(iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override { + void init(const edm::EventSetup& iSetup) override { edm::ESHandle geoHandle; iSetup.get().get(geoHandle); topology_ = std::make_unique(*geoHandle); @@ -24,7 +25,7 @@ class PFRecHitEcalEndcapNavigatorWithTime : public PFRecHitCaloNavigatorWithTime public: PFRecHitEcalEndcapNavigatorWithTime(const edm::ParameterSet& iConfig) : PFRecHitCaloNavigatorWithTime(iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override { + void init(const edm::EventSetup& iSetup) override { edm::ESHandle geoHandle; iSetup.get().get(geoHandle); topology_ = std::make_unique(*geoHandle); @@ -35,7 +36,7 @@ class PFRecHitEcalBarrelNavigator final : public PFRecHitCaloNavigator geoHandle; iSetup.get().get(geoHandle); topology_ = std::make_unique(*geoHandle); @@ -46,7 +47,7 @@ class PFRecHitEcalEndcapNavigator final : public PFRecHitCaloNavigator geoHandle; iSetup.get().get(geoHandle); topology_ = std::make_unique(*geoHandle); @@ -57,25 +58,31 @@ class PFRecHitPreshowerNavigator final : public PFRecHitCaloNavigator(); } + void init(const edm::EventSetup& iSetup) override { topology_ = std::make_unique(); } }; -class PFRecHitHCALNavigator final : public PFRecHitCaloNavigator { +class PFRecHitHCALDenseIdNavigator final : public PFHCALDenseIdNavigator { +public: + PFRecHitHCALDenseIdNavigator(const edm::ParameterSet& iConfig) : PFHCALDenseIdNavigator(iConfig) {} +}; + +class PFRecHitHCALNavigator : public PFRecHitCaloNavigator { public: PFRecHitHCALNavigator(const edm::ParameterSet& iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override { + void init(const edm::EventSetup& iSetup) override { edm::ESHandle hcalTopology; iSetup.get().get(hcalTopology); topology_.release(); topology_.reset(hcalTopology.product()); } }; + class PFRecHitHCALNavigatorWithTime : public PFRecHitCaloNavigatorWithTime { public: PFRecHitHCALNavigatorWithTime(const edm::ParameterSet& iConfig) : PFRecHitCaloNavigatorWithTime(iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override { + void init(const edm::EventSetup& iSetup) override { edm::ESHandle hcalTopology; iSetup.get().get(hcalTopology); topology_.release(); @@ -87,7 +94,7 @@ class PFRecHitCaloTowerNavigator : public PFRecHitCaloNavigator caloTowerTopology; iSetup.get().get(caloTowerTopology); topology_.release(); @@ -115,21 +122,21 @@ class PFRecHitHGCEENavigator : public PFRecHitFakeNavigator { public: PFRecHitHGCEENavigator(const edm::ParameterSet& iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override {} + void init(const edm::EventSetup& iSetup) override {} }; class PFRecHitHGCHENavigator : public PFRecHitFakeNavigator { public: PFRecHitHGCHENavigator(const edm::ParameterSet& iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override {} + void init(const edm::EventSetup& iSetup) override {} }; class PFRecHitHGCHexNavigator : public PFRecHitFakeNavigator { public: PFRecHitHGCHexNavigator(const edm::ParameterSet& iConfig) {} - void beginEvent(const edm::EventSetup& iSetup) override {} + void init(const edm::EventSetup& iSetup) override {} }; typedef HGCRecHitNavigator @@ -150,6 +157,7 @@ DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitECALNavigator, "PFRecHitECA DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitECALNavigatorWithTime, "PFRecHitECALNavigatorWithTime"); DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitCaloTowerNavigator, "PFRecHitCaloTowerNavigator"); DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitPreshowerNavigator, "PFRecHitPreshowerNavigator"); +DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitHCALDenseIdNavigator, "PFRecHitHCALDenseIdNavigator"); DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitHCALNavigator, "PFRecHitHCALNavigator"); DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitHCALNavigatorWithTime, "PFRecHitHCALNavigatorWithTime"); DEFINE_EDM_PLUGIN(PFRecHitNavigationFactory, PFRecHitHGCEENavigator, "PFRecHitHGCEENavigator"); diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc deleted file mode 100644 index 253a90576dc23..0000000000000 --- a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.cc +++ /dev/null @@ -1,789 +0,0 @@ -#include "RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.h" - -#include - -#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" - -#include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EventSetup.h" -#include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h" -#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h" -#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" -#include "DataFormats/HcalDetId/interface/HcalSubdetector.h" - -#include "DataFormats/DetId/interface/DetId.h" -#include "DataFormats/EcalDetId/interface/EEDetId.h" -#include "DataFormats/EcalDetId/interface/EBDetId.h" - -#include "FWCore/MessageLogger/interface/MessageLogger.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include "FWCore/Framework/interface/EventSetup.h" - -#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" -#include "Geometry/CaloGeometry/interface/CaloGeometry.h" -#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h" -#include "Geometry/CaloGeometry/interface/IdealZPrism.h" -#include "Geometry/Records/interface/CaloGeometryRecord.h" -#include "Geometry/Records/interface/CaloGeometryRecord.h" -#include "Geometry/Records/interface/HcalRecNumberingRecord.h" - -#include "RecoCaloTools/Navigation/interface/CaloNavigator.h" - -#include "Geometry/CaloTopology/interface/CaloTowerTopology.h" -#include "RecoCaloTools/Navigation/interface/CaloTowerNavigator.h" - -#include "DataFormats/METReco/interface/HcalCaloFlagLabels.h" -#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h" -#include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h" - -using namespace std; -using namespace edm; - -PFCTRecHitProducer::PFCTRecHitProducer(const edm::ParameterSet& iConfig) { - thresh_Barrel_ = iConfig.getParameter("thresh_Barrel"); - thresh_Endcap_ = iConfig.getParameter("thresh_Endcap"); - - thresh_HF_ = iConfig.getParameter("thresh_HF"); - navigation_HF_ = iConfig.getParameter("navigation_HF"); - weight_HFem_ = iConfig.getParameter("weight_HFem"); - weight_HFhad_ = iConfig.getParameter("weight_HFhad"); - - HCAL_Calib_ = iConfig.getParameter("HCAL_Calib"); - HF_Calib_ = iConfig.getParameter("HF_Calib"); - HCAL_Calib_29 = iConfig.getParameter("HCAL_Calib_29"); - HF_Calib_29 = iConfig.getParameter("HF_Calib_29"); - - shortFibre_Cut = iConfig.getParameter("ShortFibre_Cut"); - longFibre_Fraction = iConfig.getParameter("LongFibre_Fraction"); - - longFibre_Cut = iConfig.getParameter("LongFibre_Cut"); - shortFibre_Fraction = iConfig.getParameter("ShortFibre_Fraction"); - - applyLongShortDPG_ = iConfig.getParameter("ApplyLongShortDPG"); - - longShortFibre_Cut = iConfig.getParameter("LongShortFibre_Cut"); - minShortTiming_Cut = iConfig.getParameter("MinShortTiming_Cut"); - maxShortTiming_Cut = iConfig.getParameter("MaxShortTiming_Cut"); - minLongTiming_Cut = iConfig.getParameter("MinLongTiming_Cut"); - maxLongTiming_Cut = iConfig.getParameter("MaxLongTiming_Cut"); - - applyTimeDPG_ = iConfig.getParameter("ApplyTimeDPG"); - applyPulseDPG_ = iConfig.getParameter("ApplyPulseDPG"); - HcalMaxAllowedHFLongShortSev_ = iConfig.getParameter("HcalMaxAllowedHFLongShortSev"); - HcalMaxAllowedHFDigiTimeSev_ = iConfig.getParameter("HcalMaxAllowedHFDigiTimeSev"); - HcalMaxAllowedHFInTimeWindowSev_ = iConfig.getParameter("HcalMaxAllowedHFInTimeWindowSev"); - HcalMaxAllowedChannelStatusSev_ = iConfig.getParameter("HcalMaxAllowedChannelStatusSev"); - - ECAL_Compensate_ = iConfig.getParameter("ECAL_Compensate"); - ECAL_Threshold_ = iConfig.getParameter("ECAL_Threshold"); - ECAL_Compensation_ = iConfig.getParameter("ECAL_Compensation"); - ECAL_Dead_Code_ = iConfig.getParameter("ECAL_Dead_Code"); - - EM_Depth_ = iConfig.getParameter("EM_Depth"); - HAD_Depth_ = iConfig.getParameter("HAD_Depth"); - - //Get integer values of individual HCAL HF flags - hcalHFLongShortFlagValue_ = 1 << HcalCaloFlagLabels::HFLongShort; - hcalHFDigiTimeFlagValue_ = 1 << HcalCaloFlagLabels::HFDigiTime; - hcalHFInTimeWindowFlagValue_ = 1 << HcalCaloFlagLabels::HFInTimeWindow; - - hcalToken_ = consumes(iConfig.getParameter("hcalRecHitsHBHE")); - hfToken_ = consumes(iConfig.getParameter("hcalRecHitsHF")); - towersToken_ = consumes(iConfig.getParameter("caloTowers")); - - edm::ParameterSet navSet = iConfig.getParameter("navigator"); - navigator_ = PFRecHitNavigationFactory::get()->create(navSet.getParameter("name"), navSet); - - //--ab - produces(); - produces("Cleaned"); - produces("HFHAD").setBranchAlias("HFHADRecHits"); - produces("HFEM").setBranchAlias("HFEMRecHits"); - //--ab -} - -void PFCTRecHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { - navigator_->beginEvent(iSetup); - - edm::ESHandle geoHandle; - iSetup.get().get(geoHandle); - - // get the hcalBarrel geometry - const CaloSubdetectorGeometry* hcalBarrelGeometry = geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalBarrel); - - // get the endcap geometry - const CaloSubdetectorGeometry* hcalEndcapGeometry = geoHandle->getSubdetectorGeometry(DetId::Hcal, HcalEndcap); - - // Get Hcal Severity Level Computer, so that the severity of each rechit flag/status may be determined - edm::ESHandle hcalSevLvlComputerHndl; - iSetup.get().get(hcalSevLvlComputerHndl); - const HcalSeverityLevelComputer* hcalSevLvlComputer = hcalSevLvlComputerHndl.product(); - - // Get Hcal Topology - edm::ESHandle hcalTopology; - iSetup.get().get(hcalTopology); - theHcalTopology = hcalTopology.product(); - - auto rechits = std::make_unique>(); - auto rechitsCleaned = std::make_unique>(); - auto HFHADRecHits = std::make_unique>(); - auto HFEMRecHits = std::make_unique>(); - - edm::Handle caloTowers; - iEvent.getByToken(towersToken_, caloTowers); - edm::Handle hfHandle; - iEvent.getByToken(hfToken_, hfHandle); - - edm::Handle hbheHandle; - iEvent.getByToken(hcalToken_, hbheHandle); - // create rechits - typedef CaloTowerCollection::const_iterator ICT; - - for (auto const& ct : *caloTowers) { - // get the hadronic energy. - - // Mike: Just ask for the Hadronic part only now! - // Patrick : ARGH ! While this is ok for the HCAL, this is - // just wrong for the HF (in which em/had are artificially - // separated. - double energy = ct.hadEnergy(); - //Auguste: Photons in HF have no hadEnergy in fastsim: -> all RecHit collections are empty with photons. - double energyEM = ct.emEnergy(); // For HF ! - //so test the total energy to deal with the photons in HF: - if ((energy + energyEM) < 1e-9) - continue; - - //get the constituents of the tower - const std::vector& hits = ct.constituents(); - const std::vector& allConstituents = theTowerConstituentsMap->constituentsOf(ct.id()); - HcalDetId detid; - bool foundHCALConstituent = false; - //Loop on the calotower constituents and search for HCAL - double dead = 0.; - double alive = 0.; - for (auto hit : hits) { - if (hit.det() == DetId::Hcal) { - foundHCALConstituent = true; - detid = hit; - if (theHcalTopology->getMergePositionFlag() && detid.subdet() == HcalEndcap) { - detid = theHcalTopology->idFront(detid); - } - // An HCAL tower was found: Look for dead ECAL channels in the same CaloTower. - if (ECAL_Compensate_ && energy > ECAL_Threshold_) { - for (auto allConstituent : allConstituents) { - if (allConstituent.det() == DetId::Ecal) { - alive += 1.; - auto chIt = theEcalChStatus->find(allConstituent); - unsigned int dbStatus = chIt != theEcalChStatus->end() ? chIt->getStatusCode() : 0; - if (dbStatus > ECAL_Dead_Code_) - dead += 1.; - } - } - } - // Protection: tower 29 in HF is merged with tower 30. - // Just take the position of tower 30 in that case. - if (detid.subdet() == HcalForward && abs(detid.ieta()) == 29) - continue; - break; - } - } - - // In case of dead ECAL channel, rescale the HCAL energy... - double rescaleFactor = alive > 0. ? 1. + ECAL_Compensation_ * dead / alive : 1.; - - reco::PFRecHit* pfrh = nullptr; - reco::PFRecHit* pfrhCleaned = nullptr; - //---ab: need 2 rechits for the HF: - reco::PFRecHit* pfrhHFEM = nullptr; - reco::PFRecHit* pfrhHFHAD = nullptr; - reco::PFRecHit* pfrhHFEMCleaned = nullptr; - reco::PFRecHit* pfrhHFHADCleaned = nullptr; - reco::PFRecHit* pfrhHFEMCleaned29 = nullptr; - reco::PFRecHit* pfrhHFHADCleaned29 = nullptr; - - if (foundHCALConstituent) { - // std::cout << ", new Energy = " << energy << std::endl; - switch (detid.subdet()) { - case HcalBarrel: { - if (energy < thresh_Barrel_) - continue; - if (rescaleFactor > 1.) { - pfrhCleaned = createHcalRecHit(detid, energy, PFLayer::HCAL_BARREL1, hcalBarrelGeometry, ct.id()); - pfrhCleaned->setTime(rescaleFactor); - energy *= rescaleFactor; - } - pfrh = createHcalRecHit(detid, energy, PFLayer::HCAL_BARREL1, hcalBarrelGeometry, ct.id()); - pfrh->setTime(rescaleFactor); - } break; - case HcalEndcap: { - if (energy < thresh_Endcap_) - continue; - // Apply tower 29 calibration - if (HCAL_Calib_ && abs(detid.ieta()) == 29) - energy *= HCAL_Calib_29; - if (rescaleFactor > 1.) { - pfrhCleaned = createHcalRecHit(detid, energy, PFLayer::HCAL_ENDCAP, hcalEndcapGeometry, ct.id()); - pfrhCleaned->setTime(rescaleFactor); - energy *= rescaleFactor; - } - pfrh = createHcalRecHit(detid, energy, PFLayer::HCAL_ENDCAP, hcalEndcapGeometry, ct.id()); - pfrh->setTime(rescaleFactor); - } break; - case HcalOuter: { - } break; - case HcalForward: { - //---ab: 2 rechits for HF: - //double energyemHF = weight_HFem_*ct.emEnergy(); - //double energyhadHF = weight_HFhad_*ct.hadEnergy(); - double energyemHF = weight_HFem_ * energyEM; - double energyhadHF = weight_HFhad_ * energy; - // Some energy in the tower ! - if ((energyemHF + energyhadHF) < thresh_HF_) - continue; - - // Some cleaning in the HF - double longFibre = energyemHF + energyhadHF / 2.; - double shortFibre = energyhadHF / 2.; - int ieta = detid.ieta(); - int iphi = detid.iphi(); - HcalDetId theLongDetId(HcalForward, ieta, iphi, 1); - HcalDetId theShortDetId(HcalForward, ieta, iphi, 2); - typedef HFRecHitCollection::const_iterator iHF; - auto theLongHit = hfHandle->find(theLongDetId); - auto theShortHit = hfHandle->find(theShortDetId); - // - double theLongHitEnergy = 0.; - double theShortHitEnergy = 0.; - bool flagShortDPG = false; - bool flagLongDPG = false; - bool flagShortTimeDPG = false; - bool flagLongTimeDPG = false; - bool flagShortPulseDPG = false; - bool flagLongPulseDPG = false; - // - if (theLongHit != hfHandle->end()) { - int theLongFlag = theLongHit->flags(); - theLongHitEnergy = theLongHit->energy(); - flagLongDPG = applyLongShortDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theLongDetId, theLongFlag & hcalHFLongShortFlagValue_, 0) > - HcalMaxAllowedHFLongShortSev_); - flagLongTimeDPG = applyTimeDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theLongDetId, theLongFlag & hcalHFInTimeWindowFlagValue_, 0) > - HcalMaxAllowedHFInTimeWindowSev_); - flagLongPulseDPG = applyPulseDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theLongDetId, theLongFlag & hcalHFDigiTimeFlagValue_, 0) > - HcalMaxAllowedHFDigiTimeSev_); - - //flagLongDPG = applyLongShortDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFLongShort)==1; - //flagLongTimeDPG = applyTimeDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1; - //flagLongPulseDPG = applyPulseDPG_ && theLongHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1; - } - // - if (theShortHit != hfHandle->end()) { - int theShortFlag = theShortHit->flags(); - theShortHitEnergy = theShortHit->energy(); - flagShortDPG = applyLongShortDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theShortDetId, theShortFlag & hcalHFLongShortFlagValue_, 0) > - HcalMaxAllowedHFLongShortSev_); - flagShortTimeDPG = applyTimeDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theShortDetId, theShortFlag & hcalHFInTimeWindowFlagValue_, 0) > - HcalMaxAllowedHFInTimeWindowSev_); - flagShortPulseDPG = applyPulseDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theShortDetId, theShortFlag & hcalHFDigiTimeFlagValue_, 0) > - HcalMaxAllowedHFDigiTimeSev_); - //flagShortDPG = applyLongShortDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFLongShort)==1; - //flagShortTimeDPG = applyTimeDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1; - //flagShortPulseDPG = applyPulseDPG_ && theShortHit->flagField(HcalCaloFlagLabels::HFDigiTime)==1; - } - - // Then check the timing in short and long fibres in all other towers. - if (theShortHitEnergy > longShortFibre_Cut && - (theShortHit->time() < minShortTiming_Cut || theShortHit->time() > maxShortTiming_Cut || - flagShortTimeDPG || flagShortPulseDPG)) { - // rescaleFactor = 0. ; - pfrhHFHADCleaned = createHcalRecHit(detid, theShortHitEnergy, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - pfrhHFHADCleaned->setTime(theShortHit->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Time = " << theShortHit->time() << " " - << ". The short and long flags : " - << flagShortDPG << " " << flagLongDPG << " " - << flagShortTimeDPG << " " << flagLongTimeDPG << " " - << flagShortPulseDPG << " " << flagLongPulseDPG << " " - << ". Short fibres were cleaned." << std::endl; - */ - shortFibre -= theShortHitEnergy; - theShortHitEnergy = 0.; - } - - if (theLongHitEnergy > longShortFibre_Cut && - (theLongHit->time() < minLongTiming_Cut || theLongHit->time() > maxLongTiming_Cut || flagLongTimeDPG || - flagLongPulseDPG)) { - //rescaleFactor = 0. ; - pfrhHFEMCleaned = createHcalRecHit(detid, theLongHitEnergy, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFEMCleaned->setTime(theLongHit->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Time = " << theLongHit->time() << " " - << ". The short and long flags : " - << flagShortDPG << " " << flagLongDPG << " " - << flagShortTimeDPG << " " << flagLongTimeDPG << " " - << flagShortPulseDPG << " " << flagLongPulseDPG << " " - << ". Long fibres were cleaned." << std::endl; - */ - longFibre -= theLongHitEnergy; - theLongHitEnergy = 0.; - } - - // Some energy must be in the long fibres is there is some energy in the short fibres ! - if (theShortHitEnergy > shortFibre_Cut && - (theLongHitEnergy / theShortHitEnergy < longFibre_Fraction || flagShortDPG)) { - // Check if the long-fibre hit was not cleaned already (because hot) - // In this case don't apply the cleaning - const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId); - unsigned theStatusValue = theStatus->getValue(); - int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue); - // The channel is killed - /// if ( !theStatusValue ) - if (theSeverityLevel <= HcalMaxAllowedChannelStatusSev_) { - // rescaleFactor = 0. ; - pfrhHFHADCleaned = - createHcalRecHit(detid, theShortHitEnergy, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - pfrhHFHADCleaned->setTime(theShortHit->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Time = " << theShortHit->time() << " " - << ". The status value is " << theStatusValue - << ". The short and long flags : " - << flagShortDPG << " " << flagLongDPG << " " - << flagShortTimeDPG << " " << flagLongTimeDPG << " " - << flagShortPulseDPG << " " << flagLongPulseDPG << " " - << ". Short fibres were cleaned." << std::endl; - */ - shortFibre -= theShortHitEnergy; - theShortHitEnergy = 0.; - } - } - - if (theLongHitEnergy > longFibre_Cut && - (theShortHitEnergy / theLongHitEnergy < shortFibre_Fraction || flagLongDPG)) { - // Check if the long-fibre hit was not cleaned already (because hot) - // In this case don't apply the cleaning - const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId); - unsigned theStatusValue = theStatus->getValue(); - - int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue); - // The channel is killed - /// if ( !theStatusValue ) - if (theSeverityLevel <= HcalMaxAllowedChannelStatusSev_) { - //rescaleFactor = 0. ; - pfrhHFEMCleaned = createHcalRecHit(detid, theLongHitEnergy, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFEMCleaned->setTime(theLongHit->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". The status value is " << theStatusValue - << ". Time = " << theLongHit->time() << " " - << ". The short and long flags : " - << flagShortDPG << " " << flagLongDPG << " " - << flagShortTimeDPG << " " << flagLongTimeDPG << " " - << flagShortPulseDPG << " " << flagLongPulseDPG << " " - << ". Long fibres were cleaned." << std::endl; - */ - longFibre -= theLongHitEnergy; - theLongHitEnergy = 0.; - } - } - - // Special treatment for tower 29 - // A tower with energy only at ieta = +/- 29 is not physical -> Clean - if (abs(ieta) == 29) { - // rescaleFactor = 0. ; - // Clean long fibres - if (theLongHitEnergy > shortFibre_Cut / 2.) { - pfrhHFEMCleaned29 = - createHcalRecHit(detid, theLongHitEnergy, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFEMCleaned29->setTime(theLongHit->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Long fibres were cleaned." << std::endl; - */ - longFibre -= theLongHitEnergy; - theLongHitEnergy = 0.; - } - // Clean short fibres - if (theShortHitEnergy > shortFibre_Cut / 2.) { - pfrhHFHADCleaned29 = - createHcalRecHit(detid, theShortHitEnergy, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - pfrhHFHADCleaned29->setTime(theShortHit->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Short fibres were cleaned." << std::endl; - */ - shortFibre -= theShortHitEnergy; - theShortHitEnergy = 0.; - } - } - // Check the timing of the long and short fibre rechits - - // First, check the timing of long and short fibre in eta = 29 if tower 30. - else if (abs(ieta) == 30) { - int ieta29 = ieta > 0 ? 29 : -29; - HcalDetId theLongDetId29(HcalForward, ieta29, iphi, 1); - HcalDetId theShortDetId29(HcalForward, ieta29, iphi, 2); - auto theLongHit29 = hfHandle->find(theLongDetId29); - auto theShortHit29 = hfHandle->find(theShortDetId29); - // - double theLongHitEnergy29 = 0.; - double theShortHitEnergy29 = 0.; - bool flagShortDPG29 = false; - bool flagLongDPG29 = false; - bool flagShortTimeDPG29 = false; - bool flagLongTimeDPG29 = false; - bool flagShortPulseDPG29 = false; - bool flagLongPulseDPG29 = false; - // - if (theLongHit29 != hfHandle->end()) { - int theLongFlag29 = theLongHit29->flags(); - theLongHitEnergy29 = theLongHit29->energy(); - flagLongDPG29 = applyLongShortDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theLongDetId29, theLongFlag29 & hcalHFLongShortFlagValue_, 0) > - HcalMaxAllowedHFLongShortSev_); - flagLongTimeDPG29 = - applyTimeDPG_ && (hcalSevLvlComputer->getSeverityLevel(theLongDetId29, - theLongFlag29 & hcalHFInTimeWindowFlagValue_, - 0) > HcalMaxAllowedHFInTimeWindowSev_); - flagLongPulseDPG29 = applyPulseDPG_ && (hcalSevLvlComputer->getSeverityLevel( - theLongDetId29, theLongFlag29 & hcalHFDigiTimeFlagValue_, 0) > - HcalMaxAllowedHFDigiTimeSev_); - - //flagLongDPG29 = applyLongShortDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1; - //flagLongTimeDPG29 = applyTimeDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1; - //flagLongPulseDPG29 = applyPulseDPG_ && theLongHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1; - } - // - if (theShortHit29 != hfHandle->end()) { - int theShortFlag29 = theShortHit29->flags(); - theShortHitEnergy29 = theShortHit29->energy(); - flagShortDPG29 = - applyLongShortDPG_ && - (hcalSevLvlComputer->getSeverityLevel( - theShortDetId29, theShortFlag29 & hcalHFLongShortFlagValue_, 0) > HcalMaxAllowedHFLongShortSev_); - flagShortTimeDPG29 = - applyTimeDPG_ && (hcalSevLvlComputer->getSeverityLevel(theShortDetId29, - theShortFlag29 & hcalHFInTimeWindowFlagValue_, - 0) > HcalMaxAllowedHFInTimeWindowSev_); - flagLongPulseDPG29 = - applyPulseDPG_ && - (hcalSevLvlComputer->getSeverityLevel(theShortDetId29, theShortFlag29 & hcalHFDigiTimeFlagValue_, 0) > - HcalMaxAllowedHFDigiTimeSev_); - - //flagShortDPG29 = applyLongShortDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFLongShort)==1; - //flagShortTimeDPG29 = applyTimeDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFInTimeWindow)==1; - //flagShortPulseDPG29 = applyPulseDPG_ && theShortHit29->flagField(HcalCaloFlagLabels::HFDigiTime)==1; - } - - if (theLongHitEnergy29 > longShortFibre_Cut && - (theLongHit29->time() < minLongTiming_Cut || theLongHit29->time() > maxLongTiming_Cut || - flagLongTimeDPG29 || flagLongPulseDPG29)) { - //rescaleFactor = 0. ; - pfrhHFEMCleaned29 = - createHcalRecHit(detid, theLongHitEnergy29, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFEMCleaned29->setTime(theLongHit29->time()); - /* - std::cout << "ieta/iphi = " << ieta29 << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " - << ". Time = " << theLongHit29->time() << " " - << ". The short and long flags : " - << flagShortDPG29 << " " << flagLongDPG29 << " " - << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " " - << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " " - << ". Long fibres were cleaned." << std::endl; - */ - longFibre -= theLongHitEnergy29; - theLongHitEnergy29 = 0; - } - - if (theShortHitEnergy29 > longShortFibre_Cut && - (theShortHit29->time() < minShortTiming_Cut || theShortHit29->time() > maxShortTiming_Cut || - flagShortTimeDPG29 || flagShortPulseDPG29)) { - //rescaleFactor = 0. ; - pfrhHFHADCleaned29 = - createHcalRecHit(detid, theShortHitEnergy29, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - pfrhHFHADCleaned29->setTime(theShortHit29->time()); - /* - std::cout << "ieta/iphi = " << ieta29 << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " - << ". Time = " << theShortHit29->time() << " " - << ". The short and long flags : " - << flagShortDPG29 << " " << flagLongDPG29 << " " - << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " " - << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " " - << ". Short fibres were cleaned." << std::endl; - */ - shortFibre -= theShortHitEnergy29; - theShortHitEnergy29 = 0; - } - - // Some energy must be in the long fibres is there is some energy in the short fibres ! - if (theShortHitEnergy29 > shortFibre_Cut && - (theLongHitEnergy29 / theShortHitEnergy29 < 2. * longFibre_Fraction || flagShortDPG29)) { - // Check if the long-fibre hit was not cleaned already (because hot) - // In this case don't apply the cleaning - const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theLongDetId29); - unsigned theStatusValue = theStatus->getValue(); - - int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue); - // The channel is killed - /// if ( !theStatusValue ) - if (theSeverityLevel <= HcalMaxAllowedChannelStatusSev_) { - //rescaleFactor = 0. ; - pfrhHFHADCleaned29 = - createHcalRecHit(detid, theShortHitEnergy29, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - pfrhHFHADCleaned29->setTime(theShortHit29->time()); - /* - std::cout << "ieta/iphi = " << ieta29 << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " - << ". Time = " << theShortHit29->time() << " " - << ". The status value is " << theStatusValue - << ". The short and long flags : " - << flagShortDPG29 << " " << flagLongDPG29 << " " - << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " " - << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " " - << ". Short fibres were cleaned." << std::endl; - */ - shortFibre -= theShortHitEnergy29; - theShortHitEnergy29 = 0.; - } - } - - // Some energy must be in the short fibres is there is some energy in the long fibres ! - if (theLongHitEnergy29 > longFibre_Cut && - (theShortHitEnergy29 / theLongHitEnergy29 < shortFibre_Fraction || flagLongDPG29)) { - // Check if the long-fibre hit was not cleaned already (because hot) - // In this case don't apply the cleaning - const HcalChannelStatus* theStatus = theHcalChStatus->getValues(theShortDetId29); - unsigned theStatusValue = theStatus->getValue(); - int theSeverityLevel = hcalSevLvlComputer->getSeverityLevel(detid, 0, theStatusValue); - // The channel is killed - /// if ( !theStatusValue ) - if (theSeverityLevel <= HcalMaxAllowedChannelStatusSev_) { - //rescaleFactor = 0. ; - pfrhHFEMCleaned29 = - createHcalRecHit(detid, theLongHitEnergy29, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFEMCleaned29->setTime(theLongHit29->time()); - /* - std::cout << "ieta/iphi = " << ieta29 << " " << iphi - << ", Energy em/had/long/short = " - << energyemHF << " " << energyhadHF << " " - << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " - << ". The status value is " << theStatusValue - << ". Time = " << theLongHit29->time() << " " - << ". The short and long flags : " - << flagShortDPG29 << " " << flagLongDPG29 << " " - << flagShortTimeDPG29 << " " << flagLongTimeDPG29 << " " - << flagShortPulseDPG29 << " " << flagLongPulseDPG29 << " " - << ". Long fibres were cleaned." << std::endl; - */ - longFibre -= theLongHitEnergy29; - theLongHitEnergy29 = 0.; - } - } - - // Check that the energy in tower 29 is smaller than in tower 30 - // First in long fibres - if (theLongHitEnergy29 > std::max(theLongHitEnergy, shortFibre_Cut / 2)) { - pfrhHFEMCleaned29 = - createHcalRecHit(detid, theLongHitEnergy29, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFEMCleaned29->setTime(theLongHit29->time()); - /* - std::cout << "ieta/iphi = " << ieta29 << " " << iphi - << ", Energy L29/S29/L30/S30 = " - << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Long fibres were cleaned." << std::endl; - */ - longFibre -= theLongHitEnergy29; - theLongHitEnergy29 = 0.; - } - // Second in short fibres - if (theShortHitEnergy29 > std::max(theShortHitEnergy, shortFibre_Cut / 2.)) { - pfrhHFHADCleaned29 = - createHcalRecHit(detid, theShortHitEnergy29, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - pfrhHFHADCleaned29->setTime(theShortHit29->time()); - /* - std::cout << "ieta/iphi = " << ieta << " " << iphi - << ", Energy L29/S29/L30/S30 = " - << theLongHitEnergy29 << " " << theShortHitEnergy29 << " " - << theLongHitEnergy << " " << theShortHitEnergy << " " - << ". Short fibres were cleaned." << std::endl; - */ - shortFibre -= theShortHitEnergy29; - theShortHitEnergy29 = 0.; - } - } - - // Determine EM and HAD after cleaning of short and long fibres - energyhadHF = 2. * shortFibre; - energyemHF = longFibre - shortFibre; - - // The EM energy might be negative, as it amounts to Long - Short - // In that case, put the EM "energy" in the HAD energy - // Just to avoid systematic positive bias due to "Short" high fluctuations - if (energyemHF < thresh_HF_) { - energyhadHF += energyemHF; - energyemHF = 0.; - } - - // Apply HCAL calibration factors flor towers close to 29, if requested - if (HF_Calib_ && abs(detid.ieta()) <= 32) { - energyhadHF *= HF_Calib_29; - energyemHF *= HF_Calib_29; - } - - // Create an EM and a HAD rechit if above threshold. - if (energyemHF > thresh_HF_ || energyhadHF > thresh_HF_) { - pfrhHFEM = createHcalRecHit(detid, energyemHF, PFLayer::HF_EM, hcalEndcapGeometry, ct.id()); - pfrhHFHAD = createHcalRecHit(detid, energyhadHF, PFLayer::HF_HAD, hcalEndcapGeometry, ct.id()); - } - - } break; - default: - LogError("PFCTRecHitProducerHCAL") << "CaloTower constituent: unknown layer : " << detid.subdet() << endl; - } - - if (pfrh) { - rechits->push_back(*pfrh); - delete pfrh; - } - if (pfrhCleaned) { - rechitsCleaned->push_back(*pfrhCleaned); - delete pfrhCleaned; - } - if (pfrhHFEM) { - HFEMRecHits->push_back(*pfrhHFEM); - delete pfrhHFEM; - } - if (pfrhHFHAD) { - HFHADRecHits->push_back(*pfrhHFHAD); - delete pfrhHFHAD; - } - if (pfrhHFEMCleaned) { - rechitsCleaned->push_back(*pfrhHFEMCleaned); - delete pfrhHFEMCleaned; - } - if (pfrhHFHADCleaned) { - rechitsCleaned->push_back(*pfrhHFHADCleaned); - delete pfrhHFHADCleaned; - } - if (pfrhHFEMCleaned29) { - rechitsCleaned->push_back(*pfrhHFEMCleaned29); - delete pfrhHFEMCleaned29; - } - if (pfrhHFHADCleaned29) { - rechitsCleaned->push_back(*pfrhHFHADCleaned29); - delete pfrhHFHADCleaned29; - } - } - } - - //ok now do navigation - //create a refprod here - - edm::RefProd refProd = iEvent.getRefBeforePut(); - - for (unsigned int i = 0; i < rechits->size(); ++i) { - navigator_->associateNeighbours(rechits->at(i), rechits, refProd); - } - - if (navigation_HF_) { - edm::RefProd refProdEM = iEvent.getRefBeforePut("HFEM"); - - for (unsigned int i = 0; i < HFEMRecHits->size(); ++i) { - navigator_->associateNeighbours(HFEMRecHits->at(i), HFEMRecHits, refProdEM); - } - - edm::RefProd refProdHAD = iEvent.getRefBeforePut("HFHAD"); - - for (unsigned int i = 0; i < HFHADRecHits->size(); ++i) { - navigator_->associateNeighbours(HFHADRecHits->at(i), HFHADRecHits, refProdHAD); - } - } - - iEvent.put(std::move(rechits), ""); - iEvent.put(std::move(rechitsCleaned), "Cleaned"); - iEvent.put(std::move(HFEMRecHits), "HFEM"); - iEvent.put(std::move(HFHADRecHits), "HFHAD"); -} - -PFCTRecHitProducer::~PFCTRecHitProducer() = default; - -// ------------ method called once each job just before starting event loop ------------ -void PFCTRecHitProducer::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const EventSetup& es) { - // Get cleaned channels in the HCAL and HF - // HCAL channel status map **************************************** - edm::ESHandle hcalChStatus; - es.get().get("withTopo", hcalChStatus); - theHcalChStatus = hcalChStatus.product(); - - // Retrieve the good/bad ECAL channels from the DB - edm::ESHandle ecalChStatus; - es.get().get(ecalChStatus); - theEcalChStatus = ecalChStatus.product(); - - edm::ESHandle cttopo; - es.get().get(cttopo); - theTowerConstituentsMap = cttopo.product(); -} - -reco::PFRecHit* PFCTRecHitProducer::createHcalRecHit(const DetId& detid, - double energy, - PFLayer::Layer layer, - const CaloSubdetectorGeometry* geom, - const CaloTowerDetId& newDetId) { - std::shared_ptr thisCell = geom->getGeometry(detid); - if (!thisCell) { - edm::LogError("PFRecHitProducerHCAL") - << "warning detid " << detid.rawId() << " not found in layer " << layer << endl; - return nullptr; - } - - switch (layer) { - case PFLayer::HF_EM: - case PFLayer::HF_HAD: { - auto zp = dynamic_cast(thisCell.get()); - assert(zp); - thisCell = zp->forPF(); - } break; - default: - break; - } - - auto rh = new reco::PFRecHit(thisCell, newDetId.rawId(), layer, energy); - - return rh; -} diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.h b/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.h deleted file mode 100644 index 108a3c58187a0..0000000000000 --- a/RecoParticleFlow/PFClusterProducer/plugins/PFCTRecHitProducer.h +++ /dev/null @@ -1,125 +0,0 @@ -#ifndef RecoParticleFlow_PFClusterProducer_PFCTRecHitProducer_h_ -#define RecoParticleFlow_PFClusterProducer_PFCTRecHitProducer_h_ - -// system include files -#include -#include - -// user include files -#include "FWCore/Framework/interface/Frameworkfwd.h" -#include "FWCore/Framework/interface/stream/EDProducer.h" - -#include "FWCore/Framework/interface/Event.h" -#include "FWCore/Framework/interface/MakerMacros.h" - -#include "FWCore/ParameterSet/interface/ParameterSet.h" - -#include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h" - -#include "CondFormats/HcalObjects/interface/HcalPFCorrs.h" -#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h" -#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h" -#include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h" -#include "Geometry/CaloTopology/interface/HcalTopology.h" - -#include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h" -#include "DataFormats/HcalRecHit/interface/HFRecHit.h" -#include "DataFormats/HcalRecHit/interface/HBHERecHit.h" -#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h" -#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" - -class CaloSubdetectorTopology; -class CaloSubdetectorGeometry; -class DetId; - -class dso_hidden PFCTRecHitProducer final : public edm::stream::EDProducer<> { -public: - explicit PFCTRecHitProducer(const edm::ParameterSet&); - ~PFCTRecHitProducer() override; - - void beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& es) override; - - void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override; - - reco::PFRecHit* createHcalRecHit(const DetId& detid, - double energy, - PFLayer::Layer layer, - const CaloSubdetectorGeometry* geom, - const CaloTowerDetId& newDetId); - -protected: - double thresh_Barrel_; - double thresh_Endcap_; - const HcalChannelQuality* theHcalChStatus; - const EcalChannelStatus* theEcalChStatus; - const CaloTowerConstituentsMap* theTowerConstituentsMap; - const HcalTopology* theHcalTopology; - - // ----------access to event data - edm::EDGetTokenT hcalToken_; - edm::EDGetTokenT hfToken_; - edm::EDGetTokenT towersToken_; - - /// threshold for HF - double thresh_HF_; - - // Navigation in HF: False = no real clustering in HF; True = do clustering - bool navigation_HF_; - double weight_HFem_; - double weight_HFhad_; - - // Apply HCAL DPG rechit calibration - bool HCAL_Calib_; - bool HF_Calib_; - float HCAL_Calib_29; - float HF_Calib_29; - - // Don't allow large energy in short fibres if there is no energy in long fibres - double shortFibre_Cut; - double longFibre_Fraction; - - // Don't allow large energy in long fibres if there is no energy in short fibres - double longFibre_Cut; - double shortFibre_Fraction; - - // Also apply HCAL DPG cleaning - bool applyLongShortDPG_; - - // Don't allow too large timing excursion if energy in long/short fibres is large enough - double longShortFibre_Cut; - double minShortTiming_Cut; - double maxShortTiming_Cut; - double minLongTiming_Cut; - double maxLongTiming_Cut; - - bool applyTimeDPG_; - bool applyPulseDPG_; - int HcalMaxAllowedHFLongShortSev_; - int HcalMaxAllowedHFDigiTimeSev_; - int HcalMaxAllowedHFInTimeWindowSev_; - int HcalMaxAllowedChannelStatusSev_; - - int hcalHFLongShortFlagValue_; - int hcalHFDigiTimeFlagValue_; - int hcalHFInTimeWindowFlagValue_; - - // Compensate for dead ECAL channels - bool ECAL_Compensate_; - double ECAL_Threshold_; - double ECAL_Compensation_; - unsigned int ECAL_Dead_Code_; - - // Depth correction for EM and HAD rechits in the HF - double EM_Depth_; - double HAD_Depth_; - - int m_maxDepthHB; - int m_maxDepthHE; - - std::unique_ptr navigator_; -}; - -#include "FWCore/Framework/interface/MakerMacros.h" -DEFINE_FWK_MODULE(PFCTRecHitProducer); - -#endif diff --git a/RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducer.cc b/RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducer.cc index c6be321ec203a..8ba3648b26b34 100644 --- a/RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducer.cc +++ b/RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducer.cc @@ -21,7 +21,6 @@ PFRecHitProducer::PFRecHitProducer(const edm::ParameterSet& iConfig) { } edm::ParameterSet navSet = iConfig.getParameter("navigator"); - navigator_ = PFRecHitNavigationFactory::get()->create(navSet.getParameter("name"), navSet); } @@ -35,6 +34,7 @@ void PFRecHitProducer::beginLuminosityBlock(edm::LuminosityBlock const& iLumi, c for (const auto& creator : creators_) { creator->init(iSetup); } + navigator_->init(iSetup); } void PFRecHitProducer::endLuminosityBlock(edm::LuminosityBlock const& iLumi, const edm::EventSetup&) {} @@ -45,8 +45,6 @@ void PFRecHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup auto out = std::make_unique(); auto cleaned = std::make_unique(); - navigator_->beginEvent(iSetup); - out->reserve(localRA1.upper()); cleaned->reserve(localRA2.upper()); diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py index 596a49bd693a7..9a27ef70f53cf 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHBHE_cfi.py @@ -1,17 +1,14 @@ import FWCore.ParameterSet.Config as cms -from RecoParticleFlow.PFClusterProducer.particleFlowCaloResolution_cfi import _timeResolutionHCAL _thresholdsHB = cms.vdouble(0.8, 0.8, 0.8, 0.8) _thresholdsHE = cms.vdouble(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8) _thresholdsHBphase1 = cms.vdouble(0.1, 0.2, 0.3, 0.3) _thresholdsHEphase1 = cms.vdouble(0.1, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2) - particleFlowRecHitHBHE = cms.EDProducer("PFRecHitProducer", navigator = cms.PSet( - name = cms.string("PFRecHitHCALNavigator"), - sigmaCut = cms.double(4.0), - timeResolutionCalc = _timeResolutionHCAL + name = cms.string("PFRecHitHCALDenseIdNavigator"), + hcalEnums = cms.vint32(1,2) ), producers = cms.VPSet( cms.PSet( diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHF_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHF_cfi.py index b1d62366c097c..c9ca84b38c7ca 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHF_cfi.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHF_cfi.py @@ -1,13 +1,10 @@ import FWCore.ParameterSet.Config as cms -#HB HE HO rec hits - particleFlowRecHitHF = cms.EDProducer("PFRecHitProducer", navigator = cms.PSet( - name = cms.string("PFRecHitHCALNavigator"), - barrel = cms.PSet( ), - endcap = cms.PSet( ) + name = cms.string("PFRecHitHCALDenseIdNavigator"), + hcalEnums = cms.vint32(4) ), producers = cms.VPSet( cms.PSet( @@ -37,9 +34,9 @@ detectorEnum = cms.int32(4)) ) ) - + ) - ) + ) ) ) diff --git a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHO_cfi.py b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHO_cfi.py index b7296d1539900..8d8673ae982e8 100644 --- a/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHO_cfi.py +++ b/RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHO_cfi.py @@ -1,7 +1,9 @@ import FWCore.ParameterSet.Config as cms + particleFlowRecHitHO = cms.EDProducer("PFRecHitProducer", navigator = cms.PSet( - name = cms.string("PFRecHitHCALNavigator") + name = cms.string("PFRecHitHCALDenseIdNavigator"), + hcalEnums = cms.vint32(3) ), producers = cms.VPSet( cms.PSet(