diff --git a/DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h b/DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h new file mode 100644 index 0000000000000..5a5abd0c61898 --- /dev/null +++ b/DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h @@ -0,0 +1,107 @@ +#ifndef DataFormats_SiStripCluster_SiStripApproximateClusterCollection_v1_h +#define DataFormats_SiStripCluster_SiStripApproximateClusterCollection_v1_h + +#include + +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h" +#include +#include +/** + * This class provides a minimal interface that resembles + * edmNew::DetSetVector, but is crafted such that we are comfortable + * to provide an infinite backwards compatibility guarantee for it + * (like all RAW data). Any modifications need to be made with care. + * Please consult core software group if in doubt. +**/ +namespace v1 { + class SiStripApproximateClusterCollection { + public: + // Helper classes to make creation and iteration easier + class Filler { + public: + void push_back(v1::SiStripApproximateCluster const& cluster) { clusters_.push_back(cluster); } + + private: + friend SiStripApproximateClusterCollection; + Filler(std::vector& clusters) : clusters_(clusters) {} + + std::vector& clusters_; + }; + + class const_iterator; + class DetSet { + public: + using const_iterator = std::vector::const_iterator; + + unsigned int id() const { + return std::accumulate(coll_->detIds_.cbegin(), coll_->detIds_.cbegin() + detIndex_ + 1, 0); + } + + void move(unsigned int clusBegin) const { clusBegin_ = clusBegin; } + const_iterator begin() const { return coll_->clusters_.begin() + clusBegin_; } + const_iterator cbegin() const { return begin(); } + const_iterator end() const { return coll_->clusters_.begin() + clusEnd_; } + const_iterator cend() const { return end(); } + + private: + friend SiStripApproximateClusterCollection::const_iterator; + DetSet(SiStripApproximateClusterCollection const* coll, unsigned int detIndex) + : coll_(coll), detIndex_(detIndex), clusEnd_(coll->clusters_.size()) {} + + SiStripApproximateClusterCollection const* const coll_; + unsigned int const detIndex_; + mutable unsigned int clusBegin_ = 0; + unsigned int const clusEnd_; + }; + + class const_iterator { + public: + DetSet operator*() const { return DetSet(coll_, index_); } + + const_iterator& operator++() { + ++index_; + if (index_ == coll_->detIds_.size()) { + *this = const_iterator(); + } + return *this; + } + + const_iterator operator++(int) { + const_iterator clone = *this; + ++(*this); + return clone; + } + + bool operator==(const_iterator const& other) const { return coll_ == other.coll_ and index_ == other.index_; } + bool operator!=(const_iterator const& other) const { return not operator==(other); } + + private: + friend SiStripApproximateClusterCollection; + // default-constructed object acts as the sentinel + const_iterator() = default; + const_iterator(SiStripApproximateClusterCollection const* coll) : coll_(coll) {} + + SiStripApproximateClusterCollection const* coll_ = nullptr; + unsigned int index_ = 0; + }; + + // Actual public interface + SiStripApproximateClusterCollection() = default; + + void reserve(std::size_t dets, std::size_t clusters); + Filler beginDet(unsigned int detId); + + const_iterator begin() const { return clusters_.empty() ? end() : const_iterator(this); } + const_iterator cbegin() const { return begin(); } + const_iterator end() const { return const_iterator(); } + const_iterator cend() const { return end(); } + + private: + // The detIds_ and beginIndices_ have one element for each Det. An + // element of beginIndices_ points to the first cluster of the Det + // in clusters_. + std::vector detIds_; // DetId for the Det + std::vector clusters_; + }; +} // namespace v1 +#endif diff --git a/DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h b/DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h new file mode 100644 index 0000000000000..294cd3fbb860c --- /dev/null +++ b/DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h @@ -0,0 +1,70 @@ +#ifndef DataFormats_SiStripCluster_SiStripApproximateCluster_v1_h +#define DataFormats_SiStripCluster_SiStripApproximateCluster_v1_h + +#include "FWCore/Utilities/interface/typedefs.h" +#include "assert.h" + +class SiStripCluster; + +namespace v1 { + class SiStripApproximateCluster { + public: + SiStripApproximateCluster() {} + + explicit SiStripApproximateCluster(cms_uint16_t compBarycenter, cms_uint8_t width, cms_uint8_t compavgCharge) + : compBarycenter_(compBarycenter), width_(width), compavgCharge_(compavgCharge) {} + + explicit SiStripApproximateCluster(const SiStripCluster& cluster, + unsigned int maxNSat, + float hitPredPos, + float& previous_cluster, + unsigned int& module_length, + unsigned int& previous_module_length, + bool peakFilter); + + const cms_uint16_t compBarycenter() const { return compBarycenter_; } + + float barycenter(float previous_barycenter = 0, + unsigned int module_length = 0, + unsigned int previous_module_length = 0) const { + float barycenter; + cms_uint16_t compBarycenter = (compBarycenter_ & 0x7FFF); + if (previous_barycenter == -999) + barycenter = compBarycenter * maxBarycenter_ / maxRange_; + else { + barycenter = ((compBarycenter * maxBarycenter_ / maxRange_) - (module_length - previous_module_length)) + + previous_barycenter; + } + assert(barycenter <= maxBarycenter_ && "Returning barycenter > maxBarycenter"); + return barycenter; + } + cms_uint8_t width() const { return width_; } + float avgCharge() const { + cms_uint8_t compavgCharge = (compavgCharge_ & 0x3F); + float avgCharge_ = compavgCharge * maxavgCharge_ / maxavgChargeRange_; + assert(avgCharge_ <= maxavgCharge_ && "Returning avgCharge > maxavgCharge"); + return avgCharge_; + } + bool filter() const { return (compavgCharge_ & (1 << kfilterMask)); } + bool isSaturated() const { return (compavgCharge_ & (1 << kSaturatedMask)); } + bool peakFilter() const { return (compBarycenter_ & (1 << kpeakFilterMask)); } + + private: + cms_uint16_t compBarycenter_ = 0; + cms_uint8_t width_ = 0; + cms_uint8_t compavgCharge_ = 0; + static constexpr double maxRange_ = 32767; + static constexpr double maxBarycenter_ = 1536.; + static constexpr double maxavgChargeRange_ = 63; + static constexpr double maxavgCharge_ = 255.; + static constexpr double trimMaxADC_ = 30.; + static constexpr double trimMaxFracTotal_ = .15; + static constexpr double trimMaxFracNeigh_ = .25; + static constexpr double maxTrimmedSizeDiffNeg_ = .7; + static constexpr double maxTrimmedSizeDiffPos_ = 1.; + static constexpr int kfilterMask = 6; + static constexpr int kpeakFilterMask = 7; + static constexpr int kSaturatedMask = 15; + }; +} // namespace v1 +#endif // DataFormats_SiStripCluster_SiStripApproximateCluster_h diff --git a/DataFormats/SiStripCluster/interface/SiStripCluster.h b/DataFormats/SiStripCluster/interface/SiStripCluster.h index 83ada6f6cfc1d..245fcdaa25969 100644 --- a/DataFormats/SiStripCluster/interface/SiStripCluster.h +++ b/DataFormats/SiStripCluster/interface/SiStripCluster.h @@ -3,12 +3,11 @@ #include "DataFormats/SiStripDigi/interface/SiStripDigi.h" #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h" #include #include #include -class SiStripApproximateCluster; - class SiStripCluster { public: typedef std::vector::const_iterator SiStripDigiIter; @@ -45,6 +44,11 @@ class SiStripCluster { } SiStripCluster(const SiStripApproximateCluster cluster, const uint16_t maxStrips); + SiStripCluster(const v1::SiStripApproximateCluster cluster, + const uint16_t maxStrips, + float pc = -999, + unsigned int module_length = 0, + unsigned int previous_module_length = 0); // extend the cluster template diff --git a/DataFormats/SiStripCluster/src/SiStripApproximateClusterCollection_v1.cc b/DataFormats/SiStripCluster/src/SiStripApproximateClusterCollection_v1.cc new file mode 100644 index 0000000000000..2e8e0c4d3c185 --- /dev/null +++ b/DataFormats/SiStripCluster/src/SiStripApproximateClusterCollection_v1.cc @@ -0,0 +1,11 @@ +#include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h" +using namespace v1; +void SiStripApproximateClusterCollection::reserve(std::size_t dets, std::size_t clusters) { + detIds_.reserve(dets); + clusters_.reserve(clusters); +} + +SiStripApproximateClusterCollection::Filler SiStripApproximateClusterCollection::beginDet(unsigned int detId) { + detIds_.push_back((detIds_.empty()) ? detId : detId - (std::accumulate(detIds_.cbegin(), detIds_.cend(), 0))); + return Filler(clusters_); +} diff --git a/DataFormats/SiStripCluster/src/SiStripApproximateCluster_v1.cc b/DataFormats/SiStripCluster/src/SiStripApproximateCluster_v1.cc new file mode 100644 index 0000000000000..85187e8f303dc --- /dev/null +++ b/DataFormats/SiStripCluster/src/SiStripApproximateCluster_v1.cc @@ -0,0 +1,75 @@ +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include +#include +#include +v1::SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& cluster, + unsigned int maxNSat, + float hitPredPos, + float& previous_cluster, + unsigned int& module_length, + unsigned int& previous_module_length, + bool peakFilter) { + bool filter_, isSaturated_, peakFilter_; + if (previous_cluster == -999.) + compBarycenter_ = std::round(cluster.barycenter() * maxRange_ / maxBarycenter_); + else + compBarycenter_ = + std::round(((cluster.barycenter() - previous_cluster) + (module_length - previous_module_length)) * maxRange_ / + maxBarycenter_); + previous_cluster = barycenter(previous_cluster, module_length, previous_module_length); + assert(cluster.barycenter() <= maxBarycenter_ && "Got a barycenter > maxBarycenter"); + assert(compBarycenter_ <= maxRange_ && "Filling compBarycenter > maxRange"); + width_ = std::min(255, (int)cluster.size()); //cluster.size(); + float avgCharge_ = cluster.charge() * 1. / width_; + assert(avgCharge_ <= maxavgCharge_ && "Got a avgCharge > maxavgCharge"); + compavgCharge_ = std::round(avgCharge_ * maxavgChargeRange_ / maxavgCharge_); + assert(compavgCharge_ <= maxavgChargeRange_ && "Filling compavgCharge > maxavgChargeRange"); + filter_ = false; + isSaturated_ = false; + peakFilter_ = peakFilter; + + //mimicing the algorithm used in StripSubClusterShapeTrajectoryFilter... + //Looks for 3 adjacent saturated strips (ADC>=254) + const auto& ampls = cluster.amplitudes(); + unsigned int thisSat = (ampls[0] >= 254), maxSat = thisSat; + for (unsigned int i = 1, n = ampls.size(); i < n; ++i) { + if (ampls[i] >= 254) { + thisSat++; + } else if (thisSat > 0) { + maxSat = std::max(maxSat, thisSat); + thisSat = 0; + } + } + if (thisSat > 0) { + maxSat = std::max(maxSat, thisSat); + } + if (maxSat >= maxNSat) { + filter_ = true; + isSaturated_ = true; + } + + unsigned int hitStripsTrim = ampls.size(); + int sum = std::accumulate(ampls.begin(), ampls.end(), 0); + uint8_t trimCut = std::min(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum)); + auto begin = ampls.begin(); + auto last = ampls.end() - 1; + while (hitStripsTrim > 1 && (*begin < std::max(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) { + hitStripsTrim--; + ++begin; + } + while (hitStripsTrim > 1 && (*last < std::max(trimCut, trimMaxFracNeigh_ * (*(last - 1))))) { + hitStripsTrim--; + --last; + } + if (hitStripsTrim < std::floor(std::abs(hitPredPos) - maxTrimmedSizeDiffNeg_)) { + filter_ = false; + } else if (hitStripsTrim <= std::ceil(std::abs(hitPredPos) + maxTrimmedSizeDiffPos_)) { + filter_ = true; + } else { + filter_ = peakFilter_; + } + compavgCharge_ = (compavgCharge_ | (filter_ << kfilterMask)); + compavgCharge_ = (compavgCharge_ | (peakFilter_ << kpeakFilterMask)); + compBarycenter_ = (compBarycenter_ | (isSaturated_ << kSaturatedMask)); +} diff --git a/DataFormats/SiStripCluster/src/SiStripCluster.cc b/DataFormats/SiStripCluster/src/SiStripCluster.cc index 8586307a8e384..15517b4c0457f 100644 --- a/DataFormats/SiStripCluster/src/SiStripCluster.cc +++ b/DataFormats/SiStripCluster/src/SiStripCluster.cc @@ -39,3 +39,25 @@ SiStripCluster::SiStripCluster(const SiStripApproximateCluster cluster, const ui } firstStrip_ |= approximateMask; } + +SiStripCluster::SiStripCluster(const v1::SiStripApproximateCluster cluster, + const uint16_t maxStrips, + float p_bc, + unsigned int module_length, + unsigned int previous_module_length) + : error_x(-99999.9) { + barycenter_ = cluster.barycenter(p_bc, module_length, previous_module_length); + charge_ = cluster.width() * cluster.avgCharge(); + amplitudes_.resize(cluster.width(), cluster.avgCharge()); + filter_ = cluster.filter(); + + float halfwidth_ = 0.5f * float(cluster.width()); + + //initialize firstStrip_ + firstStrip_ = std::max(barycenter_ - halfwidth_, 0.f); + + if UNLIKELY (firstStrip_ + cluster.width() > maxStrips) { + firstStrip_ = maxStrips - cluster.width(); + } + firstStrip_ |= approximateMask; +} diff --git a/DataFormats/SiStripCluster/src/classes.h b/DataFormats/SiStripCluster/src/classes.h index 94eafdf2ecac7..b0bb9e729749c 100644 --- a/DataFormats/SiStripCluster/src/classes.h +++ b/DataFormats/SiStripCluster/src/classes.h @@ -7,7 +7,9 @@ #include "DataFormats/SiStripCluster/interface/SiStripCluster.h" #include "DataFormats/SiStripCluster/interface/SiStripClustersSOA.h" #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h" #include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h" #include "DataFormats/Common/interface/ContainerMask.h" #endif // SISTRIPCLUSTER_CLASSES_H diff --git a/DataFormats/SiStripCluster/src/classes_def.xml b/DataFormats/SiStripCluster/src/classes_def.xml index 31ff33a914296..ceeb869ab8ccf 100755 --- a/DataFormats/SiStripCluster/src/classes_def.xml +++ b/DataFormats/SiStripCluster/src/classes_def.xml @@ -59,6 +59,38 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2Clusters.cc b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2Clusters.cc index 803c8949f90c6..1d6450c49694d 100644 --- a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2Clusters.cc +++ b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2Clusters.cc @@ -1,6 +1,8 @@ #include "DataFormats/Common/interface/DetSetVectorNew.h" #include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h" #include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h" #include "DataFormats/SiStripCluster/interface/SiStripCluster.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/Frameworkfwd.h" @@ -26,35 +28,84 @@ class SiStripApprox2Clusters : public edm::global::EDProducer<> { private: edm::EDGetTokenT clusterToken_; + edm::EDGetTokenT clusterToken_v1_; edm::ESGetToken tkGeomToken_; + bool v1; }; SiStripApprox2Clusters::SiStripApprox2Clusters(const edm::ParameterSet& conf) { - clusterToken_ = consumes(conf.getParameter("inputApproxClusters")); tkGeomToken_ = esConsumes(); + v1 = conf.getParameter("v1"); + if (v1) { + clusterToken_v1_ = consumes(conf.getParameter("inputApproxClusters")); + } else { + clusterToken_ = consumes(conf.getParameter("inputApproxClusters")); + } produces>(); } void SiStripApprox2Clusters::produce(edm::StreamID id, edm::Event& event, const edm::EventSetup& iSetup) const { auto result = std::make_unique>(); - const auto& clusterCollection = event.get(clusterToken_); const auto& tkGeom = &iSetup.getData(tkGeomToken_); const auto& tkDets = tkGeom->dets(); - for (const auto& detClusters : clusterCollection) { - edmNew::DetSetVector::FastFiller ff{*result, detClusters.id()}; - unsigned int detId = detClusters.id(); + if (clusterToken_v1_.isUninitialized()) { + const auto& clusterCollection = event.get(clusterToken_); + for (const auto& detClusters : clusterCollection) { + edmNew::DetSetVector::FastFiller ff{*result, detClusters.id()}; + unsigned int detId = detClusters.id(); + + uint16_t nStrips{0}; + auto det = std::find_if(tkDets.begin(), tkDets.end(), [detId](auto& elem) -> bool { + return (elem->geographicalId().rawId() == detId); + }); + const StripTopology& p = dynamic_cast(*det)->specificTopology(); + nStrips = p.nstrips() - 1; + + for (const auto& cluster : detClusters) { + ff.push_back(SiStripCluster(cluster, nStrips)); + } + } + } else { + const auto& clusterCollection = event.get(clusterToken_v1_); + std::vector v_strip; + float previous_barycenter = -999.; + unsigned int module_length = 0; + unsigned int previous_module_length = 0; + + unsigned int clusBegin = 0; + for (const auto& detClusters : clusterCollection) { + edmNew::DetSetVector::FastFiller ff{*result, detClusters.id()}; + unsigned int detId = detClusters.id(); + + uint16_t nStrips{0}; + auto det = std::find_if(tkDets.begin(), tkDets.end(), [detId](auto& elem) -> bool { + return (elem->geographicalId().rawId() == detId); + }); + const StripTopology& p = dynamic_cast(*det)->specificTopology(); + nStrips = p.nstrips(); - uint16_t nStrips{0}; - auto det = std::find_if(tkDets.begin(), tkDets.end(), [detId](auto& elem) -> bool { - return (elem->geographicalId().rawId() == detId); - }); - const StripTopology& p = dynamic_cast(*det)->specificTopology(); - nStrips = p.nstrips() - 1; + v_strip.push_back(nStrips); + previous_module_length += (v_strip.size() < 3) ? 0 : v_strip[v_strip.size() - 3]; + module_length += (v_strip.size() < 2) ? 0 : v_strip[v_strip.size() - 2]; + bool first_cluster = true; + detClusters.move(clusBegin); - for (const auto& cluster : detClusters) { - ff.push_back(SiStripCluster(cluster, nStrips)); + for (const auto& cluster : detClusters) { + const auto convertedCluster = SiStripCluster(cluster, + nStrips - 1, + previous_barycenter, + module_length, + first_cluster ? previous_module_length : module_length); + if ((convertedCluster.barycenter()) >= nStrips) { + break; + } + previous_barycenter = convertedCluster.barycenter(); + ++clusBegin; + ff.push_back(convertedCluster); + first_cluster = false; + } } } @@ -64,6 +115,7 @@ void SiStripApprox2Clusters::produce(edm::StreamID id, edm::Event& event, const void SiStripApprox2Clusters::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("inputApproxClusters", edm::InputTag("siStripClusters")); + desc.add("v1", false); descriptions.add("SiStripApprox2Clusters", desc); } diff --git a/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters_v1.cc b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters_v1.cc new file mode 100644 index 0000000000000..5442a54cdc3a5 --- /dev/null +++ b/RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters_v1.cc @@ -0,0 +1,201 @@ +#include "CalibFormats/SiStripObjects/interface/SiStripDetInfo.h" +#include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h" +#include "CondFormats/DataRecord/interface/SiStripNoisesRcd.h" +#include "CondFormats/SiStripObjects/interface/SiStripNoises.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/Common/interface/DetSetVectorNew.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateCluster_v1.h" +#include "DataFormats/SiStripCluster/interface/SiStripApproximateClusterCollection_v1.h" +#include "DataFormats/SiStripCluster/interface/SiStripCluster.h" +#include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackBase.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/FileInPath.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "FWCore/Utilities/interface/ESInputTag.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" +#include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" +#include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h" +#include "RecoTracker/PixelLowPtUtilities/interface/SlidingPeakFinder.h" +#include "RecoTracker/Record/interface/CkfComponentsRecord.h" + +#include +#include + +class SiStripClusters2ApproxClusters_v1 : public edm::stream::EDProducer<> { +public: + explicit SiStripClusters2ApproxClusters_v1(const edm::ParameterSet& conf); + void produce(edm::Event&, const edm::EventSetup&) override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + edm::InputTag inputClusters; + edm::EDGetTokenT > clusterToken; + + unsigned int maxNSat; + static constexpr double subclusterWindow_ = .7; + static constexpr double seedCutMIPs_ = .35; + static constexpr double seedCutSN_ = 7.; + static constexpr double subclusterCutMIPs_ = .45; + static constexpr double subclusterCutSN_ = 12.; + + edm::InputTag beamSpot_; + edm::EDGetTokenT beamSpotToken_; + + edm::ESGetToken tkGeomToken_; + + edm::FileInPath fileInPath_; + SiStripDetInfo detInfo_; + + std::string csfLabel_; + edm::ESGetToken csfToken_; + + edm::ESGetToken stripNoiseToken_; + edm::ESHandle theNoise_; +}; + +SiStripClusters2ApproxClusters_v1::SiStripClusters2ApproxClusters_v1(const edm::ParameterSet& conf) { + inputClusters = conf.getParameter("inputClusters"); + maxNSat = conf.getParameter("maxSaturatedStrips"); + + clusterToken = consumes >(inputClusters); + + beamSpot_ = conf.getParameter("beamSpot"); + beamSpotToken_ = consumes(beamSpot_); + + tkGeomToken_ = esConsumes(); + + fileInPath_ = edm::FileInPath(SiStripDetInfoFileReader::kDefaultFile); + detInfo_ = SiStripDetInfoFileReader::read(fileInPath_.fullPath()); + + csfLabel_ = conf.getParameter("clusterShapeHitFilterLabel"); + csfToken_ = esConsumes(edm::ESInputTag("", csfLabel_)); + + stripNoiseToken_ = esConsumes(); + produces(); +} + +void SiStripClusters2ApproxClusters_v1::produce(edm::Event& event, edm::EventSetup const& iSetup) { + const auto& clusterCollection = event.get(clusterToken); + auto result = std::make_unique(); + result->reserve(clusterCollection.size(), clusterCollection.dataSize()); + + auto const beamSpotHandle = event.getHandle(beamSpotToken_); + auto const& bs = beamSpotHandle.isValid() ? *beamSpotHandle : reco::BeamSpot(); + if (not beamSpotHandle.isValid()) { + edm::LogError("SiStripClusters2ApproxClusters_v1") + << "didn't find a valid beamspot with label \"" << beamSpot_.encode() << "\" -> using (0,0,0)"; + } + + const auto& tkGeom = &iSetup.getData(tkGeomToken_); + const auto& theFilter = &iSetup.getData(csfToken_); + const auto& theNoise_ = &iSetup.getData(stripNoiseToken_); + + float previous_cluster = -999.; + unsigned int module_length = 0; + unsigned int previous_module_length = 0; + const auto tkDets = tkGeom->dets(); + + std::vector v_strip; + + for (const auto& detClusters : clusterCollection) { + auto ff = result->beginDet(detClusters.id()); + + unsigned int detId = detClusters.id(); + const GeomDet* det = tkGeom->idToDet(detId); + double nApvs = detInfo_.getNumberOfApvsAndStripLength(detId).first; + double stripLength = detInfo_.getNumberOfApvsAndStripLength(detId).second; + double barycenter_ypos = 0.5 * stripLength; + + const StripGeomDetUnit* stripDet = dynamic_cast(det); + float mip = 3.9 / (sistrip::MeVperADCStrip / stripDet->surface().bounds().thickness()); + + uint16_t nStrips{0}; + const auto& _detId = detId; // for the capture clause in the lambda function + auto _det = std::find_if(tkDets.begin(), tkDets.end(), [_detId](auto& elem) -> bool { + return (elem->geographicalId().rawId() == _detId); + }); + const StripTopology& p = dynamic_cast(*_det)->specificTopology(); + nStrips = p.nstrips(); + v_strip.push_back(nStrips); + + previous_module_length += (v_strip.size() < 3) ? 0 : v_strip[v_strip.size() - 3]; + module_length += (v_strip.size() < 2) ? 0 : v_strip[v_strip.size() - 2]; + assert(!detClusters.empty()); + bool first_cluster = true; + + for (const auto& cluster : detClusters) { + const LocalPoint& lp = LocalPoint(((cluster.barycenter() * 10 / (sistrip::STRIPS_PER_APV * nApvs)) - + ((stripDet->surface().bounds().width()) * 0.5f)), + barycenter_ypos - (0.5f * stripLength), + 0.); + const GlobalPoint& gpos = det->surface().toGlobal(lp); + GlobalPoint beamspot(bs.position().x(), bs.position().y(), bs.position().z()); + const GlobalVector& gdir = gpos - beamspot; + const LocalVector& ldir = det->toLocal(gdir); + + int hitStrips; + float hitPredPos; + bool usable = theFilter->getSizes(detId, cluster, lp, ldir, hitStrips, hitPredPos); + // (almost) same logic as in StripSubClusterShapeTrajectoryFilter + bool isTrivial = (std::abs(hitPredPos) < 2.f && hitStrips <= 2); + + if (!usable || isTrivial) { + ff.push_back(v1::SiStripApproximateCluster(cluster, + maxNSat, + hitPredPos, + previous_cluster, + module_length, + first_cluster ? previous_module_length : module_length, + true)); + } else { + bool peakFilter = false; + SlidingPeakFinder pf(std::max(2, std::ceil(std::abs(hitPredPos) + subclusterWindow_))); + float mipnorm = mip / std::abs(ldir.z()); + PeakFinderTest test(mipnorm, + detId, + cluster.firstStrip(), + theNoise_, + seedCutMIPs_, + seedCutSN_, + subclusterCutMIPs_, + subclusterCutSN_); + peakFilter = pf.apply(cluster.amplitudes(), test); + + ff.push_back(v1::SiStripApproximateCluster(cluster, + maxNSat, + hitPredPos, + previous_cluster, + module_length, + first_cluster ? previous_module_length : module_length, + peakFilter)); + } + first_cluster = false; + } + } + + event.put(std::move(result)); +} + +void SiStripClusters2ApproxClusters_v1::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("inputClusters", edm::InputTag("siStripClusters")); + desc.add("maxSaturatedStrips", 3); + desc.add("clusterShapeHitFilterLabel", "ClusterShapeHitFilter"); // add CSF label + desc.add("beamSpot", edm::InputTag("offlineBeamSpot")); // add BeamSpot tag + descriptions.add("SiStripClusters2ApproxClusters_v1", desc); +} + +DEFINE_FWK_MODULE(SiStripClusters2ApproxClusters_v1);