diff --git a/CommonTools/RecoAlgos/plugins/BuildFile.xml b/CommonTools/RecoAlgos/plugins/BuildFile.xml
index 14976ed78b0e4..1391d81bd6a18 100644
--- a/CommonTools/RecoAlgos/plugins/BuildFile.xml
+++ b/CommonTools/RecoAlgos/plugins/BuildFile.xml
@@ -16,4 +16,5 @@
+
diff --git a/CommonTools/RecoAlgos/plugins/HGCRecHitCollectionMerger.cc b/CommonTools/RecoAlgos/plugins/HGCRecHitCollectionMerger.cc
new file mode 100644
index 0000000000000..1d9177f6c0c55
--- /dev/null
+++ b/CommonTools/RecoAlgos/plugins/HGCRecHitCollectionMerger.cc
@@ -0,0 +1,12 @@
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "CommonTools/UtilAlgos/interface/Merger.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/Common/interface/CloneTrait.h"
+#include "DataFormats/Common/interface/RefToBaseVector.h"
+#include "DataFormats/Common/interface/CopyPolicy.h"
+
+
+typedef Merger> HGCRecHitCollectionMerger;
+
+DEFINE_FWK_MODULE(HGCRecHitCollectionMerger);
diff --git a/CommonTools/RecoAlgos/plugins/RecHitToLayerClusterAssociationProducer.cc b/CommonTools/RecoAlgos/plugins/RecHitToLayerClusterAssociationProducer.cc
new file mode 100644
index 0000000000000..2b7a3b23fc31b
--- /dev/null
+++ b/CommonTools/RecoAlgos/plugins/RecHitToLayerClusterAssociationProducer.cc
@@ -0,0 +1,127 @@
+// 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 "DataFormats/CaloRecHit/interface/CaloRecHit.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.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 "FWCore/Utilities/interface/transform.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include
+
+//
+// class decleration
+//
+typedef std::pair IdxAndFraction;
+typedef edm::AssociationMap, float>> RecHitToLayerCluster;
+
+class RecHitToLayerClusterAssociationProducer : public edm::stream::EDProducer<> {
+public:
+ explicit RecHitToLayerClusterAssociationProducer(const edm::ParameterSet&);
+ ~RecHitToLayerClusterAssociationProducer() override;
+
+private:
+ void produce(edm::Event&, const edm::EventSetup&) override;
+
+ std::vector caloRechitTags_;
+ std::vector> caloRechitCollectionTokens_;
+ edm::EDGetTokenT> layerClusterToken_;
+};
+
+RecHitToLayerClusterAssociationProducer::RecHitToLayerClusterAssociationProducer(const edm::ParameterSet& pset)
+ : caloRechitTags_(pset.getParameter>("caloRecHits")),
+ caloRechitCollectionTokens_(edm::vector_transform(
+ caloRechitTags_, [this](const edm::InputTag& tag) { return consumes(tag); })),
+ layerClusterToken_(consumes>(pset.getParameter("layerClusters"))) {
+ for (auto& tag : caloRechitTags_) {
+ const std::string& label = !tag.instance().empty() ? tag.instance() : tag.label();
+ produces>>(label + "ToBestLayerCluster");
+ produces(label + "ToLayerCluster");
+ }
+}
+
+RecHitToLayerClusterAssociationProducer::~RecHitToLayerClusterAssociationProducer() {}
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data ------------
+void RecHitToLayerClusterAssociationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+ edm::Handle> lcCollection;
+ iEvent.getByToken(layerClusterToken_, lcCollection);
+ std::unordered_map> hitDetIdToIndices;
+
+ for (size_t j = 0; j < lcCollection->size(); j++) {
+ for (auto& idAndFrac : (*lcCollection)[j].hitsAndFractions()) {
+ hitDetIdToIndices[idAndFrac.first.rawId()].push_back({j, idAndFrac.second});
+ }
+ }
+
+ for (size_t i = 0; i < caloRechitCollectionTokens_.size(); i++) {
+ std::string label = caloRechitTags_.at(i).instance();
+ if (label.empty())
+ label = caloRechitTags_.at(i).label();
+ std::vector rechitIndices;
+
+ edm::Handle caloRechitCollection;
+ iEvent.getByToken(caloRechitCollectionTokens_.at(i), caloRechitCollection);
+
+ auto assocMap = std::make_unique(caloRechitCollection, lcCollection);
+
+ for (size_t h = 0; h < caloRechitCollection->size(); h++) {
+ HGCRecHitRef caloRh(caloRechitCollection, h);
+ size_t id = caloRh->detid().rawId();
+
+ // Need to sort before inserting into AssociationMap
+ auto match = hitDetIdToIndices.find(id);
+ if (match == std::end(hitDetIdToIndices)) {
+ rechitIndices.push_back(-1);
+ continue;
+ }
+ auto& lcIdxAndFrac = match->second;
+ // Sort by energy fraction
+ std::sort(std::begin(lcIdxAndFrac), std::end(lcIdxAndFrac),
+ [](auto& a, auto& b) { return a.second > b.second; });
+
+ for (size_t m = 0; m < lcIdxAndFrac.size(); m++) {
+ float fraction = lcIdxAndFrac[m].second;
+ int lcIdx = lcIdxAndFrac[m].first;
+ // Best match is the layerCluster that carries the hit with the highest energy fraction
+ // (that is, the one responsible for the largest deposit in the detId)
+ if (m == 0)
+ rechitIndices.push_back(lcIdx);
+ edm::Ref> lc(lcCollection, lcIdx);
+ assocMap->insert(caloRh, std::make_pair(lc, fraction));
+ }
+ }
+
+ auto assoc = std::make_unique>>(lcCollection);
+ edm::Association>::Filler filler(*assoc);
+ filler.insert(caloRechitCollection, rechitIndices.begin(), rechitIndices.end());
+ filler.fill();
+ iEvent.put(std::move(assoc), label + "ToBestLayerCluster");
+ iEvent.put(std::move(assocMap), label + "ToLayerCluster");
+ }
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(RecHitToLayerClusterAssociationProducer);
+
diff --git a/CommonTools/RecoAlgos/plugins/RecHitToPFCandAssociationProducer.cc b/CommonTools/RecoAlgos/plugins/RecHitToPFCandAssociationProducer.cc
new file mode 100644
index 0000000000000..6c01bacfac3f8
--- /dev/null
+++ b/CommonTools/RecoAlgos/plugins/RecHitToPFCandAssociationProducer.cc
@@ -0,0 +1,133 @@
+// 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/TrackingAnalysis/interface/TrackingParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
+#include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
+#include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
+#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
+#include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
+
+#include "DataFormats/Common/interface/Association.h"
+
+#include "FWCore/Utilities/interface/transform.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include
+
+//
+// class decleration
+//
+typedef std::pair IdxAndFraction;
+
+class RecHitToPFCandAssociationProducer : public edm::stream::EDProducer<> {
+public:
+ explicit RecHitToPFCandAssociationProducer(const edm::ParameterSet&);
+ ~RecHitToPFCandAssociationProducer() override;
+
+private:
+ void produce(edm::Event&, const edm::EventSetup&) override;
+
+ std::vector caloRechitTags_;
+ std::vector> caloSimhitCollectionTokens_;
+ std::vector>> caloRechitCollectionTokens_;
+ edm::EDGetTokenT pfCollectionToken_;
+};
+
+RecHitToPFCandAssociationProducer::RecHitToPFCandAssociationProducer(const edm::ParameterSet& pset)
+ : caloRechitTags_(pset.getParameter>("caloRecHits")),
+ caloRechitCollectionTokens_(edm::vector_transform(
+ caloRechitTags_, [this](const edm::InputTag& tag) { return consumes>(tag); })),
+ pfCollectionToken_(consumes(pset.getParameter("pfCands"))) {
+ for (auto& tag : caloRechitTags_) {
+ const std::string& label = !tag.instance().empty() ? tag.instance() : tag.label();
+ produces>(label + "ToPFCand");
+ }
+}
+
+RecHitToPFCandAssociationProducer::~RecHitToPFCandAssociationProducer() {}
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data ------------
+void RecHitToPFCandAssociationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+ edm::Handle pfCollection;
+ iEvent.getByToken(pfCollectionToken_, pfCollection);
+ std::unordered_map hitDetIdToIndex;
+
+ for (size_t j = 0; j < pfCollection->size(); j++) {
+ const auto& pfCand = pfCollection->at(j);
+ const reco::PFCandidate::ElementsInBlocks& elements = pfCand.elementsInBlocks();
+ for (auto& element : elements) {
+ const reco::PFBlockRef blockRef = element.first;
+ if (!blockRef.isNonnull())
+ continue;
+ for (const auto& block : blockRef->elements()) {
+ // This seems to not work for PFTICL
+ //if (block.type() == reco::PFBlockElement::HGCAL) {
+ const reco::PFClusterRef cluster = block.clusterRef();
+ if (cluster.isNonnull()) {
+ const std::vector& rhf = cluster->recHitFractions();
+ for (const auto& hf : rhf) {
+ auto& hit = hf.recHitRef();
+ if (!hit)
+ throw cms::Exception("RecHitToPFCandAssociationProducer") << "Invalid RecHit ref";
+ size_t detId = hit->detId();
+ auto entry = hitDetIdToIndex.find(detId);
+ if (entry == hitDetIdToIndex.end() || entry->second.second < hf.fraction())
+ hitDetIdToIndex[detId] = {j, hf.fraction()};
+ }
+ }
+ }
+ }
+ }
+
+ for (size_t i = 0; i < caloRechitCollectionTokens_.size(); i++) {
+ std::string label = caloRechitTags_.at(i).instance();
+ if (label.empty())
+ label = caloRechitTags_.at(i).label();
+ std::vector rechitIndices;
+
+ edm::Handle> caloRechitCollection;
+ iEvent.getByToken(caloRechitCollectionTokens_.at(i), caloRechitCollection);
+
+ for (size_t h = 0; h < caloRechitCollection->size(); h++) {
+ const CaloRecHit& caloRh = caloRechitCollection->at(h);
+ size_t id = caloRh.detid().rawId();
+ int match = hitDetIdToIndex.find(id) == hitDetIdToIndex.end() ? -1 : hitDetIdToIndex.at(id).first;
+ rechitIndices.push_back(match);
+ }
+
+ auto assoc = std::make_unique>(pfCollection);
+ edm::Association::Filler filler(*assoc);
+ filler.insert(caloRechitCollection, rechitIndices.begin(), rechitIndices.end());
+ filler.fill();
+ iEvent.put(std::move(assoc), label + "ToPFCand");
+ }
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(RecHitToPFCandAssociationProducer);
diff --git a/CommonTools/RecoAlgos/plugins/SimClusterRecHitAssociationProducer.cc b/CommonTools/RecoAlgos/plugins/SimClusterRecHitAssociationProducer.cc
new file mode 100644
index 0000000000000..8af6951d69787
--- /dev/null
+++ b/CommonTools/RecoAlgos/plugins/SimClusterRecHitAssociationProducer.cc
@@ -0,0 +1,137 @@
+// 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/Common/interface/Association.h"
+#include "DataFormats/Common/interface/AssociationMap.h"
+#include "DataFormats/Common/interface/OneToManyWithQualityGeneric.h"
+#include "DataFormats/Common/interface/RefToBase.h"
+
+
+#include "FWCore/Utilities/interface/transform.h"
+#include "FWCore/Utilities/interface/EDGetToken.h"
+#include
+
+//
+// class decleration
+//
+typedef std::pair IdxAndFraction;
+typedef edm::AssociationMap> RecHitToSimCluster;
+
+class SimClusterRecHitAssociationProducer : public edm::stream::EDProducer<> {
+public:
+ explicit SimClusterRecHitAssociationProducer(const edm::ParameterSet&);
+ ~SimClusterRecHitAssociationProducer() override;
+
+private:
+ void produce(edm::Event&, const edm::EventSetup&) override;
+
+ std::vector caloRechitTags_;
+ std::vector> caloRechitCollectionTokens_;
+ edm::EDGetTokenT scCollectionToken_;
+};
+
+SimClusterRecHitAssociationProducer::SimClusterRecHitAssociationProducer(const edm::ParameterSet& pset)
+ : caloRechitTags_(pset.getParameter>("caloRecHits")),
+ caloRechitCollectionTokens_(edm::vector_transform(
+ caloRechitTags_, [this](const edm::InputTag& tag) { return consumes(tag); })),
+ scCollectionToken_(consumes(pset.getParameter("simClusters"))) {
+ for (auto& tag : caloRechitTags_) {
+ const std::string& label = !tag.instance().empty() ? tag.instance() : tag.label();
+ produces>(label + "ToBestSimClus");
+ produces(label + "ToSimClus");
+ }
+ produces>();
+}
+
+SimClusterRecHitAssociationProducer::~SimClusterRecHitAssociationProducer() {}
+
+// ------------ method called to produce the data ------------
+void SimClusterRecHitAssociationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+ auto simClusterToRecEnergy = std::make_unique>();
+ edm::Handle scCollection;
+ iEvent.getByToken(scCollectionToken_, scCollection);
+ std::unordered_map> hitDetIdToIndices;
+
+ for (size_t s = 0; s < scCollection->size(); s++) {
+ const auto& sc = scCollection->at(s);
+ // Need to initialize because not all SimClusters lead to rechits
+ (*simClusterToRecEnergy)[s] = 0.;
+ for (auto& hf : sc.hits_and_fractions()) {
+ // Can have two unique hits with the same detId
+ hitDetIdToIndices[hf.first].push_back({s, hf.second});
+ }
+ }
+
+ std::unordered_map hitDetIdToTotalRecEnergy;
+
+ for (size_t i = 0; i < caloRechitCollectionTokens_.size(); i++) {
+ std::string label = caloRechitTags_.at(i).instance();
+ if (label.empty())
+ label = caloRechitTags_.at(i).label();
+ std::vector rechitIndices;
+
+ edm::Handle caloRechitCollection;
+ iEvent.getByToken(caloRechitCollectionTokens_.at(i), caloRechitCollection);
+
+ auto assocMap = std::make_unique(caloRechitCollection, scCollection);
+
+ for (size_t h = 0; h < caloRechitCollection->size(); h++) {
+ HGCRecHitRef caloRh(caloRechitCollection, h);
+ size_t id = caloRh->detid().rawId();
+ float energy = caloRh->energy();
+ hitDetIdToTotalRecEnergy[id] += energy;
+
+ // Need to sort before inserting into AssociationMap
+ auto match = hitDetIdToIndices.find(id);
+ if (match == std::end(hitDetIdToIndices)) {
+ rechitIndices.push_back(-1);
+ continue;
+ }
+ auto& scIdxAndFrac = match->second;
+ // Sort by energy fraction
+ std::sort(std::begin(scIdxAndFrac), std::end(scIdxAndFrac),
+ [](auto& a, auto& b) { return a.second > b.second; });
+
+ for (size_t m = 0; m < scIdxAndFrac.size(); m++) {
+ float fraction = scIdxAndFrac[m].second;
+ int scIdx = scIdxAndFrac[m].first;
+ (*simClusterToRecEnergy)[scIdx] += energy * fraction;
+ // Best match is the simCluster that carries the hit with the highest energy fraction
+ // (that is, responsible for the largest deposit in the detId)
+ if (m == 0)
+ rechitIndices.push_back(scIdx);
+ SimClusterRef simclus(scCollection, scIdx);
+ assocMap->insert(caloRh, std::make_pair(simclus, fraction));
+ }
+ }
+
+ auto assoc = std::make_unique>(scCollection);
+ edm::Association::Filler filler(*assoc);
+ filler.insert(caloRechitCollection, rechitIndices.begin(), rechitIndices.end());
+ filler.fill();
+ iEvent.put(std::move(assoc), label + "ToBestSimClus");
+ iEvent.put(std::move(assocMap), label + "ToSimClus");
+ }
+ iEvent.put(std::move(simClusterToRecEnergy));
+}
+
+// define this as a plug-in
+DEFINE_FWK_MODULE(SimClusterRecHitAssociationProducer);
diff --git a/DPGAnalysis/CommonNanoAOD/BuildFile.xml b/DPGAnalysis/CommonNanoAOD/BuildFile.xml
new file mode 100644
index 0000000000000..6cc75a882ba98
--- /dev/null
+++ b/DPGAnalysis/CommonNanoAOD/BuildFile.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/DPGAnalysis/CommonNanoAOD/plugins/BuildFile.xml b/DPGAnalysis/CommonNanoAOD/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..dc57b6da53c8e
--- /dev/null
+++ b/DPGAnalysis/CommonNanoAOD/plugins/BuildFile.xml
@@ -0,0 +1,19 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/DPGAnalysis/CommonNanoAOD/plugins/HitPositionTableProducer.cc b/DPGAnalysis/CommonNanoAOD/plugins/HitPositionTableProducer.cc
new file mode 100644
index 0000000000000..f86ee30ea25d5
--- /dev/null
+++ b/DPGAnalysis/CommonNanoAOD/plugins/HitPositionTableProducer.cc
@@ -0,0 +1,119 @@
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "DataFormats/NanoAOD/interface/FlatTable.h"
+#include "DataFormats/Common/interface/View.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
+
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCScintillatorDetId.h"
+#include "DataFormats/ForwardDetId/interface/HGCSiliconDetId.h"
+
+#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
+#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
+#include "Geometry/CSCGeometry/interface/CSCGeometry.h"
+#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+#include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+
+#include
+#include
+
+template
+class HitPositionTableProducer : public edm::stream::EDProducer<> {
+public:
+ HitPositionTableProducer(edm::ParameterSet const& params)
+ : name_(params.getParameter("name")),
+ doc_(params.getParameter("doc")),
+ src_(consumes(params.getParameter("src"))),
+ cut_(params.getParameter("cut"), true) {
+ produces();
+ }
+
+ ~HitPositionTableProducer() override {}
+
+ void beginRun(const edm::Run&, const edm::EventSetup& iSetup) override {
+ // TODO: check that the geometry exists
+ iSetup.get().get(caloGeom_);
+ rhtools_.setGeometry(*caloGeom_);
+ iSetup.get().get("idealForDigi", trackGeom_);
+ // Believe this is ideal, but we're not so precise here...
+ iSetup.get().get(globalGeom_);
+ }
+
+ GlobalPoint positionFromHit(const PCaloHit& hit) {
+ DetId id = hit.id();
+ return positionFromDetId(id);
+ }
+
+ GlobalPoint positionFromHit(const CaloRecHit& hit) { return positionFromDetId(hit.detid()); }
+
+ // Should really only be used for HGCAL
+ GlobalPoint positionFromDetId(DetId id) {
+ DetId::Detector det = id.det();
+ if (det == DetId::Hcal || det == DetId::HGCalEE || det == DetId::HGCalHSi || det == DetId::HGCalHSc) {
+ return rhtools_.getPosition(id);
+ } else {
+ throw cms::Exception("HitPositionTableProducer") << "Unsupported DetId type";
+ }
+ }
+
+ GlobalPoint positionFromHit(const PSimHit& hit) {
+ auto detId = DetId(hit.detUnitId());
+ auto surface = detId.det() == DetId::Muon ? globalGeom_->idToDet(hit.detUnitId())->surface()
+ : trackGeom_->idToDet(hit.detUnitId())->surface();
+ GlobalPoint position = surface.toGlobal(hit.localPosition());
+ return position;
+ }
+
+ void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override {
+ edm::Handle objs;
+ iEvent.getByToken(src_, objs);
+
+ std::vector xvals;
+ std::vector yvals;
+ std::vector zvals;
+ for (const auto& obj : *objs) {
+ if (cut_(obj)) {
+ auto position = positionFromHit(obj);
+ xvals.emplace_back(position.x());
+ yvals.emplace_back(position.y());
+ zvals.emplace_back(position.z());
+ }
+ }
+
+ auto tab = std::make_unique(xvals.size(), name_, false, true);
+ tab->addColumn("x", xvals, "x position");
+ tab->addColumn("y", yvals, "y position");
+ tab->addColumn("z", zvals, "z position");
+
+ iEvent.put(std::move(tab));
+ }
+
+protected:
+ const std::string name_, doc_;
+ const edm::EDGetTokenT src_;
+ const StringCutObjectSelector cut_;
+ edm::ESHandle caloGeom_;
+ edm::ESHandle trackGeom_;
+ edm::ESHandle globalGeom_;
+ hgcal::RecHitTools rhtools_;
+};
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+typedef HitPositionTableProducer> PCaloHitPositionTableProducer;
+typedef HitPositionTableProducer> PSimHitPositionTableProducer;
+typedef HitPositionTableProducer HGCRecHitPositionTableProducer;
+DEFINE_FWK_MODULE(HGCRecHitPositionTableProducer);
+DEFINE_FWK_MODULE(PCaloHitPositionTableProducer);
+DEFINE_FWK_MODULE(PSimHitPositionTableProducer);
diff --git a/DPGAnalysis/HGCalNanoAOD/BuildFile.xml b/DPGAnalysis/HGCalNanoAOD/BuildFile.xml
new file mode 100644
index 0000000000000..6cc75a882ba98
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/BuildFile.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml b/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..f5ae82be55fc7
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/plugins/BuildFile.xml
@@ -0,0 +1,18 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc
new file mode 100644
index 0000000000000..636dc103d127a
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromAssociationProducer.cc
@@ -0,0 +1,39 @@
+#include "PhysicsTools/NanoAOD/interface/ObjectIndexFromAssociationProducer.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
+#include "SimDataFormats/PFAnalysis/interface/PFTruthParticle.h"
+#include "SimDataFormats/PFAnalysis/interface/PFTruthParticleFwd.h"
+#include "SimDataFormats/Track/interface/SimTrack.h"
+#include "SimDataFormats/Track/interface/SimTrackContainer.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+typedef ObjectIndexFromAssociationTableProducer
+ SimTrackToSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer
+ CaloHitToSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer, reco::PFCandidateCollection>
+ CaloRecHitToPFCandIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer, SimClusterCollection>
+ CaloRecHitToBestSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer
+ SimClusterToCaloParticleIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer
+ SimClusterToSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer, PFTruthParticleCollection>
+ CaloRecHitToPFTruthParticleIndexTableProducer;
+
+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(CaloRecHitToPFTruthParticleIndexTableProducer);
diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromOneToManyQualAssociationProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromOneToManyQualAssociationProducer.cc
new file mode 100644
index 0000000000000..2a210527206c9
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/plugins/ObjectIndexFromOneToManyQualAssociationProducer.cc
@@ -0,0 +1,21 @@
+#include "PhysicsTools/NanoAOD/interface/ObjectIndexFromOneToManyQualAssociationProducer.h"
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+typedef ObjectIndexFromOneToManyQualAssociationTableProducer, SimClusterCollection>
+ LayerClusterToSimClusterIndexTableProducer;
+typedef ObjectIndexFromOneToManyQualAssociationTableProducer
+ CaloRecHitToSimClusterIndexTableProducer;
+typedef ObjectIndexFromOneToManyQualAssociationTableProducer>
+ HGCRecHitToLayerClusterIndexTableProducer;
+typedef ObjectIndexFromOneToManyQualAssociationTableProducer
+ SimClusterToSimClustersIndexTableProducer;
+
+DEFINE_FWK_MODULE(LayerClusterToSimClusterIndexTableProducer);
+DEFINE_FWK_MODULE(HGCRecHitToLayerClusterIndexTableProducer);
+DEFINE_FWK_MODULE(CaloRecHitToSimClusterIndexTableProducer);
+DEFINE_FWK_MODULE(SimClusterToSimClustersIndexTableProducer);
diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/SimClusterRecEnergyTableProducer.cc b/DPGAnalysis/HGCalNanoAOD/plugins/SimClusterRecEnergyTableProducer.cc
new file mode 100644
index 0000000000000..295f71418841d
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/plugins/SimClusterRecEnergyTableProducer.cc
@@ -0,0 +1,7 @@
+#include "PhysicsTools/NanoAOD/interface/ObjectPropertyFromIndexMapTableProducer.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+typedef ObjectPropertyFromIndexMapTableProducer SimClusterRecEnergyTableProducer;
+DEFINE_FWK_MODULE(SimClusterRecEnergyTableProducer);
diff --git a/DPGAnalysis/HGCalNanoAOD/plugins/SimpleHGCFlatTableProducerPlugins.cc b/DPGAnalysis/HGCalNanoAOD/plugins/SimpleHGCFlatTableProducerPlugins.cc
new file mode 100644
index 0000000000000..206dc856b374a
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/plugins/SimpleHGCFlatTableProducerPlugins.cc
@@ -0,0 +1,22 @@
+#include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h"
+
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+typedef SimpleFlatTableProducer SimpleSimClusterFlatTableProducer;
+
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+typedef SimpleFlatTableProducer SimplePCaloHitFlatTableProducer;
+#include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
+typedef SimpleFlatTableProducer SimpleCaloRecHitFlatTableProducer;
+
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+typedef SimpleFlatTableProducer SimpleCaloParticleFlatTableProducer;
+
+#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
+typedef SimpleFlatTableProducer SimpleCaloClusterFlatTableProducer;
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(SimplePCaloHitFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleCaloRecHitFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleSimClusterFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleCaloParticleFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleCaloClusterFlatTableProducer);
diff --git a/DPGAnalysis/HGCalNanoAOD/python/caloParticles_cff.py b/DPGAnalysis/HGCalNanoAOD/python/caloParticles_cff.py
new file mode 100644
index 0000000000000..91617d618bded
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/caloParticles_cff.py
@@ -0,0 +1,20 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+
+caloParticleTable = cms.EDProducer("SimpleCaloParticleFlatTableProducer",
+ src = cms.InputTag("mix:MergedCaloTruth"),
+ cut = cms.string(""),
+ name = cms.string("CaloPart"),
+ doc = cms.string("CaloPart"),
+ 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,
+ simEnergy = Var('simEnergy', 'int', precision=-1, doc='Number of associated gen particles'),
+ nGenPart = Var('genParticles().size()', 'int', precision=-1, doc='Number of associated gen particles'),
+ GenPartIdx = Var('? genParticles.size() ? genParticles().at(0).key() : -1', 'int', precision=-1, doc='Number of associated gen particles'),
+ nSimHit = Var('numberOfSimHits', 'int', precision=-1, doc='Number of simhits'),
+ trackId = Var('g4Tracks().at(0).trackId', 'int', precision=-1, doc='Geant4 track ID of first track'),
+ nSimTrack = Var('g4Tracks().size', 'int', precision=-1, doc='Number of associated simtracks'),
+ )
+)
+
diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py
new file mode 100644
index 0000000000000..8fae4b33520cf
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHitSimAssociations_cff.py
@@ -0,0 +1,103 @@
+import FWCore.ParameterSet.Config as cms
+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"),
+)
+
+hgcRecHitsToMergedSimClusters = cms.EDProducer("SimClusterRecHitAssociationProducer",
+ caloRecHits = cms.VInputTag("hgcRecHits"),
+ simClusters = cms.InputTag("hgcSimTruth"),
+)
+
+hgcRecHitsToMergedDRSimClusters = cms.EDProducer("SimClusterRecHitAssociationProducer",
+ caloRecHits = cms.VInputTag("hgcRecHits"),
+ simClusters = cms.InputTag("hgcSimTruthDR"),
+)
+
+hgcRecHitsToSimClusterTable = cms.EDProducer("CaloRecHitToSimClusterIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("SimCluster"),
+ objMap = cms.InputTag("hgcRecHitsToSimClusters:hgcRecHitsToSimClus"),
+ docString = cms.string("All SimCluster responsible for sim energy in RecHit DetId (ordered by fraction of energy)")
+)
+
+hgcRecHitsToBestSimClusterTable = cms.EDProducer("CaloRecHitToBestSimClusterIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("BestSimCluster"),
+ objMap = cms.InputTag("hgcRecHitsToSimClusters:hgcRecHitsToBestSimClus"),
+ docString = cms.string("SimCluster responsible for most sim energy in RecHit DetId")
+)
+
+hgcRecHitsToMergedSimClusterTable = cms.EDProducer("CaloRecHitToSimClusterIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("MergedSimCluster"),
+ objMap = cms.InputTag("hgcRecHitsToMergedSimClusters:hgcRecHitsToSimClus"),
+ docString = cms.string("MergedSimClusters ordered by most sim energy in RecHit DetId")
+)
+
+hgcRecHitsToBestMergedSimClusterTable = cms.EDProducer("CaloRecHitToBestSimClusterIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("BestMergedSimCluster"),
+ objMap = cms.InputTag("hgcRecHitsToMergedSimClusters:hgcRecHitsToBestSimClus"),
+ docString = cms.string("MergedSimCluster responsible for most sim energy in RecHit DetId")
+)
+
+hgcRecHitsToMergedDRSimClusterTable = cms.EDProducer("CaloRecHitToSimClusterIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("MergedByDRSimCluster"),
+ objMap = cms.InputTag("hgcRecHitsToMergedDRSimClusters:hgcRecHitsToSimClus"),
+ docString = cms.string("MergedSimCluster responsible for most sim energy in RecHit DetId")
+)
+
+hgcRecHitsToBestMergedDRSimClusterTable = cms.EDProducer("CaloRecHitToBestSimClusterIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("BestMergedByDRSimCluster"),
+ objMap = cms.InputTag("hgcRecHitsToMergedDRSimClusters:hgcRecHitsToBestSimClus"),
+ docString = cms.string("MergedSimCluster (DeltaR merger) responsible for most sim energy in RecHit DetId")
+)
+
+simClusterRecEnergyTable = cms.EDProducer("SimClusterRecEnergyTableProducer",
+ src = cms.InputTag("mix:MergedCaloTruth"),
+ cut = cms.string(""),
+ objName = cms.string("SimCluster"),
+ branchName = cms.string("recEnergy"),
+ valueMap = cms.InputTag("hgcRecHitsToSimClusters"),
+ docString = cms.string("SimCluster deposited reconstructed energy associated to SimCluster")
+)
+
+mergedSimClusterRecEnergyTable = cms.EDProducer("SimClusterRecEnergyTableProducer",
+ src = cms.InputTag("hgcSimTruth"),
+ cut = cms.string(""),
+ objName = cms.string("MergedSimCluster"),
+ branchName = cms.string("recEnergy"),
+ valueMap = cms.InputTag("hgcRecHitsToMergedSimClusters"),
+ docString = cms.string("SimCluster deposited reconstructed energy associated to SimCluster")
+)
+
+hgcRecHitSimAssociationSequence = cms.Sequence(hgcRecHitsToSimClusters
+ +hgcRecHitsToMergedSimClusters
+ +hgcRecHitsToMergedDRSimClusters
+ +simClusterRecEnergyTable
+ +mergedSimClusterRecEnergyTable
+ +hgcRecHitsToSimClusterTable
+ +hgcRecHitsToBestSimClusterTable
+ +hgcRecHitsToMergedSimClusterTable
+ +hgcRecHitsToBestMergedSimClusterTable
+ +hgcRecHitsToMergedDRSimClusterTable
+ +hgcRecHitsToBestMergedDRSimClusterTable
+)
diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py
new file mode 100644
index 0000000000000..f378a107d2ecd
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/hgcRecHits_cff.py
@@ -0,0 +1,82 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import Var,P3Vars
+
+hgcRecHits = cms.EDProducer("HGCRecHitCollectionMerger",
+ src = cms.VInputTag("HGCalRecHit:HGCEERecHits",
+ "HGCalRecHit:HGCHEFRecHits", "HGCalRecHit:HGCHEBRecHits",
+ )
+)
+
+hgcRecHitsTable = cms.EDProducer("SimpleCaloRecHitFlatTableProducer",
+ src = cms.InputTag("hgcRecHits"),
+ cut = cms.string(""),
+ name = cms.string("RecHitHGC"),
+ doc = cms.string("HGCAL RecHits"),
+ 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(
+ detId = Var('detid().rawId()', 'int', precision=-1, doc='detId'),
+ energy = Var('energy', 'float', precision=14, doc='energy'),
+ time = Var('time', 'float', precision=14, doc='hit time'),
+ )
+)
+
+hgcRecHitsToPFCands = cms.EDProducer("RecHitToPFCandAssociationProducer",
+ caloRecHits = cms.VInputTag("hgcRecHits"),
+ pfCands = cms.InputTag("particleFlow"),
+)
+
+hgcRecHitsToPFCandTable = cms.EDProducer("CaloRecHitToPFCandIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("PFCand"),
+ objMap = cms.InputTag("hgcRecHitsToPFCands:hgcRecHitsToPFCand"),
+ docString = cms.string("PFCand with most associated energy in RecHit DetId")
+)
+
+hgcRecHitsToPFTICLCands = cms.EDProducer("RecHitToPFCandAssociationProducer",
+ caloRecHits = cms.VInputTag("hgcRecHits"),
+ pfCands = cms.InputTag("pfTICL"),
+)
+
+hgcRecHitsToPFTICLCandTable = cms.EDProducer("CaloRecHitToPFCandIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("PFTICLCand"),
+ objMap = cms.InputTag("hgcRecHitsToPFTICLCands:hgcRecHitsToPFCand"),
+ 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,
+ name = hgcRecHitsTable.name,
+ doc = hgcRecHitsTable.doc,
+)
+
+hgcRecHitsSequence = cms.Sequence(hgcRecHits
+ +hgcRecHitsTable
+ +hgcRecHitsToPFCands
+ +hgcRecHitsToPFCandTable
+ +hgcRecHitsToPFTICLCands
+ +hgcRecHitsToPFTICLCandTable
+ +hgcRecHitsToLayerClusters
+ +hgcRecHitsToLayerClusterTable
+ +hgcRecHitsPositionTable
+)
diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcSimHits_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcSimHits_cff.py
new file mode 100644
index 0000000000000..108c873b32d1e
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/hgcSimHits_cff.py
@@ -0,0 +1,68 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import Var,P3Vars
+
+hgcEESimHitsTable = cms.EDProducer("SimplePCaloHitFlatTableProducer",
+ src = cms.InputTag("g4SimHits:HGCHitsEE"),
+ cut = cms.string(""),
+ name = cms.string("SimHitHGCEE"),
+ doc = cms.string("Geant4 SimHits in HGCAL Electromagnetic endcap"),
+ 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(
+ detId = Var('id', 'int', precision=-1, doc='detId'),
+ energy = Var('energy', 'float', precision=14, doc='energy'),
+ trackId = Var('geantTrackId', 'int', precision=-1, doc='Geant4 track ID'),
+ )
+)
+
+hgcEEHitsToSimClusterTable = cms.EDProducer("CaloHitToSimClusterIndexTableProducer",
+ cut = hgcEESimHitsTable.cut,
+ src = hgcEESimHitsTable.src,
+ objName = hgcEESimHitsTable.name,
+ branchName = cms.string("SimCluster"),
+ objMap = cms.InputTag("mix:simHitHGCEEToSimCluster"),
+ docString = cms.string("SimCluster containing SimHit")
+)
+
+hgcHEfrontSimHitsTable = hgcEESimHitsTable.clone()
+hgcHEfrontSimHitsTable.src = "g4SimHits:HGCHitsHEfront"
+hgcHEfrontSimHitsTable.name = "SimHitHGCHEF"
+
+hgcHEfrontHitsToSimClusterTable = hgcEEHitsToSimClusterTable.clone()
+hgcHEfrontHitsToSimClusterTable.src = hgcHEfrontSimHitsTable.src
+hgcHEfrontHitsToSimClusterTable.objName = hgcHEfrontSimHitsTable.name
+hgcHEfrontHitsToSimClusterTable.objMap = "mix:simHitHGCHEfrontToSimCluster"
+
+hgcHEbackSimHitsTable = hgcEESimHitsTable.clone()
+hgcHEbackSimHitsTable.src = "g4SimHits:HGCHitsHEback"
+hgcHEbackSimHitsTable.name = "SimHitHGCHEB"
+
+hgcHEbackHitsToSimClusterTable = hgcEEHitsToSimClusterTable.clone()
+hgcHEbackHitsToSimClusterTable.src = hgcHEbackSimHitsTable.src
+hgcHEbackHitsToSimClusterTable.objName = hgcHEbackSimHitsTable.name
+hgcHEbackHitsToSimClusterTable.objMap = "mix:simHitHGCHEbackToSimCluster"
+
+hgcEESimHitsPositionTable = cms.EDProducer("PCaloHitPositionTableProducer",
+ src = hgcEESimHitsTable.src,
+ cut = hgcEESimHitsTable.cut,
+ name = hgcEESimHitsTable.name,
+ doc = hgcEESimHitsTable.doc,
+)
+
+hgcHEfrontSimHitsPositionTable = hgcEESimHitsPositionTable.clone()
+hgcHEfrontSimHitsPositionTable.name = hgcHEfrontSimHitsTable.name
+hgcHEfrontSimHitsPositionTable.src = hgcHEfrontSimHitsTable.src
+
+hgcHEbackSimHitsPositionTable = hgcEESimHitsPositionTable.clone()
+hgcHEbackSimHitsPositionTable.name = hgcHEbackSimHitsTable.name
+hgcHEbackSimHitsPositionTable.src = hgcHEbackSimHitsTable.src
+
+hgcSimHitsSequence = cms.Sequence(hgcEESimHitsTable+hgcHEbackSimHitsTable+hgcHEfrontSimHitsTable
+ +hgcEESimHitsPositionTable
+ +hgcHEfrontSimHitsPositionTable
+ +hgcHEbackSimHitsPositionTable
+ #+hgcEEHitsToSimClusterTable
+ #+hgcHEfrontHitsToSimClusterTable
+ #+hgcHEbackHitsToSimClusterTable
+ +hgcHEfrontSimHitsTable+hgcHEbackSimHitsTable)
+
diff --git a/DPGAnalysis/HGCalNanoAOD/python/hgcSimTracks_cff.py b/DPGAnalysis/HGCalNanoAOD/python/hgcSimTracks_cff.py
new file mode 100644
index 0000000000000..18be2cd601ddf
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/hgcSimTracks_cff.py
@@ -0,0 +1,46 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import Var
+
+simTrackTable = cms.EDProducer("SimpleSimTrackFlatTableProducer",
+ src = cms.InputTag("g4SimHits"),
+ cut = cms.string("abs(momentum().eta) > 1.52 || abs(getMomentumAtBoundary().eta()) > 1.52"),
+ name = cms.string("SimTrack"),
+ doc = cms.string("Geant4 sim tracks in HGCAL Electromagnetic endcap"),
+ 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(
+ crossedBoundary = Var('crossedBoundary', 'bool', doc='track crossed boundary'),
+ pdgId = Var('type', 'int', doc='pdgId (track type)'),
+ charge = Var('charge', 'int', doc='ID'),
+ trackId = Var('trackId', 'int', precision=-1, doc='ID'),
+ trackIdAtBoundary = Var('getIDAtBoundary', 'int', precision=-1, doc='ID at boundary crossing'),
+ pt = Var('momentum().pt()', 'float', precision=14, doc='pt'),
+ eta = Var('momentum().eta()', 'float', precision=14, doc='eta'),
+ phi = Var('momentum().phi()', 'float', precision=14, doc='phi'),
+ lastPos_x = Var('trackerSurfacePosition().x()', 'float', precision=14, doc='x position at HGCAL boundary'),
+ lastPos_y = Var('trackerSurfacePosition().y()', 'float', precision=14, doc='y position at HGCAL boundary'),
+ lastPos_z = Var('trackerSurfacePosition().z()', 'float', precision=14, doc='z position at HGCAL boundary'),
+ boundaryPos_x = Var('getPositionAtBoundary().x()', 'float', precision=14, doc='x position at HGCAL boundary'),
+ boundaryPos_y = Var('getPositionAtBoundary().y()', 'float', precision=14, doc='y position at HGCAL boundary'),
+ boundaryPos_z = Var('getPositionAtBoundary().z()', 'float', precision=14, doc='z position at HGCAL boundary'),
+ boundaryMomentum_pt = Var('getMomentumAtBoundary().pt()', 'float', precision=14, doc='pt at HGCAL boundary'),
+ boundaryMomentum_eta = Var('getMomentumAtBoundary().eta()', 'float', precision=14, doc='eta position at HGCAL boundary'),
+ boundaryMomentum_phi = Var('getMomentumAtBoundary().phi()', 'float', precision=14, doc='phi position at HGCAL boundary'),
+ )
+)
+
+simTrackToSimCluster = cms.EDProducer("SimTrackToSimClusterAssociationProducer",
+ simClusters = cms.InputTag("mix:MergedCaloTruth"),
+ simTracks = cms.InputTag("g4SimHits"),
+)
+
+simTrackToSimClusterTable = cms.EDProducer("SimTrackToSimClusterIndexTableProducer",
+ cut = simTrackTable.cut,
+ src = simTrackTable.src,
+ objName = simTrackTable.name,
+ branchName = cms.string("SimCluster"),
+ objMap = cms.InputTag("simTrackToSimCluster"),
+ docString = cms.string("SimCluster containing track")
+)
+
+simTrackTables = cms.Sequence(simTrackTable+simTrackToSimCluster+simTrackToSimClusterTable)
diff --git a/DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py b/DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py
new file mode 100644
index 0000000000000..555337ce94d5d
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/layerClusters_cff.py
@@ -0,0 +1,33 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import P3Vars,Var
+
+layerClusterTable = cms.EDProducer("SimpleCaloClusterFlatTableProducer",
+ src = cms.InputTag("hgcalLayerClusters"),
+ cut = cms.string(""),
+ name = cms.string("LayerCluster"),
+ doc = cms.string("LayerCluster information"),
+ 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(
+ eta = Var("eta", float,precision=12),
+ phi = Var("phi", float, precision=12),
+ energy = Var("energy", float, precision=14),
+ x = Var('position().x()', 'float', precision=14, doc='x position'),
+ y = Var('position.y()', 'float', precision=14, doc='y position'),
+ z = Var('position.z()', 'float', precision=14, doc='z position'),
+ nHits = Var('size', 'int', precision=14, doc='number of rechits'),
+ seedDetId = Var('seed().rawId()', 'int', precision=-1, doc='detId of seed hit'),
+ )
+)
+
+layerClusterToSimClusterTable = cms.EDProducer("LayerClusterToSimClusterIndexTableProducer",
+ cut = layerClusterTable.cut,
+ src = layerClusterTable.src,
+ objName = layerClusterTable.name,
+ branchName = cms.string("SimCluster"),
+ objMap = cms.InputTag("layerClusterCaloParticleAssociationProducer"),
+ docString = cms.string("Index of SimCluster matched to LayerCluster")
+)
+
+layerClusterTables = cms.Sequence(layerClusterTable+layerClusterToSimClusterTable)
+
diff --git a/DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py b/DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py
new file mode 100644
index 0000000000000..1320a72f05d38
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/nanoHGCML_cff.py
@@ -0,0 +1,58 @@
+from __future__ import print_function
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import *
+from PhysicsTools.NanoAOD.genparticles_cff import genParticleTable
+from PhysicsTools.NanoAOD.genVertex_cff import *
+from DPGAnalysis.HGCalNanoAOD.hgcSimHits_cff import *
+from DPGAnalysis.HGCalNanoAOD.hgcSimTracks_cff import *
+from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import *
+from DPGAnalysis.HGCalNanoAOD.hgcRecHitSimAssociations_cff import *
+from DPGAnalysis.HGCalNanoAOD.simClusters_cff import *
+from DPGAnalysis.HGCalNanoAOD.layerClusters_cff import *
+from DPGAnalysis.HGCalNanoAOD.caloParticles_cff import *
+from DPGAnalysis.TrackNanoAOD.trackSimHits_cff import *
+from DPGAnalysis.TrackNanoAOD.trackingParticles_cff import *
+from DPGAnalysis.TrackNanoAOD.tracks_cff import *
+from DPGAnalysis.PFNanoAOD.pfCands_cff import *
+from DPGAnalysis.PFNanoAOD.pfTruth_cff import *
+
+nanoMetadata = cms.EDProducer("UniqueStringProducer",
+ strings = cms.PSet(
+ tag = cms.string("untagged"),
+ )
+)
+
+genParticleTable.src = "genParticles"
+genParticleTable.variables = cms.PSet(genParticleTable.variables,
+ charge = CandVars.charge)
+
+nanoHGCMLSequence = cms.Sequence(nanoMetadata+genVertexTables+genParticleTable+
+ trackingParticleTable+caloParticleTable+
+ layerClusterTables+
+ simTrackTables+hgcSimHitsSequence+trackerSimHitTables+
+ simClusterTables
+)
+
+nanoHGCMLRecoSequence = cms.Sequence(hgcRecHitsSequence+
+ hgcRecHitSimAssociationSequence+
+ pfCandTable+pfTruth+
+ pfTICLCandTable+trackTables)
+
+def customizeReco(process):
+ process.nanoHGCMLSequence.insert(-1, nanoHGCMLRecoSequence)
+ return process
+
+def customizeNoMergedCaloTruth(process):
+ process.nanoHGCMLSequence.remove(simClusterTable)
+ process.nanoHGCMLSequence.remove(simClusterToCaloPartTable)
+ process.nanoHGCMLSequence.remove(hgcEEHitsToSimClusterTable)
+ process.nanoHGCMLSequence.remove(hgcHEfrontHitsToSimClusterTable)
+ process.nanoHGCMLSequence.remove(hgcHEbackHitsToSimClusterTable)
+ process.nanoHGCMLSequence.remove(simTrackToSimClusterTable)
+
+ process.nanoHGCMLSequence.remove(caloParticleTable)
+ return process
+
+def customizeMergedSimClusters(process):
+ process.nanoHGCMLSequence.insert(-1, mergedSimClusterTables)
+ return process
diff --git a/DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py b/DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py
new file mode 100644
index 0000000000000..617f40b701584
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/python/simClusters_cff.py
@@ -0,0 +1,88 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+
+simClusterTable = cms.EDProducer("SimpleSimClusterFlatTableProducer",
+ src = cms.InputTag("mix:MergedCaloTruth"),
+ cut = cms.string(""),
+ name = cms.string("SimCluster"),
+ doc = cms.string("SimCluster information"),
+ 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,
+ lastPos_x = Var('g4Tracks.at(0).trackerSurfacePosition().x()', 'float', precision=14, doc='track x final position'),
+ lastPos_y = Var('g4Tracks.at(0).trackerSurfacePosition().y()', 'float', precision=14, doc='track y final position'),
+ lastPos_z = Var('g4Tracks.at(0).trackerSurfacePosition().z()', 'float', precision=14, doc='track z final position'),
+ # NOTE: This is the cms-pepr approach, which is needed for merged simclusters
+ impactPoint_x = Var('impactPoint().x()', 'float', precision=14, doc='x position'),
+ impactPoint_y = Var('impactPoint().y()', 'float', precision=14, doc='y position'),
+ impactPoint_z = Var('impactPoint().z()', 'float', precision=14, doc='z position'),
+ impactPoint_t = Var('impactPoint().t()', 'float', precision=14, doc='Impact time'),
+ impactPoint_eta = Var('impactPoint().eta()', 'float', precision=14, doc='eta at boundary'),
+ impactPoint_phi = Var('impactPoint().phi()', 'float', precision=14, doc='phi at boundary'),
+ # For stupid reasons lost on me, the nsimhits_ variable is uninitialized, and hits_ (which are really simhits)
+ # are often referred to as rechits in the SimCluster class
+ nHits = Var('numberOfRecHits', 'int', precision=-1, doc='number of simhits'),
+ sumHitEnergy = Var('energy', 'float', precision=14, doc='total energy of simhits'),
+ boundaryPmag = Var('impactMomentum.P()', 'float', precision=14, doc='magnitude of the boundary 3-momentum vector'),
+ boundaryP4 = Var('impactMomentum.mag()', 'float', precision=14, doc='magnitude of four vector'),
+ boundaryEnergy = Var('impactMomentum.energy()', 'float', precision=14, doc='magnitude of four vector'),
+ boundaryEnergyNoMu = Var('impactMomentumNoMu.energy()', 'float', precision=14, doc='magnitude of four vector'),
+ boundaryPt = Var('impactMomentum.pt()', 'float', precision=14, doc='magnitude of four vector'),
+ trackId = Var('g4Tracks().at(0).trackId()', 'int', precision=-1, doc='Geant track id'),
+ trackIdAtBoundary = Var('g4Tracks().at(0).getIDAtBoundary()', 'int', precision=-1, doc='Track ID at boundary'),
+ hasHGCALHit = Var('hasHGCALHit', 'bool', doc='Has at least 1 simHit in HGCAL'),
+ allHitsHGCAL = Var('allHitsHGCAL', 'bool', doc='all simHits are in HGCAL'),
+ onHGCFrontFace = Var('abs(impactPoint().z()) - 322 < 1', 'bool', doc='SimCluster position is consistent with the front of the HGCAL'),
+ isTrainable = Var('numberOfRecHits > 5 && allHitsHGCAL', 'bool', doc='Should be used for training'),
+ )
+)
+
+simClusterToCaloPart = cms.EDProducer("SimClusterToCaloParticleAssociationProducer",
+ caloParticles = cms.InputTag("mix:MergedCaloTruth"),
+ simClusters = cms.InputTag("mix:MergedCaloTruth"),
+)
+
+simClusterToCaloPartTable = cms.EDProducer("SimClusterToCaloParticleIndexTableProducer",
+ cut = simClusterTable.cut,
+ src = simClusterTable.src,
+ objName = simClusterTable.name,
+ branchName = cms.string("CaloPart"),
+ objMap = cms.InputTag("simClusterToCaloPart"),
+ docString = cms.string("Index of CaloPart containing SimCluster")
+)
+
+hgcSimTruth = cms.EDProducer("simmerger",
+ MergeTheresholdsTransv = cms.vdouble(0.3,0.75,0.85), #(.3, .75, .85),
+ MergeTheresholdsLongitud = cms.vdouble(0.1,0.9) #(.1, .9)
+ )
+hgcSimTruthDR = cms.EDProducer("HGCTruthProducer")
+
+simClusterToMergedSCTable = cms.EDProducer("SimClusterToSimClusterIndexTableProducer",
+ cut = simClusterTable.cut,
+ src = simClusterTable.src,
+ objName = simClusterTable.name,
+ branchName = cms.string("MergedSimCluster"),
+ objMap = cms.InputTag("hgcSimTruth"),
+ docString = cms.string("Index of Merged SimCluster containing SimCluster")
+)
+
+mergedSimClusterTable = simClusterTable.clone()
+mergedSimClusterTable.src = "hgcSimTruth"
+mergedSimClusterTable.name = "MergedSimCluster"
+
+mergedSimClusterDRTable = simClusterTable.clone()
+mergedSimClusterDRTable.src = "hgcSimTruthDR"
+mergedSimClusterDRTable.name = "MergedByDRSimCluster"
+
+mergedToUnmergedSCTable = cms.EDProducer("SimClusterToSimClustersIndexTableProducer",
+ cut = mergedSimClusterTable.cut,
+ src = mergedSimClusterTable.src,
+ objName = mergedSimClusterTable.name,
+ branchName = cms.string("SimCluster"),
+ objMap = cms.InputTag("hgcSimTruth"),
+ docString = cms.string("Index of Merged SimCluster containing SimCluster")
+)
+
+simClusterTables = cms.Sequence(simClusterTable+simClusterToCaloPart+simClusterToCaloPartTable)
+
+mergedSimClusterTables = cms.Sequence(hgcSimTruth+mergedSimClusterTable+hgcSimTruthDR+mergedSimClusterDRTable+mergedToUnmergedSCTable+simClusterToMergedSCTable)
diff --git a/DPGAnalysis/HGCalNanoAOD/test/nanoMLGSD_cfg.py b/DPGAnalysis/HGCalNanoAOD/test/nanoMLGSD_cfg.py
new file mode 100644
index 0000000000000..425a24eab4d35
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/test/nanoMLGSD_cfg.py
@@ -0,0 +1,117 @@
+# Auto generated configuration file
+# using:
+# Revision: 1.19
+# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v
+# with command line options: step1 --filein file:test.root --fileout testNanoML.root --mc --eventcontent NANOAODSIM --datatier NANOAODSIM --conditions auto:mc --step NANO
+import FWCore.ParameterSet.Config as cms
+from FWCore.ParameterSet.VarParsing import VarParsing
+
+
+process = cms.Process('NANO')
+options = VarParsing('python')
+options.setDefault('outputFile', 'testNanoML.root')
+options.parseArguments()
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff')
+process.load('Configuration.Geometry.GeometryExtended2026D49_cff')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('DPGAnalysis.HGCalNanoAOD.nanoHGCML_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+
+process.maxEvents = cms.untracked.PSet(
+ input = cms.untracked.int32(-1),
+ output = cms.optional.untracked.allowed(cms.int32,cms.PSet)
+)
+
+# Input source
+process.source = cms.Source("PoolSource",
+ fileNames = cms.untracked.vstring(options.inputFiles),
+ secondaryFileNames = cms.untracked.vstring()
+)
+
+process.options = cms.untracked.PSet(
+ FailPath = cms.untracked.vstring(),
+ IgnoreCompletely = cms.untracked.vstring(),
+ Rethrow = cms.untracked.vstring(),
+ SkipEvent = cms.untracked.vstring(),
+ allowUnscheduled = cms.obsolete.untracked.bool,
+ canDeleteEarly = cms.untracked.vstring(),
+ emptyRunLumiMode = cms.obsolete.untracked.string,
+ eventSetup = cms.untracked.PSet(
+ forceNumberOfConcurrentIOVs = cms.untracked.PSet(
+ allowAnyLabel_=cms.required.untracked.uint32
+ ),
+ numberOfConcurrentIOVs = cms.untracked.uint32(1)
+ ),
+ fileMode = cms.untracked.string('FULLMERGE'),
+ forceEventSetupCacheClearOnNewRun = cms.untracked.bool(False),
+ makeTriggerResults = cms.obsolete.untracked.bool,
+ numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1),
+ numberOfConcurrentRuns = cms.untracked.uint32(1),
+ numberOfStreams = cms.untracked.uint32(0),
+ numberOfThreads = cms.untracked.uint32(1),
+ printDependencies = cms.untracked.bool(False),
+ sizeOfStackForThreadsInKB = cms.optional.untracked.uint32,
+ throwIfIllegalParameter = cms.untracked.bool(True),
+ wantSummary = cms.untracked.bool(False)
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+ annotation = cms.untracked.string('step1 nevts:1'),
+ name = cms.untracked.string('Applications'),
+ version = cms.untracked.string('$Revision: 1.19 $')
+)
+
+# Output definition
+
+process.NANOAODSIMoutput = cms.OutputModule("NanoAODOutputModule",
+ compressionAlgorithm = cms.untracked.string('LZMA'),
+ compressionLevel = cms.untracked.int32(9),
+ dataset = cms.untracked.PSet(
+ dataTier = cms.untracked.string('NANOAODSIM'),
+ filterName = cms.untracked.string('')
+ ),
+ fileName = cms.untracked.string(options.outputFile),
+ outputCommands = process.NANOAODSIMEventContent.outputCommands
+)
+
+process.NANOAODSIMoutput.outputCommands.remove("keep edmTriggerResults_*_*_*")
+
+# Additional output definition
+
+# Other statements
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:mc', '')
+
+# Path and EndPath definitions
+process.nanoAOD_step = cms.Path(process.nanoHGCMLSequence)
+process.endjob_step = cms.EndPath(process.endOfProcess)
+process.NANOAODSIMoutput_step = cms.EndPath(process.NANOAODSIMoutput)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.nanoAOD_step,process.endjob_step,process.NANOAODSIMoutput_step)
+from PhysicsTools.PatAlgos.tools.helpers import associatePatAlgosToolsTask
+associatePatAlgosToolsTask(process)
+
+# customisation of the process.
+from DPGAnalysis.HGCalNanoAOD.nanoHGCML_cff import customizeReco
+#process = customizeReco(process)
+
+# End of customisation functions
+
+
+# Customisation from command line
+
+# Add early deletion of temporary data products to reduce peak memory need
+from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete
+process = customiseEarlyDelete(process)
+# End adding early deletion
+
diff --git a/DPGAnalysis/HGCalNanoAOD/test/nanoML_cfg.py b/DPGAnalysis/HGCalNanoAOD/test/nanoML_cfg.py
new file mode 100644
index 0000000000000..fc1de4ff88a5d
--- /dev/null
+++ b/DPGAnalysis/HGCalNanoAOD/test/nanoML_cfg.py
@@ -0,0 +1,116 @@
+# Auto generated configuration file
+# using:
+# Revision: 1.19
+# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v
+# with command line options: step1 --filein file:test.root --fileout testNanoML.root --mc --eventcontent NANOAODSIM --datatier NANOAODSIM --conditions auto:mc --step NANO
+import FWCore.ParameterSet.Config as cms
+from FWCore.ParameterSet.VarParsing import VarParsing
+
+
+process = cms.Process('NANO')
+options = VarParsing('python')
+options.setDefault('outputFile', 'testNanoML.root')
+options.parseArguments()
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff')
+process.load('Configuration.Geometry.GeometryExtended2026D49_cff')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('DPGAnalysis.HGCalNanoAOD.nanoHGCML_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+
+process.maxEvents = cms.untracked.PSet(
+ input = cms.untracked.int32(-1),
+ output = cms.optional.untracked.allowed(cms.int32,cms.PSet)
+)
+
+# Input source
+process.source = cms.Source("PoolSource",
+ fileNames = cms.untracked.vstring(options.inputFiles),
+ secondaryFileNames = cms.untracked.vstring()
+)
+
+process.options = cms.untracked.PSet(
+ FailPath = cms.untracked.vstring(),
+ IgnoreCompletely = cms.untracked.vstring(),
+ Rethrow = cms.untracked.vstring(),
+ SkipEvent = cms.untracked.vstring(),
+ allowUnscheduled = cms.obsolete.untracked.bool,
+ canDeleteEarly = cms.untracked.vstring(),
+ emptyRunLumiMode = cms.obsolete.untracked.string,
+ eventSetup = cms.untracked.PSet(
+ forceNumberOfConcurrentIOVs = cms.untracked.PSet(
+ allowAnyLabel_=cms.required.untracked.uint32
+ ),
+ numberOfConcurrentIOVs = cms.untracked.uint32(1)
+ ),
+ fileMode = cms.untracked.string('FULLMERGE'),
+ forceEventSetupCacheClearOnNewRun = cms.untracked.bool(False),
+ makeTriggerResults = cms.obsolete.untracked.bool,
+ numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1),
+ numberOfConcurrentRuns = cms.untracked.uint32(1),
+ numberOfStreams = cms.untracked.uint32(0),
+ numberOfThreads = cms.untracked.uint32(1),
+ printDependencies = cms.untracked.bool(False),
+ sizeOfStackForThreadsInKB = cms.optional.untracked.uint32,
+ throwIfIllegalParameter = cms.untracked.bool(True),
+ wantSummary = cms.untracked.bool(False)
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+ annotation = cms.untracked.string('step1 nevts:1'),
+ name = cms.untracked.string('Applications'),
+ version = cms.untracked.string('$Revision: 1.19 $')
+)
+
+# Output definition
+
+process.NANOAODSIMoutput = cms.OutputModule("NanoAODOutputModule",
+ compressionAlgorithm = cms.untracked.string('LZMA'),
+ compressionLevel = cms.untracked.int32(9),
+ dataset = cms.untracked.PSet(
+ dataTier = cms.untracked.string('NANOAODSIM'),
+ filterName = cms.untracked.string('')
+ ),
+ fileName = cms.untracked.string(options.outputFile),
+ outputCommands = process.NANOAODSIMEventContent.outputCommands
+)
+
+process.NANOAODSIMoutput.outputCommands.remove("keep edmTriggerResults_*_*_*")
+
+# Additional output definition
+
+# Other statements
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:mc', '')
+
+# Path and EndPath definitions
+process.nanoAOD_step = cms.Path(process.nanoHGCMLSequence)
+process.endjob_step = cms.EndPath(process.endOfProcess)
+process.NANOAODSIMoutput_step = cms.EndPath(process.NANOAODSIMoutput)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.nanoAOD_step,process.endjob_step,process.NANOAODSIMoutput_step)
+from PhysicsTools.PatAlgos.tools.helpers import associatePatAlgosToolsTask
+associatePatAlgosToolsTask(process)
+
+# customisation of the process.
+from DPGAnalysis.HGCalNanoAOD.nanoHGCML_cff import customizeReco
+process = customizeReco(process)
+
+# End of customisation functions
+
+
+# Customisation from command line
+
+# Add early deletion of temporary data products to reduce peak memory need
+from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete
+process = customiseEarlyDelete(process)
+# End adding early deletion
diff --git a/DPGAnalysis/PFNanoAOD/plugins/BuildFile.xml b/DPGAnalysis/PFNanoAOD/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..41c6681105b76
--- /dev/null
+++ b/DPGAnalysis/PFNanoAOD/plugins/BuildFile.xml
@@ -0,0 +1,16 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/DPGAnalysis/PFNanoAOD/plugins/PFObjectIndexToAssociationProducer.cc b/DPGAnalysis/PFNanoAOD/plugins/PFObjectIndexToAssociationProducer.cc
new file mode 100644
index 0000000000000..c603853fadd1f
--- /dev/null
+++ b/DPGAnalysis/PFNanoAOD/plugins/PFObjectIndexToAssociationProducer.cc
@@ -0,0 +1,13 @@
+#include "PhysicsTools/NanoAOD/interface/ObjectIndexFromAssociationProducer.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "SimDataFormats/PFAnalysis/interface/PFTruthParticleFwd.h"
+#include "SimDataFormats/PFAnalysis/interface/PFTruthParticle.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+typedef ObjectIndexFromAssociationTableProducer
+ SimClusterToPFTruthParticleIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer
+ TrackingParticleToPFTruthParticleIndexTableProducer;
+DEFINE_FWK_MODULE(TrackingParticleToPFTruthParticleIndexTableProducer);
+DEFINE_FWK_MODULE(SimClusterToPFTruthParticleIndexTableProducer);
diff --git a/DPGAnalysis/PFNanoAOD/plugins/SimplePFParticleFlatTableProducerPlugins.cc b/DPGAnalysis/PFNanoAOD/plugins/SimplePFParticleFlatTableProducerPlugins.cc
new file mode 100644
index 0000000000000..6d906e216ae44
--- /dev/null
+++ b/DPGAnalysis/PFNanoAOD/plugins/SimplePFParticleFlatTableProducerPlugins.cc
@@ -0,0 +1,7 @@
+#include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h"
+
+#include "SimDataFormats/PFAnalysis/interface/PFTruthParticle.h"
+typedef SimpleFlatTableProducer SimplePFTruthParticleFlatTableProducer;
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(SimplePFTruthParticleFlatTableProducer);
diff --git a/DPGAnalysis/PFNanoAOD/python/pfCands_cff.py b/DPGAnalysis/PFNanoAOD/python/pfCands_cff.py
new file mode 100644
index 0000000000000..8b5aafaee8a1c
--- /dev/null
+++ b/DPGAnalysis/PFNanoAOD/python/pfCands_cff.py
@@ -0,0 +1,32 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+
+pfCandTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
+ src = cms.InputTag("particleFlow"),
+ cut = cms.string(""),
+ name = cms.string("PFCand"),
+ doc = cms.string("ParticleFlow 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'),
+ )
+)
+
+pfTICLCandTable = cms.EDProducer("SimpleCandidateFlatTableProducer",
+ src = cms.InputTag("pfTICL"),
+ cut = cms.string(""),
+ name = cms.string("PFTICLCand"),
+ doc = cms.string("ParticleFlow candidates"),
+ singleton = cms.bool(False), # the number of entries is variable
+ extension = cms.bool(False),
+ 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'),
+ )
+)
+
+pfCandTables = cms.Sequence(pfCandTable+pfTICLCandTable)
diff --git a/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py b/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py
new file mode 100644
index 0000000000000..f94c5a72eba58
--- /dev/null
+++ b/DPGAnalysis/PFNanoAOD/python/pfTruth_cff.py
@@ -0,0 +1,62 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+from DPGAnalysis.HGCalNanoAOD.simClusters_cff import simClusterTable
+from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import hgcRecHitsTable
+from DPGAnalysis.TrackNanoAOD.trackingParticles_cff import trackingParticleTable
+
+pfTruthParticles = cms.EDProducer("PFTruthParticleProducer",
+ trackingParticles= cms.InputTag("mix:MergedTrackTruth"),
+ caloParticles= cms.InputTag("mix:MergedCaloTruth"),
+ simClusters = cms.InputTag("mix:MergedCaloTruth"),
+ simVertices= cms.InputTag("g4SimHits"),
+ simTracks= cms.InputTag("g4SimHits"),
+ caloRecHits = cms.InputTag("hgcRecHits"),
+ rechitToSimClusAssoc = cms.InputTag("hgcRecHitsToSimClusters:hgcRecHitsToBestSimClus"),
+)
+
+pfTruthTable = cms.EDProducer("SimplePFTruthParticleFlatTableProducer",
+ src = cms.InputTag("pfTruthParticles"),
+ cut = cms.string(""),
+ name = cms.string("PFTruthPart"),
+ doc = cms.string("PFTruthParticles"),
+ 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,
+ nSimCluster = Var("nSimCluster", 'int', precision=-1, doc='Number of associated SimClusters'),
+ nTrackingPart = Var("nTrackingParticle", 'int', precision=-1, doc='Number of associated SimClusters'),
+ )
+)
+
+simClusterToPFTruthTable = cms.EDProducer("SimClusterToPFTruthParticleIndexTableProducer",
+ cut = cms.string(""),
+ src = cms.InputTag("mix:MergedCaloTruth"),
+ objName = simClusterTable.name,
+ branchName = pfTruthTable.name,
+ objMap = cms.InputTag("pfTruthParticles:simClusToPFTruth"),
+ docString = cms.string("PFTruth particle to which the SimCluster is associated")
+)
+
+trackingPartToPFTruthTable = cms.EDProducer("TrackingParticleToPFTruthParticleIndexTableProducer",
+ cut = cms.string(""),
+ src = cms.InputTag("mix:MergedTrackTruth"),
+ objName = trackingParticleTable.name,
+ branchName = pfTruthTable.name,
+ objMap = cms.InputTag("pfTruthParticles:trackingPartToPFTruth"),
+ docString = cms.string("PFTruth particle to which the TrackingPart is associated")
+)
+
+recHitToPFTruthTable = cms.EDProducer("CaloRecHitToPFTruthParticleIndexTableProducer",
+ cut = hgcRecHitsTable.cut,
+ src = hgcRecHitsTable.src,
+ objName = hgcRecHitsTable.name,
+ branchName = cms.string("BestPFTruthPart"),
+ objMap = cms.InputTag("pfTruthParticles:caloRecHitToPFTruth"),
+ docString = cms.string("PFTruthParticle to which the RecHit is matched (via the sim cluster)")
+)
+
+pfTruth = cms.Sequence(pfTruthParticles+pfTruthTable
+ +simClusterToPFTruthTable
+ +trackingPartToPFTruthTable
+ +recHitToPFTruthTable
+ )
+
diff --git a/DPGAnalysis/TrackNanoAOD/BuildFile.xml b/DPGAnalysis/TrackNanoAOD/BuildFile.xml
new file mode 100644
index 0000000000000..6cc75a882ba98
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/BuildFile.xml
@@ -0,0 +1,6 @@
+
+
+
+
+
+
diff --git a/DPGAnalysis/TrackNanoAOD/plugins/BuildFile.xml b/DPGAnalysis/TrackNanoAOD/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..eba77cbba8b36
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/plugins/BuildFile.xml
@@ -0,0 +1,18 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/DPGAnalysis/TrackNanoAOD/plugins/SimpleTrackFlatTableProducerPlugins.cc b/DPGAnalysis/TrackNanoAOD/plugins/SimpleTrackFlatTableProducerPlugins.cc
new file mode 100644
index 0000000000000..b23e57e8bebaf
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/plugins/SimpleTrackFlatTableProducerPlugins.cc
@@ -0,0 +1,18 @@
+#include "PhysicsTools/NanoAOD/interface/SimpleFlatTableProducer.h"
+
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+typedef SimpleFlatTableProducer SimplePSimHitFlatTableProducer;
+
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
+typedef SimpleFlatTableProducer SimpleTrackingParticleFlatTableProducer;
+
+#include "SimDataFormats/Track/interface/SimTrack.h"
+typedef SimpleFlatTableProducer SimpleSimTrackFlatTableProducer;
+#include "DataFormats/TrackReco/interface/Track.h"
+typedef SimpleFlatTableProducer SimpleTrackFlatTableProducer;
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+DEFINE_FWK_MODULE(SimplePSimHitFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleSimTrackFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleTrackingParticleFlatTableProducer);
+DEFINE_FWK_MODULE(SimpleTrackFlatTableProducer);
diff --git a/DPGAnalysis/TrackNanoAOD/plugins/TrackPositionAtHGCALTableProducer.cc b/DPGAnalysis/TrackNanoAOD/plugins/TrackPositionAtHGCALTableProducer.cc
new file mode 100644
index 0000000000000..603a6a22b5f23
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/plugins/TrackPositionAtHGCALTableProducer.cc
@@ -0,0 +1,65 @@
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+#include "DataFormats/NanoAOD/interface/FlatTable.h"
+#include "DataFormats/Common/interface/View.h"
+#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
+#include "RecoHGCal/GraphReco/interface/HGCalTrackPropagator.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include
+
+class TrackPositionAtHGCALTableProducer : public edm::stream::EDProducer<> {
+public:
+ TrackPositionAtHGCALTableProducer(edm::ParameterSet const& params)
+ : name_(params.getParameter("name")),
+ src_(consumes(params.getParameter("src"))),
+ cut_(params.getParameter("cut"), true) {
+ produces();
+ }
+
+ ~TrackPositionAtHGCALTableProducer() override {}
+
+ void beginRun(const edm::Run&, const edm::EventSetup& iSetup) override {
+ trackprop_.getEventSetup(iSetup);
+ }
+
+ void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override {
+ edm::Handle objs;
+ iEvent.getByToken(src_, objs);
+
+ std::vector xvals;
+ std::vector yvals;
+ std::vector zvals;
+ std::vector phis;
+ std::vector etas;
+ for (const auto& obj : *objs) {
+ if (cut_(obj)) {
+ auto prop = trackprop_.propagateObject(obj);
+ xvals.emplace_back(prop.pos.x());
+ yvals.emplace_back(prop.pos.y());
+ zvals.emplace_back(prop.pos.z());
+ etas.emplace_back(prop.pos.eta());
+ phis.emplace_back(prop.pos.phi());
+ }
+ }
+
+ auto tab = std::make_unique(xvals.size(), name_, false, true);
+ tab->addColumn("HGCFront_x", xvals, "x position at front face of HGCAL");
+ tab->addColumn("HGCFront_y", yvals, "y position at front face of HGCAL");
+ tab->addColumn("HGCFront_z", zvals, "z position at front face of HGCAL");
+ tab->addColumn("HGCFront_eta", etas, "eta position at front face of HGCAL");
+ tab->addColumn("HGCFront_phi", phis, "phi position at front face of HGCAL");
+
+ iEvent.put(std::move(tab));
+ }
+
+protected:
+ const std::string name_;
+ const edm::EDGetTokenT src_;
+ const StringCutObjectSelector cut_;
+ HGCalTrackPropagator trackprop_;
+};
+
+DEFINE_FWK_MODULE(TrackPositionAtHGCALTableProducer);
diff --git a/DPGAnalysis/TrackNanoAOD/python/trackSimHits_cff.py b/DPGAnalysis/TrackNanoAOD/python/trackSimHits_cff.py
new file mode 100644
index 0000000000000..2067b9a50ca79
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/python/trackSimHits_cff.py
@@ -0,0 +1,53 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+
+trackerHitsPixelEndcapLowTofTable = cms.EDProducer("SimplePSimHitFlatTableProducer",
+ src = cms.InputTag("g4SimHits:TrackerHitsPixelEndcapLowTof"),
+ cut = cms.string(""),
+ name = cms.string("SimHitPixelECLowTof"),
+ doc = cms.string("Geant4 SimHits in tracker endcap"),
+ 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(
+ detId = Var('detUnitId', 'int', precision=-1, doc='detId'),
+ energy = Var('energyLoss', 'float', precision=14, doc='energy'),
+ pMag = Var('pabs', 'float', precision=14, doc='magnitude of momentum'),
+ trackId = Var('trackId', 'int', precision=-1, doc='Geant4 track ID'),
+ pdgId = Var('particleType', 'int', doc='PDG ID of associated track'),
+ )
+)
+
+trackerHitsPixelEndcapLowTofPositionTable = cms.EDProducer("PSimHitPositionTableProducer",
+ src = trackerHitsPixelEndcapLowTofTable.src,
+ cut = trackerHitsPixelEndcapLowTofTable.cut,
+ name = trackerHitsPixelEndcapLowTofTable.name,
+ doc = trackerHitsPixelEndcapLowTofTable.doc,
+)
+
+trackerHitsPixelBarrelLowTofTable = trackerHitsPixelEndcapLowTofTable.clone()
+trackerHitsPixelBarrelLowTofTable.src = "g4SimHits:TrackerHitsPixelBarrelLowTof"
+trackerHitsPixelBarrelLowTofTable.name = "SimHitPixelLowTof"
+trackerHitsPixelBarrelLowTofTable.doc = "Geant4 SimHits in pixel barrel"
+
+trackerHitsPixelBarrelLowTofPositionTable = cms.EDProducer("PSimHitPositionTableProducer",
+ src = trackerHitsPixelBarrelLowTofTable.src,
+ cut = trackerHitsPixelBarrelLowTofTable.cut,
+ name = trackerHitsPixelBarrelLowTofTable.name,
+ doc = trackerHitsPixelBarrelLowTofTable.doc,
+)
+
+muonCSCHitsTable = trackerHitsPixelEndcapLowTofTable.clone()
+muonCSCHitsTable.src = "g4SimHits:MuonCSCHits"
+muonCSCHitsTable.name = "SimHitMuonCSC"
+muonCSCHitsTable.doc = "Geant4 SimHits in Muon CSCs"
+
+muonCSCHitsPositionTable = cms.EDProducer("PSimHitPositionTableProducer",
+ src = muonCSCHitsTable.src,
+ cut = muonCSCHitsTable.cut,
+ name = muonCSCHitsTable.name,
+ doc = muonCSCHitsTable.doc,
+)
+
+trackerSimHitTables = cms.Sequence(trackerHitsPixelEndcapLowTofTable+trackerHitsPixelEndcapLowTofPositionTable+
+ trackerHitsPixelBarrelLowTofTable+trackerHitsPixelBarrelLowTofPositionTable+
+ muonCSCHitsTable+muonCSCHitsPositionTable)
diff --git a/DPGAnalysis/TrackNanoAOD/python/trackingParticles_cff.py b/DPGAnalysis/TrackNanoAOD/python/trackingParticles_cff.py
new file mode 100644
index 0000000000000..c2c7e4f4d72e0
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/python/trackingParticles_cff.py
@@ -0,0 +1,28 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import CandVars,Var
+
+trackingParticleTable = cms.EDProducer("SimpleTrackingParticleFlatTableProducer",
+ src = cms.InputTag("mix:MergedTrackTruth"),
+ cut = cms.string(""),
+ name = cms.string("TrackingPart"),
+ doc = cms.string("TrackingPart"),
+ 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,
+ nGenPart = Var('genParticles().size()', 'int', precision=-1, doc='Number of associated gen particles'),
+ GenPartIdx = Var('? genParticles.size() ? genParticles().at(0).key() : -1', 'int', precision=-1, doc='Number of associated gen particles'),
+ trackId = Var('g4Tracks.at(0).trackId', 'int', precision=-1, doc='Geant4 track ID of first track'),
+ nSimTrack = Var('g4Tracks().size', 'int', precision=-1, doc='Number of associated simtracks'),
+ Vtx_x = Var('vx()', 'float', precision=14, doc='parent vertex x pos'),
+ Vtx_y = Var('vy()', 'float', precision=14, doc='parent vertex y pos'),
+ Vtx_z = Var('vz()', 'float', precision=14, doc='parent vertex z pos'),
+ Vtx_t = Var('parentVertex().position().T()', 'float', precision=14, doc='parent vertex time'),
+ nDecayVtx = Var('decayVertices().size()', 'int', precision=-1, doc='number of decay vertices'),
+ DecayVtx_y = Var('? decayVertices().size() > 0 ? decayVertices().at(0).position().x : 10000', 'float', precision=14, doc='x position of first decay vertex'),
+ DecayVtx_x = Var('? decayVertices().size() > 0 ? decayVertices().at(0).position().y : 10000', 'float', precision=14, doc='y position of first decay vertex'),
+ DecayVtx_z = Var('? decayVertices().size() > 0 ? decayVertices().at(0).position().z : 10000', 'float', precision=14, doc='z position of first decay vertex'),
+ DecayVtx_t = Var('? decayVertices().size() > 0 ? decayVertices().at(0).position().t : 10000', 'float', precision=14, doc='time of first decay vertex'),
+ )
+)
+
+
diff --git a/DPGAnalysis/TrackNanoAOD/python/tracks_cff.py b/DPGAnalysis/TrackNanoAOD/python/tracks_cff.py
new file mode 100644
index 0000000000000..819f0afab71c0
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/python/tracks_cff.py
@@ -0,0 +1,54 @@
+import FWCore.ParameterSet.Config as cms
+from PhysicsTools.NanoAOD.common_cff import P3Vars,Var
+
+generalTrackTable = cms.EDProducer("SimpleTrackFlatTableProducer",
+ src = cms.InputTag("generalTracks"),
+ cut = cms.string(""),
+ name = cms.string("Track"),
+ doc = cms.string("reco::Track"),
+ 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(P3Vars,
+ charge = Var("charge", int, doc="electric charge"),
+ normChiSq = Var("normalizedChi2", float, precision=14, doc="Chi^2/ndof"),
+ numberOfValidHits = Var('numberOfValidHits()', 'int', precision=-1, doc='Number of valid hits in track'),
+ numberOfLostHits = Var('numberOfLostHits()', 'int', precision=-1, doc='Number of lost hits in track'),
+ Vtx_x = Var('vx()', 'float', precision=14, doc='parent vertex x pos'),
+ Vtx_y = Var('vy()', 'float', precision=14, doc='parent vertex y pos'),
+ Vtx_z = Var('vz()', 'float', precision=14, doc='parent vertex z pos'),
+ Vtx_t = Var('t0', 'float', precision=14, doc='parent vertex time'),
+ # Be careful here, this isn't really a decay vertex
+ DecayVtx_y = Var('outerPosition().x()', 'float', precision=14, doc='x position of first decay vertex'),
+ DecayVtx_x = Var('outerPosition().y()', 'float', precision=14, doc='y position of first decay vertex'),
+ DecayVtx_z = Var('outerPosition().z()', 'float', precision=14, doc='z position of first decay vertex'),
+ DecayVtx_t = Var('0', 'float', precision=14, doc='DUMMY VALUE! for time of first decay vertex'),
+ )
+)
+
+generalTrackHGCPositionTable = cms.EDProducer("TrackPositionAtHGCALTableProducer",
+ src = generalTrackTable.src,
+ name = generalTrackTable.name,
+ cut = generalTrackTable.cut,
+)
+
+trackConversionsTable = generalTrackTable.clone()
+trackConversionsTable.src = "conversionStepTracks"
+trackConversionsTable.name = "TrackConv"
+
+trackDisplacedTable = cms.EDProducer("SimpleTrackFlatTableProducer",
+ src = cms.InputTag("displacedTracks"),
+ cut = cms.string(""),
+ name = cms.string("TrackDisp"),
+ doc = cms.string("reco::Track"),
+ 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(P3Vars,
+ charge = Var("charge", int, doc="electric charge"),
+ Vtx_x = Var('vx()', 'float', precision=14, doc='parent vertex x pos'),
+ Vtx_y = Var('vy()', 'float', precision=14, doc='parent vertex y pos'),
+ Vtx_z = Var('vz()', 'float', precision=14, doc='parent vertex z pos'),
+ Vtx_t = Var('t0', 'float', precision=14, doc='parent vertex time'),
+ )
+)
+
+trackTables = cms.Sequence(generalTrackTable+generalTrackHGCPositionTable+trackConversionsTable+trackDisplacedTable)
diff --git a/DPGAnalysis/TrackNanoAOD/test/nanoMLGSD_cfg.py b/DPGAnalysis/TrackNanoAOD/test/nanoMLGSD_cfg.py
new file mode 100644
index 0000000000000..a696831220612
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/test/nanoMLGSD_cfg.py
@@ -0,0 +1,118 @@
+# Auto generated configuration file
+# using:
+# Revision: 1.19
+# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v
+# with command line options: step1 --filein file:test.root --fileout testNanoML.root --mc --eventcontent NANOAODSIM --datatier NANOAODSIM --conditions auto:mc --step NANO
+import FWCore.ParameterSet.Config as cms
+from FWCore.ParameterSet.VarParsing import VarParsing
+
+
+process = cms.Process('NANO')
+options = VarParsing('python')
+options.setDefault('outputFile', 'testNanoML.root')
+options.parseArguments()
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff')
+process.load('Configuration.Geometry.GeometryExtended2026D49_cff')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('PhysicsTools.NanoAOD.nanoHGCML_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+
+process.maxEvents = cms.untracked.PSet(
+ input = cms.untracked.int32(-1),
+ output = cms.optional.untracked.allowed(cms.int32,cms.PSet)
+)
+
+# Input source
+process.source = cms.Source("PoolSource",
+ fileNames = cms.untracked.vstring(options.inputFiles),
+ secondaryFileNames = cms.untracked.vstring()
+)
+
+process.options = cms.untracked.PSet(
+ FailPath = cms.untracked.vstring(),
+ IgnoreCompletely = cms.untracked.vstring(),
+ Rethrow = cms.untracked.vstring(),
+ SkipEvent = cms.untracked.vstring(),
+ allowUnscheduled = cms.obsolete.untracked.bool,
+ canDeleteEarly = cms.untracked.vstring(),
+ emptyRunLumiMode = cms.obsolete.untracked.string,
+ eventSetup = cms.untracked.PSet(
+ forceNumberOfConcurrentIOVs = cms.untracked.PSet(
+ allowAnyLabel_=cms.required.untracked.uint32
+ ),
+ numberOfConcurrentIOVs = cms.untracked.uint32(1)
+ ),
+ fileMode = cms.untracked.string('FULLMERGE'),
+ forceEventSetupCacheClearOnNewRun = cms.untracked.bool(False),
+ makeTriggerResults = cms.obsolete.untracked.bool,
+ numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1),
+ numberOfConcurrentRuns = cms.untracked.uint32(1),
+ numberOfStreams = cms.untracked.uint32(0),
+ numberOfThreads = cms.untracked.uint32(1),
+ printDependencies = cms.untracked.bool(False),
+ sizeOfStackForThreadsInKB = cms.optional.untracked.uint32,
+ throwIfIllegalParameter = cms.untracked.bool(True),
+ wantSummary = cms.untracked.bool(False)
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+ annotation = cms.untracked.string('step1 nevts:1'),
+ name = cms.untracked.string('Applications'),
+ version = cms.untracked.string('$Revision: 1.19 $')
+)
+
+# Output definition
+
+process.NANOAODSIMoutput = cms.OutputModule("NanoAODOutputModule",
+ compressionAlgorithm = cms.untracked.string('LZMA'),
+ compressionLevel = cms.untracked.int32(9),
+ dataset = cms.untracked.PSet(
+ dataTier = cms.untracked.string('NANOAODSIM'),
+ filterName = cms.untracked.string('')
+ ),
+ fileName = cms.untracked.string(options.outputFile),
+ outputCommands = process.NANOAODSIMEventContent.outputCommands
+)
+
+process.NANOAODSIMoutput.outputCommands.remove("keep edmTriggerResults_*_*_*")
+
+# Additional output definition
+
+# Other statements
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:mc', '')
+
+# Path and EndPath definitions
+process.nanoAOD_step = cms.Path(process.nanoHGCMLSequence)
+process.endjob_step = cms.EndPath(process.endOfProcess)
+process.NANOAODSIMoutput_step = cms.EndPath(process.NANOAODSIMoutput)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.nanoAOD_step,process.endjob_step,process.NANOAODSIMoutput_step)
+from PhysicsTools.PatAlgos.tools.helpers import associatePatAlgosToolsTask
+associatePatAlgosToolsTask(process)
+
+# customisation of the process.
+
+# End of customisation functions
+from PhysicsTools.NanoAOD.nanoHGCML_cff import customizeNoMergedCaloTruth,customizeMergedSimClusters
+# Uncomment if you didn't schedule SimClusters/CaloParticles
+# process = customizeNoMergedCaloTruth(process)
+# merged simclusters (turn off if you aren't running through PEPR)
+process = customizeMergedSimClusters(process)
+
+# Customisation from command line
+
+# Add early deletion of temporary data products to reduce peak memory need
+from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete
+process = customiseEarlyDelete(process)
+# End adding early deletion
diff --git a/DPGAnalysis/TrackNanoAOD/test/nanoML_cfg.py b/DPGAnalysis/TrackNanoAOD/test/nanoML_cfg.py
new file mode 100644
index 0000000000000..c3e61dbb2b302
--- /dev/null
+++ b/DPGAnalysis/TrackNanoAOD/test/nanoML_cfg.py
@@ -0,0 +1,120 @@
+# Auto generated configuration file
+# using:
+# Revision: 1.19
+# Source: /local/reps/CMSSW/CMSSW/Configuration/Applications/python/ConfigBuilder.py,v
+# with command line options: step1 --filein file:test.root --fileout testNanoML.root --mc --eventcontent NANOAODSIM --datatier NANOAODSIM --conditions auto:mc --step NANO
+import FWCore.ParameterSet.Config as cms
+from FWCore.ParameterSet.VarParsing import VarParsing
+
+
+process = cms.Process('NANO')
+options = VarParsing('python')
+options.setDefault('outputFile', 'testNanoML.root')
+options.parseArguments()
+
+# import of standard configurations
+process.load('Configuration.StandardSequences.Services_cff')
+process.load('SimGeneral.HepPDTESSource.pythiapdt_cfi')
+process.load('FWCore.MessageService.MessageLogger_cfi')
+process.load('Configuration.EventContent.EventContent_cff')
+process.load('SimGeneral.MixingModule.mixNoPU_cfi')
+process.load('Configuration.Geometry.GeometryExtended2026D49Reco_cff')
+process.load('Configuration.Geometry.GeometryExtended2026D49_cff')
+process.load('Configuration.StandardSequences.MagneticField_cff')
+process.load('PhysicsTools.NanoAOD.nanoHGCML_cff')
+process.load('Configuration.StandardSequences.EndOfProcess_cff')
+process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff')
+
+process.maxEvents = cms.untracked.PSet(
+ input = cms.untracked.int32(-1),
+ output = cms.optional.untracked.allowed(cms.int32,cms.PSet)
+)
+
+# Input source
+process.source = cms.Source("PoolSource",
+ fileNames = cms.untracked.vstring(options.inputFiles),
+ secondaryFileNames = cms.untracked.vstring()
+)
+
+process.options = cms.untracked.PSet(
+ FailPath = cms.untracked.vstring(),
+ IgnoreCompletely = cms.untracked.vstring(),
+ Rethrow = cms.untracked.vstring(),
+ SkipEvent = cms.untracked.vstring(),
+ allowUnscheduled = cms.obsolete.untracked.bool,
+ canDeleteEarly = cms.untracked.vstring(),
+ emptyRunLumiMode = cms.obsolete.untracked.string,
+ eventSetup = cms.untracked.PSet(
+ forceNumberOfConcurrentIOVs = cms.untracked.PSet(
+ allowAnyLabel_=cms.required.untracked.uint32
+ ),
+ numberOfConcurrentIOVs = cms.untracked.uint32(1)
+ ),
+ fileMode = cms.untracked.string('FULLMERGE'),
+ forceEventSetupCacheClearOnNewRun = cms.untracked.bool(False),
+ makeTriggerResults = cms.obsolete.untracked.bool,
+ numberOfConcurrentLuminosityBlocks = cms.untracked.uint32(1),
+ numberOfConcurrentRuns = cms.untracked.uint32(1),
+ numberOfStreams = cms.untracked.uint32(0),
+ numberOfThreads = cms.untracked.uint32(1),
+ printDependencies = cms.untracked.bool(False),
+ sizeOfStackForThreadsInKB = cms.optional.untracked.uint32,
+ throwIfIllegalParameter = cms.untracked.bool(True),
+ wantSummary = cms.untracked.bool(False)
+)
+
+# Production Info
+process.configurationMetadata = cms.untracked.PSet(
+ annotation = cms.untracked.string('step1 nevts:1'),
+ name = cms.untracked.string('Applications'),
+ version = cms.untracked.string('$Revision: 1.19 $')
+)
+
+# Output definition
+
+process.NANOAODSIMoutput = cms.OutputModule("NanoAODOutputModule",
+ compressionAlgorithm = cms.untracked.string('LZMA'),
+ compressionLevel = cms.untracked.int32(9),
+ dataset = cms.untracked.PSet(
+ dataTier = cms.untracked.string('NANOAODSIM'),
+ filterName = cms.untracked.string('')
+ ),
+ fileName = cms.untracked.string(options.outputFile),
+ outputCommands = process.NANOAODSIMEventContent.outputCommands
+)
+
+process.NANOAODSIMoutput.outputCommands.remove("keep edmTriggerResults_*_*_*")
+
+# Additional output definition
+
+# Other statements
+from Configuration.AlCa.GlobalTag import GlobalTag
+process.GlobalTag = GlobalTag(process.GlobalTag, 'auto:mc', '')
+
+# Path and EndPath definitions
+process.nanoAOD_step = cms.Path(process.nanoHGCMLSequence)
+process.endjob_step = cms.EndPath(process.endOfProcess)
+process.NANOAODSIMoutput_step = cms.EndPath(process.NANOAODSIMoutput)
+
+# Schedule definition
+process.schedule = cms.Schedule(process.nanoAOD_step,process.endjob_step,process.NANOAODSIMoutput_step)
+from PhysicsTools.PatAlgos.tools.helpers import associatePatAlgosToolsTask
+associatePatAlgosToolsTask(process)
+
+# customisation of the process.
+from PhysicsTools.NanoAOD.nanoHGCML_cff import customizeReco,customizeMergedSimClusters
+process = customizeReco(process)
+# Uncomment if you didn't schedule SimClusters/CaloParticles
+# process = customizeNoMergedCaloTruth(process)
+# merged simclusters (turn off if you aren't running through PEPR)
+process = customizeMergedSimClusters(process)
+
+# End of customisation functions
+
+
+# Customisation from command line
+
+# Add early deletion of temporary data products to reduce peak memory need
+from Configuration.StandardSequences.earlyDeleteSettings_cff import customiseEarlyDelete
+process = customiseEarlyDelete(process)
+# End adding early deletion
diff --git a/HGCSimTruth/HGCSimTruth/.DS_Store b/HGCSimTruth/HGCSimTruth/.DS_Store
new file mode 100644
index 0000000000000..d4c3b3b44335f
Binary files /dev/null and b/HGCSimTruth/HGCSimTruth/.DS_Store differ
diff --git a/HGCSimTruth/HGCSimTruth/._.DS_Store b/HGCSimTruth/HGCSimTruth/._.DS_Store
new file mode 100644
index 0000000000000..a5b28df1cbc6e
Binary files /dev/null and b/HGCSimTruth/HGCSimTruth/._.DS_Store differ
diff --git a/HGCSimTruth/HGCSimTruth/BuildFile.xml b/HGCSimTruth/HGCSimTruth/BuildFile.xml
new file mode 100644
index 0000000000000..4383c46093a58
--- /dev/null
+++ b/HGCSimTruth/HGCSimTruth/BuildFile.xml
@@ -0,0 +1,25 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/HGCSimTruth/HGCSimTruth/interface/SimClusterTools.h b/HGCSimTruth/HGCSimTruth/interface/SimClusterTools.h
new file mode 100644
index 0000000000000..589d534a1d6da
--- /dev/null
+++ b/HGCSimTruth/HGCSimTruth/interface/SimClusterTools.h
@@ -0,0 +1,42 @@
+/*
+ * SimClusterTools.h
+ *
+ * Created on: 16 Dec 2019
+ * Author: jkiesele
+ */
+
+#ifndef PRODUCTION_HGCALSIM_CMSSW_JKIESELE_CMSSW_11_0_0_PRE9_SRC_RECOHGCAL_GRAPHRECO_INTERFACE_SIMCLUSTERTOOLS_H_
+#define PRODUCTION_HGCALSIM_CMSSW_JKIESELE_CMSSW_11_0_0_PRE9_SRC_RECOHGCAL_GRAPHRECO_INTERFACE_SIMCLUSTERTOOLS_H_
+
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/Math/interface/deltaR.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include
+#include
+
+class SimClusterTools {
+public:
+ SimClusterTools() : recHitTools_(0), detid_to_rh_index_(0), allrechits_(0), considerhitdistance_(0.1) {}
+
+ void recalculatePosition(SimCluster& cluster, const double& assigned_time) const;
+
+ //has at least one hit in the hgcal
+ bool isHGCal(const SimCluster& cluster)const;
+
+ void setRechitTools(const hgcal::RecHitTools& rht) { recHitTools_ = &rht; }
+
+ void setRechitVector(const std::vector& rhv) { allrechits_ = &rhv; }
+
+ void setIndexMap(const std::unordered_map& idxmap) { detid_to_rh_index_ = &idxmap; }
+
+private:
+ const hgcal::RecHitTools* recHitTools_;
+ const std::unordered_map* detid_to_rh_index_;
+ const std::vector* allrechits_;
+ float considerhitdistance_;
+};
+
+#endif /* PRODUCTION_HGCALSIM_CMSSW_JKIESELE_CMSSW_11_0_0_PRE9_SRC_RECOHGCAL_GRAPHRECO_INTERFACE_SIMCLUSTERTOOLS_H_ */
diff --git a/HGCSimTruth/HGCSimTruth/plugins/BuildFile.xml b/HGCSimTruth/HGCSimTruth/plugins/BuildFile.xml
new file mode 100644
index 0000000000000..0b4c0b3575559
--- /dev/null
+++ b/HGCSimTruth/HGCSimTruth/plugins/BuildFile.xml
@@ -0,0 +1,25 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/HGCSimTruth/HGCSimTruth/plugins/HGCSimTruth.cc b/HGCSimTruth/HGCSimTruth/plugins/HGCSimTruth.cc
new file mode 100644
index 0000000000000..455ad8064821c
--- /dev/null
+++ b/HGCSimTruth/HGCSimTruth/plugins/HGCSimTruth.cc
@@ -0,0 +1,722 @@
+/*
+ * CMSSW module to create the trees used for the HH2bbtautau analysis.
+ * Disclainer: This prototype is not meant for creating an official PR. Major refactoring ahead.
+ */
+
+#include
+
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/ProducerBase.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
+#include "DataFormats/HcalDetId/interface/HcalDetId.h"
+#include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
+#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/Math/interface/deltaR.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
+#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
+
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
+#include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
+#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
+#include "SimDataFormats/Vertex/interface/SimVertex.h"
+
+#include "SimGeneral/MixingModule/interface/DigiAccumulatorMixMod.h"
+#include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
+#include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
+#include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
+
+#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
+#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
+#include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
+#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
+#include "Geometry/Records/interface/CaloGeometryRecord.h"
+#include "DataFormats/Common/interface/Association.h"
+#include "FWCore/Utilities/interface/transform.h"
+
+#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
+
+#include "HGCSimTruth/HGCSimTruth/interface/SimClusterTools.h"
+#include "DataFormats/Common/interface/AssociationMap.h"
+#include "DataFormats/Common/interface/OneToManyWithQualityGeneric.h"
+
+#include "TMath.h"
+#include "TCanvas.h"
+#include "TH1F.h"
+
+float getShowerRadius(float f, float r1, float f1) {
+ float a = r1 / (sqrt(2) * TMath::ErfInverse(f1 - 0.00001));
+ return a * sqrt(2) * TMath::ErfInverse(f - 0.00001);
+}
+
+float getShowerRadius(float f, float r1, float f1, float r2, float f2) {
+ float d1 = getShowerRadius(f, r1, f1);
+ float d2 = getShowerRadius(f, r2, f2);
+ return (d1 + d2) / 2.f;
+}
+
+float getShowerRadius(float f, float r1, float f1, float r2, float f2, float r3, float f3) {
+ float d1 = getShowerRadius(f, r1, f1);
+ float d2 = getShowerRadius(f, r2, f2);
+ float d3 = getShowerRadius(f, r3, f3);
+ return (d1 + d2 + d3) / 3.f;
+}
+
+int commonAncestorPdgId(const std::vector>& v1, const std::vector>& v2) {
+ // integer pairs represent the vertex id in the graph and the pdg id of the corresponding particle
+ // to get the pdg id of the first common ancestor, walk through v1 and check if any vertex id is
+ // contained in v2, and if so, return the pdg id (should be the same for the matching elements)
+ for (const auto& p1 : v1) {
+ for (const auto& p2 : v2) {
+ if (p1.first == p2.first) {
+ return p1.second;
+ }
+ }
+ }
+ return 0;
+}
+
+struct ChainIndex {
+ ChainIndex(int iSC, ChainIndex* prev, ChainIndex* next) : iSC(iSC), prev(prev), next(next) {}
+
+ ChainIndex(int iSC) : ChainIndex(iSC, nullptr, nullptr) {}
+
+ ChainIndex(const ChainIndex& o) : ChainIndex(o.iSC, o.prev, o.next) {}
+
+ bool operator==(const ChainIndex& o) const { return o.iSC == iSC; }
+
+ int iSC;
+ ChainIndex* prev;
+ ChainIndex* next;
+};
+
+typedef std::pair IdxAndFraction;
+typedef edm::AssociationMap> SimClusterToSimClusters;
+
+class HGCTruthProducer : public edm::stream::EDProducer<> {
+public:
+ static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+ explicit HGCTruthProducer(const edm::ParameterSet&);
+ ~HGCTruthProducer();
+
+ virtual void beginStream(edm::StreamID) override;
+ virtual void endStream() override;
+ virtual void produce(edm::Event&, const edm::EventSetup&) override;
+
+ void determineSimClusterGroups(const SimClusterCollection&,
+ std::vector>&,
+ std::unique_ptr>&,
+ std::unique_ptr>&) const;
+
+ bool checkSimClusterMerging(int,
+ int,
+ const SimClusterCollection&,
+ const std::unique_ptr>&,
+ const std::unique_ptr>&) const;
+
+ void mergeSimClusters(const SimClusterCollection&,
+ const std::vector>&,
+ std::unique_ptr&,
+ const std::vector&,
+ const std::unordered_map&) const;
+
+private:
+ float showerContainment_;
+ float minEta_;
+ float maxEta_;
+ int minHitsRadius_;
+ float maxDeltaEtaGroup_;
+ bool verbose_;
+
+ edm::EDGetTokenT> caloParticleToken_;
+ edm::EDGetTokenT> simClusterToken_;
+ std::vector> recHitTokens_;
+ //std::vector simHitTags_
+ //std::vector> simHitTokens_;
+
+ hgcal::RecHitTools recHitTools_;
+};
+
+void HGCTruthProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+ edm::ParameterSetDescription desc;
+
+ desc.add("minEta", 1.52);
+ desc.add("maxEta", 3.00);
+ desc.add("minHitsRadius", 5);
+ desc.add("maxDeltaEtaGroup", 0.15);
+ desc.addUntracked("verbose", false);
+ desc.add("caloParticleCollection", edm::InputTag("mix", "MergedCaloTruth"));
+ desc.add("simClusterCollection", edm::InputTag("mix", "MergedCaloTruth"));
+ desc.add>("recHitCollections",
+ {
+ edm::InputTag("HGCalRecHit", "HGCEERecHits"),
+ edm::InputTag("HGCalRecHit", "HGCHEFRecHits"),
+ edm::InputTag("HGCalRecHit", "HGCHEBRecHits"),
+ });
+ //desc.add>("simHitCollections",
+ // {
+ // edm::InputTag("g4SimHits", "HGCHitsEE"),
+ // edm::InputTag("g4SimHits", "HGCHitsHEfront"),
+ // edm::InputTag("g4SimHits", "HGCHitsHEback"),
+ // });
+
+
+ descriptions.add("hgcTruthProducer", desc);
+}
+
+HGCTruthProducer::HGCTruthProducer(const edm::ParameterSet& params)
+ : showerContainment_(0.6827),
+ minEta_(params.getParameter("minEta")),
+ maxEta_(params.getParameter("maxEta")),
+ minHitsRadius_(params.getParameter("minHitsRadius")),
+ maxDeltaEtaGroup_(params.getParameter("maxDeltaEtaGroup")),
+ verbose_(params.getUntrackedParameter("verbose")),
+ caloParticleToken_(
+ consumes>(params.getParameter("caloParticleCollection"))),
+ simClusterToken_(consumes>(params.getParameter("simClusterCollection"))),
+ recHitTokens_(edm::vector_transform(params.getParameter>("recHitCollections"),
+ [this](const edm::InputTag& tag) {return consumes(tag); })) {
+ //simHitTags_(pset.getParameter>("simHitCollections")),
+ //simHitTokens_(edm::vector_transform(caloSimHitTags_,
+ // [this](const edm::InputTag& tag) {return consumes>(tag); })) {
+ if (verbose_) {
+ std::cout << "running TreeWriter in verbose mode" << std::endl;
+ }
+
+ // define products
+ produces>(); // SimClusters
+ produces>(); // radii
+ produces>(); // rechit-base four-momenta of merged clusters
+
+ //for (auto& tag : simHitTags_) {
+ // std::string label = tag.instance();
+ // produces>(label+"ToSimClus");
+ //}
+ produces>();
+ produces();
+}
+
+HGCTruthProducer::~HGCTruthProducer() {}
+
+void HGCTruthProducer::beginStream(edm::StreamID) {}
+
+void HGCTruthProducer::endStream() {}
+
+void HGCTruthProducer::produce(edm::Event& event, const edm::EventSetup& setup) {
+ std::cout << "HGCTruthProducer::produce" << std::endl;
+ edm::ESHandle geom;
+ setup.get().get(geom);
+ recHitTools_.setGeometry(*geom);
+
+ // create unique pointers for output products
+ std::unique_ptr> mergedSimClusters = std::make_unique>();
+ std::unique_ptr> radii = std::make_unique>();
+ std::unique_ptr> vectors =
+ std::make_unique>();
+
+ // read original SimClusters
+ edm::Handle> simClusterHandle;
+ event.getByToken(simClusterToken_, simClusterHandle);
+ SimClusterCollection simClusters(*simClusterHandle);
+
+ std::vector allrechits;
+ std::unordered_map detid_to_rh_index;
+ size_t rhindex = 0;
+ for (auto& token : recHitTokens_) {
+ for (const auto& rh : event.get(token)) {
+ detid_to_rh_index[rh.detid()] = rhindex;
+ rhindex++;
+ allrechits.push_back(&rh);
+ }
+ }
+
+ // determine clustering groups
+ std::vector> groups;
+ determineSimClusterGroups(simClusters, groups, radii, vectors);
+
+ mergeSimClusters(simClusters, groups, mergedSimClusters, allrechits, detid_to_rh_index);
+
+ std::cout << "initial simclusters " << simClusters.size() << " merged: " << mergedSimClusters->size()
+ << std::endl; //DEBUG Jan
+
+ // save outputs
+ const auto& mergedSCHandle = event.put(std::move(mergedSimClusters));
+ event.put(std::move(radii));
+ event.put(std::move(vectors));
+
+ auto assocMap = std::make_unique(mergedSCHandle, simClusterHandle);
+
+ std::vector mergedIndices(simClusters.size(), 0);
+ for (size_t i = 0; i < groups.size(); i++) {
+ SimClusterRef msc(mergedSCHandle, i);
+ float mergedE = msc->impactMomentum().energy();
+ for (auto idx : groups.at(i)) {
+ mergedIndices.at(idx) = i;
+ SimClusterRef sc(simClusterHandle, idx);
+ float scE = sc->impactMomentum().energy();
+ assocMap->insert(msc, std::make_pair(sc, mergedE/scE));
+ }
+ }
+ event.put(std::move(assocMap));
+
+ // do the actual merging
+ auto assoc = std::make_unique>(mergedSCHandle);
+ edm::Association::Filler filler(*assoc);
+ filler.insert(simClusterHandle, mergedIndices.begin(), mergedIndices.end());
+ filler.fill();
+ event.put(std::move(assoc));
+ //std::unordered_map hitDetIdToIndex;
+ //for (size_t s = 0; s < mergedSimClusters->size(); s++) {
+ // const auto& sc = mergedSimClusters->at(s);
+ // for (auto& hf : sc.hits_and_fractions()) {
+ // auto entry = hitDetIdToIndex.find(hf.first);
+ // // Update SimCluster assigment if detId has been found in no other SCs or if
+ // // SC has greater fraction of energy in DetId than the SC already found
+ // if (entry == hitDetIdToIndex.end() || entry->second.second < hf.second)
+ // hitDetIdToIndex[hf.first] = {s, hf.second};
+ // }
+ //}
+
+}
+
+void HGCTruthProducer::determineSimClusterGroups(const SimClusterCollection& simClusters,
+ std::vector>& groups,
+ std::unique_ptr>& radii,
+ std::unique_ptr>& vectors) const {
+ // Note: In many places, the implementation is based on indices for performance purposes. To
+ // simplify the distinction, they are referred to as {i,j}Name, such as iSC for sim clusters or iH
+ // for hits. Sometimes, lists (vectors) of these indices are used. The index to those values are
+ // (e.g.) iiSC or iiH.
+ int nSimClusters = (int)simClusters.size();
+
+ // clear the radii vector and reserve space
+ radii->clear();
+ radii->resize(nSimClusters);
+ vectors->clear();
+ vectors->resize(nSimClusters);
+
+ // for the moment, only SimClusters in the HGCal are considered
+ // clusters in the HCal barrel are copied unchanged
+ std::vector scIndices;
+ for (int iSC = 0; iSC < nSimClusters; iSC++) {
+ /*
+ * Note: we cannot exclude simclusters by eta here, because some of them have wrong information.
+ * If we remove those, their hits will be wrongly labelled as noise later.
+ */
+
+ // auto absEta = fabs(simClusters[iSC].impactPoint().eta() );
+ //if (absEta == absEta && absEta >= minEta_ && absEta <= maxEta_) {
+ // cluster located in HGCal, store the index for treatment downstream
+ scIndices.push_back(iSC);
+ //} else {
+ // cluster located in HCal, only store a non-physical radius
+ (*radii)[iSC] = -1.;
+ // }
+ }
+
+ // the number of sim clusters in the HGCal
+ int nHGCalSimClusters = (int)scIndices.size();
+
+ // determine simCluster radii
+ for (const int& iSC : scIndices) {
+ // create four-vectors for each hit and the summed shower vector
+ auto atthispointthis_is_hitsandfractions_hitsAndEnergies = simClusters[iSC].hits_and_fractions();
+ int nHits = (int)atthispointthis_is_hitsandfractions_hitsAndEnergies.size();
+ math::XYZTLorentzVectorD showerVector;
+ std::vector hitVectors;
+ hitVectors.resize(atthispointthis_is_hitsandfractions_hitsAndEnergies.size());
+ std::vector hitVectorIndices;
+ hitVectorIndices.resize(atthispointthis_is_hitsandfractions_hitsAndEnergies.size());
+ for (int iH = 0; iH < nHits; iH++) {
+ // get the hit position and create the four-vector, assuming no mass
+ GlobalPoint position = recHitTools_.getPosition(atthispointthis_is_hitsandfractions_hitsAndEnergies[iH].first);
+ float energy = atthispointthis_is_hitsandfractions_hitsAndEnergies[iH].second;
+ float scale = energy / position.mag();
+ hitVectors[iH].SetPxPyPzE(position.x() * scale, position.y() * scale, position.z() * scale, energy);
+ showerVector += hitVectors[iH];
+ hitVectorIndices[iH] = iH;
+ }
+ //overwrites whatever was there before
+ //Jan: hack in the propagated vector
+ showerVector = simClusters[iSC].impactPoint();
+
+ // sort the hit vector indices according to distance to the shower vector, closest first
+ // notice the change from iH to iiH as from now on, the hitVectorIndices hold indices referring
+ // to the actual index in hitVectors
+ sort(
+ hitVectorIndices.begin(), hitVectorIndices.end(), [&showerVector, &hitVectors](const int& iiH, const int& jjH) {
+ return deltaR(showerVector, hitVectors[iiH]) < deltaR(showerVector, hitVectors[jjH]);
+ });
+
+ // loop through hit vectors in order of the sorted indices until the target energy is reached
+ float radius = deltaR(showerVector, hitVectors[hitVectorIndices[nHits - 1]]);
+ if (nHits >= minHitsRadius_) {
+ float sumEnergy = 0.;
+
+ for (int iiH = 0; iiH < nHits; iiH++) {
+ int iH = hitVectorIndices[iiH];
+ sumEnergy += hitVectors[iH].E();
+ // define the radius as the average of two or three guidance points close to the requested
+ // shower containment fraction
+ if (sumEnergy / showerVector.E() > showerContainment_) {
+ if (iiH >= 1 && iiH <= nHits - 2) {
+ float r1 = deltaR(showerVector, hitVectors[iH]);
+ float f1 = sumEnergy / showerVector.E();
+ float r2 = deltaR(showerVector, hitVectors[hitVectorIndices[iiH - 1]]);
+ float r3 = deltaR(showerVector, hitVectors[hitVectorIndices[iiH + 1]]);
+ float f2 = (sumEnergy - hitVectors[hitVectorIndices[iiH - 1]].E()) / showerVector.E();
+ float f3 = (sumEnergy + hitVectors[hitVectorIndices[iiH + 1]].E()) / showerVector.E();
+ radius = getShowerRadius(showerContainment_, r1, f1, r2, f2, r3, f3);
+ } else if (iiH >= 1) {
+ float r1 = deltaR(showerVector, hitVectors[iH]);
+ float f1 = sumEnergy / showerVector.E();
+ float r2 = deltaR(showerVector, hitVectors[hitVectorIndices[iiH - 1]]);
+ float f2 = (sumEnergy - hitVectors[hitVectorIndices[iiH - 1]].E()) / showerVector.E();
+ radius = getShowerRadius(showerContainment_, r1, f1, r2, f2);
+ } else if (iiH <= nHits - 2) {
+ float r1 = deltaR(showerVector, hitVectors[iH]);
+ float f1 = sumEnergy / showerVector.E();
+ float r3 = deltaR(showerVector, hitVectors[hitVectorIndices[iiH + 1]]);
+ float f3 = (sumEnergy + hitVectors[hitVectorIndices[iiH + 1]].E()) / showerVector.E();
+ radius = getShowerRadius(showerContainment_, r1, f1, r3, f3);
+ }
+ }
+ }
+ }
+
+ (*vectors)[iSC] = showerVector;
+ (*radii)[iSC] = radius;
+ }
+
+ // sort the cluster indices by eta
+ sort(scIndices.begin(), scIndices.end(), [&simClusters](const int& iSC, const int& jSC) {
+ return simClusters[iSC].impactPoint().eta() < simClusters[jSC].impactPoint().eta();
+ });
+
+ // define a vector of chained indices for easy lookup of closeby clusters
+ // this might look trivial, but the next/prev pointers are changed dynamically downstream to
+ // describe / close holes in the chain during cluster merging
+ std::vector chain;
+ for (int iiSC = 0; iiSC < nHGCalSimClusters; iiSC++) {
+ // add a new chain element
+ int iSC = scIndices[iiSC];
+ chain.push_back(new ChainIndex(iSC));
+
+ // connect elements
+ if (iiSC > 0) {
+ chain[iiSC]->prev = chain[iiSC - 1];
+ chain[iiSC - 1]->next = chain[iiSC];
+ }
+ }
+
+ while (chain.size() > 0) {
+ // define the merging group, seeded by the current front-most element of the chain
+ // (fancy talk for the sim cluster with the smallest eta that hasn't been checked yet)
+ std::vector chainIndices = {chain[0]};
+ std::vector group = {chain[0]->iSC};
+ int iG = 0;
+
+ while (iG < (int)chainIndices.size()) {
+ // look for nearby clusters of the sim cluster in the group referred to by iG
+ // this approach mimics recursion as the group might be extended within this while loop
+ ChainIndex* curr = chainIndices[iG];
+
+ // check clusters with smaller eta
+ ChainIndex* prev = curr->prev;
+ while (prev != nullptr) {
+ // do nothing when the previous element is already in the group
+ if (std::find(group.begin(), group.end(), prev->iSC) == group.end()) {
+ // stop when the difference in eta is too high already
+ const auto& currCluster = simClusters[curr->iSC];
+ const auto& prevCluster = simClusters[prev->iSC];
+ auto dEta =
+ fabs(currCluster.impactPoint().eta() - prevCluster.impactPoint().eta()); // fabs actually not required
+ if (dEta > maxDeltaEtaGroup_) {
+ break;
+ }
+ // check if the two clusters should be merged
+ // TODO: shouldn't the check consider all elements already in the group? ordering issue ahead?
+ if (checkSimClusterMerging(curr->iSC, prev->iSC, simClusters, radii, vectors)) {
+ // add the element to the group
+ chainIndices.push_back(prev);
+ group.push_back(prev->iSC);
+ }
+ }
+
+ // go backward
+ prev = prev->prev;
+ }
+
+ // check clusters with higher eta
+ ChainIndex* next = curr->next;
+ while (next != nullptr) {
+ // do nothing when the next element is already in the group
+ if (std::find(group.begin(), group.end(), next->iSC) == group.end()) {
+ // stop when the difference in eta is too high already
+ const auto& currCluster = simClusters[curr->iSC];
+ const auto& nextCluster = simClusters[next->iSC];
+ auto dEta =
+ fabs(nextCluster.impactPoint().eta() - currCluster.impactPoint().eta()); // fabs actually not required
+ if (dEta > maxDeltaEtaGroup_) {
+ break;
+ }
+ // check if the two clusters should be merged
+ // TODO: see not above
+ if (checkSimClusterMerging(next->iSC, curr->iSC, simClusters, radii, vectors)) {
+ // add the element to the group
+ chainIndices.push_back(next);
+ group.push_back(next->iSC);
+ }
+ }
+
+ // go forward
+ next = next->next;
+ }
+
+ // remove the current element and reconnect the chain
+ if (curr->prev != nullptr) {
+ if (curr->next != nullptr) {
+ curr->prev->next = curr->next;
+ curr->next->prev = curr->prev;
+ } else {
+ curr->prev->next = nullptr;
+ }
+ } else if (curr->next != nullptr) {
+ curr->next->prev = nullptr;
+ }
+
+ std::vector::iterator it = std::find(chain.begin(), chain.end(), curr);
+ delete *it;
+ chain.erase(it);
+
+ // choose the next element in the group for the next iteration
+ iG++;
+ }
+
+ // store the actual group of sim cluster indices
+ groups.push_back(group);
+ }
+
+ // group debugging
+ std::cout << "created " << groups.size() << " groups" << std::endl;
+ double sumSize = 0.;
+ for (const auto& g : groups) {
+ sumSize += g.size();
+ }
+ std::cout << "average size: " << (sumSize / groups.size()) << std::endl;
+
+ // radius debugging
+ // {
+ // TCanvas* canvas = new TCanvas();
+ // TH1F* hist = new TH1F(("h" + std::to_string(firstI)).c_str(), ";Radius;# Entries", 51, -0.02, 2.02);
+
+ // int nRadii = 0;
+ // double sumRadius = 0.;
+ // for (const float &r : radii) {
+ // // if (isinf(r)) std::cout << "got inf radius!" << std::endl;
+ // hist->Fill(r);
+ // if (r > 0) {
+ // sumRadius += r;
+ // nRadii++;
+ // }
+ // }
+ // std::cout << "sum: " << sumRadius << std::endl;
+ // std::cout << "n radii: " << nRadii << std::endl;
+ // double meanRadius = sumRadius / nRadii;
+ // std::cout << "mean: " << meanRadius << std::endl;
+ // double sumDiffRadius2 = 0.;
+ // for (const float &r : radii) {
+ // if (r > 0) {
+ // double diffRadius = r - meanRadius;
+ // sumDiffRadius2 += diffRadius * diffRadius;
+ // }
+ // }
+ // double varianceRadius = sumDiffRadius2 / (nRadii - 1);
+ // double stdRadius = sqrt(varianceRadius);
+ // std::cout << "radius mean: " << meanRadius << ", std: " << stdRadius << std::endl;
+
+ // hist->Draw();
+ // canvas->SaveAs(("hist_" + std::to_string(firstI) + ".pdf").c_str());
+
+ // delete hist;
+ // delete canvas;
+ // }
+}
+
+bool HGCTruthProducer::checkSimClusterMerging(
+ int iSC,
+ int jSC,
+ const SimClusterCollection& simClusters,
+ const std::unique_ptr>& radii,
+ const std::unique_ptr>& vectors) const {
+ // get some values
+ auto dR = deltaR(vectors->at(iSC), vectors->at(jSC));
+
+ //Jan
+ return dR < 0.016;
+
+ auto radius1 = radii->at(iSC);
+ auto radius2 = radii->at(jSC);
+
+ // define a combined radius using error propagation rules (as the pull is defined statistics-like)
+ auto combinedRadius = pow(radius1 * radius1 + radius2 * radius2, 0.5);
+
+ // when both clusters have no radius, i.e., they are coming from single hits, make a quick
+ // decision solely based on dR (FREE PARAMETER)
+ if (combinedRadius == 0) {
+ return dR < 0.005;
+ }
+
+ // if at least one clusters has a non-zero radius, define a pull
+ auto pull = dR / combinedRadius;
+
+ // use a simple merging criterion based on the pull (FREE PARAMETER)
+ return pull < 0.5;
+
+ // TODO: exploit information about common ancestors (e.g. useful for e/gamma, pi0, ...)
+}
+
+void HGCTruthProducer::mergeSimClusters(const SimClusterCollection& simClusters,
+ const std::vector>& groups,
+ std::unique_ptr& mergedSimClusters,
+ const std::vector& rechits,
+ const std::unordered_map& rh_detid_to_idx) const {
+ for (std::vector group : groups) {
+ if (group.size() == 0) {
+ continue;
+ }
+
+ //DEBUG >>>
+
+ auto groupcp=group;
+ auto it = std::unique (groupcp.begin(),groupcp.end());
+ groupcp.resize(std::distance(groupcp.begin(),it));
+
+ if(groupcp.size() != group.size())
+ throw std::runtime_error("not unique group");
+
+ // <<<< DEBUG end
+
+ // sort group elements by sim cluster energies
+ sort(group.begin(), group.end(), [&simClusters](const int& iSC, const int& jSC) {
+ return simClusters[iSC].energy() > simClusters[jSC].energy();
+ });
+
+ // get a vector of all sim tracks, also store the total energy
+ ROOT::Math::LorentzVector> combinedmomentum ;
+ ROOT::Math::LorentzVector> combinedimpact;
+ std::vector mergedSCTracks;
+
+ //std::cout << "grouping: " << std::endl;
+ SimCluster newsc;
+ for (const int& iSC : group) {
+ // there is currently only one track per sim clusters
+ newsc += simClusters[iSC];
+ auto& scTracks = simClusters[iSC].g4Tracks();
+ mergedSCTracks.insert(mergedSCTracks.end(), scTracks.begin(), scTracks.end());
+ }
+
+ // TODO: determine the pdg id using infos about common ancestors
+ // For now just take the PDG of the highest energy track
+ std::sort(mergedSCTracks.begin(), mergedSCTracks.end(), [](auto& t1, auto& t2)
+ { return t1.momentum().E() > t2.momentum().E(); });
+
+ int pdgId = mergedSCTracks.front().type();
+
+ ROOT::Math::LorentzVector> combinedmomentumD(
+ combinedmomentum.x(), combinedmomentum.y(), combinedmomentum.z(), combinedmomentum.t());
+ //still needs to be added
+ newsc.setPdgId(pdgId);
+
+ //std::cout << "ESum group: " << esum << " vs " << combinedmomentum.E() << std::endl;
+
+ // create the cluster
+ mergedSimClusters->emplace_back(newsc);
+ auto& cluster = mergedSimClusters->back();
+
+ // fill hits and energy fractions
+ // since we merge, we have to care about duplicate hits and this might become drastically slow
+ // therefore, use the default addRecHitAndFraction method for the first cluster, and then
+ // loop over the remaining ones to invoke the slower but necessary addDuplicateRecHitAndFraction
+
+ for (const auto& hf : simClusters[group[0]].hits_and_fractions()) {
+ cluster.addRecHitAndFraction(hf.first, hf.second);
+ }
+ for (int iiSC = 1; iiSC < (int)group.size(); iiSC++) {
+ int iSC = group[iiSC];
+ for (const auto& hf : simClusters[iSC].hits_and_fractions()) {
+ cluster.addDuplicateRecHitAndFraction(hf.first, hf.second);
+ }
+ }
+
+ //DEBUG >>>>
+ cluster.simEnergy();
+
+
+ }
+}
+
+DEFINE_FWK_MODULE(HGCTruthProducer);
+
+// unused helper methods, to be removed
+// bool vectorsIntersect(const std::vector &v1, const std::vector &v2) {
+// // for large vectors, it would be useful to loop through sorted indices and stop early
+// // when no match of a number of v1 is possible with any element in v2, but when comparing
+// // vectors of incoming vertices, i.e., parent particles, those vectors are rather small
+// for (const int &i1 : v1) {
+// for (const int &i2 : v2) {
+// if (i1 == i2) {
+// return true;
+// }
+// }
+// }
+// return false;
+
+// // second possible approach based on set intersections, but this one compares all elements
+// // without stopping early
+// // std::vector result;
+// // std::set_intersection(v1.begin(), v1.end(), v2.begin(), v2.end(), back_inserter(result));
+// // return result.size() > 0;
+// }
+//
+// float getShowerRadius(float f, float r1, float f1, float r2, float f2) {
+// float b = log(log(1.f - f2) / log(1.f - f1)) / log(r2 / r1);
+// float a = - log(1.f - f1) / pow(r1, b);
+// float res = pow(-log(1.f - f) / a, 1.f / b);
+// if (isinf(res)) {
+// std::cout << "produced inf" << std::endl;
+// std::cout << "f : " << f << std::endl;
+// std::cout << "r1: " << r1 << std::endl;
+// std::cout << "f1: " << f1 << std::endl;
+// std::cout << "r2: " << r2 << std::endl;
+// std::cout << "f2: " << f2 << std::endl;
+// std::cout << "--------------" << std::endl;
+// }
+// return pow(-log(1.f - f) / a, 1.f / b);
+// }
+//
+// float getShowerRadius(float f, float r1, float f1, float r2, float f2, float r3, float f3) {
+// float d12 = getShowerRadius(f, r1, f1, r2, f2);
+// float d13 = getShowerRadius(f, r1, f1, r3, f3);
+// float d23 = getShowerRadius(f, r2, f2, r3, f3);
+// return (d12 + d13 + d23) / 3.f;
+// }
diff --git a/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc b/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc
new file mode 100644
index 0000000000000..c7b1b0e71e0ad
--- /dev/null
+++ b/HGCSimTruth/HGCSimTruth/plugins/SimClusterTreeMerger.cc
@@ -0,0 +1,1293 @@
+#include
+#include
+#include
+#include
+#include
+#include
+#include