New HeavyIon ClusterCompatibility object and producer for beam scraping removal

DataFormats/HeavyIonEvent/interface/ClusterCompatibility.h

@@ -0,0 +1,45 @@
+#ifndef DataFormats_ClusterCompatibility_h
+#define DataFormats_ClusterCompatibility_h
+
+#include <vector>
+
+namespace reco { class ClusterCompatibility {
+public:
+
+  ClusterCompatibility();
+  virtual ~ClusterCompatibility();
+
+  /// Number of valid pixel clusters
+  int nValidPixelHits() const { return nValidPixelHits_; }
+
+  /// Number of vertex-position hypotheses tested
+  int size() const { return z0_.size(); }
+
+  /// Vertex z position for the i-th vertex-position hypothesis
+  float z0(int i) const { return z0_[i]; }
+
+  /// Number of compatible non-edge pixel-barrel clusters 
+  /// for the i-th vertex-position hypothesis
+  int nHit(int i) const { return nHit_[i]; }
+
+  /// Sum of the difference between the expected and actual 
+  /// width of all compatible non-edge pixel-barrel clusters 
+  /// for the i-th vertex-position hypothesis
+  float chi(int i) const { return chi_[i]; }
+
+  void append(float, int, float);
+  void setNValidPixelHits(int nPxl) { nValidPixelHits_ = nPxl; }
+
+protected:
+
+
+  int nValidPixelHits_;
+
+  std::vector<float> z0_;
+  std::vector<int> nHit_;
+  std::vector<float> chi_;
+
+};
+
+}
+#endif

DataFormats/HeavyIonEvent/src/ClusterCompatibility.cc

@@ -0,0 +1,18 @@
+#include "DataFormats/HeavyIonEvent/interface/ClusterCompatibility.h"
+using namespace reco;
+
+ClusterCompatibility::ClusterCompatibility():
+  nValidPixelHits_(0),
+  z0_(),
+  nHit_(),
+  chi_()
+{}
+
+ClusterCompatibility::~ClusterCompatibility() {}
+
+void
+ClusterCompatibility::append(float z0, int nHit, float chi) {
+  z0_.push_back(z0);
+  nHit_.push_back(nHit);
+  chi_.push_back(chi);
+}

DataFormats/HeavyIonEvent/src/classes.h

@@ -1,4 +1,5 @@
 #include "DataFormats/HeavyIonEvent/interface/Centrality.h"
+#include "DataFormats/HeavyIonEvent/interface/ClusterCompatibility.h"
 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
 #include "DataFormats/HeavyIonEvent/interface/HeavyIon.h"
 #include "DataFormats/HeavyIonEvent/interface/VoronoiBackground.h"
@@ -20,6 +21,9 @@ namespace DataFormats_HeavyIonEvent {
      reco::EvtPlaneCollection evcol;
      edm::Wrapper<reco::EvtPlaneCollection> wevcol;
 
+     reco::ClusterCompatibility clus_compat;
+     edm::Wrapper<reco::ClusterCompatibility> w_clus_compat;
+
      reco::VoronoiBackground vor;
      edm::Wrapper<reco::VoronoiBackground> wvor;
 

DataFormats/HeavyIonEvent/src/classes_def.xml

@@ -21,4 +21,6 @@
 <class name="edm::Wrapper<edm::ValueMap<reco::VoronoiBackground> >" />
 <class name="std::vector<reco::VoronoiBackground>" />
 <class name="edm::Wrapper<std::vector<reco::VoronoiBackground> >" />
+<class name="reco::ClusterCompatibility" />
+<class name="edm::Wrapper<reco::ClusterCompatibility>"/>
 </lcgdict>

RecoHI/Configuration/python/Reconstruction_HI_cff.py

@@ -22,6 +22,7 @@
 
 # Heavy Ion Event Characterization
 from RecoHI.HiCentralityAlgos.HiCentrality_cfi import *
+from RecoHI.HiCentralityAlgos.HiClusterCompatibility_cfi import *
 from RecoHI.HiEvtPlaneAlgos.HiEvtPlane_cfi import *
 
 # HCAL noise producer
@@ -38,6 +39,7 @@
                               * hiEgammaSequence
                               * hiParticleFlowReco
                               * hiCentrality
+                              * hiClusterCompatibility
                               * hiEvtPlane
                               * hcalnoise
                               )
