diff --git a/PhysicsTools/NanoAOD/plugins/BuildFile.xml b/PhysicsTools/NanoAOD/plugins/BuildFile.xml index b66dea70dbb0e..53f4751cb37fe 100644 --- a/PhysicsTools/NanoAOD/plugins/BuildFile.xml +++ b/PhysicsTools/NanoAOD/plugins/BuildFile.xml @@ -2,6 +2,8 @@ + + @@ -11,6 +13,8 @@ + + diff --git a/PhysicsTools/NanoAOD/plugins/ObjectIndexFromAssociationProducer.cc b/PhysicsTools/NanoAOD/plugins/ObjectIndexFromAssociationProducer.cc new file mode 100644 index 0000000000000..b28c4aee93aa8 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/ObjectIndexFromAssociationProducer.cc @@ -0,0 +1,72 @@ +#include "FWCore/Framework/interface/global/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/CaloAnalysis/interface/SimCluster.h" +#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h" +#include "SimDataFormats/Track/interface/SimTrack.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h" +#include "DataFormats/Common/interface/Association.h" +#include "CommonTools/Utils/interface/StringCutObjectSelector.h" + +#include +#include + +template +class ObjectIndexFromAssociationTableProducer : public edm::global::EDProducer<> { +public: + ObjectIndexFromAssociationTableProducer(edm::ParameterSet const& params) + : objName_(params.getParameter("objName")), + branchName_(params.getParameter("branchName")), + doc_(params.getParameter("docString")), + src_(consumes(params.getParameter("src"))), + objMap_(consumes>(params.getParameter("objMap"))), + cut_(params.getParameter("cut"), true) { + produces(); + } + + ~ObjectIndexFromAssociationTableProducer() override {} + + void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override { + edm::Handle objs; + iEvent.getByToken(src_, objs); + + edm::Handle> assoc; + iEvent.getByToken(objMap_, assoc); + + std::vector keys; + for (unsigned int i = 0; i < objs->size(); ++i) { + edm::Ref tk(objs, i); + if (cut_(*tk)) { + edm::Ref match = (*assoc)[tk]; + int key = match.isNonnull() ? match.key() : -1; + keys.emplace_back(key); + } + } + + auto tab = std::make_unique(keys.size(), objName_, false, true); + tab->addColumn(branchName_ + "Idx", keys, doc_); + + iEvent.put(std::move(tab)); + } + +protected: + const std::string objName_, branchName_, doc_; + const edm::EDGetTokenT src_; + const edm::EDGetTokenT> objMap_; + const StringCutObjectSelector cut_; +}; + +#include "FWCore/Framework/interface/MakerMacros.h" +typedef ObjectIndexFromAssociationTableProducer SimTrackToSimClusterIndexTableProducer; +typedef ObjectIndexFromAssociationTableProducer CaloHitToSimClusterIndexTableProducer; +typedef ObjectIndexFromAssociationTableProducer SimClusterToCaloParticleIndexTableProducer; +DEFINE_FWK_MODULE(SimTrackToSimClusterIndexTableProducer); +DEFINE_FWK_MODULE(CaloHitToSimClusterIndexTableProducer); +DEFINE_FWK_MODULE(SimClusterToCaloParticleIndexTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/PositionFromDetIDTableProducer.cc b/PhysicsTools/NanoAOD/plugins/PositionFromDetIDTableProducer.cc new file mode 100644 index 0000000000000..94ae11fb8fd86 --- /dev/null +++ b/PhysicsTools/NanoAOD/plugins/PositionFromDetIDTableProducer.cc @@ -0,0 +1,112 @@ +#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 +#include +#include +#include +#include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h" +#include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" + +#include +#include + +template +class PositionFromDetIDTableProducer : public edm::stream::EDProducer<> { +public: + PositionFromDetIDTableProducer(edm::ParameterSet const& params) + : name_(params.getParameter("name")), + doc_(params.getParameter("doc")), + src_(consumes(params.getParameter("src"))), + cut_(params.getParameter("cut"), true) { + produces(); + } + + ~PositionFromDetIDTableProducer() override {} + + void beginRun(const edm::Run&, const edm::EventSetup& iSetup) { + // TODO: check that the geometry exists + iSetup.get().get(caloGeom_); + iSetup.get().get(trackGeom_); + } + + GlobalPoint positionFromHit(const PCaloHit& hit) { + DetId id = hit.id(); + DetId::Detector det = id.det(); + int subid = (det == DetId::HGCalEE || det == DetId::HGCalHSi || det == DetId::HGCalHSc) + ? ForwardSubdetector::ForwardEmpty + : id.subdetId(); + auto geom = caloGeom_->getSubdetectorGeometry(det, subid); + + GlobalPoint position; + if (id.det() == DetId::Hcal) { + position = geom->getGeometry(id)->getPosition(); + } else if (id.det() == DetId::HGCalEE || id.det() == DetId::HGCalHSi || id.det() == DetId::HGCalHSc) { + auto hg = static_cast(geom); + position = hg->getPosition(id); + } + else { + throw cms::Exception("PositionFromDetIDTableProducer") << "Unsupported DetId type"; + } + + return position; + } + + GlobalPoint positionFromHit(const PSimHit& hit) { + auto surface = trackGeom_->idToDet(hit.detUnitId())->surface(); + //LocalPoint localPos = surface.position(); + 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_; + +}; + +#include "FWCore/Framework/interface/MakerMacros.h" +typedef PositionFromDetIDTableProducer> PCaloHitPositionFromDetIDTableProducer; +typedef PositionFromDetIDTableProducer> PSimHitPositionFromDetIDTableProducer; +DEFINE_FWK_MODULE(PCaloHitPositionFromDetIDTableProducer); +DEFINE_FWK_MODULE(PSimHitPositionFromDetIDTableProducer); diff --git a/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc b/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc index 1207a565c8a40..8426b258fd718 100644 --- a/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc +++ b/PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc @@ -3,6 +3,27 @@ #include "DataFormats/Candidate/interface/Candidate.h" typedef SimpleFlatTableProducer SimpleCandidateFlatTableProducer; +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +typedef SimpleFlatTableProducer SimpleSimClusterFlatTableProducer; + +#include "SimDataFormats/CaloHit/interface/PCaloHit.h" +typedef SimpleFlatTableProducer SimplePCaloHitFlatTableProducer; + +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +typedef SimpleFlatTableProducer SimplePSimHitFlatTableProducer; + +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +typedef SimpleFlatTableProducer SimpleTrackingParticleFlatTableProducer; +#include "SimDataFormats/PFAnalysis/interface/PFParticle.h" +typedef SimpleFlatTableProducer SimplePFParticleFlatTableProducer; +#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h" +typedef SimpleFlatTableProducer SimpleCaloParticleFlatTableProducer; +#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" +typedef SimpleFlatTableProducer SimplePFCandidateFlatTableProducer; + +#include "SimDataFormats/Track/interface/SimTrack.h" +typedef SimpleFlatTableProducer SimpleSimTrackFlatTableProducer; + #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" typedef EventSingletonSimpleFlatTableProducer SimpleGenEventFlatTableProducer; @@ -19,9 +40,21 @@ typedef SimpleFlatTableProducer SimpleLocalTrackFlatTablePr typedef EventSingletonSimpleFlatTableProducer SimpleXYZPointFlatTableProducer; #include "FWCore/Framework/interface/MakerMacros.h" +DEFINE_FWK_MODULE(SimplePCaloHitFlatTableProducer); +DEFINE_FWK_MODULE(SimplePSimHitFlatTableProducer); +DEFINE_FWK_MODULE(SimpleSimTrackFlatTableProducer); +DEFINE_FWK_MODULE(SimpleSimClusterFlatTableProducer); +DEFINE_FWK_MODULE(SimpleTrackingParticleFlatTableProducer); +DEFINE_FWK_MODULE(SimplePFParticleFlatTableProducer); +DEFINE_FWK_MODULE(SimplePFCandidateFlatTableProducer); +DEFINE_FWK_MODULE(SimpleCaloParticleFlatTableProducer); DEFINE_FWK_MODULE(SimpleCandidateFlatTableProducer); DEFINE_FWK_MODULE(SimpleGenEventFlatTableProducer); DEFINE_FWK_MODULE(SimpleHTXSFlatTableProducer); +<<<<<<< HEAD DEFINE_FWK_MODULE(SimpleProtonTrackFlatTableProducer); DEFINE_FWK_MODULE(SimpleLocalTrackFlatTableProducer); -DEFINE_FWK_MODULE(SimpleXYZPointFlatTableProducer); \ No newline at end of file +DEFINE_FWK_MODULE(SimpleXYZPointFlatTableProducer); +======= +DEFINE_FWK_MODULE(SimpleXYZPointFlatTableProducer); +>>>>>>> Add PFParticle outline and ntuples diff --git a/PhysicsTools/NanoAOD/python/caloParticles_cff.py b/PhysicsTools/NanoAOD/python/caloParticles_cff.py new file mode 100644 index 0000000000000..91617d618bded --- /dev/null +++ b/PhysicsTools/NanoAOD/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/PhysicsTools/NanoAOD/python/hgcSimHits_cff.py b/PhysicsTools/NanoAOD/python/hgcSimHits_cff.py new file mode 100644 index 0000000000000..d45e32626c810 --- /dev/null +++ b/PhysicsTools/NanoAOD/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'), + fineTrackId = Var('geantFineTrackId', 'int', precision=-1, doc='granular 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 = "SimHitHGCHEfront" + +hgcHEfrontHitsToSimClusterTable = hgcEEHitsToSimClusterTable.clone() +hgcHEfrontHitsToSimClusterTable.src = hgcHEfrontSimHitsTable.src +hgcHEfrontHitsToSimClusterTable.objName = hgcHEfrontSimHitsTable.name +hgcHEfrontHitsToSimClusterTable.objMap = "mix:simHitHGCHEfrontToSimCluster" + +hgcHEbackSimHitsTable = hgcEESimHitsTable.clone() +hgcHEbackSimHitsTable.src = "g4SimHits:HGCHitsHEback" +hgcHEbackSimHitsTable.name = "SimHitHGCHEback" + +hgcHEbackHitsToSimClusterTable = hgcEEHitsToSimClusterTable.clone() +hgcHEbackHitsToSimClusterTable.src = hgcHEbackSimHitsTable.src +hgcHEbackHitsToSimClusterTable.objName = hgcHEbackSimHitsTable.name +hgcHEbackHitsToSimClusterTable.objMap = "mix:simHitHGCHEbackToSimCluster" + +hgcEESimHitsPositionTable = cms.EDProducer("PCaloHitPositionFromDetIDTableProducer", + 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/PhysicsTools/NanoAOD/python/hgcSimTracks_cff.py b/PhysicsTools/NanoAOD/python/hgcSimTracks_cff.py new file mode 100644 index 0000000000000..8a9870c44bb0c --- /dev/null +++ b/PhysicsTools/NanoAOD/python/hgcSimTracks_cff.py @@ -0,0 +1,41 @@ +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'), + ) +) + +simTrackToSimClusterTable = cms.EDProducer("SimTrackToSimClusterIndexTableProducer", + cut = simTrackTable.cut, + src = simTrackTable.src, + objName = simTrackTable.name, + branchName = cms.string("SimCluster"), + objMap = cms.InputTag("mix:simTrackToSimCluster"), + docString = cms.string("SimCluster containing track") +) + +simTrackTables = cms.Sequence(simTrackTable+simTrackToSimClusterTable) diff --git a/PhysicsTools/NanoAOD/python/nanoHGCML_cff.py b/PhysicsTools/NanoAOD/python/nanoHGCML_cff.py new file mode 100644 index 0000000000000..ef6f44f464fc9 --- /dev/null +++ b/PhysicsTools/NanoAOD/python/nanoHGCML_cff.py @@ -0,0 +1,27 @@ +from __future__ import print_function +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import * +from hgcSimHits_cff import * +from hgcSimTracks_cff import * +from trackSimHits_cff import * +from simClusters_cff import * +from caloParticles_cff import * +from trackingParticles_cff import * +from genparticles_cff import genParticleTable +from genVertex_cff import * +from pfCands_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+simClusterTables+ \ + #pfCandTable+pfTICLCandTable \ + simTrackTables+hgcSimHitsSequence+trackerSimHitTables) diff --git a/PhysicsTools/NanoAOD/python/pfCands_cff.py b/PhysicsTools/NanoAOD/python/pfCands_cff.py new file mode 100644 index 0000000000000..389aa257e7283 --- /dev/null +++ b/PhysicsTools/NanoAOD/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("SimplePFCandidateFlatTableProducer", + 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("SimplePFCandidateFlatTableProducer", + 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/PhysicsTools/NanoAOD/python/simClusters_cff.py b/PhysicsTools/NanoAOD/python/simClusters_cff.py new file mode 100644 index 0000000000000..2d49f5f3ed0c4 --- /dev/null +++ b/PhysicsTools/NanoAOD/python/simClusters_cff.py @@ -0,0 +1,33 @@ +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, + y = Var('impactPoint().y()', 'float', precision=14, doc='y position'), + z = Var('impactPoint().z()', 'float', precision=14, doc='z position'), + impactPointX = Var('impactPoint().x()', 'float', precision=14, doc='x position'), + impactPointY = Var('impactPoint().y()', 'float', precision=14, doc='y position'), + impactPointZ = Var('impactPoint().z()', 'float', precision=14, doc='z position'), + nSimHits = Var('numberOfSimHits', 'int', precision=-1, doc='total energy of simhits'), + simEnergy = Var('simEnergy', 'float', precision=14, doc='total energy of simhits'), + trackId = Var('g4Tracks().at(0).trackId()', 'int', precision=10, doc='Geant track id'), + trackIdAtBoundary = Var('g4Tracks().at(0).getIDAtBoundary()', 'int', precision=-1, doc='Track ID at boundary'), + ) +) + +simClusterToCaloPartTable = cms.EDProducer("SimClusterToCaloParticleIndexTableProducer", + cut = simClusterTable.cut, + src = simClusterTable.src, + objName = simClusterTable.name, + branchName = cms.string("CaloPart"), + objMap = cms.InputTag("mix:simClusterToCaloParticle"), + docString = cms.string("Index of CaloPart containing SimCluster") +) + +simClusterTables = cms.Sequence(simClusterTable+simClusterToCaloPartTable) diff --git a/PhysicsTools/NanoAOD/python/trackSimHits_cff.py b/PhysicsTools/NanoAOD/python/trackSimHits_cff.py new file mode 100644 index 0000000000000..c30c5f58066cc --- /dev/null +++ b/PhysicsTools/NanoAOD/python/trackSimHits_cff.py @@ -0,0 +1,55 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.NanoAOD.common_cff import CandVars,Var +#TrackerHitsTECLowTof +#TrackerHitsPixelEndcapHighTof +#TrackerHitsPixelEndcapLowTof + +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'), + 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("PSimHitPositionFromDetIDTableProducer", + 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("PSimHitPositionFromDetIDTableProducer", + 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("PSimHitPositionFromDetIDTableProducer", + src = muonCSCHitsTable.src, + cut = muonCSCHitsTable.cut, + name = muonCSCHitsTable.name, + doc = muonCSCHitsTable.doc, +) + +trackerSimHitTables = cms.Sequence(trackerHitsPixelEndcapLowTofTable+trackerHitsPixelEndcapLowTofPositionTable+ + trackerHitsPixelBarrelLowTofTable+trackerHitsPixelBarrelLowTofPositionTable+ + muonCSCHitsTable+muonCSCHitsPositionTable) diff --git a/PhysicsTools/NanoAOD/python/trackingParticles_cff.py b/PhysicsTools/NanoAOD/python/trackingParticles_cff.py new file mode 100644 index 0000000000000..c2c7e4f4d72e0 --- /dev/null +++ b/PhysicsTools/NanoAOD/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/SimDataFormats/PFAnalysis/BuildFile.xml b/SimDataFormats/PFAnalysis/BuildFile.xml new file mode 100644 index 0000000000000..f51d00ade2d32 --- /dev/null +++ b/SimDataFormats/PFAnalysis/BuildFile.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/SimDataFormats/PFAnalysis/interface/PFParticle.h b/SimDataFormats/PFAnalysis/interface/PFParticle.h new file mode 100644 index 0000000000000..802c50af32f0e --- /dev/null +++ b/SimDataFormats/PFAnalysis/interface/PFParticle.h @@ -0,0 +1,164 @@ +#ifndef SimDataFormats_PFParticle_h +#define SimDataFormats_PFParticle_h + +#include +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +#include "DataFormats/Math/interface/Point3D.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/Math/interface/LorentzVector.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + +// +// Forward declarations +// +class TrackingVertex; +class SimTrack; +class EncodedEventId; + +class PFParticle { + friend std::ostream& operator<<(std::ostream& s, PFParticle const& tp); + +public: + typedef int Charge; ///< electric charge type + typedef math::XYZTLorentzVectorD LorentzVector; ///< Lorentz vector + typedef math::PtEtaPhiMLorentzVector PolarLorentzVector; ///< Lorentz vector + typedef math::XYZPointD Point; ///< point in the space + typedef math::XYZVectorD Vector; ///< point in the space + + /// reference to reco::GenParticle + typedef reco::GenParticleRefVector::iterator genp_iterator; + typedef std::vector::const_iterator g4t_iterator; + + /** @brief Default constructor. Note that the object will be useless until it is provided + * with a SimTrack and parent TrackingVertex. + * + * Most of the methods assume there is a SimTrack and parent TrackingVertex set, so will either + * crash or give undefined results if this isn't true. This constructor should only be used to + * create a placeholder until setParentVertex() and addG4Track() can be called. + */ + PFParticle(); + + PFParticle(const TrackingParticleRefVector& trackingParticles, const SimClusterRefVector& simClusters); + + // destructor + ~PFParticle(); + void setTrackingParticles(const TrackingParticleRefVector& refs); + void setSimClusters(const SimClusterRefVector& refs); + void addSimCluster(const SimClusterRef& sc); + void addTrackingParticle(const TrackingParticleRef& tp); + + /** @brief PDG ID. + * + * Returns the PDG ID of the first associated gen particle. If there are no gen particles associated + * then it returns type() from the first SimTrack. */ + int pdgId() const { + if (genParticles_.empty()) + return g4Tracks_[0].type(); + else + return (*genParticles_.begin())->pdgId(); + } + + /** @brief Signal source, crossing number. + * + * Note this is taken from the first SimTrack only, but there shouldn't be any SimTracks from different + * crossings in the PFParticle. */ + EncodedEventId eventId() const { return g4Tracks_[0].eventId(); } + + // Setters for G4 and reco::GenParticle + void addGenParticle(const reco::GenParticleRef& ); + void addG4Track(const SimTrack& t); + /// iterators + genp_iterator genParticle_begin() const; + genp_iterator genParticle_end() const; + g4t_iterator g4Track_begin() const; + g4t_iterator g4Track_end() const; + + // Getters for Embd and Sim Tracks + const reco::GenParticleRefVector& genParticles() const { return genParticles_; } + const std::vector& g4Tracks() const { return g4Tracks_; } + + /// @brief Electric charge. Note this is taken from the first SimTrack only. + float charge() const { return g4Tracks_[0].charge(); } + /// Gives charge in unit of quark charge (should be 3 times "charge()") + int threeCharge() const { return lrintf(3.f * charge()); } + + /// @brief Four-momentum Lorentz vector. Note this is taken from the first SimTrack only. + const LorentzVector& p4() const { return g4Tracks_[0].momentum(); } + + /// @brief spatial momentum vector + Vector momentum() const { return p4().Vect(); } + + /// @brief Vector to boost to the particle centre of mass frame. + Vector boostToCM() const { return p4().BoostToCM(); } + + /// @brief Magnitude of momentum vector. Note this is taken from the first SimTrack only. + double p() const { return p4().P(); } + + /// @brief Energy. Note this is taken from the first SimTrack only. + double energy() const { return p4().E(); } + + /// @brief Transverse energy. Note this is taken from the first SimTrack only. + double et() const { return p4().Et(); } + + /// @brief Mass. Note this is taken from the first SimTrack only. + double mass() const { return p4().M(); } + + /// @brief Mass squared. Note this is taken from the first SimTrack only. + double massSqr() const { return pow(mass(), 2); } + + /// @brief Transverse mass. Note this is taken from the first SimTrack only. + double mt() const { return p4().Mt(); } + + /// @brief Transverse mass squared. Note this is taken from the first SimTrack only. + double mtSqr() const { return p4().Mt2(); } + + /// @brief x coordinate of momentum vector. Note this is taken from the first SimTrack only. + double px() const { return p4().Px(); } + + /// @brief y coordinate of momentum vector. Note this is taken from the first SimTrack only. + double py() const { return p4().Py(); } + + /// @brief z coordinate of momentum vector. Note this is taken from the first SimTrack only. + double pz() const { return p4().Pz(); } + + /// @brief Transverse momentum. Note this is taken from the first SimTrack only. + double pt() const { return p4().Pt(); } + + /// @brief Momentum azimuthal angle. Note this is taken from the first SimTrack only. + double phi() const { return p4().Phi(); } + + /// @brief Momentum polar angle. Note this is taken from the first SimTrack only. + double theta() const { return p4().Theta(); } + + /// @brief Momentum pseudorapidity. Note this is taken from the first SimTrack only. + double eta() const { return p4().Eta(); } + + /// @brief Rapidity. Note this is taken from the first SimTrack only. + double rapidity() const { return p4().Rapidity(); } + + /// @brief Same as rapidity(). + double y() const { return rapidity(); } + + /** @brief Status word. + * + * Returns status() from the first gen particle, or -99 if there are no gen particles attached. */ + int status() const { return genParticles_.empty() ? -99 : (*genParticles_[0]).status(); } + + static const unsigned int longLivedTag; ///< long lived flag + + /// is long lived? + bool longLived() const { return status() & longLivedTag; } + +private: + /// references to G4 and reco::GenParticle tracks + std::vector g4Tracks_; + reco::GenParticleRefVector genParticles_; + SimClusterRefVector simClusters_; + TrackingParticleRefVector trackingParticles_; +}; + +#endif // SimDataFormats_PFParticle_H + diff --git a/SimDataFormats/PFAnalysis/interface/PFParticleFwd.h b/SimDataFormats/PFAnalysis/interface/PFParticleFwd.h new file mode 100644 index 0000000000000..4e85443dbfa2b --- /dev/null +++ b/SimDataFormats/PFAnalysis/interface/PFParticleFwd.h @@ -0,0 +1,16 @@ +#ifndef CaloAnalysis_PFParticleFwd_h +#define CaloAnalysis_PFParticleFwd_h +#include "DataFormats/Common/interface/Ref.h" +#include "DataFormats/Common/interface/RefProd.h" +#include "DataFormats/Common/interface/RefVector.h" +#include + +class PFParticle; +typedef std::vector PFParticleCollection; +typedef edm::Ref PFParticleRef; +typedef edm::RefVector PFParticleRefVector; +typedef edm::RefProd PFParticleRefProd; +typedef edm::RefVector PFParticleContainer; + +#endif + diff --git a/SimDataFormats/PFAnalysis/src/PFParticle.cc b/SimDataFormats/PFAnalysis/src/PFParticle.cc new file mode 100644 index 0000000000000..de0d8127e51e3 --- /dev/null +++ b/SimDataFormats/PFAnalysis/src/PFParticle.cc @@ -0,0 +1,36 @@ +#include "SimDataFormats/PFAnalysis/interface/PFParticle.h" + +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" + +#include + +PFParticle::PFParticle() { + // No operation +} + +PFParticle::~PFParticle() {} + +PFParticle::PFParticle(const TrackingParticleRefVector& trackingParticles, const SimClusterRefVector& simClusters) { + setTrackingParticles(trackingParticles); + setSimClusters(simClusters); +} + +void PFParticle::setTrackingParticles(const TrackingParticleRefVector& refs) { trackingParticles_ = refs; } + +void PFParticle::setSimClusters(const SimClusterRefVector& refs) { simClusters_ = refs; } + +void PFParticle::addSimCluster(const SimClusterRef& sc) { simClusters_.push_back(sc); } + +void PFParticle::addTrackingParticle(const TrackingParticleRef& tp) { trackingParticles_.push_back(tp); } + +void PFParticle::addGenParticle(const reco::GenParticleRef& gp) { genParticles_.push_back(gp); } + +void PFParticle::addG4Track(const SimTrack& t) { g4Tracks_.push_back(t); } + +PFParticle::genp_iterator PFParticle::genParticle_begin() const { return genParticles_.begin(); } + +PFParticle::genp_iterator PFParticle::genParticle_end() const { return genParticles_.end(); } + +PFParticle::g4t_iterator PFParticle::g4Track_begin() const { return g4Tracks_.begin(); } + +PFParticle::g4t_iterator PFParticle::g4Track_end() const { return g4Tracks_.end(); } diff --git a/SimDataFormats/PFAnalysis/src/classes.h b/SimDataFormats/PFAnalysis/src/classes.h new file mode 100644 index 0000000000000..3c0ea47e67ece --- /dev/null +++ b/SimDataFormats/PFAnalysis/src/classes.h @@ -0,0 +1,8 @@ +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" +#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h" +#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h" +#include "SimDataFormats/PFAnalysis/interface/PFParticle.h" +#include "SimDataFormats/PFAnalysis/interface/PFParticleFwd.h" +#include "DataFormats/Common/interface/Wrapper.h" + diff --git a/SimDataFormats/PFAnalysis/src/classes_def.xml b/SimDataFormats/PFAnalysis/src/classes_def.xml new file mode 100644 index 0000000000000..044b8d921029d --- /dev/null +++ b/SimDataFormats/PFAnalysis/src/classes_def.xml @@ -0,0 +1,12 @@ + + + + + + + + + + + + diff --git a/SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc b/SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc index 4f13f5c49f5a7..6389f8e1cc5ea 100644 --- a/SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc +++ b/SimGeneral/CaloAnalysis/plugins/CaloTruthAccumulator.cc @@ -147,6 +147,8 @@ class CaloTruthAccumulator : public DigiAccumulatorMixMod { std::unordered_map> &simTrackDetIdEnergyMap, const T &event, const edm::EventSetup &setup); + std::unique_ptr> fillSimHitAssociation( + edm::Handle>& inColl, edm::OrphanHandle& matchColl); const std::string messageCategory_; @@ -183,6 +185,10 @@ class CaloTruthAccumulator : public DigiAccumulatorMixMod { struct OutputCollections { std::unique_ptr pSimClusters; std::unique_ptr pCaloParticles; + std::unique_ptr> pTrackSCAssoc; + std::unique_ptr> pHGCEEHitSCAssoc; + std::unique_ptr> pHGCHEfrontHitSCAssoc; + std::unique_ptr> pHGCHEbackHitSCAssoc; }; struct calo_particles { @@ -209,6 +215,8 @@ class CaloTruthAccumulator : public DigiAccumulatorMixMod { // geometry type (0 pre-TDR; 1 TDR) int geometryType_; bool doHGCAL; + std::vector trackToSimClusterIndices_; + std::unordered_map caloHitToSimClusterIdx_; }; /* Graph utility functions */ @@ -291,11 +299,13 @@ namespace { CaloTruthAccumulator::calo_particles &caloParticles, std::unordered_multimap &simHitBarcodeToIndex, std::unordered_map> &simTrackDetIdEnergyMap, + std::unordered_map &caloHitToSimClusterIdxMap, Selector selector) : output_(output), caloParticles_(caloParticles), simHitBarcodeToIndex_(simHitBarcodeToIndex), simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap), + caloHitIdToSimClusterIdxMap_(caloHitToSimClusterIdxMap), selector_(selector) {} template void discover_vertex(Vertex u, const Graph &g) { @@ -314,6 +324,7 @@ namespace { auto &simcluster = output_.pSimClusters->back(); std::unordered_map acc_energy; for (auto const &hit_and_energy : simTrackDetIdEnergyMap_[trackIdx]) { + caloHitIdToSimClusterIdxMap_[hit_and_energy.first] = output_.pSimClusters->size()-1; acc_energy[hit_and_energy.first] += hit_and_energy.second; } for (auto const &hit_and_energy : acc_energy) { @@ -353,6 +364,7 @@ namespace { CaloTruthAccumulator::calo_particles &caloParticles_; std::unordered_multimap &simHitBarcodeToIndex_; std::unordered_map> &simTrackDetIdEnergyMap_; + std::unordered_map &caloHitIdToSimClusterIdxMap_; Selector selector_; }; } // namespace @@ -379,6 +391,16 @@ CaloTruthAccumulator::CaloTruthAccumulator(const edm::ParameterSet &config, producesCollector.produces>>("MergedCaloTruth"); } + if (true) { + producesCollector.produces>("simClusterToCaloParticle"); + producesCollector.produces>("simTrackToSimCluster"); + producesCollector.produces>("simHitHGCEEToSimCluster"); + producesCollector.produces>("simHitHGCHEfrontToSimCluster"); + producesCollector.produces>("simHitHGCHEbackToSimCluster"); + } + + iC.consumes>(simTrackLabel_); + iC.consumes>(simTrackLabel_); iC.consumes>(simVertexLabel_); iC.consumes>(genParticleLabel_); @@ -443,6 +465,8 @@ void CaloTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSe output_.pCaloParticles = std::make_unique(); m_detIdToTotalSimEnergy.clear(); + caloHitToSimClusterIdx_.clear(); + trackToSimClusterIndices_.clear(); } /** Create handle to edm::HepMCProduct here because event.getByLabel with @@ -511,6 +535,34 @@ void CaloTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup cons // save the SimCluster orphan handle so we can fill the calo particles auto scHandle = event.put(std::move(output_.pSimClusters), "MergedCaloTruth"); + output_.pTrackSCAssoc = std::make_unique>(scHandle); + edm::Association::Filler filler(*output_.pTrackSCAssoc); + filler.insert(hSimTracks, trackToSimClusterIndices_.begin(), trackToSimClusterIndices_.end()); + filler.fill(); + event.put(std::move(output_.pTrackSCAssoc), "simTrackToSimCluster"); + + for (auto const &collectionTag : collectionTags_) { + edm::Handle> hSimHits; + const bool isHGCal = (collectionTag.instance().find("HGCHits") != std::string::npos); + if (!isHGCal) + continue; + event.getByLabel(collectionTag, hSimHits); + + if (collectionTag.instance().find("HEfront") != std::string::npos) { + output_.pHGCHEfrontHitSCAssoc = fillSimHitAssociation(hSimHits, scHandle); + } + else if (collectionTag.instance().find("HEback") != std::string::npos) { + output_.pHGCHEbackHitSCAssoc = fillSimHitAssociation(hSimHits, scHandle); + } + else + output_.pHGCEEHitSCAssoc = fillSimHitAssociation(hSimHits, scHandle); + + } + event.put(std::move(output_.pHGCEEHitSCAssoc), "simHitHGCEEToSimCluster"); + event.put(std::move(output_.pHGCHEfrontHitSCAssoc), "simHitHGCHEfrontToSimCluster"); + event.put(std::move(output_.pHGCHEbackHitSCAssoc), "simHitHGCHEbackToSimCluster"); + + std::unordered_map simClusIdxToCaloParticleIdxMap_; // now fill the calo particles for (unsigned i = 0; i < output_.pCaloParticles->size(); ++i) { @@ -521,12 +573,25 @@ void CaloTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup cons } } - event.put(std::move(output_.pCaloParticles), "MergedCaloTruth"); + auto cpHandle = event.put(std::move(output_.pCaloParticles), "MergedCaloTruth"); calo_particles().swap(m_caloParticles); std::unordered_map().swap(m_detIdToTotalSimEnergy); std::unordered_multimap().swap(m_simHitBarcodeToIndex); + + std::vector caloPartIdx; + for (size_t i = 0; i < scHandle->size(); i++) { + int matchIdx = -1; + if (simClusIdxToCaloParticleIdxMap_.find(i) != simClusIdxToCaloParticleIdxMap_.end()) + matchIdx = simClusIdxToCaloParticleIdxMap_[i]; + caloPartIdx.emplace_back(matchIdx); + } + auto assocMap = std::make_unique>(cpHandle); + edm::Association::Filler cpfiller(*assocMap); + cpfiller.insert(scHandle, caloPartIdx.begin(), caloPartIdx.end()); + cpfiller.fill(); + event.put(std::move(assocMap), "simClusterToCaloParticle"); } template @@ -652,6 +717,7 @@ void CaloTruthAccumulator::accumulateEvent(const T &event, m_caloParticles, m_simHitBarcodeToIndex, simTrackDetIdEnergyMap, + caloHitToSimClusterIdx_, [&](EdgeProperty &edge_property) -> bool { // Apply selection on SimTracks in order to promote them to be // CaloParticles. The function returns TRUE if the particle satisfies @@ -663,6 +729,19 @@ void CaloTruthAccumulator::accumulateEvent(const T &event, }); depth_first_search(decay, visitor(caloParticleCreator)); + std::unordered_map trackIdxToSimClusterIdx; + for (size_t i = 0; i < output_.pSimClusters->size(); i ++) { + auto& sc = output_.pSimClusters->at(i); + for (auto tk : sc.g4Tracks()) { + trackIdxToSimClusterIdx[trackid_to_track_index.at(tk.trackId())] = i; + } + } + + for (size_t i = 0; i < hSimTracks->size(); i++) { + int idx = trackIdxToSimClusterIdx.find(i) != trackIdxToSimClusterIdx.end() ? trackIdxToSimClusterIdx[i] : -1; + trackToSimClusterIndices_.emplace_back(idx); + } + #if DEBUG boost::write_graphviz(std::cout, decay, @@ -671,6 +750,23 @@ void CaloTruthAccumulator::accumulateEvent(const T &event, #endif } +std::unique_ptr> CaloTruthAccumulator::fillSimHitAssociation( + edm::Handle>& inColl, edm::OrphanHandle& matchColl) { + std::vector matchIndices(inColl->size(), -1); + + for (size_t i = 0; i < inColl->size(); i++) { + auto& hit = inColl->at(i); + if (caloHitToSimClusterIdx_.find(hit.id()) != caloHitToSimClusterIdx_.end()) + matchIndices[i] = caloHitToSimClusterIdx_[hit.id()]; + } + + auto assocMap = std::make_unique>(matchColl); + edm::Association::Filler filler(*assocMap); + filler.insert(inColl, matchIndices.begin(), matchIndices.end()); + filler.fill(); + return assocMap; +} + template void CaloTruthAccumulator::fillSimHits(std::vector> &returnValue, std::unordered_map> &simTrackDetIdEnergyMap,