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/LayerClusterToTICLCandAssociationProducer.cc b/CommonTools/RecoAlgos/plugins/LayerClusterToTICLCandAssociationProducer.cc new file mode 100644 index 0000000000000..d4d9e50178a62 --- /dev/null +++ b/CommonTools/RecoAlgos/plugins/LayerClusterToTICLCandAssociationProducer.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 > LayerClusterToTICLCands; +*/ + +class LayerClusterToTICLCandidateAssociationProducer : public edm::stream::EDProducer<> { +public: + explicit LayerClusterToTICLCandidateAssociationProducer(const edm::ParameterSet&); + ~LayerClusterToTICLCandidateAssociationProducer() override; + +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + edm::EDGetTokenT layerClustersToken_; + edm::EDGetTokenT ticlCandToken_; + edm::EDGetTokenT recHitToken_; + hgcal::RecHitTools recHitTools_; +}; + +LayerClusterToTICLCandidateAssociationProducer::LayerClusterToTICLCandidateAssociationProducer(const edm::ParameterSet& pset) + : layerClustersToken_(consumes>(pset.getParameter("layerClusters"))), + ticlCandToken_(consumes(pset.getParameter("ticlCandidates"))), + recHitToken_(consumes(pset.getParameter("caloRecHits"))){ + + produces>(); +} + +LayerClusterToTICLCandidateAssociationProducer::~LayerClusterToTICLCandidateAssociationProducer() {} + +// ------------ method called to produce the data ------------ +void LayerClusterToTICLCandidateAssociationProducer::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 << "LayerClusterToTICLCandidateAssociationProducer 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(LayerClusterToTICLCandidateAssociationProducer); diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py index 8fae4b33520cf..ee810f647d93d 100644 --- a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py +++ b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py @@ -2,6 +2,13 @@ from PhysicsTools.NanoAOD.common_cff import Var,P3Vars from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import * + +hgcRecHitsToTiclCands = cms.EDProducer("LayerClusterToTICLCandidateAssociationProducer", + caloRecHits = cms.InputTag("hgcRecHits"), + layerClusters = cms.InputTag("hgcalLayerClusters"), + ticlCandidates = cms.InputTag("ticlTrackstersMerge") +) + hgcRecHitsToSimClusters = cms.EDProducer("SimClusterRecHitAssociationProducer", caloRecHits = cms.VInputTag("hgcRecHits"), simClusters = cms.InputTag("mix:MergedCaloTruth"), @@ -92,6 +99,7 @@ hgcRecHitSimAssociationSequence = cms.Sequence(hgcRecHitsToSimClusters +hgcRecHitsToMergedSimClusters +hgcRecHitsToMergedDRSimClusters + +hgcRecHitsToTiclCands +simClusterRecEnergyTable +mergedSimClusterRecEnergyTable +hgcRecHitsToSimClusterTable diff --git a/IOMC/ParticleGuns/src/FlatEtaRangeNoTrackerGunProducer.cc b/IOMC/ParticleGuns/src/FlatEtaRangeNoTrackerGunProducer.cc index 9722ffc5683ca..f00a523710357 100644 --- a/IOMC/ParticleGuns/src/FlatEtaRangeNoTrackerGunProducer.cc +++ b/IOMC/ParticleGuns/src/FlatEtaRangeNoTrackerGunProducer.cc @@ -102,6 +102,7 @@ void edm::FlatEtaRangeNoTrackerGunProducer::produce(edm::Event& event, const edm } int particle_counter=0; + int unsuccessful_trials=0; std::vector previous_impacts; // shoot particles @@ -155,17 +156,21 @@ void edm::FlatEtaRangeNoTrackerGunProducer::produce(edm::Event& event, const edm bool next=false; for(const auto& previ: previous_impacts){ math::XYZPoint this_im; - if(reco::deltaR(previ.eta(),previ.phi(),vertexPos.eta(),vertexPos.phi()) < minDistDR_){ + if(reco::deltaR2(previ.eta(),previ.phi(),vertexPos.eta(),vertexPos.phi()) < minDistDR_*minDistDR_){ next=true; break; } } if(next){ i--; + unsuccessful_trials++; if(debug_) LogDebug("FlatEtaRangeNoTrackerGunProducer") << " : skipping too close particle" << std::endl; + if(unsuccessful_trials>20) + break; //stop here continue; } + unsuccessful_trials=0;//reset previous_impacts.push_back(math::XYZPoint(vertexPos.x(),vertexPos.y(),vertexPos.z())); }