From bee9abcda9563bddb553ce1fa923c18740769b26 Mon Sep 17 00:00:00 2001 From: Igor Date: Mon, 18 May 2015 18:33:23 +0200 Subject: [PATCH] Backport of the HCAL negative energy filter code --- CondCore/HcalPlugins/src/plugin.cc | 4 + CondCore/Utilities/src/CondDBFetch.cc | 1 + CondCore/Utilities/src/CondDBImport.cc | 1 + CondCore/Utilities/src/CondFormats.h | 1 + .../interface/HBHENegativeEFilterRcd.h | 26 +++ .../DataRecord/src/HBHENegativeEFilterRcd.cc | 15 ++ .../interface/HBHENegativeEFilter.h | 87 +++++++++ .../HcalObjects/src/HBHENegativeEFilter.cc | 133 +++++++++++++ .../src/T_EventSetup_HBHENegativeEFilter.cc | 4 + CondFormats/HcalObjects/src/classes.h | 4 + CondFormats/HcalObjects/src/classes_def.xml | 2 + CondFormats/HcalObjects/src/headers.h | 1 + .../plugins/HBHENegativeEFilterDBModules.cc | 13 ++ .../test/HBHENegativeEFilterDBReader_cfg.py | 29 +++ .../test/HBHENegativeEFilterDBWriter_cfg.py | 51 +++++ .../HcalRecAlgos/interface/HBHENegativeFlag.h | 54 ++---- .../HcalRecAlgos/src/HBHENegativeFlag.cc | 181 +----------------- .../HcalRecAlgos/src/HBHETimingShapedFlag.cc | 4 + .../python/HcalHitReconstructor_hbhe_cfi.py | 11 +- .../src/HcalHitReconstructor.cc | 74 +++---- .../src/HcalHitReconstructor.h | 2 - 21 files changed, 438 insertions(+), 260 deletions(-) create mode 100644 CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h create mode 100644 CondFormats/DataRecord/src/HBHENegativeEFilterRcd.cc create mode 100644 CondFormats/HcalObjects/interface/HBHENegativeEFilter.h create mode 100644 CondFormats/HcalObjects/src/HBHENegativeEFilter.cc create mode 100644 CondFormats/HcalObjects/src/T_EventSetup_HBHENegativeEFilter.cc create mode 100644 CondTools/Hcal/plugins/HBHENegativeEFilterDBModules.cc create mode 100644 CondTools/Hcal/test/HBHENegativeEFilterDBReader_cfg.py create mode 100644 CondTools/Hcal/test/HBHENegativeEFilterDBWriter_cfg.py diff --git a/CondCore/HcalPlugins/src/plugin.cc b/CondCore/HcalPlugins/src/plugin.cc index 419a3413f12d0..265cfbaf4faab 100644 --- a/CondCore/HcalPlugins/src/plugin.cc +++ b/CondCore/HcalPlugins/src/plugin.cc @@ -22,6 +22,9 @@ #include "CondFormats/DataRecord/interface/HcalInterpolatedPulseCollRcd.h" #include "CondFormats/HcalObjects/interface/HcalInterpolatedPulseColl.h" +#include "CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h" +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" + // #include "CondCore/CondDB/interface/Serialization.h" @@ -66,3 +69,4 @@ REGISTER_PLUGIN(HcalOOTPileupCorrectionRcd,OOTPileupCorrectionColl); REGISTER_PLUGIN(HcalOOTPileupCompatibilityRcd,OOTPileupCorrectionBuffer); REGISTER_PLUGIN(HcalOOTPileupCorrectionMapCollRcd,OOTPileupCorrectionMapColl); REGISTER_PLUGIN(HcalInterpolatedPulseCollRcd,HcalInterpolatedPulseColl); +REGISTER_PLUGIN(HBHENegativeEFilterRcd,HBHENegativeEFilter); diff --git a/CondCore/Utilities/src/CondDBFetch.cc b/CondCore/Utilities/src/CondDBFetch.cc index b9b235e8f0fa8..f4eade5d7983a 100644 --- a/CondCore/Utilities/src/CondDBFetch.cc +++ b/CondCore/Utilities/src/CondDBFetch.cc @@ -133,6 +133,7 @@ namespace cond { FETCH_PAYLOAD_CASE( FileBlob ) FETCH_PAYLOAD_CASE( GBRForest ) FETCH_PAYLOAD_CASE( GBRForestD ) + FETCH_PAYLOAD_CASE( HBHENegativeEFilter ) FETCH_PAYLOAD_CASE( HcalChannelQuality ) FETCH_PAYLOAD_CASE( HcalCholeskyMatrices ) FETCH_PAYLOAD_CASE( HcalElectronicsMap ) diff --git a/CondCore/Utilities/src/CondDBImport.cc b/CondCore/Utilities/src/CondDBImport.cc index 9a77e4d11592c..3b482bbb74e9b 100644 --- a/CondCore/Utilities/src/CondDBImport.cc +++ b/CondCore/Utilities/src/CondDBImport.cc @@ -148,6 +148,7 @@ namespace cond { IMPORT_PAYLOAD_CASE( FileBlob ) IMPORT_PAYLOAD_CASE( GBRForest ) IMPORT_PAYLOAD_CASE( GBRForestD ) + IMPORT_PAYLOAD_CASE( HBHENegativeEFilter ) IMPORT_PAYLOAD_CASE( HcalChannelQuality ) IMPORT_PAYLOAD_CASE( HcalCholeskyMatrices ) IMPORT_PAYLOAD_CASE( HcalDcsValues ) diff --git a/CondCore/Utilities/src/CondFormats.h b/CondCore/Utilities/src/CondFormats.h index 275c331a6bee8..8d9054f31ff80 100644 --- a/CondCore/Utilities/src/CondFormats.h +++ b/CondCore/Utilities/src/CondFormats.h @@ -65,6 +65,7 @@ #include "CondFormats/HcalObjects/interface/OOTPileupCorrectionBuffer.h" #include "CondFormats/HcalObjects/interface/StorableDoubleMap.h" #include "CondFormats/HcalObjects/interface/HcalInterpolatedPulseColl.h" +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h" #include "CondFormats/JetMETObjects/interface/QGLikelihoodObject.h" #include "CondFormats/L1TObjects/interface/L1CaloEcalScale.h" diff --git a/CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h b/CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h new file mode 100644 index 0000000000000..4583f4e9dbe8f --- /dev/null +++ b/CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h @@ -0,0 +1,26 @@ +#ifndef CondFormats_HBHENegativeEFilterRcd_h +#define CondFormats_HBHENegativeEFilterRcd_h +// -*- C++ -*- +// +// Package: CondFormats/DataRecord +// Class : HBHENegativeEFilterRcd +// +/**\class HBHENegativeEFilterRcd HBHENegativeEFilterRcd.h CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h + + Description: record for storing HCAL negative energy filter data + + Usage: + + +*/ +// +// Author: Igor Volobouev +// Created: Fri Mar 20 17:05:13 CDT 2015 +// + +#include "FWCore/Framework/interface/EventSetupRecordImplementation.h" + +class HBHENegativeEFilterRcd : public edm::eventsetup::EventSetupRecordImplementation {}; + +#endif // CondFormats_HBHENegativeEFilterRcd_h + diff --git a/CondFormats/DataRecord/src/HBHENegativeEFilterRcd.cc b/CondFormats/DataRecord/src/HBHENegativeEFilterRcd.cc new file mode 100644 index 0000000000000..db7e711dcf939 --- /dev/null +++ b/CondFormats/DataRecord/src/HBHENegativeEFilterRcd.cc @@ -0,0 +1,15 @@ +// -*- C++ -*- +// +// Package: CondFormats/DataRecord +// Class : HBHENegativeEFilterRcd +// +// Implementation: +// [Notes on implementation] +// +// Author: Igor Volobouev +// Created: Fri Mar 20 17:06:03 CDT 2015 + +#include "CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h" +#include "FWCore/Framework/interface/eventsetuprecord_registration_macro.h" + +EVENTSETUP_RECORD_REG(HBHENegativeEFilterRcd); diff --git a/CondFormats/HcalObjects/interface/HBHENegativeEFilter.h b/CondFormats/HcalObjects/interface/HBHENegativeEFilter.h new file mode 100644 index 0000000000000..f0a3c1bd9f9c9 --- /dev/null +++ b/CondFormats/HcalObjects/interface/HBHENegativeEFilter.h @@ -0,0 +1,87 @@ +#ifndef CondFormats_HcalObjects_HBHENegativeEFilter_h_ +#define CondFormats_HcalObjects_HBHENegativeEFilter_h_ + +#include +#include +#include "FWCore/Utilities/interface/Exception.h" + +#include "boost/cstdint.hpp" +#include "boost/serialization/utility.hpp" +#include "boost/serialization/access.hpp" +#include "boost/serialization/split_member.hpp" + +#include "DataFormats/HcalDetId/interface/HcalDetId.h" +#include "CondFormats/HcalObjects/interface/PiecewiseScalingPolynomial.h" + +class HBHENegativeEFilter +{ +public: + inline HBHENegativeEFilter() : minCharge_(0.), tFirst_(0), tLast_(0) {} + + // If the vector of cuts is empty, the filter will be disabled + HBHENegativeEFilter(const std::vector& a1vec, + const std::vector& a2vec, + const std::vector& iEtaLimits, + const std::vector >& cut, + double minCharge, unsigned firstTimeSlice, + unsigned lastTimeSlice); + + // Does the sequence of time slices pass the filter? + bool checkPassFilter(const HcalDetId& id, + const double* ts, unsigned lenTS) const; + + // Examing various filter data elements + inline const PiecewiseScalingPolynomial& getA1(const HcalDetId& id) const + {return a1v_.at(getEtaIndex(id));} + inline const PiecewiseScalingPolynomial& getA2(const HcalDetId& id) const + {return a2v_.at(getEtaIndex(id));} + inline const std::vector& getEtaLimits() const + {return iEtaLimits_;} + inline const std::vector >& getCut() const + {return cut_;} + inline double getMinCharge() const {return minCharge_;} + inline unsigned getFirstTimeSlice() const {return tFirst_;} + inline unsigned getLastTimeSlice() const {return tLast_;} + inline bool isEnabled() const {return !cut_.empty();} + + // Comparison operators + bool operator==(const HBHENegativeEFilter& r) const; + inline bool operator!=(const HBHENegativeEFilter& r) const + {return !(*this == r);} + +private: + unsigned getEtaIndex(const HcalDetId& id) const; + bool validate() const; + + std::vector a1v_; + std::vector a2v_; + std::vector iEtaLimits_; + std::vector > cut_; + double minCharge_; + uint32_t tFirst_; + uint32_t tLast_; + + friend class boost::serialization::access; + + template + inline void save(Archive & ar, const unsigned /* version */) const + { + if (!validate()) throw cms::Exception( + "In HBHENegativeEFilter::save: invalid data"); + ar & a1v_ & a2v_ & iEtaLimits_ & cut_ & minCharge_ & tFirst_ & tLast_; + } + + template + inline void load(Archive & ar, const unsigned /* version */) + { + ar & a1v_ & a2v_ & iEtaLimits_ & cut_ & minCharge_ & tFirst_ & tLast_; + if (!validate()) throw cms::Exception( + "In HBHENegativeEFilter::load: invalid data"); + } + + BOOST_SERIALIZATION_SPLIT_MEMBER() +}; + +BOOST_CLASS_VERSION(HBHENegativeEFilter, 1) + +#endif // CondFormats_HcalObjects_HBHENegativeEFilter_h_ diff --git a/CondFormats/HcalObjects/src/HBHENegativeEFilter.cc b/CondFormats/HcalObjects/src/HBHENegativeEFilter.cc new file mode 100644 index 0000000000000..e8d079357d5a3 --- /dev/null +++ b/CondFormats/HcalObjects/src/HBHENegativeEFilter.cc @@ -0,0 +1,133 @@ +#include +#include + +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" + +HBHENegativeEFilter::HBHENegativeEFilter( + const std::vector& a1vec, + const std::vector& a2vec, + const std::vector& iEtaLimits, + const std::vector >& cut, + const double minCharge, + const unsigned firstTimeSlice, + const unsigned lastTimeSlice) + : a1v_(a1vec), + a2v_(a2vec), + iEtaLimits_(iEtaLimits), + cut_(cut), + minCharge_(minCharge), + tFirst_(firstTimeSlice), + tLast_(lastTimeSlice) +{ + if (!validate()) throw cms::Exception( + "Invalid HBHENegativeEFilter constructor arguments"); +} + +bool HBHENegativeEFilter::validate() const +{ + if (cut_.empty()) + return true; + + const std::size_t nLimits(iEtaLimits_.size()); + if (nLimits >= static_cast(UINT_MAX - 1U)) + return false; + for (std::size_t i=1; i= static_cast(UINT_MAX - 1U)) + return false; + for (std::size_t i=1; i= minCharge_) + { + // Figure out the cut value for this charge + const std::pair* cut = &cut_[0]; + double cutValue = cut[0].second; + if (sz > 1U) + { + // First point larger than charge + unsigned largerPoint = 0; + for (; cut[largerPoint].first <= chargeInWindow; ++largerPoint) {} + + // Constant extrapolation beyond min and max coords + if (largerPoint >= sz) + cutValue = cut[sz - 1U].second; + else if (largerPoint) + { + const double slope = (cut[largerPoint].second - cut[largerPoint-1U].second)/ + (cut[largerPoint].first - cut[largerPoint-1U].first); + cutValue = cut[largerPoint-1U].second + slope* + (chargeInWindow - cut[largerPoint-1U].first); + } + } + + // Compare the modified time slices with the cut + const unsigned itaIdx = getEtaIndex(id); + const PiecewiseScalingPolynomial& a1(a1v_[itaIdx]); + const PiecewiseScalingPolynomial& a2(a2v_[itaIdx]); + + for (unsigned i=tFirst_; i<=tLast_ && i= cutValue; + } + } + } + return passes; +} diff --git a/CondFormats/HcalObjects/src/T_EventSetup_HBHENegativeEFilter.cc b/CondFormats/HcalObjects/src/T_EventSetup_HBHENegativeEFilter.cc new file mode 100644 index 0000000000000..3f9c4aa4990b4 --- /dev/null +++ b/CondFormats/HcalObjects/src/T_EventSetup_HBHENegativeEFilter.cc @@ -0,0 +1,4 @@ +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" +#include "FWCore/Utilities/interface/typelookup.h" + +TYPELOOKUP_DATA_REG(HBHENegativeEFilter); diff --git a/CondFormats/HcalObjects/src/classes.h b/CondFormats/HcalObjects/src/classes.h index fae17a2305dcd..383f782af1aa7 100644 --- a/CondFormats/HcalObjects/src/classes.h +++ b/CondFormats/HcalObjects/src/classes.h @@ -102,6 +102,10 @@ namespace CondFormats_HcalObjects { std::vector myHcalInterpolatedPulseVec; HBHEChannelGroups myHBHEChannelGroups; HcalInterpolatedPulseColl myHcalInterpolatedPulseColl; + + // HBHE negative energy filter + std::vector myPiecewiseScalingPolynomialVec; + HBHENegativeEFilter myHBHENegativeEFilter; }; } diff --git a/CondFormats/HcalObjects/src/classes_def.xml b/CondFormats/HcalObjects/src/classes_def.xml index 4b5e9ce4501bd..c736cd423ce59 100644 --- a/CondFormats/HcalObjects/src/classes_def.xml +++ b/CondFormats/HcalObjects/src/classes_def.xml @@ -415,4 +415,6 @@ + + diff --git a/CondFormats/HcalObjects/src/headers.h b/CondFormats/HcalObjects/src/headers.h index 8187025a3bd46..043c9f328032d 100644 --- a/CondFormats/HcalObjects/src/headers.h +++ b/CondFormats/HcalObjects/src/headers.h @@ -11,3 +11,4 @@ #include "CondFormats/HcalObjects/interface/HcalInterpolatedPulse.h" #include "CondFormats/HcalObjects/interface/HcalInterpolatedPulseColl.h" #include "CondFormats/HcalObjects/interface/HBHEChannelGroups.h" +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" diff --git a/CondTools/Hcal/plugins/HBHENegativeEFilterDBModules.cc b/CondTools/Hcal/plugins/HBHENegativeEFilterDBModules.cc new file mode 100644 index 0000000000000..fac4b8cae4ce2 --- /dev/null +++ b/CondTools/Hcal/plugins/HBHENegativeEFilterDBModules.cc @@ -0,0 +1,13 @@ +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" +#include "CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h" + +#include "CondTools/Hcal/interface/BoostIODBWriter.h" +#include "CondTools/Hcal/interface/BoostIODBReader.h" + +#include "FWCore/Framework/interface/MakerMacros.h" + +typedef BoostIODBWriter HBHENegativeEFilterDBWriter; +typedef BoostIODBReader HBHENegativeEFilterDBReader; + +DEFINE_FWK_MODULE(HBHENegativeEFilterDBWriter); +DEFINE_FWK_MODULE(HBHENegativeEFilterDBReader); diff --git a/CondTools/Hcal/test/HBHENegativeEFilterDBReader_cfg.py b/CondTools/Hcal/test/HBHENegativeEFilterDBReader_cfg.py new file mode 100644 index 0000000000000..8895e025b6b6d --- /dev/null +++ b/CondTools/Hcal/test/HBHENegativeEFilterDBReader_cfg.py @@ -0,0 +1,29 @@ +database = "sqlite_file:HBHENegativeEFilter_V00_data.db" +tag = "HBHENegativeEFilter_V00_data" +outputfile = "dbread.bbin" + +import FWCore.ParameterSet.Config as cms + +process = cms.Process('HBHENegativeEFilterDBRead') + +process.source = cms.Source('EmptySource') +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1)) + +process.load("CondCore.CondDB.CondDB_cfi") +process.CondDB.connect = database +# process.CondDB.dbFormat = cms.untracked.int32(1) + +process.PoolDBESSource = cms.ESSource("PoolDBESSource", + process.CondDB, + toGet = cms.VPSet(cms.PSet( + record = cms.string("HBHENegativeEFilterRcd"), + tag = cms.string(tag) + )) +) + +process.dumper = cms.EDAnalyzer( + 'HBHENegativeEFilterDBReader', + outputFile = cms.string(outputfile) +) + +process.p = cms.Path(process.dumper) diff --git a/CondTools/Hcal/test/HBHENegativeEFilterDBWriter_cfg.py b/CondTools/Hcal/test/HBHENegativeEFilterDBWriter_cfg.py new file mode 100644 index 0000000000000..a6cb9eca5ccd9 --- /dev/null +++ b/CondTools/Hcal/test/HBHENegativeEFilterDBWriter_cfg.py @@ -0,0 +1,51 @@ +database = "sqlite_file:HBHENegativeEFilter_V00_data.db" +tag = "HBHENegativeEFilter_V00_data" +inputfile = "HBHENegativeEFilter_V00_data.bbin" + +import FWCore.ParameterSet.Config as cms + +process = cms.Process('HBHENegativeEFilterDBWrite') + +process.source = cms.Source('EmptyIOVSource', + lastValue = cms.uint64(1), + timetype = cms.string('runnumber'), + firstValue = cms.uint64(1), + interval = cms.uint64(1) +) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(1)) + +process.load("CondCore.CondDB.CondDB_cfi") +process.CondDB.connect = database +# process.CondDB.dbFormat = cms.untracked.int32(1) + +# Data is tagged in the database by the "tag" parameter specified in the +# PoolDBOutputService configuration. We then check if the tag already exists. +# +# -- If the tag does not exist, a new interval of validity (IOV) for this tag +# is created, valid till "end of time". +# +# -- If the tag already exists: the IOV of the previous data is stopped at +# "current time" and we register new data valid from now on (currentTime +# is the time of the current event!). +# +# The "record" parameter should be the same in the PoolDBOutputService +# configuration and in the module which writes the object. It is basically +# used in order to just associate the record with the tag. +# +process.PoolDBOutputService = cms.Service( + "PoolDBOutputService", + process.CondDB, + timetype = cms.untracked.string('runnumber'), + toPut = cms.VPSet(cms.PSet( + record = cms.string("HBHENegativeEFilterRcd"), + tag = cms.string(tag) + )) +) + +process.filereader = cms.EDAnalyzer( + 'HBHENegativeEFilterDBWriter', + inputFile = cms.string(inputfile), + record = cms.string("HBHENegativeEFilterRcd") +) + +process.p = cms.Path(process.filereader) diff --git a/RecoLocalCalo/HcalRecAlgos/interface/HBHENegativeFlag.h b/RecoLocalCalo/HcalRecAlgos/interface/HBHENegativeFlag.h index c63a3e524c2a7..d00b05314e0a5 100755 --- a/RecoLocalCalo/HcalRecAlgos/interface/HBHENegativeFlag.h +++ b/RecoLocalCalo/HcalRecAlgos/interface/HBHENegativeFlag.h @@ -1,53 +1,27 @@ -//--------------------------------------------------------------------------- #ifndef HBHENegativeFlag_H #define HBHENegativeFlag_H + //--------------------------------------------------------------------------- // Negative filter algorithms for HBHE noise flagging -// -// Original Author: Yi Chen (Caltech), (1)3364 (Aug. 21, 2014) -//--------------------------------------------------------------------------- -#include -#include -#include -#include "boost/shared_ptr.hpp" //--------------------------------------------------------------------------- + #include "DataFormats/HcalDigi/interface/HBHEDataFrame.h" -#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" #include "DataFormats/HcalRecHit/interface/HBHERecHit.h" -#include "DataFormats/METReco/interface/HcalCaloFlagLabels.h" #include "CalibFormats/HcalObjects/interface/HcalCalibrations.h" #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h" -#include "CondFormats/HcalObjects/interface/AbsOOTPileupCorrection.h" -//--------------------------------------------------------------------------- -class HBHENegativeFlagSetter; -//--------------------------------------------------------------------------- +#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h" + class HBHENegativeFlagSetter { - public: - HBHENegativeFlagSetter(); - HBHENegativeFlagSetter(double minimumChargeThreshold, - double tS4TS5ChargeThreshold, - int first, int last, - std::vector threshold, - std::vector cut); - ~HBHENegativeFlagSetter(); - void setPulseShapeFlags(HBHERecHit& hbhe, const HBHEDataFrame &digi, - const HcalCoder &coder, const HcalCalibrations &calib); - void setHBHEPileupCorrection(boost::shared_ptr corr); - void setBXInfo(const BunchXParameter *info, unsigned length); - private: - double mMinimumChargeThreshold; - double mTS4TS5ChargeThreshold; - boost::shared_ptr hbhePileupCorr_; - int mFirst; - int mLast; - const BunchXParameter *mBunchCrossingInfo; - unsigned mLengthBunchCrossingInfo; - std::vector > mCut; - private: - bool checkPassFilter(double charge, double discriminant, std::vector > &cuts, - int side); +public: + inline HBHENegativeFlagSetter() : filter_(0) {} + + inline void configFilter(const HBHENegativeEFilter* f) {filter_ = f;} + + void setPulseShapeFlags(HBHERecHit& hbhe, const HBHEDataFrame &digi, + const HcalCoder &coder, const HcalCalibrations &calib); +private: + const HBHENegativeEFilter* filter_; }; -//--------------------------------------------------------------------------- -#endif +#endif diff --git a/RecoLocalCalo/HcalRecAlgos/src/HBHENegativeFlag.cc b/RecoLocalCalo/HcalRecAlgos/src/HBHENegativeFlag.cc index 92ba1972a5413..941efead647c5 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/HBHENegativeFlag.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/HBHENegativeFlag.cc @@ -1,98 +1,14 @@ -//--------------------------------------------------------------------------- -#include -#include -#include -#include -#include -#include - -//--------------------------------------------------------------------------- #include "RecoLocalCalo/HcalRecAlgos/interface/HBHENegativeFlag.h" +#include "CalibFormats/CaloObjects/interface/CaloSamples.h" #include "DataFormats/METReco/interface/HcalCaloFlagLabels.h" -#include "FWCore/Utilities/interface/Exception.h" -#include "FWCore/MessageLogger/interface/MessageLogger.h" - -#include "DataFormats/HcalDigi/interface/HBHEDataFrame.h" -#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h" -#include "DataFormats/HcalRecHit/interface/HBHERecHit.h" - -#include "CalibFormats/HcalObjects/interface/HcalCalibrations.h" -#include "CalibFormats/HcalObjects/interface/HcalCoderDb.h" -#include "CalibCalorimetry/HcalAlgos/interface/HcalPulseShapes.h" - -#include "CondFormats/DataRecord/interface/HcalOOTPileupCorrectionRcd.h" -#include "CondFormats/HcalObjects/interface/OOTPileupCorrectionColl.h" -#include "CondFormats/HcalObjects/interface/OOTPileupCorrData.h" - -#include "Geometry/CaloTopology/interface/HcalTopology.h" -#include "Geometry/Records/interface/IdealGeometryRecord.h" -//--------------------------------------------------------------------------- -HBHENegativeFlagSetter::HBHENegativeFlagSetter() -{ - // - // Argumentless constructor; should not be used - // - // If arguments not properly specified for the constructor, I don't think - // we'd trust the flagging algorithm. - // Set the minimum charge threshold large enough so that nothing will be flagged. - // - - mMinimumChargeThreshold = 99999999; - mTS4TS5ChargeThreshold = 99999999; - mFirst = 4; - mLast = 6; - mBunchCrossingInfo = NULL; - mLengthBunchCrossingInfo = 0; -} -//--------------------------------------------------------------------------- -HBHENegativeFlagSetter::HBHENegativeFlagSetter(double minimumChargeThreshold, - double tS4TS5ChargeThreshold, int first, int last, - std::vector threshold, std::vector cut) -{ - // - // The constructor that should be used - // - // Copies various thresholds and limits and parameters to the class for future use. - // - - mMinimumChargeThreshold = minimumChargeThreshold; - mTS4TS5ChargeThreshold = tS4TS5ChargeThreshold; - mFirst = first; - mLast = last; - - if (mFirst < 2 || mLast < mFirst) - throw cms::Exception("Invalid mFirst, mLast specification"); - - mBunchCrossingInfo = NULL; - mLengthBunchCrossingInfo = 0; - - for(int i = 0; i < (int)std::min(threshold.size(), cut.size()); i++) - mCut.push_back(std::pair(threshold[i], cut[i])); - std::sort(mCut.begin(), mCut.end()); -} -//--------------------------------------------------------------------------- -HBHENegativeFlagSetter::~HBHENegativeFlagSetter() +void HBHENegativeFlagSetter::setPulseShapeFlags( + HBHERecHit &hbhe, + const HBHEDataFrame &digi, + const HcalCoder &coder, + const HcalCalibrations &calib) { - // Dummy destructor - there's nothing to destruct by hand here -} -//--------------------------------------------------------------------------- -void HBHENegativeFlagSetter::setPulseShapeFlags(HBHERecHit &hbhe, - const HBHEDataFrame &digi, - const HcalCoder &coder, - const HcalCalibrations &calib) -{ - // - // Decide if a digi/pulse is good or bad using negative energy discriminants - // - // SetPulseShapeFlags determines the total charge in the digi. - // If the charge is above the minimum threshold, the code then - // runs the flag-setting algorithms to determine whether the - // flags should be set. - // - - const OOTPileupCorrData* corrObj = dynamic_cast(hbhePileupCorr_.get()); - if (corrObj) + if (filter_) { CaloSamples cs; coder.adc2fC(digi,cs); @@ -105,89 +21,8 @@ void HBHENegativeFlagSetter::setPulseShapeFlags(HBHERecHit &hbhe, ts[i] = cs[i] - calib.pedestal(capid); } - double ChargeInWindow = 0.0; - for(int i = mFirst; i <= mLast && i < CaloSamples::MAXSAMPLES; i++) - ChargeInWindow += ts[i]; - if(ChargeInWindow < mMinimumChargeThreshold) - return; - - const OOTPileupCorrDataFcn& fcn = corrObj->getCorrectionByID(hbhe.id()); - const PiecewiseScalingPolynomial& a1 = fcn.getA1(); - const PiecewiseScalingPolynomial& a2 = fcn.getA2(); - - bool passes = true; - for (int i = mFirst; i <= mLast && passes; i++) - { - const double ecorr = ts[i] - a1(ts[i-1]) - a2(ts[i-2]); - passes = checkPassFilter(ChargeInWindow, ecorr, mCut, -1); - } + const bool passes = filter_->checkPassFilter(hbhe.id(), &ts[0], nRead); if (!passes) hbhe.setFlagField(1, HcalCaloFlagLabels::HBHENegativeNoise); } } -//--------------------------------------------------------------------------- -void HBHENegativeFlagSetter::setHBHEPileupCorrection(boost::shared_ptr corr) -{ - hbhePileupCorr_ = corr; -} -//--------------------------------------------------------------------------- -void HBHENegativeFlagSetter::setBXInfo(const BunchXParameter *info, unsigned length) -{ - mBunchCrossingInfo = info; - mLengthBunchCrossingInfo = length; -} -//--------------------------------------------------------------------------- -bool HBHENegativeFlagSetter::checkPassFilter(double charge, - double discriminant, - std::vector > &cuts, - int side) -{ - // - // Checks whether Discriminant value passes Cuts for the specified Charge. True if pulse is good. - // - // The "Cuts" pairs are assumed to be sorted in terms of size from small to large, - // where each "pair" = (Charge, Discriminant) - // "Side" is either positive or negative, which determines whether to discard the pulse if discriminant - // is greater or smaller than the cut value - // - - if(cuts.size() == 0) // safety check that there are some cuts defined - return true; - - if(charge <= cuts[0].first) // too small to cut on - return true; - - int indexLargerThanCharge = -1; // find the range it is falling in - for(int i = 1; i < (int)cuts.size(); i++) - { - if(cuts[i].first > charge) - { - indexLargerThanCharge = i; - break; - } - } - - double limit = 1000000; - - if(indexLargerThanCharge == -1) // if charge is greater than the last entry, assume flat line - limit = cuts[cuts.size()-1].second; - else // otherwise, do a linear interpolation to find the cut position - { - double C1 = cuts[indexLargerThanCharge].first; - double C2 = cuts[indexLargerThanCharge-1].first; - double L1 = cuts[indexLargerThanCharge].second; - double L2 = cuts[indexLargerThanCharge-1].second; - - limit = (charge - C1) / (C2 - C1) * (L2 - L1) + L1; - } - - if(side > 0 && discriminant > limit) - return false; - if(side < 0 && discriminant < limit) - return false; - - return true; -} -//--------------------------------------------------------------------------- - - diff --git a/RecoLocalCalo/HcalRecAlgos/src/HBHETimingShapedFlag.cc b/RecoLocalCalo/HcalRecAlgos/src/HBHETimingShapedFlag.cc index 72b213250db55..c81e36a612dad 100644 --- a/RecoLocalCalo/HcalRecAlgos/src/HBHETimingShapedFlag.cc +++ b/RecoLocalCalo/HcalRecAlgos/src/HBHETimingShapedFlag.cc @@ -14,6 +14,10 @@ HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter() { } +HBHETimingShapedFlagSetter::~HBHETimingShapedFlagSetter() +{ +} + HBHETimingShapedFlagSetter::HBHETimingShapedFlagSetter(const std::vector& v_userEnvelope) { makeTfilterEnvelope(v_userEnvelope); diff --git a/RecoLocalCalo/HcalRecProducers/python/HcalHitReconstructor_hbhe_cfi.py b/RecoLocalCalo/HcalRecProducers/python/HcalHitReconstructor_hbhe_cfi.py index 7b29db8042a9d..c3ef1202f030e 100644 --- a/RecoLocalCalo/HcalRecProducers/python/HcalHitReconstructor_hbhe_cfi.py +++ b/RecoLocalCalo/HcalRecProducers/python/HcalHitReconstructor_hbhe_cfi.py @@ -31,7 +31,9 @@ setTimingShapedCutsFlags = cms.bool(True), setTimingTrustFlags = cms.bool(False), # timing flags currently only implemented for HF setPulseShapeFlags = cms.bool(True), - setNegativeFlags = cms.bool(True), # only in HBHE + + # Disable negative energy filter pending db support + setNegativeFlags = cms.bool(False), flagParameters= cms.PSet(nominalPedestal=cms.double(3.0), #fC hitEnergyMinimum=cms.double(1.0), #GeV @@ -77,13 +79,6 @@ UseDualFit = cms.bool(True), TriangleIgnoreSlow = cms.bool(False)), - negativeParameters = cms.PSet(MinimumChargeThreshold = cms.double(20), - TS4TS5ChargeThreshold = cms.double(70), - First = cms.int32(4), - Last = cms.int32(6), - Threshold = cms.vdouble(100, 120, 160, 200, 300, 500, 1.0e4), - Cut = cms.vdouble(-50, -100, -100, -100, -100, -100, -1.0e6)), - # shaped cut parameters are triples of (energy, low time threshold, high time threshold) values. # The low and high thresholds must straddle zero (i.e., low<0, high>0); use win_offset to shift. # win_gain is applied to both threshold values before win_offset. diff --git a/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc b/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc index 7987f17d0004e..198b03c73502a 100644 --- a/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc +++ b/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.cc @@ -13,6 +13,7 @@ #include "Geometry/Records/interface/IdealGeometryRecord.h" #include "CondFormats/DataRecord/interface/HcalOOTPileupCorrectionRcd.h" #include "CondFormats/DataRecord/interface/HcalOOTPileupCompatibilityRcd.h" +#include "CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h" #include "CondFormats/HcalObjects/interface/OOTPileupCorrectionColl.h" #include "CondFormats/HcalObjects/interface/OOTPileupCorrData.h" #include @@ -45,7 +46,6 @@ HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf): mcOOTCorrectionName_(""), mcOOTCorrectionCategory_("MC"), setPileupCorrection_(0), - setPileupCorrectionForNegative_(0), paramTS(0), puCorrMethod_(conf.existsAs("puCorrMethod") ? conf.getParameter("puCorrMethod") : 0), cntprtCorrMethod_(0), @@ -97,7 +97,6 @@ HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf): setPileupCorrection_ = 0; if(puCorrMethod_ == 1) setPileupCorrection_ = &HcalSimpleRecAlgo::setHBHEPileupCorrection; - setPileupCorrectionForNegative_ = &HBHENegativeFlagSetter::setHBHEPileupCorrection; bool timingShapedCutsFlags = conf.getParameter("setTimingShapedCutsFlags"); if (timingShapedCutsFlags) @@ -159,17 +158,8 @@ HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf): psPulseShape.getParameter("TriangleIgnoreSlow")); } // if (setPulseShapeFlags_) if (setNegativeFlags_) - { - const edm::ParameterSet &psNegative = conf.getParameter("negativeParameters"); - hbheNegativeFlagSetter_ = new HBHENegativeFlagSetter( - psNegative.getParameter("MinimumChargeThreshold"), - psNegative.getParameter("TS4TS5ChargeThreshold"), - psNegative.getParameter("First"), - psNegative.getParameter("Last"), - psNegative.getParameter >("Threshold"), - psNegative.getParameter >("Cut")); - } - + hbheNegativeFlagSetter_ = new HBHENegativeFlagSetter(); + produces(); } else if (!strcasecmp(subd.c_str(),"HO")) { subdet_=HcalOuter; @@ -253,10 +243,7 @@ HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf): if (conf.existsAs("mcOOTCorrectionCategory")) mcOOTCorrectionCategory_ = conf.getParameter("mcOOTCorrectionCategory"); if (dataOOTCorrectionName_.empty() && mcOOTCorrectionName_.empty()) - { setPileupCorrection_ = 0; - setPileupCorrectionForNegative_ = 0; - } reco_.setpuCorrMethod(puCorrMethod_); if(puCorrMethod_ == 2) { @@ -287,11 +274,18 @@ HcalHitReconstructor::HcalHitReconstructor(edm::ParameterSet const& conf): HcalHitReconstructor::~HcalHitReconstructor() { delete hbheFlagSetter_; - delete hfdigibit_; delete hbheHSCPFlagSetter_; delete hbhePulseShapeFlagSetter_; + delete hbheNegativeFlagSetter_; + delete hbheTimingShapedFlagSetter_; + delete hfdigibit_; + delete hfS9S1_; + delete hfS8S1_; delete hfPET_; + delete saturationFlagSetter_; + delete HFTimingTrustFlagSetter_; + delete paramTS; } @@ -385,24 +379,30 @@ void HcalHitReconstructor::produce(edm::Event& e, const edm::EventSetup& eventSe if( testMethod1Ptr ) isMethod1Set = true; (reco_.*setPileupCorrection_)(pileupCorrections->get(corrName_, cat_)); } + } - if(setPileupCorrectionForNegative_ && hbheNegativeFlagSetter_) - (hbheNegativeFlagSetter_->*setPileupCorrectionForNegative_)(pileupCorrections->get(corrName_, cat_)); - } -// Only for HBHE - if( subdet_ == HcalBarrel ){ - if( !cntprtCorrMethod_ ){ - cntprtCorrMethod_++; - if( puCorrMethod_ == 2 ) LogTrace("HcalPUcorrMethod") << "Using Hcal OOTPU method 2" << std::endl; - else if( puCorrMethod_ == 1 ){ - if( isMethod1Set ) LogTrace("HcalPUcorrMethod") << "Using Hcal OOTPU method 1" << std::endl; - else edm::LogWarning("HcalPUcorrMethod") <<"puCorrMethod_ set to be 1 but method 1 is NOT activated (method 0 used instead)!\n" - <<"Please check GlobalTag usage or method 1 seperately disabled by dataOOTCorrectionName & mcOOTCorrectionName?" << std::endl; - }else if (puCorrMethod_ == 3) { + // Configure the negative energy filter + edm::ESHandle negEhandle; + if (hbheNegativeFlagSetter_) + { + eventSetup.get().get(negEhandle); + hbheNegativeFlagSetter_->configFilter(negEhandle.product()); + } + + // Only for HBHE + if( subdet_ == HcalBarrel ) { + if( !cntprtCorrMethod_ ) { + cntprtCorrMethod_++; + if( puCorrMethod_ == 2 ) LogTrace("HcalPUcorrMethod") << "Using Hcal OOTPU method 2" << std::endl; + else if( puCorrMethod_ == 1 ){ + if( isMethod1Set ) LogTrace("HcalPUcorrMethod") << "Using Hcal OOTPU method 1" << std::endl; + else edm::LogWarning("HcalPUcorrMethod") <<"puCorrMethod_ set to be 1 but method 1 is NOT activated (method 0 used instead)!\n" + <<"Please check GlobalTag usage or method 1 separately disabled by dataOOTCorrectionName & mcOOTCorrectionName?" << std::endl; + } else if (puCorrMethod_ == 3) { LogTrace("HcalPUcorrMethod") << "Using Hcal Deterministic Fit Method!" << std::endl; - } else LogTrace("HcalPUcorrMethod") << "Using Hcal OOTPU method 0" << std::endl; - } - } + } else LogTrace("HcalPUcorrMethod") << "Using Hcal OOTPU method 0" << std::endl; + } + } // GET THE BEAM CROSSING INFO HERE, WHEN WE UNDERSTAND HOW THINGS WORK. // Then, call "setBXInfo" method of the reco_ object. @@ -532,11 +532,11 @@ void HcalHitReconstructor::produce(edm::Event& e, const edm::EventSetup& eventSe hbheTimingShapedFlagSetter_->SetTimingShapedFlags(rec->back()); if (setNoiseFlags_) hbheFlagSetter_->SetFlagsFromDigi(&(*topo),rec->back(),*i,coder,calibrations,first,toadd); - if (setPulseShapeFlags_ == true) + if (setPulseShapeFlags_) hbhePulseShapeFlagSetter_->SetPulseShapeFlags(rec->back(), *i, coder, calibrations); - if (setNegativeFlags_ == true) - hbheNegativeFlagSetter_->setPulseShapeFlags(rec->back(), *i, coder, calibrations); - if (setSaturationFlags_) + if (setNegativeFlags_) + hbheNegativeFlagSetter_->setPulseShapeFlags(rec->back(), *i, coder, calibrations); + if (setSaturationFlags_) saturationFlagSetter_->setSaturationFlag(rec->back(),*i); if (correctTiming_) HcalTimingCorrector::Correct(rec->back(), *i, favorite_capid); diff --git a/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.h b/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.h index ddedb21e1986b..71f40ed029851 100644 --- a/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.h +++ b/RecoLocalCalo/HcalRecProducers/src/HcalHitReconstructor.h @@ -54,7 +54,6 @@ class HcalTopology; private: typedef void (HcalSimpleRecAlgo::*SetCorrectionFcn)(boost::shared_ptr); - typedef void (HBHENegativeFlagSetter::*SetCorrectionFcnForNegative)(boost::shared_ptr); HcalSimpleRecAlgo reco_; HcalADCSaturationFlag* saturationFlagSetter_; @@ -106,7 +105,6 @@ class HcalTopology; std::string mcOOTCorrectionName_; std::string mcOOTCorrectionCategory_; SetCorrectionFcn setPileupCorrection_; - SetCorrectionFcnForNegative setPileupCorrectionForNegative_; HcalRecoParams* paramTS; // firstSample & sampleToAdd from DB std::unique_ptr HFDigiTimeParams; // HF DigiTime parameters