Adding filter flag to SiStripApproximateClusters

Configuration/StandardSequences/python/DigiToRaw_Repack_cff.py

@@ -72,13 +72,17 @@
                                         cms.InputTag('siStripZeroSuppressionHLT','ScopeMode')),
                                 )
 
-from RecoLocalTracker.SiStripClusterizer.SiStripClusters2ApproxClusters_cff import hltSiStripClusters2ApproxClusters
+from RecoLocalTracker.SiStripClusterizer.SiStripClusters2ApproxClusters_cff import * 
 
 from EventFilter.Utilities.EvFFEDExcluder_cfi import EvFFEDExcluder as _EvFFEDExcluder
 rawPrimeDataRepacker = _EvFFEDExcluder.clone(
     src = 'rawDataCollector',
     fedsToExclude = [foo for foo in range(50, 490)]
 )
 
-DigiToApproxClusterRawTask = cms.Task(siStripDigisHLT,siStripZeroSuppressionHLT,siStripClustersHLT,hltSiStripClusters2ApproxClusters,rawPrimeDataRepacker)
+hltScalersRawToDigi =  cms.EDProducer( "ScalersRawToDigi",
+   scalersInputTag = cms.InputTag( "rawDataRepacker" )
+)
+
+DigiToApproxClusterRawTask = cms.Task(siStripDigisHLT,siStripZeroSuppressionHLT,hltScalersRawToDigi,hltBeamSpotProducer,siStripClustersHLT,hltSiStripClusters2ApproxClusters,rawPrimeDataRepacker)
 DigiToApproxClusterRaw = cms.Sequence(DigiToApproxClusterRawTask)

DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h

