Fix merge mistakes

Remove PFParticle from ntuple defines (it's a WIP)

Restore pepr-specific sim cluster stuff

Fix pepr simcluster merging

Add deposited energy producer

CommonTools/RecoAlgos/plugins/SimHitRecHitAssocitionProducer.cc → CommonTools/RecoAlgos/plugins/SimClusterRecHitAssocitionProducer.cc

@@ -32,10 +32,10 @@
 //
 typedef std::pair<size_t, float> IdxAndFraction;
 
-class SimHitRecHitAssociationProducer : public edm::stream::EDProducer<> {
+class SimClusterRecHitAssociationProducer : public edm::stream::EDProducer<> {
 public:
-  explicit SimHitRecHitAssociationProducer(const edm::ParameterSet &);
-  ~SimHitRecHitAssociationProducer() override;
+  explicit SimClusterRecHitAssociationProducer(const edm::ParameterSet &);
+  ~SimClusterRecHitAssociationProducer() override;
 
 private:
   virtual void produce(edm::Event&, const edm::EventSetup&) override;
@@ -47,7 +47,7 @@ class SimHitRecHitAssociationProducer : public edm::stream::EDProducer<> {
   edm::EDGetTokenT<SimClusterCollection> scCollectionToken_;
 };
 
-SimHitRecHitAssociationProducer::SimHitRecHitAssociationProducer(const edm::ParameterSet &pset)
+SimClusterRecHitAssociationProducer::SimClusterRecHitAssociationProducer(const edm::ParameterSet &pset)
     : //caloSimhitCollectionTokens_(edm::vector_transform(pset.getParameter<std::vector<edm::InputTag>>("caloSimHits"),
       //  [this](const edm::InputTag& tag) {return mayConsume<edm::PCaloHitContainer>(tag); })),
     //trackSimhitCollectionTokens_(edm::vector_transform(pset.getParameter<edm::InputTag>("trackSimHits"),
@@ -61,44 +61,68 @@ SimHitRecHitAssociationProducer::SimHitRecHitAssociationProducer(const edm::Para
     std::string label = tag.instance();
     produces<edm::Association<SimClusterCollection>>(label+"ToSimClus");
   }
+  produces<std::unordered_map<int, float>>();
 }
 
-SimHitRecHitAssociationProducer::~SimHitRecHitAssociationProducer() {}
+SimClusterRecHitAssociationProducer::~SimClusterRecHitAssociationProducer() {}
 
 //
 // member functions
 //
 
 // ------------ method called to produce the data  ------------
-void SimHitRecHitAssociationProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
+void SimClusterRecHitAssociationProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
+  auto simClusterToRecEnergy = std::make_unique<std::unordered_map<int, float>>();
   edm::Handle<SimClusterCollection> scCollection;
   iEvent.getByToken(scCollectionToken_, scCollection);
   std::unordered_map<size_t, IdxAndFraction> hitDetIdToIndex;
+  std::unordered_map<size_t, float> hitDetIdToTotalSimFrac;
 
   for (size_t s = 0; s < scCollection->size(); s++) {
     const auto& sc = scCollection->at(s);
+    (*simClusterToRecEnergy)[s] = 0.;
     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)
+        if (entry == hitDetIdToIndex.end()) {
+            hitDetIdToTotalSimFrac[hf.first] = hf.second;
             hitDetIdToIndex[hf.first] = {s, hf.second};
+        }
+        else {
+            hitDetIdToTotalSimFrac[hf.first] += hf.second;
+            if (entry->second.second < hf.second)
+                hitDetIdToIndex[hf.first] = {s, hf.second};
+        }
     }
   }
 
+  std::unordered_map<int, float> hitDetIdToTotalRecEnergy;
+
   for (size_t i = 0; i < caloRechitCollectionTokens_.size(); i++) {
     std::string label = caloRechitTags_.at(i).instance();
     std::vector<size_t> rechitIndices;
 
     edm::Handle<edm::View<CaloRecHit>> caloRechitCollection;
     iEvent.getByToken(caloRechitCollectionTokens_.at(i), caloRechitCollection);
-    //edm::Handle<edm::PCaloHitContainer> caloSimhitCollection;
-    //iEvent.getByToken(caloSimhitCollectionTokens_.at(i), caloSimhitCollection);
 
     for (size_t h = 0; h < caloRechitCollection->size(); h++) {
         const CaloRecHit& caloRh = caloRechitCollection->at(h);
         size_t id = caloRh.detid().rawId();
+        auto entry = hitDetIdToTotalRecEnergy.find(id);
+        float energy = caloRh.energy();
+        if (entry == hitDetIdToTotalRecEnergy.end())
+            hitDetIdToTotalRecEnergy[id] = energy;
+        else
+            hitDetIdToTotalRecEnergy.at(id) += energy;
+
         int match = hitDetIdToIndex.find(id) == hitDetIdToIndex.end() ? -1 : hitDetIdToIndex.at(id).first;
+        float fraction = match != -1 ? hitDetIdToTotalSimFrac.at(id) : 1.;
+        if (simClusterToRecEnergy->find(match) == simClusterToRecEnergy->end())
+            (*simClusterToRecEnergy)[match] = energy*fraction;
+        else
+            simClusterToRecEnergy->at(match) += energy*fraction;
+
         rechitIndices.push_back(match);
     }
 
@@ -108,8 +132,9 @@ void SimHitRecHitAssociationProducer::produce(edm::Event &iEvent, const edm::Eve
     filler.fill();
     iEvent.put(std::move(assoc), label+"ToSimClus");
   }
+  iEvent.put(std::move(simClusterToRecEnergy));
 }
 
 // define this as a plug-in
