SimTauCPLink - Adding a link between GenTaus and CaloParticles

RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py

@@ -35,6 +35,7 @@
       'keep *_ticlSimTracksters_*_*',
       'keep *_ticlSimTICLCandidates_*_*',
       'keep *_ticlSimTrackstersFromCP_*_*',
+      'keep *_SimTau*_*_*'
       )
     )
 TICL_FEVT.outputCommands.extend(TICL_RECO.outputCommands)

SimCalorimetry/HGCalAssociatorProducers/plugins/SimTauCPLinkProducer.cc

@@ -0,0 +1,133 @@
+// -*- C++ -*-
+//
+//
+// Authors:  Marco Rovere, Andreas Gruber
+//         Created:  Mon, 16 Oct 2023 14:24:35 GMT
+//
+//
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/global/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+#include "DataFormats/HGCalReco/interface/TICLCandidate.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h"
+
+class SimTauProducer : public edm::global::EDProducer<> {
+public:
+  explicit SimTauProducer(const edm::ParameterSet&);
+  ~SimTauProducer() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+private:
+  void produce(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;
+
+  void buildSimTau(SimTauCPLink&,
+                   uint8_t,
+                   int,
+                   const reco::GenParticle&,
+                   int,
+                   edm::Handle<std::vector<CaloParticle>>,
+                   const std::vector<int>&) const;
+  // ----------member data ---------------------------
+  const edm::EDGetTokenT<std::vector<CaloParticle>> caloParticles_token_;
+  const edm::EDGetTokenT<std::vector<reco::GenParticle>> genParticles_token_;
+  const edm::EDGetTokenT<std::vector<int>> genBarcodes_token_;
+};
+
+SimTauProducer::SimTauProducer(const edm::ParameterSet& iConfig)
+    : caloParticles_token_(consumes<std::vector<CaloParticle>>(iConfig.getParameter<edm::InputTag>("caloParticles"))),
+      genParticles_token_(consumes<reco::GenParticleCollection>(iConfig.getParameter<edm::InputTag>("genParticles"))),
+      genBarcodes_token_(consumes<std::vector<int>>(iConfig.getParameter<edm::InputTag>("genBarcodes"))) {
+  produces<std::vector<SimTauCPLink>>();
+}
+
+void SimTauProducer::buildSimTau(SimTauCPLink& t,
+                                 uint8_t generation,
+                                 int resonance_idx,
+                                 const reco::GenParticle& gen_particle,
+                                 int gen_particle_key,
+                                 edm::Handle<std::vector<CaloParticle>> calo_particle_h,
+                                 const std::vector<int>& gen_particle_barcodes) const {
+  const auto& caloPartVec = *calo_particle_h;
+  auto& daughters = gen_particle.daughterRefVector();
+  bool is_leaf = (daughters.empty());
+
+  if (is_leaf) {
+    LogDebug("SimTauProducer").format(" TO BE SAVED {} ", resonance_idx);
+    auto const& gen_particle_barcode = gen_particle_barcodes[gen_particle_key];
+    auto const& found_in_caloparticles = std::find_if(caloPartVec.begin(), caloPartVec.end(), [&](const auto& p) {
+      return p.g4Tracks()[0].genpartIndex() == gen_particle_barcode;
+    });
+    if (found_in_caloparticles != caloPartVec.end()) {
+      auto calo_particle_idx = (found_in_caloparticles - caloPartVec.begin());
+      t.calo_particle_leaves.push_back(CaloParticleRef(calo_particle_h, calo_particle_idx));
+      t.leaves.push_back(
+          {gen_particle.pdgId(), resonance_idx, (int)t.calo_particle_leaves.size() - 1, gen_particle_key});
+      LogDebug("SimTauProducer").format(" CP {} {}", calo_particle_idx, caloPartVec[calo_particle_idx]);
+    } else {
+      t.leaves.push_back({gen_particle.pdgId(), resonance_idx, -1, gen_particle_key});
+    }
+    return;
+  } else if (generation != 0) {
+    t.resonances.push_back({gen_particle.pdgId(), resonance_idx});
+    resonance_idx = t.resonances.size() - 1;
+    LogDebug("SimTauProducer").format(" RESONANCE/INTERMEDIATE {} ", resonance_idx);
+  }
+
+  ++generation;
+  for (auto daughter = daughters.begin(); daughter != daughters.end(); ++daughter) {
+    int gen_particle_key = (*daughter).key();
+    LogDebug("SimTauProducer").format(" gen {} {} {} ", generation, gen_particle_key, (*daughter)->pdgId());
+    buildSimTau(t, generation, resonance_idx, *(*daughter), gen_particle_key, calo_particle_h, gen_particle_barcodes);
+  }
+}
+
+void SimTauProducer::produce(edm::StreamID, edm::Event& iEvent, edm::EventSetup const& iSetup) const {
+  using namespace edm;
+
+  auto caloParticles_h = iEvent.getHandle(caloParticles_token_);
+  const auto& genParticles = iEvent.get(genParticles_token_);
+  const auto& genBarcodes = iEvent.get(genBarcodes_token_);
+
+  auto tauDecayVec = std::make_unique<std::vector<SimTauCPLink>>();
+  for (auto const& g : genParticles) {
+    auto const& flags = g.statusFlags();
+    if (std::abs(g.pdgId()) == 15 and flags.isPrompt() and flags.isDecayedLeptonHadron()) {
+      SimTauCPLink t;
+      buildSimTau(t, 0, -1, g, -1, caloParticles_h, genBarcodes);
+      t.decayMode = t.buildDecayModes();
+#ifdef EDM_ML_DEBUG
+      t.dumpFullDecay();
+      t.dump();
+#endif
+      (*tauDecayVec).push_back(t);
+    }
+  }
+  iEvent.put(std::move(tauDecayVec));
+}
+
+void SimTauProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("caloParticles", edm::InputTag("mix", "MergedCaloTruth"));
+  desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles"));
+  desc.add<edm::InputTag>("genBarcodes", edm::InputTag("genParticles"));
+  descriptions.add("SimTauProducer", desc);
+}
+
+DEFINE_FWK_MODULE(SimTauProducer);

SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h

@@ -0,0 +1,157 @@
+#ifndef SimDataFormats_SimTauCPLink_h
+#define SimDataFormats_SimTauCPLink_h
+
+#include "DataFormats/HGCalReco/interface/TICLCandidate.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+class SimTauCPLink {
+public:
+  SimTauCPLink(){};
+  ~SimTauCPLink(){};
+  struct DecayNav {
+    int pdgId_;
+    int resonance_idx_;
+    int calo_particle_idx_;
+    int gen_particle_idx_;
+    int pdgId() const { return pdgId_; }
+    int resonance_idx() const { return resonance_idx_; }
+    int calo_particle_idx() const { return calo_particle_idx_; }
+    int gen_particle_idx() const { return gen_particle_idx_; }
+
+  private:
+  };
+
+  std::vector<std::pair<int, int>> resonances;
+  std::vector<DecayNav> leaves;
+  CaloParticleRefVector calo_particle_leaves;
+  int decayMode;
+
+  enum decayModes {
+    kNull = -1,
+    kOneProng0PiZero,
+    kOneProng1PiZero,
+    kOneProng2PiZero,
+    kOneProng3PiZero,
+    kOneProngNPiZero,
+    kTwoProng0PiZero,
+    kTwoProng1PiZero,
+    kTwoProng2PiZero,
+    kTwoProng3PiZero,
+    kTwoProngNPiZero,
+    kThreeProng0PiZero,
+    kThreeProng1PiZero,
+    kThreeProng2PiZero,
+    kThreeProng3PiZero,
+    kThreeProngNPiZero,
+    kRareDecayMode,
+    kElectron,
+    kMuon
+  };
+
+  void dump(void) const {
+    for (auto const &l : leaves) {
+      LogDebug("SimTauProducer")
+          .format(
+              "L {} {} CP: {} GenP idx: {}", l.pdgId(), l.resonance_idx(), l.calo_particle_idx(), l.gen_particle_idx());
+    }
+    for (auto const &r : resonances) {
+      LogDebug("SimTauProducer").format("R {} {}", r.first, r.second);
+    }
+  }
+
+  void dumpDecay(const std::pair<int, int> &entry) const {
+    if (entry.second == -1) {  // No intermediate mother.
+      LogDebug("SimTauProducer").format("{} {}", entry.first, entry.second);
+    } else {
+      LogDebug("SimTauProducer").format("{} {} coming from: ", entry.first, entry.second);
+      auto const &mother = resonances[entry.second];
+      dumpDecay(mother);
+    }
+  }
+
+  void dumpFullDecay(void) const {
+    for (auto const &leaf : leaves) {
+      dumpDecay({leaf.pdgId(), leaf.resonance_idx()});
+    }
+  }
+
+  int buildDecayModes() {
+    int numElectrons = 0;
+    int numMuons = 0;
+    int numHadrons = 0;
+    int numPhotons = 0;
+    for (auto leaf : leaves) {
+      int pdg_id = abs(leaf.pdgId());
+      switch (pdg_id) {
+        case 22:
+          numPhotons++;
+          break;
+        case 11:
+          numElectrons++;
+          break;
+        case 13:
+          numMuons++;
+          break;
+        case 16:
+          break;
+        default:
+          numHadrons++;
+      }
+    }
+
+    if (numElectrons == 1)
+      return kElectron;
+    else if (numMuons == 1)
+      return kMuon;
+    switch (numHadrons) {
+      case 1:
+        switch (numPhotons) {
+          case 0:
+            return kOneProng0PiZero;
+          case 2:
+            return kOneProng1PiZero;
+          case 4:
+            return kOneProng2PiZero;
+          case 6:
+            return kOneProng3PiZero;
+          default:
+            return kOneProngNPiZero;
+        }
+      case 2:
+        switch (numPhotons) {
+          case 0:
+            return kTwoProng0PiZero;
+          case 2:
+            return kTwoProng1PiZero;
+          case 4:
+            return kTwoProng2PiZero;
+          case 6:
+            return kTwoProng3PiZero;
+          default:
+            return kTwoProngNPiZero;
+        }
+      case 3:
+        switch (numPhotons) {
+          case 0:
+            return kThreeProng0PiZero;
+          case 2:
+            return kThreeProng1PiZero;
+          case 4:
+            return kThreeProng2PiZero;
+          case 6:
+            return kThreeProng3PiZero;
+          default:
+            return kThreeProngNPiZero;
+        }
+      default:
+        return kRareDecayMode;
+    }
+  }
+
+private:
+};
+
+#endif  //SimTauCPLink

