Add TICL candidates and matching to rechits

CommonTools/RecoAlgos/BuildFile.xml

@@ -15,6 +15,7 @@
 <use name="Geometry/TrackerGeometryBuilder"/>
 <use name="SimDataFormats/TrackingAnalysis"/>
 <use name="fastjet"/>
+<use name="RecoLocalCalo/HGCalRecAlgos" />
 <export>
   <lib name="1"/>
 </export>

CommonTools/RecoAlgos/plugins/RecHitToTICLCandAssociationProducer.cc

@@ -0,0 +1,180 @@
+//
+// system include files
+#include <memory>
+#include <string>
+
+// user include files
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/Framework/interface/ESHandle.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/HGCalReco/interface/TICLCandidate.h"
+
+#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+
+#include "DataFormats/Common/interface/Association.h"
+#include "DataFormats/Common/interface/AssociationMap.h"
+#include "DataFormats/Common/interface/OneToManyWithQualityGeneric.h"
+#include "DataFormats/Common/interface/RefToBase.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+
+#include "DataFormats/Math/interface/deltaR.h"
+
+#include "FWCore/Utilities/interface/transform.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include <set>
+
+//helper
+#include <numeric>
+#include <algorithm>
+
+template <typename T>
+std::vector<size_t> argsort(const std::vector<T> &v) {
+    using namespace std;
+    vector<size_t> idx(v.size());
+    iota(idx.begin(), idx.end(), 0);
+    stable_sort(idx.begin(), idx.end(),
+            [&v](size_t i1, size_t i2) {return v[i1] < v[i2];});
+    return idx;
+}
+
+
+//
+// class declaration
+//
+typedef std::vector<TICLCandidate> TICLCandidateCollection;
+typedef edm::Ref<TICLCandidateCollection> TICLCandidateRef;
+
+typedef edm::Ref<reco::CaloClusterCollection> CaloClusterRef;
+
+/* doesn't compile
+typedef edm::AssociationMap<edm::OneToManyWithQualityGeneric<
+    HGCRecHitCollection, TICLCandidate, float> > RecHitToTICLCands;
+typedef edm::AssociationMap<edm::OneToManyWithQualityGeneric<
+    reco::CaloCluster, TICLCandidate, float> > RecHitToTICLCands;
+*/
+
+class RecHitToTICLCandidateAssociationProducer : public edm::stream::EDProducer<> {
+public:
+  explicit RecHitToTICLCandidateAssociationProducer(const edm::ParameterSet&);
+  ~RecHitToTICLCandidateAssociationProducer() override;
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+
+  edm::EDGetTokenT<reco::CaloClusterCollection> layerClustersToken_;
+  edm::EDGetTokenT<TICLCandidateCollection> ticlCandToken_;
+  edm::EDGetTokenT<HGCRecHitCollection> recHitToken_;
+  hgcal::RecHitTools recHitTools_;
+};
+
+RecHitToTICLCandidateAssociationProducer::RecHitToTICLCandidateAssociationProducer(const edm::ParameterSet& pset)
+    : layerClustersToken_(consumes<std::vector<reco::CaloCluster>>(pset.getParameter<edm::InputTag>("layerClusters"))),
+      ticlCandToken_(consumes<TICLCandidateCollection>(pset.getParameter<edm::InputTag>("ticlCandidates"))),
+      recHitToken_(consumes<HGCRecHitCollection>(pset.getParameter<edm::InputTag>("caloRecHits"))){
+
+    produces<edm::Association<TICLCandidateCollection>>();
+}
+
+RecHitToTICLCandidateAssociationProducer::~RecHitToTICLCandidateAssociationProducer() {}
+
+// ------------ method called to produce the data  ------------
+void RecHitToTICLCandidateAssociationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::ESHandle<CaloGeometry> geom;
+  iSetup.get<CaloGeometryRecord>().get(geom);
+  recHitTools_.setGeometry(*geom);
+
+  edm::Handle<TICLCandidateCollection> ticlCands;
+  iEvent.getByToken(ticlCandToken_, ticlCands);
+
+  edm::Handle<HGCRecHitCollection> recHits;
+  iEvent.getByToken(recHitToken_, recHits);
+
+  //convenient access. indices will refer to this later
+  std::vector<const HGCRecHit*> allrechits;
+
+  //for convenience
+  std::unordered_map<DetId, size_t> detid_to_rh_index;
+
+  //hit energies (copied) that will be depleted during association
+  size_t rhindex = 0;
+  for (const auto& rh : *recHits) {
+    detid_to_rh_index[rh.detid()] = rhindex;
+    rhindex++;
+    allrechits.push_back(&rh);
+  }
+
+  edm::Handle<std::vector<reco::CaloCluster>> layerClusterCollection;
+  iEvent.getByToken(layerClustersToken_, layerClusterCollection);
+
+  //assigns a ticl candidate to each rechit
+  std::vector<int> hittotc(allrechits.size(), -1);
+
+  for (size_t tcidx = 0; tcidx < ticlCands->size(); tcidx++) {
+    TICLCandidateRef ticlCand(ticlCands, tcidx);
+    for (const auto& trackster : ticlCand->tracksters()) {
+      for (int lcidx : trackster->vertices()) {
+        //this is the simple part, use equal sharing
+        float targetfrac = 1. / trackster->vertex_multiplicity(lcidx);
+        const auto& lc = layerClusterCollection->at(lcidx);
+
+        float targetEnergy = targetfrac * lc.energy();
+
+        auto hafv = lc.hitsAndFractions();
+
+        //calculate hit distances^2
+        std::vector<float> distances;
+        for (const auto& hwf : hafv) {
+          auto pos = recHitTools_.getPosition(hwf.first);
+          float dist =
+              reco::deltaR2(pos.eta(), pos.phi(), trackster->barycenter().eta(), trackster->barycenter().phi());
+          distances.push_back(dist);
+        }
+
+        //re-order by distance^2
+        auto distidx = argsort(distances);
+
+        //deplete hit energies starting from closest until target energy reached
+        for (const auto& hidx : distidx) {
+          //float dist = distances.at(hidx);//doesn'tcidx matter here
+          auto rhdetid = hafv.at(hidx).first;
+          size_t rhidx = detid_to_rh_index[rhdetid];
+          float frachitenergy = allrechits.at(rhidx)->energy() * hafv.at(hidx).second;
+
+          //if there is still energy left in the hit, even if it's just a bit use it
+          if (hittotc.at(rhidx)<0) {//not assigned yet
+            //less target energy remaining
+            targetEnergy -= frachitenergy;
+            //associate
+            hittotc.at(rhidx) = tcidx;
+          }
+          if (targetEnergy <= 0)  //done with layer cluster
+            break;
+        }
+      }
+    }
+  }
+
+  std::cout << "RecHitToTICLCandidateAssociationProducer done" << std::endl;
+
+  auto assoc = std::make_unique<edm::Association<TICLCandidateCollection>>(ticlCands);
+  edm::Association<TICLCandidateCollection>::Filler filler(*assoc);
+  filler.insert(recHits, hittotc.begin(), hittotc.end());
+  filler.fill();
+  iEvent.put(std::move(assoc));
+
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(RecHitToTICLCandidateAssociationProducer);

DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc

@@ -11,6 +11,7 @@
 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
 #include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/HGCalReco/interface/TICLCandidate.h"
 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
 #include "FWCore/Framework/interface/MakerMacros.h"
@@ -23,17 +24,23 @@ typedef ObjectIndexFromAssociationTableProducer<edm::View<CaloRecHit>, reco::PFC
     CaloRecHitToPFCandIndexTableProducer;
 typedef ObjectIndexFromAssociationTableProducer<edm::View<CaloRecHit>, SimClusterCollection>
     CaloRecHitToBestSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer<edm::View<CaloRecHit>, reco::CaloClusterCollection>
+    CaloRecHitToBestLayerClusterIndexTableProducer;
 typedef ObjectIndexFromAssociationTableProducer<SimClusterCollection, CaloParticleCollection>
     SimClusterToCaloParticleIndexTableProducer;
 typedef ObjectIndexFromAssociationTableProducer<SimClusterCollection, SimClusterCollection>
     SimClusterToSimClusterIndexTableProducer;
 typedef ObjectIndexFromAssociationTableProducer<edm::View<CaloRecHit>, PFTruthParticleCollection>
     CaloRecHitToPFTruthParticleIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer<edm::View<CaloRecHit>, std::vector<TICLCandidate>>
+    RecHitToTICLCandidateIndexTableProducer;
 
 DEFINE_FWK_MODULE(SimTrackToSimClusterIndexTableProducer);
 DEFINE_FWK_MODULE(CaloHitToSimClusterIndexTableProducer);
 DEFINE_FWK_MODULE(SimClusterToCaloParticleIndexTableProducer);
 DEFINE_FWK_MODULE(SimClusterToSimClusterIndexTableProducer);
 DEFINE_FWK_MODULE(CaloRecHitToPFCandIndexTableProducer);
 DEFINE_FWK_MODULE(CaloRecHitToBestSimClusterIndexTableProducer);
