Add PFParticle outline and ntuples

PhysicsTools/NanoAOD/plugins/BuildFile.xml

@@ -2,6 +2,8 @@
 <use name="FWCore/ParameterSet"/>
 <use name="FWCore/ServiceRegistry"/>
 <use name="FWCore/Utilities"/>
+<use name="Geometry/HGCalGeometry"/>
+<use name="Geometry/CaloGeometry"/>
 <use name="DataFormats/Candidate"/>
 <use name="DataFormats/PatCandidates"/>
 <use name="CommonTools/Egamma"/>
@@ -11,6 +13,8 @@
 <use name="RecoVertex/VertexTools"/>
 <use name="RecoVertex/VertexPrimitives"/>
 <use name="DataFormats/L1TGlobal"/>
+<use name="SimDataFormats/CaloAnalysis"/>
+<use name="SimDataFormats/Track"/>
 <use name="IOPool/Provenance"/>
 <use name="DQMServices/Core"/>
 <use name="CondFormats/BTauObjects"/>

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 <vector>
+#include <iostream>
+
+template <typename T, typename M>
+class ObjectIndexFromAssociationTableProducer : public edm::global::EDProducer<> {
+public:
+  ObjectIndexFromAssociationTableProducer(edm::ParameterSet const& params)
+      : objName_(params.getParameter<std::string>("objName")),
+        branchName_(params.getParameter<std::string>("branchName")),
+        doc_(params.getParameter<std::string>("docString")),
+        src_(consumes<T>(params.getParameter<edm::InputTag>("src"))),
+        objMap_(consumes<edm::Association<M>>(params.getParameter<edm::InputTag>("objMap"))),
+        cut_(params.getParameter<std::string>("cut"), true) {
+    produces<nanoaod::FlatTable>();
+  }
+
+  ~ObjectIndexFromAssociationTableProducer() override {}
+
+  void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
+    edm::Handle<T> objs;
+    iEvent.getByToken(src_, objs);
+
+    edm::Handle<edm::Association<M>> assoc;
+    iEvent.getByToken(objMap_, assoc);
+
+    std::vector<int> keys;
+    for (unsigned int i = 0; i < objs->size(); ++i) {
+        edm::Ref<T> tk(objs, i);
+      if (cut_(*tk)) {
+        edm::Ref<M> match = (*assoc)[tk];
+        int key = match.isNonnull() ? match.key() : -1;
+        keys.emplace_back(key);
+      }
+    }
+
+    auto tab = std::make_unique<nanoaod::FlatTable>(keys.size(), objName_, false, true);
+    tab->addColumn<int>(branchName_ + "Idx", keys, doc_);
+
+    iEvent.put(std::move(tab));
+  }
+
+protected:
+  const std::string objName_, branchName_, doc_;
+  const edm::EDGetTokenT<T> src_;
+  const edm::EDGetTokenT<edm::Association<M>> objMap_;
+  const StringCutObjectSelector<typename T::value_type> cut_;
+};
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+typedef ObjectIndexFromAssociationTableProducer<edm::SimTrackContainer, SimClusterCollection> SimTrackToSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer<edm::PCaloHitContainer, SimClusterCollection> CaloHitToSimClusterIndexTableProducer;
+typedef ObjectIndexFromAssociationTableProducer<SimClusterCollection, CaloParticleCollection> SimClusterToCaloParticleIndexTableProducer;
+DEFINE_FWK_MODULE(SimTrackToSimClusterIndexTableProducer);
+DEFINE_FWK_MODULE(CaloHitToSimClusterIndexTableProducer);
+DEFINE_FWK_MODULE(SimClusterToCaloParticleIndexTableProducer);

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 <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/Records/interface/GlobalTrackingGeometryRecord.h"
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+
+#include <vector>
+#include <iostream>
+
+template <typename T>
+class PositionFromDetIDTableProducer : public edm::stream::EDProducer<> {
+public:
+  PositionFromDetIDTableProducer(edm::ParameterSet const& params)
+      : name_(params.getParameter<std::string>("name")),
+        doc_(params.getParameter<std::string>("doc")),
+        src_(consumes<T>(params.getParameter<edm::InputTag>("src"))),
+        cut_(params.getParameter<std::string>("cut"), true) {
+    produces<nanoaod::FlatTable>();
+  }
+
+  ~PositionFromDetIDTableProducer() override {}
+
+  void beginRun(const edm::Run&, const edm::EventSetup& iSetup) {
+      // TODO: check that the geometry exists
+      iSetup.get<CaloGeometryRecord>().get(caloGeom_);
+      iSetup.get<GlobalTrackingGeometryRecord>().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<const HGCalGeometry*>(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<T> objs;
+    iEvent.getByToken(src_, objs);
+
+    std::vector<float> xvals;
+    std::vector<float> yvals;
+    std::vector<float> 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<nanoaod::FlatTable>(xvals.size(), name_, false, true);
+    tab->addColumn<float>("x", xvals, "x position");
+    tab->addColumn<float>("y", yvals, "y position");
+    tab->addColumn<float>("z", zvals, "z position");
+
+    iEvent.put(std::move(tab));
+  }
+
+protected:
+  const std::string name_, doc_;
+  const edm::EDGetTokenT<T> src_;
+  const StringCutObjectSelector<typename T::value_type> cut_;
+  edm::ESHandle<CaloGeometry> caloGeom_;
+  edm::ESHandle<GlobalTrackingGeometry> trackGeom_;
+
+};
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+typedef PositionFromDetIDTableProducer<std::vector<PCaloHit>> PCaloHitPositionFromDetIDTableProducer;
+typedef PositionFromDetIDTableProducer<std::vector<PSimHit>> PSimHitPositionFromDetIDTableProducer;
+DEFINE_FWK_MODULE(PCaloHitPositionFromDetIDTableProducer);
+DEFINE_FWK_MODULE(PSimHitPositionFromDetIDTableProducer);

PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc

@@ -3,6 +3,27 @@
 #include "DataFormats/Candidate/interface/Candidate.h"
 typedef SimpleFlatTableProducer<reco::Candidate> SimpleCandidateFlatTableProducer;
 
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+typedef SimpleFlatTableProducer<SimCluster> SimpleSimClusterFlatTableProducer;
+
+#include "SimDataFormats/CaloHit/interface/PCaloHit.h"
+typedef SimpleFlatTableProducer<PCaloHit> SimplePCaloHitFlatTableProducer;
+
+#include "SimDataFormats/TrackingHit/interface/PSimHit.h"
+typedef SimpleFlatTableProducer<PSimHit> SimplePSimHitFlatTableProducer;
+
+#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
+typedef SimpleFlatTableProducer<TrackingParticle> SimpleTrackingParticleFlatTableProducer;
+#include "SimDataFormats/PFAnalysis/interface/PFParticle.h"
+typedef SimpleFlatTableProducer<PFParticle> SimplePFParticleFlatTableProducer;
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+typedef SimpleFlatTableProducer<CaloParticle> SimpleCaloParticleFlatTableProducer;
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
+typedef SimpleFlatTableProducer<reco::PFCandidate> SimplePFCandidateFlatTableProducer;
+
+#include "SimDataFormats/Track/interface/SimTrack.h"
+typedef SimpleFlatTableProducer<SimTrack> SimpleSimTrackFlatTableProducer;
+
 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
 typedef EventSingletonSimpleFlatTableProducer<GenEventInfoProduct> SimpleGenEventFlatTableProducer;
 
@@ -19,9 +40,21 @@ typedef SimpleFlatTableProducer<CTPPSLocalTrackLite> SimpleLocalTrackFlatTablePr
 typedef EventSingletonSimpleFlatTableProducer<math::XYZPointF> 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);
+DEFINE_FWK_MODULE(SimpleXYZPointFlatTableProducer);
+=======
+DEFINE_FWK_MODULE(SimpleXYZPointFlatTableProducer);
+>>>>>>> Add PFParticle outline and ntuples

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'),
+    )
+)
+

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)

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)