Merge pull request #8447 from igv4321/negative-energy-filter-fix

Database object for the HCAL negative energy filter

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);

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:
+    <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<HBHENegativeEFilterRcd> {};
+
+#endif // CondFormats_HBHENegativeEFilterRcd_h
+

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);

CondFormats/HcalObjects/interface/HBHENegativeEFilter.h

@@ -0,0 +1,87 @@
+#ifndef CondFormats_HcalObjects_HBHENegativeEFilter_h_
+#define CondFormats_HcalObjects_HBHENegativeEFilter_h_
+
+#include <vector>
+#include <utility>
+#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<PiecewiseScalingPolynomial>& a1vec,
+                        const std::vector<PiecewiseScalingPolynomial>& a2vec,
+                        const std::vector<uint32_t>& iEtaLimits,
+                        const std::vector<std::pair<double,double> >& 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<uint32_t>& getEtaLimits() const
+        {return iEtaLimits_;}
+    inline const std::vector<std::pair<double,double> >& 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<PiecewiseScalingPolynomial> a1v_;
+    std::vector<PiecewiseScalingPolynomial> a2v_;
+    std::vector<uint32_t> iEtaLimits_;
+    std::vector<std::pair<double,double> > cut_;
+    double minCharge_;
+    uint32_t tFirst_;
+    uint32_t tLast_;
+
+    friend class boost::serialization::access;
+
+    template<class Archive>
+    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<class Archive>
+    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_

CondFormats/HcalObjects/src/HBHENegativeEFilter.cc

@@ -0,0 +1,133 @@
+#include <cmath>
+#include <climits>
+
+#include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h"
+
+HBHENegativeEFilter::HBHENegativeEFilter(
+    const std::vector<PiecewiseScalingPolynomial>& a1vec,
+    const std::vector<PiecewiseScalingPolynomial>& a2vec,
+    const std::vector<uint32_t>& iEtaLimits,
+    const std::vector<std::pair<double,double> >& 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<std::size_t>(UINT_MAX - 1U))
+        return false;
+    for (std::size_t i=1; i<nLimits; ++i)
+        if (!(iEtaLimits_[i-1] < iEtaLimits_[i]))
+            return false;
+
+    if (a1v_.size() != nLimits + 1)
+        return false;
+    if (a2v_.size() != nLimits + 1)
+        return false;
+
+    const std::size_t sz = cut_.size();
+    if (sz >= static_cast<std::size_t>(UINT_MAX - 1U))
+        return false;
+    for (std::size_t i=1; i<sz; ++i)
+        if (!(cut_[i-1U].first < cut_[i].first))
+            return false;
+
+    if (tFirst_ < 2U)
+        return false;
+    if (!(tFirst_ <= tLast_))
+        return false;
+
+    return true;
+}
+
+bool HBHENegativeEFilter::operator==(const HBHENegativeEFilter& r) const
+{
+    if (cut_.empty() && r.cut_.empty())
+        return true;
+    else
+        return a1v_ == r.a1v_ &&
+               a2v_ == r.a2v_ &&
+               iEtaLimits_ == r.iEtaLimits_ &&
+               cut_ == r.cut_ &&
+               minCharge_ == r.minCharge_ &&
+               tFirst_ == r.tFirst_ &&
+               tLast_ == r.tLast_;
+}
+
+unsigned HBHENegativeEFilter::getEtaIndex(const HcalDetId& id) const
+{
+    const unsigned nLimits = iEtaLimits_.size();
+    unsigned which(0U);
+    if (nLimits)
+    {
+        const uint32_t uEta = std::abs(id.ieta());
+        const uint32_t* limits(&iEtaLimits_[0]);
+        for (; which<nLimits; ++which)
+            if (uEta < limits[which])
+                break;
+    }
+    return which;
+}
+
+bool HBHENegativeEFilter::checkPassFilter(const HcalDetId& id,
+                                          const double* ts, const unsigned lenTS) const
+{
+    bool passes = true;
+    const unsigned sz = cut_.size();
+    if (sz)
+    {
+        double chargeInWindow = 0.0;
+        for (unsigned i=tFirst_; i<=tLast_ && i<lenTS; ++i)
+            chargeInWindow += ts[i];
+        if (chargeInWindow >= minCharge_)
+        {
+            // Figure out the cut value for this charge
+            const std::pair<double,double>* 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<lenTS && passes; ++i)
+            {
+                const double ecorr = ts[i] - a1(ts[i-1U]) - a2(ts[i-2U]);
+                passes = ecorr >= cutValue;
+            }
+        }
+    }
+    return passes;
+}

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);

CondFormats/HcalObjects/src/classes.h

@@ -102,6 +102,10 @@ namespace CondFormats_HcalObjects {
     std::vector<HcalInterpolatedPulse> myHcalInterpolatedPulseVec;
     HBHEChannelGroups myHBHEChannelGroups;
     HcalInterpolatedPulseColl myHcalInterpolatedPulseColl;
+
+    // HBHE negative energy filter
+    std::vector<PiecewiseScalingPolynomial> myPiecewiseScalingPolynomialVec;
+    HBHENegativeEFilter myHBHENegativeEFilter;
   };
 }
 

CondFormats/HcalObjects/src/classes_def.xml

@@ -415,4 +415,6 @@
   <field name="group_" mapping="blob" />
  </class>
  <class name="HcalInterpolatedPulseColl" class_version="1"/>
+ <class name="std::vector<PiecewiseScalingPolynomial>"/>
+ <class name="HBHENegativeEFilter" class_version="1"/>
 </lcgdict>   

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"

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<HBHENegativeEFilter> HBHENegativeEFilterDBWriter;
+typedef BoostIODBReader<HBHENegativeEFilter,HBHENegativeEFilterRcd> HBHENegativeEFilterDBReader;
+
+DEFINE_FWK_MODULE(HBHENegativeEFilterDBWriter);
+DEFINE_FWK_MODULE(HBHENegativeEFilterDBReader);

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)