+DEFINE_FWK_MODULE(CaloRecHitToBestLayerClusterIndexTableProducer);
 DEFINE_FWK_MODULE(CaloRecHitToPFTruthParticleIndexTableProducer);
+DEFINE_FWK_MODULE(RecHitToTICLCandidateIndexTableProducer);

DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py

@@ -2,6 +2,7 @@
 from PhysicsTools.NanoAOD.common_cff import Var,P3Vars
 from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import *
 
+
 hgcRecHitsToSimClusters = cms.EDProducer("SimClusterRecHitAssociationProducer",
     caloRecHits = cms.VInputTag("hgcRecHits"),
     simClusters = cms.InputTag("mix:MergedCaloTruth"),

DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py

@@ -49,20 +49,6 @@
     docString = cms.string("PFTICLCand with most associated energy in RecHit DetId")
 )
 
-hgcRecHitsToLayerClusters = cms.EDProducer("RecHitToLayerClusterAssociationProducer",
-    caloRecHits = cms.VInputTag("hgcRecHits"),
-    layerClusters = cms.InputTag("hgcalLayerClusters"),
-)
-
-hgcRecHitsToLayerClusterTable = cms.EDProducer("HGCRecHitToLayerClusterIndexTableProducer",
-    cut = hgcRecHitsTable.cut,
-    src = hgcRecHitsTable.src,
-    objName = hgcRecHitsTable.name,
-    branchName = cms.string("LayerCluster"),
-    objMap = cms.InputTag("hgcRecHitsToLayerClusters:hgcRecHitsToLayerCluster"),
-    docString = cms.string("LayerCluster assigned largest RecHit fraction")
-)
-
 hgcRecHitsPositionTable = cms.EDProducer("HGCRecHitPositionTableProducer",
     src = hgcRecHitsTable.src,
     cut = hgcRecHitsTable.cut, 
@@ -76,7 +62,5 @@
                 +hgcRecHitsToPFCandTable
                 +hgcRecHitsToPFTICLCands
                 +hgcRecHitsToPFTICLCandTable
-                +hgcRecHitsToLayerClusters 
-                +hgcRecHitsToLayerClusterTable
                 +hgcRecHitsPositionTable
 )

DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py

@@ -1,5 +1,6 @@
 import FWCore.ParameterSet.Config as cms
 from PhysicsTools.NanoAOD.common_cff import P3Vars,Var
+from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import hgcRecHitsTable
 
 layerClusterTable = cms.EDProducer("SimpleCaloClusterFlatTableProducer",
     src = cms.InputTag("hgcalLayerClusters"),
@@ -29,5 +30,31 @@
     docString = cms.string("Index of SimCluster matched to LayerCluster")
 )
 
-layerClusterTables = cms.Sequence(layerClusterTable+layerClusterToSimClusterTable)
+hgcRecHitsToLayerClusters = cms.EDProducer("RecHitToLayerClusterAssociationProducer",
+    caloRecHits = cms.VInputTag("hgcRecHits"),
+    layerClusters = cms.InputTag("hgcalLayerClusters"),
+)
+
+hgcRecHitsToLayerClusterTable = cms.EDProducer("HGCRecHitToLayerClusterIndexTableProducer",
+    cut = hgcRecHitsTable.cut,
+    src = hgcRecHitsTable.src,
+    objName = hgcRecHitsTable.name,
+    branchName = cms.string("LayerCluster"),
+    objMap = cms.InputTag("hgcRecHitsToLayerClusters:hgcRecHitsToLayerCluster"),
+    docString = cms.string("LayerCluster assigned largest RecHit fraction")
+)
+
+# For now there is only ever a 1 to 1 rechit to lc match
+hgcRecHitsToBestLayerClusterTable = cms.EDProducer("CaloRecHitToBestLayerClusterIndexTableProducer",
+    cut = hgcRecHitsTable.cut,
+    src = hgcRecHitsTable.src,
+    objName = hgcRecHitsTable.name,
+    branchName = layerClusterTable.name,
+    objMap = cms.InputTag("hgcRecHitsToLayerClusters:hgcRecHitsToBestLayerCluster"),
+    docString = cms.string("LayerCluster to which the RecHit is associated")
+)
+
+
+layerClusterTables = cms.Sequence(layerClusterTable+layerClusterToSimClusterTable+
+    hgcRecHitsToLayerClusters+hgcRecHitsToLayerClusterTable+hgcRecHitsToBestLayerClusterTable)
 

DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py