-DEFINE_FWK_MODULE(SimHitRecHitAssociationProducer);
+DEFINE_FWK_MODULE(SimClusterRecHitAssociationProducer);
 

IOMC/ParticleGuns/src/SealModule.cc

@@ -23,7 +23,6 @@
 #include "IOMC/ParticleGuns/interface/RandomtXiGunProducer.h"
 #include "IOMC/ParticleGuns/interface/RandomMultiParticlePGunProducer.h"
 #include "IOMC/ParticleGuns/interface/RandomXiThetaGunProducer.h"
-#include "IOMC/ParticleGuns/interface/FlatEtaRangeGunProducer.h"
 // particle gun prototypes
 //
 
@@ -68,5 +67,3 @@ using edm::RandomMultiParticlePGunProducer;
 DEFINE_FWK_MODULE(RandomMultiParticlePGunProducer);
 using edm::RandomXiThetaGunProducer;
 DEFINE_FWK_MODULE(RandomXiThetaGunProducer);
-using edm::FlatEtaRangeGunProducer;
-DEFINE_FWK_MODULE(FlatEtaRangeGunProducer);

PhysicsTools/NanoAOD/plugins/ObjectIndexFromAssociationProducer.cc

@@ -45,7 +45,7 @@ class ObjectIndexFromAssociationTableProducer : public edm::global::EDProducer<>
 
     std::vector<int> keys;
     for (unsigned int i = 0; i < objs->size(); ++i) {
-        edm::Ref<T> tk(objs, i);
+      edm::Ref<T> tk(objs, i);
       if (cut_(*tk)) {
         edm::Ref<M> match = (*assoc)[tk];
         int key = match.isNonnull() ? match.key() : -1;

PhysicsTools/NanoAOD/plugins/ObjectPropertyFromIndexMapProducer.cc

@@ -0,0 +1,76 @@
+#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 "DataFormats/Common/interface/Association.h"
+#include "CommonTools/Utils/interface/StringCutObjectSelector.h"
+
+#include <unordered_map>
+#include <vector>
+#include <iostream>
+
+template <typename T, typename M>
+class ObjectPropertyFromIndexMapTableProducer : public edm::global::EDProducer<> {
+public:
+  ObjectPropertyFromIndexMapTableProducer(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"))),
+        mapToken_(consumes<std::unordered_map<int, M>>(params.getParameter<edm::InputTag>("valueMap"))),
+        cut_(params.getParameter<std::string>("cut"), true) {
+    produces<nanoaod::FlatTable>();
+  }
+
+  ~ObjectPropertyFromIndexMapTableProducer() override {}
+  //
+  // Because I'm not sure if this can be templated, overload instead...
+  std::unique_ptr<nanoaod::FlatTable> fillTable(const std::vector<float>& values, const std::string& objName) const {
+    auto tab = std::make_unique<nanoaod::FlatTable>(values.size(), objName, false, true);
+    tab->addColumn<float>(branchName_, values, doc_);
+    return tab;
+  }
+
+  // Because I'm not sure if this can be templated, overload instead...
+  std::unique_ptr<nanoaod::FlatTable> fillTable(const std::vector<int>& values, const std::string& objName) const {
+    auto tab = std::make_unique<nanoaod::FlatTable>(values.size(), objName, false, true);
+    tab->addColumn<int>(branchName_, values, doc_);
+    return tab;
+  }
+
+  void produce(edm::StreamID id, edm::Event& iEvent, const edm::EventSetup& iSetup) const override {
+    edm::Handle<T> objs;
+    iEvent.getByToken(src_, objs);
+
+    edm::Handle<std::unordered_map<int, M>> valueMap;
+    iEvent.getByToken(mapToken_, valueMap);
+
+    std::vector<M> values;
+    for (unsigned int i = 0; i < objs->size(); ++i) {
+      edm::Ref<T> obj(objs, i);
+      if (cut_(*obj)) {
+        if (valueMap->find(i) == valueMap->end()) 
+            throw cms::Exception("ObjectPropertyFromIndexMapTableProducer") << "No entry in value map for candidate " << i;
+        values.emplace_back(valueMap->at(i));
+      }
+    }
+
+    auto tab = fillTable(values, objName_);
+    iEvent.put(std::move(tab));
+  }
+
+protected:
+  const std::string objName_, branchName_, doc_;
+  const edm::EDGetTokenT<T> src_;
+  const edm::EDGetTokenT<std::unordered_map<int, M>> mapToken_;
+  const StringCutObjectSelector<typename T::value_type> cut_;
+};
+
+#include "FWCore/Framework/interface/MakerMacros.h"
+typedef ObjectPropertyFromIndexMapTableProducer<SimClusterCollection, float> SimClusterRecEnergyTableProducer;
+DEFINE_FWK_MODULE(SimClusterRecEnergyTableProducer);

PhysicsTools/NanoAOD/plugins/SimpleFlatTableProducerPlugins.cc

@@ -15,8 +15,6 @@ typedef SimpleFlatTableProducer<CaloRecHit> SimpleCaloRecHitFlatTableProducer;
 
 #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"
@@ -50,7 +48,6 @@ DEFINE_FWK_MODULE(SimpleSimTrackFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleSimClusterFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleTrackingParticleFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleTrackFlatTableProducer);
-DEFINE_FWK_MODULE(SimplePFParticleFlatTableProducer);
 DEFINE_FWK_MODULE(SimplePFCandidateFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleCaloParticleFlatTableProducer);
 DEFINE_FWK_MODULE(SimpleCandidateFlatTableProducer);

PhysicsTools/NanoAOD/python/hgcRecHits_cff.py

@@ -15,7 +15,7 @@
     )
 )
 
-hgcRecHitsToSimClusters = cms.EDProducer("SimHitRecHitAssociationProducer",
+hgcRecHitsToSimClusters = cms.EDProducer("SimClusterRecHitAssociationProducer",
     caloRecHits = cms.VInputTag("HGCalRecHit:HGCEERecHits",
         "HGCalRecHit:HGCHEFRecHits", "HGCalRecHit:HGCHEBRecHits",
     ),
@@ -31,6 +31,16 @@
     docString = cms.string("SimCluster responsible for most sim energy in RecHit DetId")
 )
 
+# TODO: Pair with simclusters
+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")
+)
+
 hgcRecHitsToPFCands = cms.EDProducer("RecHitToPFCandAssociationProducer",
     caloRecHits = cms.VInputTag("HGCalRecHit:HGCEERecHits",
         "HGCalRecHit:HGCHEFRecHits", "HGCalRecHit:HGCHEBRecHits",
@@ -92,6 +102,7 @@
 
 hgcRecHitsSequence = cms.Sequence(hgcEERecHitsTable+hgcHEbackRecHitsTable+hgcHEfrontRecHitsTable
                 +hgcRecHitsToSimClusters
+                +simClusterRecEnergyTable 
                 +hgcRecHitsToPFCands
                 +hgcEERecHitsToPFCandTable+hgcHEfrontRecHitsToPFCandTable+hgcHEbackRecHitsToPFCandTable
                 +hgcEERecHitsPositionTable

PhysicsTools/NanoAOD/python/simClusters_cff.py

@@ -6,17 +6,12 @@
     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
+    singleton = cms.bool(False),
+    extension = cms.bool(False), 
     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_x = Var('g4Tracks().at(0).getPositionAtBoundary().x()', 'float', precision=14, doc='x position'),
         impactPoint_y = Var('g4Tracks().at(0).getPositionAtBoundary().y()', 'float', precision=14, doc='y position'),
         impactPoint_z = Var('g4Tracks().at(0).getPositionAtBoundary().z()', 'float', precision=14, doc='z position'),
@@ -25,10 +20,10 @@
         # 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'),
-        #boundaryPt = Var('impactMomentum.pt()', 'float', precision=14, doc='magnitude of four vector'),
+        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'),
+        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'),
     )
@@ -43,21 +38,4 @@
     docString = cms.string("Index of CaloPart containing SimCluster")
 )
 
-hgcSimTruth = 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"
-
 simClusterTables = cms.Sequence(simClusterTable+simClusterToCaloPartTable)
-
-mergedSimClusterTables = cms.Sequence(hgcSimTruth+mergedSimClusterTable+simClusterToMergedSCTable)

PhysicsTools/NanoAOD/python/trackSimHits_cff.py

@@ -1,8 +1,5 @@
 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"),

SimDataFormats/Associations/src/classes_def.xml

@@ -38,6 +38,9 @@
   <class name="map<unsigned int,vector<pair<unsigned int,pair<float,float> > > >" />
   <class name="vector<pair<unsigned int,pair<float,float> > >" />
 
+  <class name="edm::Wrapper<std::unordered_map<int,float>>"/>
+  <class name="std::unordered_map<int,float>"/>
+
   <class name="hgcal::SimToRecoCollectionWithSimClusters" >
    <field name="transientMap_" transient="true"/>
   </class>