SimDataFormats/CaloAnalysis/src/classes.h

@@ -10,3 +10,5 @@
 #include "SimDataFormats/CaloAnalysis/interface/MtdSimLayerClusterFwd.h"
 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTrackster.h"
 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTracksterFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h"
+#include "DataFormats/Common/interface/Wrapper.h"

SimDataFormats/CaloAnalysis/src/classes_def.xml

@@ -52,6 +52,12 @@
   <class name="MtdSimTracksterRefProd"/>
   <class name="MtdSimTracksterContainer"/>
 
+  <class name="SimTauCPLink" ClassVersion="3">
+   <version ClassVersion="3" checksum="213448491"/>
+  </class>
+  <class name="std::vector<SimTauCPLink>"/>
+  <class name="edm::Wrapper<std::vector<SimTauCPLink> >" />
+
   <class name="std::pair<uint64_t,float>"/>
   <class name="std::vector<std::pair<uint64_t,float>>"/>
   <class name="edm::Wrapper<std::vector<std::pair<uint64_t,float>> >"/>

SimGeneral/Debugging/plugins/BuildFile.xml

@@ -9,4 +9,6 @@
 <use name="Geometry/HcalTowerAlgo"/>
 <use name="Geometry/HGCalGeometry"/>
 <use name="Geometry/Records"/>
+<use name="CommonTools/UtilAlgos"/>
+<use name="FWCore/ServiceRegistry"/>
 <flags EDM_PLUGIN="1"/>

SimGeneral/Debugging/plugins/SimTauAnalyzer.cc

@@ -0,0 +1,65 @@
+// -*- C++ -*-
+//
+//
+// Original Author:  Andreas Gruber
+//         Created:  Mon, 16 Oct 2023 14:24:35 GMT
+//
+//
+
+// system include files
+#include <memory>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/ConsumesCollector.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "TH1.h"
+
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimTauCPLink.h"
+
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "CommonTools/UtilAlgos/interface/TFileService.h"
+
+class SimTauAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
+public:
+  explicit SimTauAnalyzer(const edm::ParameterSet&);
+  ~SimTauAnalyzer() override = default;
+
+private:
+  void analyze(const edm::Event&, const edm::EventSetup&) override;
+
+  const edm::EDGetTokenT<std::vector<SimTauCPLink>> simTau_token_;
+  TH1D* DM_histo;
+};
+
+SimTauAnalyzer::SimTauAnalyzer(const edm::ParameterSet& iConfig)
+    : simTau_token_(consumes<std::vector<SimTauCPLink>>(iConfig.getParameter<edm::InputTag>("simTau"))) {
+  usesResource(TFileService::kSharedResource);
+  edm::Service<TFileService> fs;
+  DM_histo = fs->make<TH1D>("DM_histo", "DM_histo", 20, -1, 19);
+}
+
+// ------------ method called for each event  ------------
+void SimTauAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::Handle<std::vector<SimTauCPLink>> simTau_h;
+  iEvent.getByToken(simTau_token_, simTau_h);
+
+  const auto& simTaus = *simTau_h;
+
+  for (auto const& simTau : simTaus) {
+#ifdef EDM_ML_DEBUG
+    simTau.dumpFullDecay();
+    simTau.dump();
+#endif
+    DM_histo->Fill(simTau.decayMode);
+  }
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(SimTauAnalyzer);