Merge pull request #6349 from wmtford/mergedCPE_730pr2

Flag merged clusters for relaxed CPE errors

DataFormats/SiStripCluster/interface/SiStripCluster.h

@@ -12,6 +12,9 @@ class SiStripCluster  {
   typedef std::vector<SiStripDigi>::const_iterator   SiStripDigiIter;
   typedef std::pair<SiStripDigiIter,SiStripDigiIter>   SiStripDigiRange;
 
+  static const uint16_t stripIndexMask = 0x7FFF;  // The first strip index is in the low 15 bits of firstStrip_
+  static const uint16_t mergedValueMask = 0x8000;  // The merged state is given by the high bit of firstStrip_
+
   /** Construct from a range of digis that form a cluster and from 
    *  a DetID. The range is assumed to be non-empty.
    */
@@ -31,9 +34,16 @@ class SiStripCluster  {
   // errors to the corresponding rechit.
   error_x(-99999.9){}
 
-  /** The number of the first strip in the cluster
+  template<typename Iter>
+  SiStripCluster(const uint16_t& firstStrip, Iter begin, Iter end, bool merged):
+	 amplitudes_(begin,end), firstStrip_(firstStrip), error_x(-99999.9) {
+	   if (merged) firstStrip_ |= mergedValueMask;  // if this is a candidate merged cluster
+	 }
+
+  /** The number of the first strip in the cluster.
+   *  The high bit of firstStrip_ indicates whether the cluster is a candidate for being merged.
    */
-  uint16_t firstStrip() const {return firstStrip_;}
+  uint16_t firstStrip() const {return firstStrip_ & stripIndexMask;}
 
   /** The amplitudes of the strips forming the cluster.
    *  The amplitudes are on consecutive strips; if a strip is missing
@@ -58,6 +68,12 @@ class SiStripCluster  {
    */
    int  charge() const { return std::accumulate(amplitudes().begin(), amplitudes().end(), int(0)); }
 
+  /** Test (set) the merged status of the cluster
+   *
+   */
+  bool isMerged() const {return (firstStrip_ & mergedValueMask) != 0;}
+  void setMerged(bool mergedState) {mergedState ? firstStrip_ |= mergedValueMask : firstStrip_ &= stripIndexMask;}
+
   float getSplitClusterError () const    {  return error_x;  }
   void  setSplitClusterError ( float errx ) { error_x = errx; }
 

DataFormats/SiStripCluster/src/SiStripCluster.cc

@@ -36,6 +36,7 @@ float SiStripCluster::barycenter() const{
   }
 
   // strip centers are offcet by half pitch w.r.t. strip numbers,
-  // so one has to add 0.5 to get the correct barycenter position
-  return float(firstStrip_) + float(sumx) / float(suma) + 0.5f;
+  // so one has to add 0.5 to get the correct barycenter position.
+  // Need to mask off the high bit of firstStrip_, which contains the merged status.
+  return float((firstStrip_ & stripIndexMask)) + float(sumx) / float(suma) + 0.5f;
 }

RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusterizer.cc

@@ -4,13 +4,19 @@
 #include "FWCore/Framework/interface/Event.h"
 #include "FWCore/Utilities/interface/transform.h"
 #include "boost/foreach.hpp"
+#include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
+#include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
 
 SiStripClusterizer::
 SiStripClusterizer(const edm::ParameterSet& conf) 
-  : inputTags( conf.getParameter<std::vector<edm::InputTag> >("DigiProducersList") ),
+  : confClusterizer_(conf.getParameter<edm::ParameterSet>("Clusterizer")),
+    inputTags( conf.getParameter<std::vector<edm::InputTag> >("DigiProducersList") ),
     algorithm( StripClusterizerAlgorithmFactory::create(conf.getParameter<edm::ParameterSet>("Clusterizer")) ) {
   produces< edmNew::DetSetVector<SiStripCluster> > ();
   inputTokens = edm::vector_transform( inputTags, [this](edm::InputTag const & tag) { return consumes< edm::DetSetVector<SiStripDigi> >(tag);} );
+  doRefineCluster_ = confClusterizer_.existsAs<bool>("doRefineCluster") ? confClusterizer_.getParameter<bool>("doRefineCluster") : false;
+  occupancyThreshold_ = confClusterizer_.existsAs<double>("occupancyThreshold") ? confClusterizer_.getParameter<double>("occupancyThreshold") : 0.05;
+  widthThreshold_ = confClusterizer_.existsAs<unsigned>("widthThreshold") ? confClusterizer_.getParameter<unsigned>("widthThreshold") : 4;
 }
 
 void SiStripClusterizer::
@@ -25,7 +31,10 @@ produce(edm::Event& event, const edm::EventSetup& es)  {
   algorithm->initialize(es);  
 
   BOOST_FOREACH( const edm::EDGetTokenT< edm::DetSetVector<SiStripDigi> >& token, inputTokens) {
-    if(      findInput( token, inputOld, event) ) algorithm->clusterize(*inputOld, *output); 
+    if(      findInput( token, inputOld, event) ) {
+      algorithm->clusterize(*inputOld, *output);
+      if (doRefineCluster_) refineCluster(inputOld, output);
+    } 
 //     else if( findInput( tag, inputNew, event) ) algorithm->clusterize(*inputNew, *output);
     else edm::LogError("Input Not Found") << "[SiStripClusterizer::produce] ";// << tag;
   }
@@ -43,3 +52,38 @@ findInput(const edm::EDGetTokenT<T>& tag, edm::Handle<T>& handle, const edm::Eve
     e.getByToken( tag, handle);
     return handle.isValid();
 }
+
+void SiStripClusterizer::
+refineCluster(const edm::Handle< edm::DetSetVector<SiStripDigi> >& input,
+	      std::auto_ptr< edmNew::DetSetVector<SiStripCluster> >& output) {
+  if (input->size() == 0) return;
+
+  // Flag merge-prone clusters for relaxed CPE errors
+  // Criterion is sensor occupancy and cluster width exceeding thresholds
+
+  for (edmNew::DetSetVector<SiStripCluster>::const_iterator det=output->begin(); det!=output->end(); det++) {
+    uint32_t detId = det->id();
+    // Find the number of good strips in this sensor
+    int nchannideal = SiStripDetCabling_->nApvPairs(detId) * 2 * 128;
+    int nchannreal = 0;
+    for(int strip = 0; strip < nchannideal; ++strip)
+      if(!quality_->IsStripBad(detId,strip)) ++nchannreal;
+
+    edm::DetSetVector<SiStripDigi>::const_iterator digis = input->find(detId);
+    if (digis != input->end()) {
+      int ndigi = digis->size();
+      for (edmNew::DetSet<SiStripCluster>::iterator clust = det->begin(); clust != det->end(); clust++) {
+	if (ndigi > occupancyThreshold_*nchannreal && clust->amplitudes().size() >= widthThreshold_) clust->setMerged(true);
+	else clust->setMerged(false);
+      }
+      // std::cout << "Sensor:strips_occStrips_clust " << nchannreal << " " << ndigi << " " << det->size() << std::endl;
+    }
+  }  // traverse sensors
+}
+
+void SiStripClusterizer::beginRun(const edm::Run& run, const edm::EventSetup& es ) {
+  if (doRefineCluster_) {
+    es.get<SiStripDetCablingRcd>().get( SiStripDetCabling_);
+    es.get<SiStripQualityRcd>().get("", quality_);
+  }
+}

RecoLocalTracker/SiStripClusterizer/plugins/SiStripClusterizer.h

@@ -6,6 +6,8 @@
 #include "RecoLocalTracker/SiStripClusterizer/interface/StripClusterizerAlgorithm.h"
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "FWCore/Utilities/interface/InputTag.h"
+#include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
+#include "CalibTracker/Records/interface/SiStripQualityRcd.h"
 
 #include <vector>
 #include <memory>
@@ -15,17 +17,26 @@ class SiStripClusterizer : public edm::stream::EDProducer<>  {
 public:
 
   explicit SiStripClusterizer(const edm::ParameterSet& conf);
+  void beginRun(const edm::Run& run, const edm::EventSetup& es) override;
   virtual void produce(edm::Event&, const edm::EventSetup&);
 
 private:
 
+  edm::ParameterSet confClusterizer_;
+  bool doRefineCluster_;
+  float occupancyThreshold_;
+  unsigned widthThreshold_;
   template<class T> bool findInput(const edm::EDGetTokenT<T>&, edm::Handle<T>&, const edm::Event&);
   const std::vector<edm::InputTag> inputTags;
   std::auto_ptr<StripClusterizerAlgorithm> algorithm;
+  void refineCluster(const edm::Handle< edm::DetSetVector<SiStripDigi> >& input,
+		     std::auto_ptr< edmNew::DetSetVector<SiStripCluster> >& output);
   typedef edm::EDGetTokenT< edm::DetSetVector<SiStripDigi> > token_t;
   typedef std::vector<token_t> token_v;
   token_v inputTokens;
 
+  edm::ESHandle<SiStripDetCabling> SiStripDetCabling_;
+  edm::ESHandle<SiStripQuality> quality_;
 };
 
 #endif

RecoLocalTracker/SiStripClusterizer/python/DefaultClusterizer_cff.py

@@ -10,5 +10,8 @@
     MaxAdjacentBad = cms.uint32(0),
     QualityLabel = cms.string(""),
     RemoveApvShots     = cms.bool(True),
-    minGoodCharge = cms.double(-2069)
+    minGoodCharge = cms.double(-2069),
+    doRefineCluster = cms.bool(False),
+    occupancyThreshold = cms.double(0.05),
+    widthThreshold = cms.uint32(4)
     )

RecoLocalTracker/SiStripRecHitConverter/src/StripCPEfromTrackAngle.cc

@@ -38,7 +38,7 @@ localParameters( const SiStripCluster& cluster, const GeomDetUnit& det, const Lo
 
   const unsigned N = cluster.amplitudes().size();
   const float fullProjection = p.coveredStrips( track+p.drift, ltp.position());
-  const float uerr2 = useLegacyError ? legacyStripErrorSquared(N,std::abs(fullProjection)) : stripErrorSquared( N, std::abs(fullProjection),ssdid.subDetector() );
+  const float uerr2 = useLegacyError || cluster.isMerged() ? legacyStripErrorSquared(N,std::abs(fullProjection)) : stripErrorSquared( N, std::abs(fullProjection),ssdid.subDetector() );
   const float strip = cluster.barycenter() -  0.5f*(1.f-p.backplanecorrection) * fullProjection
     + 0.5f*p.coveredStrips(track, ltp.position());