@@ -51,6 +53,7 @@
                                               * hiEgammaSequence
                                               * hiParticleFlowReco
                                               * hiCentrality
+                                              * hiClusterCompatibility
                                               * hiEvtPlane
                                               * hcalnoise
                                               )

RecoHI/HiCentralityAlgos/plugins/ClusterCompatibilityProducer.cc

@@ -0,0 +1,179 @@
+//
+// Derived from HLTrigger/special/src/HLTPixelClusterShapeFilter.cc
+// at version 7_5_0_pre3
+//
+// Original Author (of Derivative Producer):  Eric Appelt
+//         Created:  Mon Apr 27, 2015
+
+#include <iostream>
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/Common/interface/Handle.h"
+#include "DataFormats/GeometryVector/interface/LocalPoint.h"
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
+#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
+#include "DataFormats/HeavyIonEvent/interface/ClusterCompatibility.h"
+
+#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
+#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
+#include "Geometry/CommonTopologies/interface/PixelTopology.h"
+#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
+#include "Geometry/CommonDetUnit/interface/GeomDet.h"
+
+
+//
+// class declaration
+//
+
+class ClusterCompatibilityProducer : public edm::stream::EDProducer<> {
+
+  public:
+    explicit ClusterCompatibilityProducer(const edm::ParameterSet&);
+    ~ClusterCompatibilityProducer();
+
+    virtual void produce(edm::Event&, const edm::EventSetup&) override;
+
+  private:
+
+    edm::EDGetTokenT<SiPixelRecHitCollection> inputToken_;
+    edm::InputTag       inputTag_;      // input tag identifying product containing pixel clusters
+    double              minZ_;          // beginning z-vertex position
+    double              maxZ_;          // end z-vertex position
+    double              zStep_;         // size of steps in z-vertex test
+
+    struct VertexHit
+    {
+      float z;
+      float r;
+      float w;
+    };
+
+    struct ContainedHits
+    {
+      float z0;
+      int nHit;
+      float chi;
+    };
+
+    ContainedHits getContainedHits(const std::vector<VertexHit> &hits, double z0) const;
+
+};
+
+ClusterCompatibilityProducer::ClusterCompatibilityProducer(const edm::ParameterSet& config):
+  inputTag_     (config.getParameter<edm::InputTag>("inputTag")),
+  minZ_         (config.getParameter<double>("minZ")),
+  maxZ_         (config.getParameter<double>("maxZ")),
+  zStep_        (config.getParameter<double>("zStep"))
+{
+  inputToken_ = consumes<SiPixelRecHitCollection>(inputTag_);
+  LogDebug("") << "Using the " << inputTag_ << " input collection";
+  produces<reco::ClusterCompatibility>();
+}
+
+ClusterCompatibilityProducer::~ClusterCompatibilityProducer() {}
+
+void
+ClusterCompatibilityProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
+{
+  std::auto_ptr<reco::ClusterCompatibility> creco(new reco::ClusterCompatibility());
+
+  // get hold of products from Event
+  edm::Handle<SiPixelRecHitCollection> hRecHits;
+  iEvent.getByToken(inputToken_, hRecHits);
+
+  // get tracker geometry
+  if (hRecHits.isValid()) {
+    edm::ESHandle<TrackerGeometry> trackerHandle;
+    iSetup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
+    const TrackerGeometry *tgeo = trackerHandle.product();
+    const SiPixelRecHitCollection *hits = hRecHits.product();
+
+    // loop over pixel rechits
+    int nPxlHits=0;
+    std::vector<VertexHit> vhits;
+    for(SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(),
+          end = hits->data().end(); hit != end; ++hit) {
+      if (!hit->isValid())
+        continue;
+      ++nPxlHits;
+      DetId id(hit->geographicalId());
+      if(id.subdetId() != int(PixelSubdetector::PixelBarrel))
+        continue;
+      const PixelGeomDetUnit *pgdu = static_cast<const PixelGeomDetUnit*>(tgeo->idToDet(id));
+      const PixelTopology *pixTopo = &(pgdu->specificTopology());
+      std::vector<SiPixelCluster::Pixel> pixels(hit->cluster()->pixels());
+      bool pixelOnEdge = false;
+      for(std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin();
+          pixel != pixels.end(); ++pixel) {
+        int pixelX = pixel->x;
+        int pixelY = pixel->y;
+        if(pixTopo->isItEdgePixelInX(pixelX) || pixTopo->isItEdgePixelInY(pixelY)) {
+          pixelOnEdge = true;
+          break;
+        }
+      }
+      if (pixelOnEdge)
+        continue;
+
+      LocalPoint lpos = LocalPoint(hit->localPosition().x(),
+                                   hit->localPosition().y(),
+                                   hit->localPosition().z());
+      GlobalPoint gpos = pgdu->toGlobal(lpos);
+      VertexHit vh;
+      vh.z = gpos.z();
+      vh.r = gpos.perp();
+      vh.w = hit->cluster()->sizeY();
+      vhits.push_back(vh);
+    }
+
+    creco->setNValidPixelHits(nPxlHits);
+
+    // append cluster compatibility  for each z-position
+    for(double z0 = minZ_; z0 <= maxZ_; z0 += zStep_)
+    { 
+      ContainedHits c = getContainedHits(vhits, z0);
+      creco->append(c.z0, c.nHit, c.chi);
+    }
+
+  }
+  iEvent.put(creco);
+
+}
+
+
+ClusterCompatibilityProducer::ContainedHits ClusterCompatibilityProducer::getContainedHits(const std::vector<VertexHit> &hits, double z0) const
+{
+
+  // Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
+  int n = 0;
+  double chi = 0.;
+
+  for(std::vector<VertexHit>::const_iterator hit = hits.begin(); hit!= hits.end(); hit++) {
+    // the calculation of the predicted cluster width p was 
+    // marked 'FIXME' in the HLTPixelClusterShapeFilter. It should
+    // be revisited but is retained as it was for compatibility with the 
+    // older filter.
+    double p = 2 * fabs(hit->z - z0)/hit->r + 0.5; 
+    if(fabs(p - hit->w) <= 1.) {                   
+      chi += fabs(p - hit->w);
+      n++;
+    }
+  }
+  ClusterCompatibilityProducer::ContainedHits output;
+  output.z0 = z0;
+  output.nHit = n;
+  output.chi = chi;
+  return output;
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(ClusterCompatibilityProducer);

RecoHI/HiCentralityAlgos/python/HiClusterCompatibility_cfi.py

@@ -0,0 +1,8 @@
+import FWCore.ParameterSet.Config as cms
+
+hiClusterCompatibility = cms.EDProducer("ClusterCompatibilityProducer",
+   inputTag      = cms.InputTag( "siPixelRecHits" ),
+   minZ          = cms.double(-40.0),
+   maxZ          = cms.double(40.05),
+   zStep         = cms.double(0.2)
+)

RecoHI/HiCentralityAlgos/python/RecoHiCentrality_EventContent_cff.py

@@ -1,13 +1,16 @@
 import FWCore.ParameterSet.Config as cms
 
 RecoHiCentralityFEVT = cms.PSet(
-    outputCommands = cms.untracked.vstring('keep recoCentrality*_hiCentrality_*_*')
+    outputCommands = cms.untracked.vstring('keep recoCentrality*_hiCentrality_*_*',
+                                           'keep recoClusterCompatibility*_hiClusterCompatibility_*_*')
     )
 
 RecoHiCentralityRECO = cms.PSet(
-    outputCommands = cms.untracked.vstring('keep recoCentrality*_hiCentrality_*_*')
+    outputCommands = cms.untracked.vstring('keep recoCentrality*_hiCentrality_*_*',
+                                           'keep recoClusterCompatibility*_hiClusterCompatibility_*_*')
     )
 
 RecoHiCentralityAOD = cms.PSet(
-    outputCommands = cms.untracked.vstring('keep recoCentrality*_hiCentrality_*_*')
+    outputCommands = cms.untracked.vstring('keep recoCentrality*_hiCentrality_*_*',
+                                           'keep recoClusterCompatibility*_hiClusterCompatibility_*_*')
     )