@@ -8,27 +8,38 @@ class SiStripApproximateCluster {
 public:
   SiStripApproximateCluster() {}
 
-  explicit SiStripApproximateCluster(cms_uint16_t barycenter,
-                                     cms_uint8_t width,
-                                     cms_uint8_t avgCharge,
-                                     bool isSaturated) {
+  explicit SiStripApproximateCluster(
+      cms_uint16_t barycenter, cms_uint8_t width, cms_uint8_t avgCharge, bool filter, bool isSaturated) {
     barycenter_ = barycenter;
     width_ = width;
     avgCharge_ = avgCharge;
+    filter_ = filter;
     isSaturated_ = isSaturated;
   }
 
-  explicit SiStripApproximateCluster(const SiStripCluster& cluster, unsigned int maxNSat);
+  explicit SiStripApproximateCluster(const SiStripCluster& cluster,
+                                     unsigned int maxNSat,
+                                     float hitPredPos,
+                                     bool peakFilter);
 
   cms_uint16_t barycenter() const { return barycenter_; }
   cms_uint8_t width() const { return width_; }
   cms_uint8_t avgCharge() const { return avgCharge_; }
+  bool filter() const { return filter_; }
   bool isSaturated() const { return isSaturated_; }
+  bool peakFilter() const { return peakFilter_; }
 
 private:
   cms_uint16_t barycenter_ = 0;
   cms_uint8_t width_ = 0;
   cms_uint8_t avgCharge_ = 0;
+  bool filter_ = false;
   bool isSaturated_ = false;
+  bool peakFilter_ = false;
+  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.;
 };
 #endif  // DataFormats_SiStripCluster_SiStripApproximateCluster_h

DataFormats/SiStripCluster/interface/SiStripCluster.h

@@ -82,6 +82,10 @@ class SiStripCluster {
    */
   int charge() const;
 
+  bool filter() const;
+
+  bool isFromApprox() const;
+
   /** Test (set) the merged status of the cluster
    *
    */
@@ -99,6 +103,7 @@ class SiStripCluster {
   //these are used if amplitude information is not available (using approximate cluster constructor)
   float barycenter_ = 0;
   int charge_ = 0;
+  bool filter_ = false;
 
   // [email protected], 01/05/12
   // Add cluster errors to be used by rechits from split clusters.

DataFormats/SiStripCluster/src/SiStripApproximateCluster.cc

@@ -3,11 +3,16 @@
 #include <algorithm>
 #include <cmath>
 
-SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& cluster, unsigned int maxNSat) {
+SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& cluster,
+                                                     unsigned int maxNSat,
+                                                     float hitPredPos,
+                                                     bool peakFilter) {
   barycenter_ = std::round(cluster.barycenter() * 10);
   width_ = cluster.size();
   avgCharge_ = cluster.charge() / cluster.size();
+  filter_ = false;
   isSaturated_ = false;
+  peakFilter_ = peakFilter;
 
   //mimicing the algorithm used in StripSubClusterShapeTrajectoryFilter...
   //Looks for 3 adjacent saturated strips (ADC>=254)
@@ -25,6 +30,28 @@ SiStripApproximateCluster::SiStripApproximateCluster(const SiStripCluster& clust
     maxSat = std::max<int>(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<uint8_t>(trimMaxADC_, std::floor(trimMaxFracTotal_ * sum));
+  auto begin = ampls.begin();
+  auto last = ampls.end() - 1;
+  while (hitStripsTrim > 1 && (*begin < std::max<uint8_t>(trimCut, trimMaxFracNeigh_ * (*(begin + 1))))) {
+    hitStripsTrim--;
+    ++begin;
+  }
+  while (hitStripsTrim > 1 && (*last < std::max<uint8_t>(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_;
+  }
 }

DataFormats/SiStripCluster/src/SiStripCluster.cc

@@ -26,6 +26,7 @@ SiStripCluster::SiStripCluster(const SiStripApproximateCluster cluster, const ui
   barycenter_ = cluster.barycenter() / 10.0;
   charge_ = cluster.width() * cluster.avgCharge();
   amplitudes_.resize(cluster.width(), cluster.avgCharge());
+  filter_ = cluster.filter();
 
   float halfwidth_ = 0.5f * float(cluster.width());
 
@@ -60,3 +61,10 @@ float SiStripCluster::barycenter() const {
   // Need to mask off the high bit of firstStrip_, which contains the merged status.
   return float((firstStrip_ & stripIndexMask)) + float(sumx) / float(suma) + 0.5f;
 }
+bool SiStripCluster::filter() const {
+  if (barycenter_ > 0)
+    return filter_;
+  return false;
+}
+
+bool SiStripCluster::isFromApprox() const { return (barycenter_ > 0); }

DataFormats/SiStripCluster/src/classes_def.xml

@@ -1,6 +1,7 @@
 <lcgdict>
 
- <class name="SiStripCluster" ClassVersion="12">
+ <class name="SiStripCluster" ClassVersion="13">
+  <version ClassVersion="13" checksum="1374720584"/>
   <version ClassVersion="12" checksum="2984011925"/>
   <version ClassVersion="11" checksum="3702468681"/>
   <version ClassVersion="10" checksum="3791198690"/>
@@ -24,7 +25,8 @@
  <class name="edm::Wrapper<edmNew::DetSetVector<edm::Ref<edmNew::DetSetVector<SiStripCluster>,SiStripCluster,edmNew::DetSetVector<SiStripCluster>::FindForDetSetVector> > >" />
 
 
- <class name="SiStripApproximateCluster" ClassVersion="5">
+ <class name="SiStripApproximateCluster" ClassVersion="6">
+  <version ClassVersion="6" checksum="132211472"/>
   <version ClassVersion="5" checksum="3495825183"/>
   <version ClassVersion="4" checksum="2854791577"/>
   <version ClassVersion="3" checksum="2041370183"/>

DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h

@@ -44,6 +44,7 @@ namespace sistrip {
   static const uint16_t STRIPS_PER_FEDCH = STRIPS_PER_APV * APVS_PER_FEDCH;
   static const uint16_t STRIPS_PER_FEUNIT = STRIPS_PER_FEDCH * FEDCH_PER_FEUNIT;  // 3072
   static const uint16_t STRIPS_PER_FED = STRIPS_PER_FEUNIT * FEUNITS_PER_FED;     // 24576
+  static constexpr float MeVperADCStrip = 9.5665E-4;
 
   // -------------------- FED buffers --------------------
 

HLTrigger/Configuration/python/customizeHLTforCMSSW.py

@@ -226,6 +226,16 @@ def customizeHLTfor41495(process):
 
     return process
 
+def customizeHLTfor41815(process):
+    # use hlt online BeamSpot for SiStripClusters2ApproxClusters
+    for producer in producers_by_type(process, 'SiStripClusters2ApproxClusters'):
+        producer.beamSpot = cms.InputTag('hltOnlineBeamSpot')
+
+    if hasattr(process, 'HLT_HIRandom_v4'):
+        getattr(process,'HLT_HIRandom_v4').insert(2,process.HLTBeamSpot)
+
+    return process
+
 # CMSSW version specific customizations
 def customizeHLTforCMSSW(process, menuType="GRun"):
 
@@ -236,5 +246,6 @@ def customizeHLTforCMSSW(process, menuType="GRun"):
 
     process = customizeHLTfor41058(process)
     process = customizeHLTfor41495(process)
+    process = customizeHLTfor41815(process)
 
     return process

RecoLocalTracker/SiStripClusterizer/plugins/BuildFile.xml

@@ -6,5 +6,8 @@
   <use name="HeterogeneousCore/CUDAUtilities"/>
   <use name="CUDADataFormats/SiStripCluster"/>
   <use name="cuda"/>
+  <use name="RecoTracker/PixelLowPtUtilities"/>
+  <use name="CalibFormats/SiStripObjects"/>
+  <use name="CalibTracker/SiStripCommon"/>
   <flags EDM_PLUGIN="1"/>
 </library>

RecoLocalTracker/SiStripClusterizer/plugins/SiStripApprox2ApproxClusters.cc

@@ -59,6 +59,7 @@ void SiStripApprox2ApproxClusters::produce(edm::Event& event, edm::EventSetup co
       float barycenter = cluster.barycenter();
       uint8_t width = cluster.width();
       float avgCharge = cluster.avgCharge();
+      bool filter = cluster.filter();
       bool isSaturated = cluster.isSaturated();
 
       switch (approxVersion) {
@@ -86,7 +87,7 @@ void SiStripApprox2ApproxClusters::produce(edm::Event& event, edm::EventSetup co
           break;
       }
 
-      ff.push_back(SiStripApproximateCluster(barycenter, width, avgCharge, isSaturated));
+      ff.push_back(SiStripApproximateCluster(barycenter, width, avgCharge, filter, isSaturated));
     }
   }
 

RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusters2ApproxClusters.cc

@@ -1,16 +1,33 @@
-
-
-#include "FWCore/Framework/interface/MakerMacros.h"
+#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.h"
+#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackBase.h"
+#include "DataFormats/SiStripCommon/interface/ConstantsForHardwareSystems.h"
 #include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
 #include "FWCore/Framework/interface/stream/EDProducer.h"
-#include "FWCore/ParameterSet/interface/ParameterSet.h"
-#include "FWCore/Utilities/interface/InputTag.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 "DataFormats/SiStripCluster/interface/SiStripApproximateCluster.h"
-#include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
-#include "DataFormats/Common/interface/DetSetVectorNew.h"
-#include "DataFormats/Common/interface/DetSetVector.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 <vector>
 #include <memory>
@@ -27,25 +44,105 @@ class SiStripClusters2ApproxClusters : public edm::stream::EDProducer<> {
   edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > 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<reco::BeamSpot> beamSpotToken_;
+
+  edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken_;
+
+  edm::FileInPath fileInPath_;
+  SiStripDetInfo detInfo_;
+
+  std::string csfLabel_;
+  edm::ESGetToken<ClusterShapeHitFilter, CkfComponentsRecord> csfToken_;
+
+  edm::ESGetToken<SiStripNoises, SiStripNoisesRcd> stripNoiseToken_;
+  edm::ESHandle<SiStripNoises> theNoise_;
 };
 
 SiStripClusters2ApproxClusters::SiStripClusters2ApproxClusters(const edm::ParameterSet& conf) {
   inputClusters = conf.getParameter<edm::InputTag>("inputClusters");
   maxNSat = conf.getParameter<unsigned int>("maxSaturatedStrips");
 
   clusterToken = consumes<edmNew::DetSetVector<SiStripCluster> >(inputClusters);
+
+  beamSpot_ = conf.getParameter<edm::InputTag>("beamSpot");
+  beamSpotToken_ = consumes<reco::BeamSpot>(beamSpot_);
+
+  tkGeomToken_ = esConsumes();
+
+  fileInPath_ = edm::FileInPath(SiStripDetInfoFileReader::kDefaultFile);
+  detInfo_ = SiStripDetInfoFileReader::read(fileInPath_.fullPath());
+
+  csfLabel_ = conf.getParameter<std::string>("clusterShapeHitFilterLabel");
+  csfToken_ = esConsumes(edm::ESInputTag("", csfLabel_));
+
+  stripNoiseToken_ = esConsumes();
+
   produces<edmNew::DetSetVector<SiStripApproximateCluster> >();
 }
 
-void SiStripClusters2ApproxClusters::produce(edm::Event& event, edm::EventSetup const&) {
+void SiStripClusters2ApproxClusters::produce(edm::Event& event, edm::EventSetup const& iSetup) {
   auto result = std::make_unique<edmNew::DetSetVector<SiStripApproximateCluster> >();
   const auto& clusterCollection = event.get(clusterToken);
 
+  auto const beamSpotHandle = event.getHandle(beamSpotToken_);
+  auto const& bs = beamSpotHandle.isValid() ? *beamSpotHandle : reco::BeamSpot();
+  if (not beamSpotHandle.isValid()) {
+    edm::LogError("SiStripClusters2ApproxClusters")
+        << "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_);
+
   for (const auto& detClusters : clusterCollection) {
     edmNew::DetSetVector<SiStripApproximateCluster>::FastFiller ff{*result, detClusters.id()};
 
-    for (const auto& cluster : detClusters)
-      ff.push_back(SiStripApproximateCluster(cluster, maxNSat));
+    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<const StripGeomDetUnit*>(det);
+    float mip = 3.9 / (sistrip::MeVperADCStrip / stripDet->surface().bounds().thickness());
+
+    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;
+      theFilter->getSizes(detId, cluster, lp, ldir, hitStrips, hitPredPos);
+
+      bool peakFilter = false;
+      SlidingPeakFinder pf(std::max<int>(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(SiStripApproximateCluster(cluster, maxNSat, hitPredPos, peakFilter));
+    }
   }
 
   event.put(std::move(result));
@@ -55,6 +152,8 @@ void SiStripClusters2ApproxClusters::fillDescriptions(edm::ConfigurationDescript
   edm::ParameterSetDescription desc;
   desc.add<edm::InputTag>("inputClusters", edm::InputTag("siStripClusters"));
   desc.add<unsigned int>("maxSaturatedStrips", 3);
+  desc.add<std::string>("clusterShapeHitFilterLabel", "ClusterShapeHitFilter");  // add CSF label
+  desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot"));         // add BeamSpot tag
   descriptions.add("SiStripClusters2ApproxClusters", desc);
 }