Merge pull request #5627 from lgray/HGC-Linking

HGC Linking Incrementals

RecoParticleFlow/PFClusterProducer/interface/PFRecHitCaloNavigator.h

@@ -41,84 +41,113 @@ class PFRecHitCaloNavigator : public PFRecHitNavigatorBase {
 
  virtual ~PFRecHitCaloNavigator() { if(!ownsTopo) { topology_.release(); } }
 
-  void associateNeighbours(reco::PFRecHit& hit,std::auto_ptr<reco::PFRecHitCollection>& hits,edm::RefProd<reco::PFRecHitCollection>& refProd) {
-    constexpr unsigned short halfDim = 1; // range in each dimension
-    constexpr unsigned short dimSize = 3;
-    constexpr CaloDirection directions[6] = { NORTH, SOUTH, EAST, WEST, UP, DOWN }; // max available directions
-
-    DET detid( hit.detId() );
-    if( detid.null() ) return;
-    // this is a ripoff of getWindow() from CaloSubDetTopology
-    NeighbourInfo cellsInWindow;
-    DET tmpId(detid);
-    NeighbourInfo fringe;
-    fringe.emplace_back(tmpId,std::make_tuple(0,0,0));
-
-    std::vector<std::pair<CellInfo,Coordinate> > visited_cells;
-    visited_cells.resize(std::pow(dimSize,DIM));
-
-    while (fringe.size() > 0) {
-      NeighbourInfo::value_type cur = fringe.back();
-      fringe.pop_back();
-      // check all 2*DIM neighbours (in case of 2D only check NSEW)
-      for (unsigned dirnum = 0; dirnum < 2*DIM; ++dirnum) {
-	Coordinate neighbour = getNeighbourIndex(cur.second,directions[dirnum]);
-	//If outside the window range
-	if ( std::get<0>(neighbour) < -halfDim ||
-	     std::get<0>(neighbour) > halfDim ||
-	     std::get<1>(neighbour) < -halfDim ||
-	     std::get<1>(neighbour) > halfDim ||
-	     std::get<2>(neighbour) < -halfDim ||
-	     std::get<2>(neighbour) > halfDim )
-	  continue;
-
-	//Found integer index in the matrix
-	unsigned int_index = ( std::get<0>(neighbour) + halfDim +
-			       dimSize * (std::get<1>(neighbour) + halfDim ) +
-			       (DIM==3)*dimSize*dimSize * ( std::get<2>(neighbour) + halfDim ) );	
-	assert(int_index < visited_cells.size());
-
-	// check whether we have seen this neighbour already
-	if (visited_cells[int_index].first.visited)
-	  // we have seen this one already
-	  continue;
-
-	// a new cell, get the DetId of the neighbour, mark it
-	// as visited and add it to the fringe
-	visited_cells[int_index].first.visited = true;
-	std::vector<DetId> neighbourCells = getNeighbours(cur.first,directions[dirnum]);
-
-	if ( neighbourCells.size() == 1 ) {
-	  visited_cells[int_index].first.cell = neighbourCells[0];
-	  visited_cells[int_index].second = neighbour;
-	} else if ( neighbourCells.size() == 0 ) {
-	  visited_cells[int_index].first.cell = DetId(0);
-	  visited_cells[int_index].second = neighbour;
-	} else {
-	  throw cms::Exception("getWindowError") << "Not supported subdetector for getWindow method";
-	}
-
-	if (!visited_cells[int_index].first.cell.null())
-	  fringe.emplace_back(visited_cells[int_index].first.cell,neighbour);
-
-      } // loop over all possible directions
-    } // while some cells are left on the fringe
-
-    for (unsigned int i=0; i<visited_cells.size(); i++) {
-      if (!visited_cells[i].first.cell.null()) {
-	const DetId& id = visited_cells[i].first.cell;
-	const Coordinate& coord = visited_cells[i].second;
-	if( std::get<0>(coord) != 0 ||
-	    std::get<1>(coord) != 0 ||
-	    std::get<2>(coord) != 0 ) {
-	  associateNeighbour( id,hit,hits,refProd,
-			      std::get<0>(coord),
-			      std::get<1>(coord),
-			      std::get<2>(coord) );
-	}
-      }
-    } 
-  }
+ void associateNeighbours(reco::PFRecHit& hit, 
+			  std::auto_ptr<reco::PFRecHitCollection>& hits,
+			  const DetIdToHitIdx& hitmap,
+			  edm::RefProd<reco::PFRecHitCollection>& refProd) override {
+   auto visited_cells = std::move( getWindow( hit ) );
+   for (unsigned int i=0; i<visited_cells.size(); i++) {
+     if (!visited_cells[i].first.cell.null()) {
+       const DetId& id = visited_cells[i].first.cell;
+       const Coordinate& coord = visited_cells[i].second;
+       if( std::get<0>(coord) != 0 ||
+	   std::get<1>(coord) != 0 ||
+	   std::get<2>(coord) != 0 ) {
+	 associateNeighbour( id,hit,hits,hitmap,refProd,
+			     std::get<0>(coord),
+			     std::get<1>(coord),
+			     std::get<2>(coord) );
+       }
+     }
+   } 
+ }
+
+ void associateNeighbours(reco::PFRecHit& hit,
+			  std::auto_ptr<reco::PFRecHitCollection>& hits,
+			  edm::RefProd<reco::PFRecHitCollection>& refProd) {
+   auto visited_cells = std::move( getWindow( hit ) );
+   for (unsigned int i=0; i<visited_cells.size(); i++) {
+     if (!visited_cells[i].first.cell.null()) {
+       const DetId& id = visited_cells[i].first.cell;
+       const Coordinate& coord = visited_cells[i].second;
+       if( std::get<0>(coord) != 0 ||
+	   std::get<1>(coord) != 0 ||
+	   std::get<2>(coord) != 0 ) {
+	 associateNeighbour( id,hit,hits,refProd,
+			     std::get<0>(coord),
+			     std::get<1>(coord),
+			     std::get<2>(coord) );
+       }
+     }
+   }
+ }
+
+
+ std::vector<std::pair<CellInfo,Coordinate> > getWindow(reco::PFRecHit& hit) {
+   constexpr unsigned short halfDim = 1; // range in each dimension
+   constexpr unsigned short dimSize = 3;
+   constexpr CaloDirection directions[6] = { NORTH, SOUTH, EAST, WEST, UP, DOWN }; // max available directions
+
+   DET detid( hit.detId() );
+   if( detid.null() ) return std::vector<std::pair<CellInfo,Coordinate> >();
+   // this is a ripoff of getWindow() from CaloSubDetTopology
+   NeighbourInfo cellsInWindow;
+   DET tmpId(detid);
+   NeighbourInfo fringe;
+   fringe.emplace_back(tmpId,std::make_tuple(0,0,0));
+
+   std::vector<std::pair<CellInfo,Coordinate> > visited_cells;
+   visited_cells.resize(std::pow(dimSize,DIM));
+
+   while (fringe.size() > 0) {
+     NeighbourInfo::value_type cur = fringe.back();
+     fringe.pop_back();
+     // check all 2*DIM neighbours (in case of 2D only check NSEW)
+     for (unsigned dirnum = 0; dirnum < 2*DIM; ++dirnum) {
+       Coordinate neighbour = getNeighbourIndex(cur.second,directions[dirnum]);
+       //If outside the window range
+       if ( std::get<0>(neighbour) < -halfDim ||
+	    std::get<0>(neighbour) > halfDim ||
+	    std::get<1>(neighbour) < -halfDim ||
+	    std::get<1>(neighbour) > halfDim ||
+	    std::get<2>(neighbour) < -halfDim ||
+	    std::get<2>(neighbour) > halfDim )
+	 continue;
+
+       //Found integer index in the matrix
+       unsigned int_index = ( std::get<0>(neighbour) + halfDim +
+			      dimSize * (std::get<1>(neighbour) + halfDim ) +
+			      (DIM==3)*dimSize*dimSize * ( std::get<2>(neighbour) + halfDim ) );	
+       assert(int_index < visited_cells.size());
+
+       // check whether we have seen this neighbour already
+       if (visited_cells[int_index].first.visited)
+	 // we have seen this one already
+	 continue;
+
+       // a new cell, get the DetId of the neighbour, mark it
+       // as visited and add it to the fringe
+       visited_cells[int_index].first.visited = true;
+       std::vector<DetId> neighbourCells = getNeighbours(cur.first,directions[dirnum]);
+
+       if ( neighbourCells.size() == 1 ) {
+	 visited_cells[int_index].first.cell = neighbourCells[0];
+	 visited_cells[int_index].second = neighbour;
+       } else if ( neighbourCells.size() == 0 ) {
+	 visited_cells[int_index].first.cell = DetId(0);
+	 visited_cells[int_index].second = neighbour;
+       } else {
+	 throw cms::Exception("getWindowError") << "Not supported subdetector for getWindow method";
+       }
+
+       if (!visited_cells[int_index].first.cell.null())
+	 fringe.emplace_back(visited_cells[int_index].first.cell,neighbour);
+
+     } // loop over all possible directions
+   } // while some cells are left on the fringe
+
+   return visited_cells;
+ }
 
  protected:
   std::unique_ptr<const TOPO> topology_;

RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h

@@ -21,20 +21,27 @@
 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
 #include "Geometry/Records/interface/CaloGeometryRecord.h"
 
+#include <unordered_map>
 
 class PFRecHitNavigatorBase {
  public:
+  typedef std::unordered_map<unsigned,unsigned> DetIdToHitIdx;
+
   PFRecHitNavigatorBase() {}
   PFRecHitNavigatorBase(const edm::ParameterSet& iConfig) {}
 
   virtual ~PFRecHitNavigatorBase() {}
 
   virtual void beginEvent(const edm::EventSetup&)=0;
   virtual void associateNeighbours(reco::PFRecHit&,std::auto_ptr<reco::PFRecHitCollection>&,edm::RefProd<reco::PFRecHitCollection>&)=0;
+  virtual void associateNeighbours(reco::PFRecHit&,
+				   std::auto_ptr<reco::PFRecHitCollection>&,
+				   const DetIdToHitIdx&,
+				   edm::RefProd<reco::PFRecHitCollection>&) {};
 
 
  protected:
-
+  // assumes sorted
   void associateNeighbour(const DetId& id, reco::PFRecHit& hit,std::auto_ptr<reco::PFRecHitCollection>& hits,edm::RefProd<reco::PFRecHitCollection>& refProd,short eta, short phi,short depth) {
     const reco::PFRecHit temp(id,PFLayer::NONE,0.0,math::XYZPoint(0,0,0),math::XYZVector(0,0,0),std::vector<math::XYZPoint>());
     auto found_hit = std::lower_bound(hits->begin(),hits->end(),
@@ -47,7 +54,20 @@ class PFRecHitNavigatorBase {
       hit.addNeighbour(eta,phi,depth,reco::PFRecHitRef(refProd,std::distance(hits->begin(),found_hit)));
     }    
   }
-
+  // map to indices is provided
+  void associateNeighbour(const DetId& id, 
+			  reco::PFRecHit& hit,
+			  std::auto_ptr<reco::PFRecHitCollection>& hits,
+			  const DetIdToHitIdx& hitmap,
+			  edm::RefProd<reco::PFRecHitCollection>& refProd,
+			  short eta, short phi,short depth) {
+    auto found_hit = hitmap.find(id.rawId());
+    if( found_hit != hitmap.end() ) {
+      auto begin = hits->begin();
+      auto hit_pos = begin + found_hit->second;
+      hit.addNeighbour(eta,phi,depth,reco::PFRecHitRef(refProd,std::distance(begin,hit_pos)));
+    }    
+  }
 
 };
 

RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducer.cc

@@ -7,7 +7,8 @@ namespace {
   }
 }
 
- PFRecHitProducer:: PFRecHitProducer(const edm::ParameterSet& iConfig)
+PFRecHitProducer:: PFRecHitProducer(const edm::ParameterSet& iConfig):
+  _useHitMap(iConfig.getUntrackedParameter<bool>("useHitMap",false))
 {
 
   produces<reco::PFRecHitCollection>();
@@ -52,14 +53,26 @@ void
      creators_.at(i)->importRecHits(out,cleaned,iEvent,iSetup);
    }
 
-   std::sort(out->begin(),out->end(),sortByDetId);
+   if(!_useHitMap) std::sort(out->begin(),out->end(),sortByDetId);
+
+   reco::PFRecHitCollection& outprod = *out;
+   PFRecHitNavigatorBase::DetIdToHitIdx hitmap(outprod.size());
+   if( _useHitMap ) {
+     for( unsigned i = 0 ; i < outprod.size(); ++i ) {
+       hitmap[outprod[i].detId()] = i;
+     }
+   }
 
    //create a refprod here
    edm::RefProd<reco::PFRecHitCollection> refProd = 
      iEvent.getRefBeforePut<reco::PFRecHitCollection>();
 
-   for( unsigned int i=0;i<out->size();++i) {
-     navigator_->associateNeighbours(out->at(i),out,refProd);
+   for( unsigned int i=0;i<outprod.size();++i) {
+     if( _useHitMap ) {
+       navigator_->associateNeighbours(outprod[i],out,hitmap,refProd);
+     } else {
+       navigator_->associateNeighbours(outprod[i],out,refProd);
+     }
    }
 
    iEvent.put(out,"");

RecoParticleFlow/PFClusterProducer/plugins/PFRecHitProducer.h

@@ -26,6 +26,7 @@ class PFRecHitProducer : public edm::EDProducer {
       static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
 
    private:
+      const bool _useHitMap;
       virtual void produce(edm::Event&, const edm::EventSetup&) override;
       std::vector<std::unique_ptr<PFRecHitCreatorBase> > creators_;
       std::unique_ptr<PFRecHitNavigatorBase> navigator_;

RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCEE_cfi.py

@@ -26,8 +26,8 @@
 _manqiArborClusterizer_HGCEE = cms.PSet(
     algoName = cms.string("SimpleArborClusterizer"), 
     # use basic pad sizes in HGCEE
-    cellSize = cms.double(10.0),
-    layerThickness = cms.double(16.0),
+    cellSize = cms.double(16.0),
+    layerThickness = cms.double(18.0),
     distSeedForMerge = cms.double(20.0),
     killNoiseClusters = cms.bool(True),
     maxNoiseClusterSize = cms.uint32(3),

RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCHEB_cfi.py

@@ -23,8 +23,8 @@
 _manqiArborClusterizer_HGCHEB = cms.PSet(
     algoName = cms.string("SimpleArborClusterizer"), 
     # use basic pad sizes in HGCEE
-    cellSize = cms.double(30.0),
-    layerThickness = cms.double(55.0),
+    cellSize = cms.double(35.0),
+    layerThickness = cms.double(65.0),
     distSeedForMerge = cms.double(20.0),
     killNoiseClusters = cms.bool(True),
     maxNoiseClusterSize = cms.uint32(3),

RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGCHEF_cfi.py

@@ -23,9 +23,9 @@
 _manqiArborClusterizer_HGCHEF = cms.PSet(
     algoName = cms.string("SimpleArborClusterizer"), 
     # use basic pad sizes in HGCEE
-    cellSize = cms.double(10.0),
-    layerThickness = cms.double(45.0),
-    distSeedForMerge = cms.double(20.0),
+    cellSize = cms.double(15.0),
+    layerThickness = cms.double(55.0),
+    distSeedForMerge = cms.double(30.0),
     killNoiseClusters = cms.bool(True),
     maxNoiseClusterSize = cms.uint32(3),
     thresholdsByDetector = cms.VPSet( )

RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCEE_cfi.py

@@ -1,6 +1,7 @@
 import FWCore.ParameterSet.Config as cms
 
 particleFlowRecHitHGCEE = cms.EDProducer("PFRecHitProducer",
+    useHitMap = cms.untracked.bool(True),
     navigator = cms.PSet(
         name = cms.string("PFRecHitHGCEENavigator"),
         topologySource = cms.string("HGCalEESensitive"),

RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCHEB_cfi.py

@@ -1,6 +1,7 @@
 import FWCore.ParameterSet.Config as cms
 
 particleFlowRecHitHGCHEB = cms.EDProducer("PFRecHitProducer",
+    useHitMap = cms.untracked.bool(True),
     navigator = cms.PSet(
         name = cms.string("PFRecHitHGCHENavigator"),
         topologySource = cms.string("HGCalHEScintillatorSensitive"),

RecoParticleFlow/PFClusterProducer/python/particleFlowRecHitHGCHEF_cfi.py

@@ -1,6 +1,7 @@
 import FWCore.ParameterSet.Config as cms
 
 particleFlowRecHitHGCHEF = cms.EDProducer("PFRecHitProducer",
+    useHitMap = cms.untracked.bool(True),
     navigator = cms.PSet(
         name = cms.string("PFRecHitHGCHENavigator"),
         topologySource = cms.string("HGCalHESiliconSensitive"),

RecoParticleFlow/PFClusterProducer/src/Arbor.cc

@@ -286,7 +286,7 @@ void BuildInitLink()
 	  const auto& PosB = cleanedHits[found[j0].data];
 	  PosDiffAB = PosA - PosB;
 
-	  if( PosDiffAB.Mag2() < InitLinkThreshold2 ) // || ( PosDiffAB.Mag() < 1.6*InitLinkThreshold && PosDiffAB.Dot(PosB) < 0.9*PosDiffAB.Mag()*PosB.Mag() )  )	//Distance threshold to be optimized - should also depends on Geometry
+	  if( std::abs(PosDiffAB.Z()) > 1e-3 && PosDiffAB.Mag2() < InitLinkThreshold2 ) // || ( PosDiffAB.Mag() < 1.6*InitLinkThreshold && PosDiffAB.Dot(PosB) < 0.9*PosDiffAB.Mag()*PosB.Mag() )  )	//Distance threshold to be optimized - should also depends on Geometry
 	    {
 	      std::pair<int, int> a_Link;
 	      if( PosA.Mag2() > PosB.Mag2() )
@@ -461,7 +461,8 @@ void LinkIteration()	//Energy corrections, semi-local correction
       PosB = cleanedHits[found[j1].data];
       DiffPosAB = PosB - PosA; 
 
-      if( DiffPosAB.Mag2() < IterLinkThreshold2 && 
+      if( std::abs(DiffPosAB.Z()) > 1e-3 &&
+	  DiffPosAB.Mag2() < IterLinkThreshold2 && 
 	  DiffPosAB.Mag2() > InitLinkThreshold2 && 
 	  DiffPosAB.Angle(RefDir[i1]) < 0.8 ) {