diff --git a/CommonTools/RecoAlgos/BuildFile.xml b/CommonTools/RecoAlgos/BuildFile.xml index 736edd2d2c27e..94a12a477243d 100644 --- a/CommonTools/RecoAlgos/BuildFile.xml +++ b/CommonTools/RecoAlgos/BuildFile.xml @@ -15,6 +15,7 @@ + diff --git a/CommonTools/RecoAlgos/plugins/RecHitToTICLCandAssociationProducer.cc b/CommonTools/RecoAlgos/plugins/RecHitToTICLCandAssociationProducer.cc new file mode 100644 index 0000000000000..341b329100ca3 --- /dev/null +++ b/CommonTools/RecoAlgos/plugins/RecHitToTICLCandAssociationProducer.cc @@ -0,0 +1,180 @@ +// +// system include files +#include +#include + +// 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 + +//helper +#include +#include + +template +std::vector argsort(const std::vector &v) { + using namespace std; + vector 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 TICLCandidateCollection; +typedef edm::Ref TICLCandidateRef; + +typedef edm::Ref CaloClusterRef; + +/* doesn't compile +typedef edm::AssociationMap > RecHitToTICLCands; +typedef edm::AssociationMap > 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 layerClustersToken_; + edm::EDGetTokenT ticlCandToken_; + edm::EDGetTokenT recHitToken_; + hgcal::RecHitTools recHitTools_; +}; + +RecHitToTICLCandidateAssociationProducer::RecHitToTICLCandidateAssociationProducer(const edm::ParameterSet& pset) + : layerClustersToken_(consumes>(pset.getParameter("layerClusters"))), + ticlCandToken_(consumes(pset.getParameter("ticlCandidates"))), + recHitToken_(consumes(pset.getParameter("caloRecHits"))){ + + produces>(); +} + +RecHitToTICLCandidateAssociationProducer::~RecHitToTICLCandidateAssociationProducer() {} + +// ------------ method called to produce the data ------------ +void RecHitToTICLCandidateAssociationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + edm::ESHandle geom; + iSetup.get().get(geom); + recHitTools_.setGeometry(*geom); + + edm::Handle ticlCands; + iEvent.getByToken(ticlCandToken_, ticlCands); + + edm::Handle recHits; + iEvent.getByToken(recHitToken_, recHits); + + //convenient access. indices will refer to this later + std::vector allrechits; + + //for convenience + std::unordered_map 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> layerClusterCollection; + iEvent.getByToken(layerClustersToken_, layerClusterCollection); + + //assigns a ticl candidate to each rechit + std::vector 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 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>(ticlCands); + edm::Association::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); diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc index 636dc103d127a..bdb1e69d7bbf8 100644 --- a/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc +++ b/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,12 +24,16 @@ typedef ObjectIndexFromAssociationTableProducer, reco::PFC CaloRecHitToPFCandIndexTableProducer; typedef ObjectIndexFromAssociationTableProducer, SimClusterCollection> CaloRecHitToBestSimClusterIndexTableProducer; +typedef ObjectIndexFromAssociationTableProducer, reco::CaloClusterCollection> + CaloRecHitToBestLayerClusterIndexTableProducer; typedef ObjectIndexFromAssociationTableProducer SimClusterToCaloParticleIndexTableProducer; typedef ObjectIndexFromAssociationTableProducer SimClusterToSimClusterIndexTableProducer; typedef ObjectIndexFromAssociationTableProducer, PFTruthParticleCollection> CaloRecHitToPFTruthParticleIndexTableProducer; +typedef ObjectIndexFromAssociationTableProducer, std::vector> + RecHitToTICLCandidateIndexTableProducer; DEFINE_FWK_MODULE(SimTrackToSimClusterIndexTableProducer); DEFINE_FWK_MODULE(CaloHitToSimClusterIndexTableProducer); @@ -36,4 +41,6 @@ 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); diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py index 8fae4b33520cf..740e6457c691d 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py +++ b/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"), diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py index f378a107d2ecd..d8307e35f1c27 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py +++ b/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 ) diff --git a/DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py b/DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py index 555337ce94d5d..a1746eb9c7e3e 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py +++ b/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) diff --git a/DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py b/DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py index 88b81e0dd03db..d6c183c81c60e 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py +++ b/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) diff --git a/DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py b/DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py index 0f2ae22766640..338be8c6a7dc4 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py +++ b/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) ) diff --git a/DPGAnalysis/HGCalNanoAOD/python/ticlCandidates_cff.py b/DPGAnalysis/HGCalNanoAOD/python/ticlCandidates_cff.py new file mode 100644 index 0000000000000..3620e9e62a72d --- /dev/null +++ b/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) + diff --git a/DataFormats/HGCalReco/src/classes.h b/DataFormats/HGCalReco/src/classes.h index d871bfb485a71..9ff7e83b6cbed 100644 --- a/DataFormats/HGCalReco/src/classes.h +++ b/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" diff --git a/DataFormats/HGCalReco/src/classes_def.xml b/DataFormats/HGCalReco/src/classes_def.xml index 9ac5b29b4f2e2..7098a1a8e5243 100644 --- a/DataFormats/HGCalReco/src/classes_def.xml +++ b/DataFormats/HGCalReco/src/classes_def.xml @@ -54,4 +54,8 @@ + + + + diff --git a/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc b/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc index ec393a5611ec0..5c5126177a7bb 100644 --- a/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc +++ b/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc @@ -917,7 +917,6 @@ bool merge_leafparent(Node* leafparent){ // Parent itself can be mergeable, if it has hits and is not a root if (leafparent->hasParent() && leafparent->hadHits()) mergeable.push_back(leafparent); // Start merging - std::cout << "while start" << std::endl; while(true){ bool didUpdateThisIteration = false; float currentMinDistanceSq = FLT_MAX; @@ -968,7 +967,6 @@ bool merge_leafparent(Node* leafparent){ } - std::cout << "while ends" << std::endl; // Make a string representation of the mergeable nodes for debugging std::string mergeableStr = ""; if(mergeable.size()){