@@ -10,6 +10,7 @@
 from DPGAnalysis.HGCalNanoAOD.simClusters_cff import *
 from DPGAnalysis.HGCalNanoAOD.layerClusters_cff import *
 from DPGAnalysis.HGCalNanoAOD.caloParticles_cff import *
+from DPGAnalysis.HGCalNanoAOD.ticlCandidates_cff import *
 from DPGAnalysis.TrackNanoAOD.trackSimHits_cff import *
 from DPGAnalysis.TrackNanoAOD.trackingParticles_cff import *
 from DPGAnalysis.TrackNanoAOD.tracks_cff import *
@@ -28,7 +29,6 @@
 
 nanoHGCMLSequence = cms.Sequence(nanoMetadata+
     genVertexTable+genVertexT0Table+genParticleTable+
-    layerClusterTables+
     simTrackTables+hgcSimHitsSequence+trackerSimHitTables+
     simClusterTables+
     trackingParticleTables+
@@ -37,8 +37,10 @@
 
 nanoHGCMLRecoSequence = cms.Sequence(hgcRecHitsSequence+
         hgcRecHitSimAssociationSequence+
+        layerClusterTables+
         pfCandTable+pfTruth+
         pfTICLCandTable+
+        ticlTables+
         trackTables+
         trackSCAssocTable)
 

DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py

@@ -55,11 +55,11 @@
 )
 
 hgcSimTruth = cms.EDProducer("simmerger",
-    useNLayers = cms.int32(1),
-    searchRadiusScale = cms.double(2.),
+    useNLayers = cms.int32(3),
+    searchRadiusScale = cms.double(2.0),
     clusterRadiusScale = cms.double(1.5),
-    mergeRadiusScale = cms.double(1.5),
-    energyContainment = cms.double(0.5),
+    mergeRadiusScale = cms.double(1.0),
+    energyContainment = cms.double(0.7),
     relOverlapDistance = cms.double(0.9)
                            )
 

DPGAnalysis/HGCalNanoAOD/python/ticlCandidates_cff.py

@@ -0,0 +1,36 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import hgcRecHitsTable
+
+ticlTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
+    src = cms.InputTag("ticlTrackstersMerge"),
+    cut = cms.string(""), 
+    name = cms.string("TICLCand"),
+    doc  = cms.string("TICL Candidates"),
+    singleton = cms.bool(False), # the number of entries is variable
+    extension = cms.bool(False), # this is the main table for the muons
+    variables = cms.PSet(CandVars,
+        Vtx_x = Var('vertex().x()', 'float', precision=14, doc='vertex x pos'),
+        Vtx_y = Var('vertex().y()', 'float', precision=14, doc='vertex y pos'),
+        Vtx_z = Var('vertex().z()', 'float', precision=14, doc='vertex z pos'),
+    )
+)
+
+hgcRecHitsToTiclCands = cms.EDProducer("RecHitToTICLCandidateAssociationProducer",
+    caloRecHits = hgcRecHitsTable.src,
+    layerClusters = cms.InputTag("hgcalLayerClusters"),
+    ticlCandidates = ticlTable.src,
+)
+
+hgcRecHitToTiclTable = cms.EDProducer("RecHitToTICLCandidateIndexTableProducer",
+    cut = hgcRecHitsTable.cut,
+    src = hgcRecHitsTable.src,
+    objName = hgcRecHitsTable.name,
+    branchName = ticlTable.name,
+    objMap = cms.InputTag("hgcRecHitsToTiclCands"),
+    docString = cms.string("Index of TICLCand containing RecHit")
+)
+
+
+ticlTables = cms.Sequence(ticlTable+hgcRecHitsToTiclCands+hgcRecHitToTiclTable)
+

DataFormats/HGCalReco/src/classes.h

@@ -5,3 +5,4 @@
 #include "DataFormats/HGCalReco/interface/TICLSeedingRegion.h"
 #include "DataFormats/HGCalReco/interface/TICLCandidate.h"
 #include "DataFormats/Common/interface/Wrapper.h"
+#include "DataFormats/Common/interface/Association.h"

DataFormats/HGCalReco/src/classes_def.xml

@@ -54,4 +54,8 @@
   <class name="std::vector<TICLCandidate>" />
   <class name="edm::Wrapper<TICLCandidate>" />
   <class name="edm::Wrapper<std::vector<TICLCandidate> >" />
+  <class name="edm::RefProd<vector<TICLCandidate> >"/>
+  <class name="edm::Wrapper<edm::RefProd<vector<TICLCandidate> >>"/>
+  <class name="edm::Association<std::vector<TICLCandidate>>"/>
+  <class name="edm::Wrapper<edm::Association<std::vector<TICLCandidate>> >"/>
 </lcgdict>