CaloAnalysis/plugins/CaloTruthAccumulator.cc

0001 #define DEBUG false
0002
0003 #if DEBUG
0004 #pragma GCC diagnostic pop
0005 #endif
0006
0007 #include <iterator>
0008 #include <memory>
0009
0010 #include <numeric>  // for std::accumulate
0011 #include <unordered_map>
0012
0013 #include "FWCore/Framework/interface/ConsumesCollector.h"
0014 #include "FWCore/Framework/interface/ESWatcher.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/ProducesCollector.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/Utilities/interface/InputTag.h"
0021
0022 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0023 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0024 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0025 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0026 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0027 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0028
0029 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0030 #include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
0031 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0032 #include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
0033 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0034 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0035 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0036 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0037
0038 #include "SimGeneral/MixingModule/interface/DecayGraph.h"
0039 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixMod.h"
0040 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
0041 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0042 #include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
0043
0044 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0045 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0046 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0047 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0048 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0049
0050 namespace {
0051   using Index_t = unsigned;
0052   using Barcode_t = int;
0053   const std::string messageCategoryGraph_("CaloTruthAccumulatorGraphProducer");
0054 }  // namespace
0055
0056 class CaloTruthAccumulator : public DigiAccumulatorMixMod {
0057 public:
0058   explicit CaloTruthAccumulator(const edm::ParameterSet &config, edm::ProducesCollector, edm::ConsumesCollector &iC);
0059
0060 private:
0061   void initializeEvent(const edm::Event &event, const edm::EventSetup &setup) override;
0062   void accumulate(const edm::Event &event, const edm::EventSetup &setup) override;
0063   void accumulate(const PileUpEventPrincipal &event, const edm::EventSetup &setup, edm::StreamID const &) override;
0064   void finalizeEvent(edm::Event &event, const edm::EventSetup &setup) override;
0065
0066   /** @brief Both forms of accumulate() delegate to this templated method. */
0067   template <class T>
0068   void accumulateEvent(const T &event,
0069                        const edm::EventSetup &setup,
0070                        const edm::Handle<edm::HepMCProduct> &hepMCproduct);
0071
0072   /** @brief Fills the supplied vector with pointers to the SimHits, checking
0073    * for bad modules if required */
0074   template <class T>
0075   void fillSimHits(std::vector<std::pair<DetId, const PCaloHit *>> &returnValue,
0076                    std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0077                    const T &event,
0078                    const edm::EventSetup &setup);
0079
0080   const std::string messageCategory_;
0081
0082   std::unordered_map<Index_t, float> m_detIdToTotalSimEnergy;  // keep track of cell normalizations
0083   std::unordered_multimap<Barcode_t, Index_t> m_simHitBarcodeToIndex;
0084
0085   /** The maximum bunch crossing BEFORE the signal crossing to create
0086       TrackinParticles for. Use positive values. If set to zero no
0087       previous bunches are added and only in-time, signal and after bunches
0088       (defined by maximumSubsequentBunchCrossing_) are used.
0089   */
0090   const unsigned int maximumPreviousBunchCrossing_;
0091   /** The maximum bunch crossing AFTER the signal crossing to create
0092       TrackinParticles for. E.g. if set to zero only
0093       uses the signal and in time pileup (and previous bunches defined by the
0094       maximumPreviousBunchCrossing_ parameter).
0095   */
0096   const unsigned int maximumSubsequentBunchCrossing_;
0097
0098   const edm::InputTag simTrackLabel_;
0099   const edm::InputTag simVertexLabel_;
0100   edm::Handle<std::vector<SimTrack>> hSimTracks;
0101   edm::Handle<std::vector<SimVertex>> hSimVertices;
0102
0103   std::vector<edm::InputTag> collectionTags_;
0104   edm::InputTag genParticleLabel_;
0105   /// Needed to add HepMC::GenVertex to SimVertex
0106   edm::InputTag hepMCproductLabel_;
0107   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0108   edm::ESWatcher<CaloGeometryRecord> geomWatcher_;
0109
0110   const double minEnergy_, maxPseudoRapidity_;
0111   const bool premixStage1_;
0112
0113 public:
0114   struct OutputCollections {
0115     std::unique_ptr<SimClusterCollection> pSimClusters;
0116     std::unique_ptr<CaloParticleCollection> pCaloParticles;
0117   };
0118
0119   struct calo_particles {
0120     std::vector<uint32_t> sc_start_;
0121     std::vector<uint32_t> sc_stop_;
0122
0123     void swap(calo_particles &oth) {
0124       sc_start_.swap(oth.sc_start_);
0125       sc_stop_.swap(oth.sc_stop_);
0126     }
0127
0128     void clear() {
0129       sc_start_.clear();
0130       sc_stop_.clear();
0131     }
0132   };
0133
0134 private:
0135   const HGCalTopology *hgtopo_[3] = {nullptr, nullptr, nullptr};
0136   const HGCalDDDConstants *hgddd_[3] = {nullptr, nullptr, nullptr};
0137   const HcalDDDRecConstants *hcddd_ = nullptr;
0138   OutputCollections output_;
0139   calo_particles m_caloParticles;
0140   // geometry type (0 pre-TDR; 1 TDR)
0141   int geometryType_;
0142   bool doHGCAL;
0143 };
0144
0145 /* Graph utility functions */
0146
0147 namespace {
0148   class CaloParticle_dfs_visitor : public boost::default_dfs_visitor {
0149   public:
0150     CaloParticle_dfs_visitor(CaloTruthAccumulator::OutputCollections &output,
0151                              CaloTruthAccumulator::calo_particles &caloParticles,
0152                              std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex,
0153                              std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0154                              Selector selector)
0155         : output_(output),
0156           caloParticles_(caloParticles),
0157           simHitBarcodeToIndex_(simHitBarcodeToIndex),
0158           simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap),
0159           selector_(selector) {}
0160     template <typename Vertex, typename Graph>
0161     void discover_vertex(Vertex u, const Graph &g) {
0162       // If we reach the vertex 0, it means that we are backtracking with respect
0163       // to the first generation of stable particles: simply return;
0164       //    if (u == 0) return;
0165       print_vertex(u, g);
0166       auto const vertex_property = get(vertex_name, g, u);
0167       if (!vertex_property.simTrack)
0168         return;
0169       auto trackIdx = vertex_property.simTrack->trackId();
0170       IfLogDebug(DEBUG, messageCategoryGraph_)
0171           << " Found " << simHitBarcodeToIndex_.count(trackIdx) << " associated simHits" << std::endl;
0172       if (simHitBarcodeToIndex_.count(trackIdx)) {
0173         output_.pSimClusters->emplace_back(*vertex_property.simTrack);
0174         auto &simcluster = output_.pSimClusters->back();
0175         std::unordered_map<uint32_t, float> acc_energy;
0176         for (auto const &hit_and_energy : simTrackDetIdEnergyMap_[trackIdx]) {
0177           acc_energy[hit_and_energy.first] += hit_and_energy.second;
0178         }
0179         for (auto const &hit_and_energy : acc_energy) {
0180           simcluster.addRecHitAndFraction(hit_and_energy.first, hit_and_energy.second);
0181         }
0182       }
0183     }
0184     template <typename Edge, typename Graph>
0185     void examine_edge(Edge e, const Graph &g) {
0186       auto src = source(e, g);
0187       auto vertex_property = get(vertex_name, g, src);
0188       if (src == 0 or (vertex_property.simTrack == nullptr)) {
0189         auto edge_property = get(edge_weight, g, e);
0190         IfLogDebug(DEBUG, messageCategoryGraph_) << "Considering CaloParticle: " << edge_property.simTrack->trackId();
0191         if (selector_(edge_property)) {
0192           IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
0193           output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
0194           caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
0195         }
0196       }
0197     }
0198
0199     template <typename Edge, typename Graph>
0200     void finish_edge(Edge e, const Graph &g) {
0201       auto src = source(e, g);
0202       auto vertex_property = get(vertex_name, g, src);
0203       if (src == 0 or (vertex_property.simTrack == nullptr)) {
0204         auto edge_property = get(edge_weight, g, e);
0205         if (selector_(edge_property)) {
0206           caloParticles_.sc_stop_.push_back(output_.pSimClusters->size());
0207         }
0208       }
0209     }
0210
0211   private:
0212     CaloTruthAccumulator::OutputCollections &output_;
0213     CaloTruthAccumulator::calo_particles &caloParticles_;
0214     std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex_;
0215     std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap_;
0216     Selector selector_;
0217   };
0218 }  // namespace
0219
0220 CaloTruthAccumulator::CaloTruthAccumulator(const edm::ParameterSet &config,
0221                                            edm::ProducesCollector producesCollector,
0222                                            edm::ConsumesCollector &iC)
0223     : messageCategory_("CaloTruthAccumulator"),
0224       maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
0225       maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
0226       simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
0227       simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
0228       collectionTags_(),
0229       genParticleLabel_(config.getParameter<edm::InputTag>("genParticleCollection")),
0230       hepMCproductLabel_(config.getParameter<edm::InputTag>("HepMCProductLabel")),
0231       geomToken_(iC.esConsumes()),
0232       minEnergy_(config.getParameter<double>("MinEnergy")),
0233       maxPseudoRapidity_(config.getParameter<double>("MaxPseudoRapidity")),
0234       premixStage1_(config.getParameter<bool>("premixStage1")),
0235       geometryType_(-1),
0236       doHGCAL(config.getParameter<bool>("doHGCAL")) {
0237   producesCollector.produces<SimClusterCollection>("MergedCaloTruth");
0238   producesCollector.produces<CaloParticleCollection>("MergedCaloTruth");
0239   if (premixStage1_) {
0240     producesCollector.produces<std::vector<std::pair<unsigned int, float>>>("MergedCaloTruth");
0241   }
0242
0243   iC.consumes<std::vector<SimTrack>>(simTrackLabel_);
0244   iC.consumes<std::vector<SimVertex>>(simVertexLabel_);
0245   iC.consumes<std::vector<reco::GenParticle>>(genParticleLabel_);
0246   iC.consumes<std::vector<int>>(genParticleLabel_);
0247   iC.consumes<std::vector<int>>(hepMCproductLabel_);
0248
0249   // Fill the collection tags
0250   const edm::ParameterSet &simHitCollectionConfig = config.getParameterSet("simHitCollections");
0251   std::vector<std::string> parameterNames = simHitCollectionConfig.getParameterNames();
0252
0253   for (auto const &parameterName : parameterNames) {
0254     std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
0255     collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
0256   }
0257
0258   for (auto const &collectionTag : collectionTags_) {
0259     iC.consumes<std::vector<PCaloHit>>(collectionTag);
0260   }
0261 }
0262
0263 void CaloTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
0264   output_.pSimClusters = std::make_unique<SimClusterCollection>();
0265   output_.pCaloParticles = std::make_unique<CaloParticleCollection>();
0266
0267   m_detIdToTotalSimEnergy.clear();
0268
0269   if (geomWatcher_.check(setup)) {
0270     auto const &geom = setup.getData(geomToken_);
0271     const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
0272     const HcalGeometry *bhgeom = nullptr;
0273     bhgeom = static_cast<const HcalGeometry *>(geom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0274
0275     if (doHGCAL) {
0276       eegeom = static_cast<const HGCalGeometry *>(
0277           geom.getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
0278       // check if it's the new geometry
0279       if (eegeom) {
0280         geometryType_ = 1;
0281         fhgeom = static_cast<const HGCalGeometry *>(
0282             geom.getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
0283         bhgeomnew = static_cast<const HGCalGeometry *>(
0284             geom.getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0285       } else {
0286         geometryType_ = 0;
0287         eegeom = static_cast<const HGCalGeometry *>(geom.getSubdetectorGeometry(DetId::Forward, HGCEE));
0288         fhgeom = static_cast<const HGCalGeometry *>(geom.getSubdetectorGeometry(DetId::Forward, HGCHEF));
0289         bhgeom = static_cast<const HcalGeometry *>(geom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0290       }
0291       hgtopo_[0] = &(eegeom->topology());
0292       hgtopo_[1] = &(fhgeom->topology());
0293       if (bhgeomnew)
0294         hgtopo_[2] = &(bhgeomnew->topology());
0295
0296       for (unsigned i = 0; i < 3; ++i) {
0297         if (hgtopo_[i])
0298           hgddd_[i] = &(hgtopo_[i]->dddConstants());
0299       }
0300     }
0301
0302     if (bhgeom) {
0303       hcddd_ = bhgeom->topology().dddConstants();
0304     }
0305   }
0306 }
0307
0308 /** Create handle to edm::HepMCProduct here because event.getByLabel with
0309     edm::HepMCProduct only works for edm::Event but not for
0310     PileUpEventPrincipal; PileUpEventPrincipal::getByLabel tries to call
0311     T::value_type and T::iterator (where T is the type of the object one wants
0312     to get a handle to) which is only implemented for container-like objects
0313     like std::vector but not for edm::HepMCProduct!
0314 */
0315 void CaloTruthAccumulator::accumulate(edm::Event const &event, edm::EventSetup const &setup) {
0316   edm::Handle<edm::HepMCProduct> hepmc;
0317   event.getByLabel(hepMCproductLabel_, hepmc);
0318
0319   edm::LogInfo(messageCategory_) << " CaloTruthAccumulator::accumulate (signal)";
0320   accumulateEvent(event, setup, hepmc);
0321 }
0322
0323 void CaloTruthAccumulator::accumulate(PileUpEventPrincipal const &event,
0324                                       edm::EventSetup const &setup,
0325                                       edm::StreamID const &) {
0326   if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
0327       event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
0328     // simply create empty handle as we do not have a HepMCProduct in PU anyway
0329     edm::Handle<edm::HepMCProduct> hepmc;
0330     edm::LogInfo(messageCategory_) << " CaloTruthAccumulator::accumulate (pileup) bunchCrossing="
0331                                    << event.bunchCrossing();
0332     accumulateEvent(event, setup, hepmc);
0333   } else {
0334     edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
0335   }
0336 }
0337
0338 void CaloTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const &setup) {
0339   edm::LogInfo(messageCategory_) << "Adding " << output_.pSimClusters->size() << " SimParticles and "
0340                                  << output_.pCaloParticles->size() << " CaloParticles to the event.";
0341
0342   // We need to normalize the hits and energies into hits and fractions (since
0343   // we have looped over all pileup events)
0344   // For premixing stage1 we keep the energies, they will be normalized to
0345   // fractions in stage2
0346
0347   if (premixStage1_) {
0348     auto totalEnergies = std::make_unique<std::vector<std::pair<unsigned int, float>>>();
0349     totalEnergies->reserve(m_detIdToTotalSimEnergy.size());
0350     std::copy(m_detIdToTotalSimEnergy.begin(), m_detIdToTotalSimEnergy.end(), std::back_inserter(*totalEnergies));
0351     std::sort(totalEnergies->begin(), totalEnergies->end());
0352     event.put(std::move(totalEnergies), "MergedCaloTruth");
0353   } else {
0354     for (auto &sc : *(output_.pSimClusters)) {
0355       auto hitsAndEnergies = sc.hits_and_fractions();
0356       sc.clearHitsAndFractions();
0357       sc.clearHitsEnergy();
0358       for (auto &hAndE : hitsAndEnergies) {
0359         const float totalenergy = m_detIdToTotalSimEnergy[hAndE.first];
0360         float fraction = 0.;
0361         if (totalenergy > 0)
0362           fraction = hAndE.second / totalenergy;
0363         else
0364           edm::LogWarning(messageCategory_)
0365               << "TotalSimEnergy for hit " << hAndE.first << " is 0! The fraction for this hit cannot be computed.";
0366         sc.addRecHitAndFraction(hAndE.first, fraction);
0367         sc.addHitEnergy(hAndE.second);
0368       }
0369     }
0370   }
0371
0372   // save the SimCluster orphan handle so we can fill the calo particles
0373   auto scHandle = event.put(std::move(output_.pSimClusters), "MergedCaloTruth");
0374
0375   // now fill the calo particles
0376   for (unsigned i = 0; i < output_.pCaloParticles->size(); ++i) {
0377     auto &cp = (*output_.pCaloParticles)[i];
0378     for (unsigned j = m_caloParticles.sc_start_[i]; j < m_caloParticles.sc_stop_[i]; ++j) {
0379       edm::Ref<SimClusterCollection> ref(scHandle, j);
0380       cp.addSimCluster(ref);
0381     }
0382   }
0383
0384   event.put(std::move(output_.pCaloParticles), "MergedCaloTruth");
0385
0386   calo_particles().swap(m_caloParticles);
0387
0388   std::unordered_map<Index_t, float>().swap(m_detIdToTotalSimEnergy);
0389   std::unordered_multimap<Barcode_t, Index_t>().swap(m_simHitBarcodeToIndex);
0390 }
0391
0392 template <class T>
0393 void CaloTruthAccumulator::accumulateEvent(const T &event,
0394                                            const edm::EventSetup &setup,
0395                                            const edm::Handle<edm::HepMCProduct> &hepMCproduct) {
0396   edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
0397   edm::Handle<std::vector<int>> hGenParticleIndices;
0398
0399   event.getByLabel(simTrackLabel_, hSimTracks);
0400   event.getByLabel(simVertexLabel_, hSimVertices);
0401
0402   event.getByLabel(genParticleLabel_, hGenParticles);
0403   event.getByLabel(genParticleLabel_, hGenParticleIndices);
0404
0405   std::vector<std::pair<DetId, const PCaloHit *>> simHitPointers;
0406   std::unordered_map<int, std::map<int, float>> simTrackDetIdEnergyMap;
0407   fillSimHits(simHitPointers, simTrackDetIdEnergyMap, event, setup);
0408
0409   // Clear maps from previous event fill them for this one
0410   m_simHitBarcodeToIndex.clear();
0411   for (unsigned int i = 0; i < simHitPointers.size(); ++i) {
0412     m_simHitBarcodeToIndex.emplace(simHitPointers[i].second->geantTrackId(), i);
0413   }
0414
0415   auto const &tracks = *hSimTracks;
0416   auto const &vertices = *hSimVertices;
0417   std::unordered_map<int, int> trackid_to_track_index;
0418   DecayChain decay;
0419   int idx = 0;
0420
0421   IfLogDebug(DEBUG, messageCategory_) << " TRACKS" << std::endl;
0422   for (auto const &t : tracks) {
0423     IfLogDebug(DEBUG, messageCategory_) << " " << idx << "\t" << t.trackId() << "\t" << t << std::endl;
0424     trackid_to_track_index[t.trackId()] = idx;
0425     idx++;
0426   }
0427
0428   /**
0429   Build the main decay graph and assign the SimTrack to each edge. The graph
0430   built here will only contain the particles that have a decay vertex
0431   associated to them. In order to recover also the particles that will not
0432   decay, we need to keep track of the SimTrack used here and add, a-posteriori,
0433   the ones not used, associating a ghost vertex (starting from the highest
0434   simulated vertex number), in order to build the edge and identify them
0435   immediately as stable (i.e. not decayed).
0436
0437   To take into account the multi-bremsstrahlung effects in which a single
0438   particle is emitting photons in different vertices **keeping the same
0439   track index**, we also collapsed those vertices into 1 unique vertex. The
0440   other approach of fully representing the decay chain keeping the same
0441   track index would have the problem of over-counting the contributions of
0442   that track, especially in terms of hits.
0443
0444   The 2 auxiliary vectors are structured as follow:
0445
0446   1. used_sim_tracks is a vector that has the same size as the overall
0447      number of simulated tracks. The associated integer is the vertexId of
0448      the **decaying vertex for that track**.
0449   2. collapsed_vertices is a vector that has the same size as the overall
0450      number of simulated vertices. The vector's index is the vertexId
0451      itself, the associated value is the vertexId of the vertex on which
0452      this should collapse.
0453   */
0454   idx = 0;
0455   std::vector<int> used_sim_tracks(tracks.size(), 0);
0456   std::vector<int> collapsed_vertices(vertices.size(), 0);
0457   IfLogDebug(DEBUG, messageCategory_) << " VERTICES" << std::endl;
0458   for (auto const &v : vertices) {
0459     IfLogDebug(DEBUG, messageCategory_) << " " << idx++ << "\t" << v << std::endl;
0460     if (v.parentIndex() != -1) {
0461       auto trk_idx = trackid_to_track_index[v.parentIndex()];
0462       auto origin_vtx = tracks[trk_idx].vertIndex();
0463       if (used_sim_tracks[trk_idx]) {
0464         // collapse the vertex into the original first vertex we saw associated
0465         // to this track. Omit adding the edge in order to avoid double
0466         // counting of the very same particles  and its associated hits.
0467         collapsed_vertices[v.vertexId()] = used_sim_tracks[trk_idx];
0468         continue;
0469       }
0470       // Perform the actual vertex collapsing, if needed.
0471       if (collapsed_vertices[origin_vtx])
0472         origin_vtx = collapsed_vertices[origin_vtx];
0473       add_edge(origin_vtx,
0474                v.vertexId(),
0475                EdgeProperty(&tracks[trk_idx], simTrackDetIdEnergyMap[v.parentIndex()].size(), 0),
0476                decay);
0477       used_sim_tracks[trk_idx] = v.vertexId();
0478     }
0479   }
0480   // Build the motherParticle property to each vertex
0481   auto const &vertexMothersProp = get(vertex_name, decay);
0482   // Now recover the particles that did not decay. Append them with an index
0483   // bigger than the size of the generated vertices.
0484   int offset = vertices.size();
0485   for (size_t i = 0; i < tracks.size(); ++i) {
0486     if (!used_sim_tracks[i]) {
0487       auto origin_vtx = tracks[i].vertIndex();
0488       // Perform the actual vertex collapsing, if needed.
0489       if (collapsed_vertices[origin_vtx])
0490         origin_vtx = collapsed_vertices[origin_vtx];
0491       add_edge(
0492           origin_vtx, offset, EdgeProperty(&tracks[i], simTrackDetIdEnergyMap[tracks[i].trackId()].size(), 0), decay);
0493       // The properties for "fake" vertices associated to stable particles have
0494       // to be set inside this loop, since they do not belong to the vertices
0495       // collection and would be skipped by that loop (coming next)
0496       put(vertexMothersProp, offset, VertexProperty(&tracks[i], 0));
0497       offset++;
0498     }
0499   }
0500   for (auto const &v : vertices) {
0501     if (v.parentIndex() != -1) {
0502       // Skip collapsed_vertices
0503       if (collapsed_vertices[v.vertexId()])
0504         continue;
0505       put(vertexMothersProp, v.vertexId(), VertexProperty(&tracks[trackid_to_track_index[v.parentIndex()]], 0));
0506     }
0507   }
0508   SimHitsAccumulator_dfs_visitor vis;
0509   depth_first_search(decay, visitor(vis));
0510   CaloParticle_dfs_visitor caloParticleCreator(
0511       output_,
0512       m_caloParticles,
0513       m_simHitBarcodeToIndex,
0514       simTrackDetIdEnergyMap,
0515       [&](EdgeProperty &edge_property) -> bool {
0516         // Apply selection on SimTracks in order to promote them to be
0517         // CaloParticles. The function returns TRUE if the particle satisfies
0518         // the selection, FALSE otherwise. Therefore the correct logic to select
0519         // the particle is to ask for TRUE as return value.
0520         return (edge_property.cumulative_simHits != 0 and !edge_property.simTrack->noGenpart() and
0521                 edge_property.simTrack->momentum().E() > minEnergy_ and
0522                 std::abs(edge_property.simTrack->momentum().Eta()) < maxPseudoRapidity_);
0523       });
0524   depth_first_search(decay, visitor(caloParticleCreator));
0525
0526 #if DEBUG
0527   boost::write_graphviz(std::cout,
0528                         decay,
0529                         make_label_writer(make_transform_value_property_map(&graphviz_vertex, get(vertex_name, decay))),
0530                         make_label_writer(make_transform_value_property_map(&graphviz_edge, get(edge_weight, decay))));
0531 #endif
0532 }
0533
0534 template <class T>
0535 void CaloTruthAccumulator::fillSimHits(std::vector<std::pair<DetId, const PCaloHit *>> &returnValue,
0536                                        std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0537                                        const T &event,
0538                                        const edm::EventSetup &setup) {
0539   for (auto const &collectionTag : collectionTags_) {
0540     edm::Handle<std::vector<PCaloHit>> hSimHits;
0541     const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
0542     event.getByLabel(collectionTag, hSimHits);
0543
0544     for (auto const &simHit : *hSimHits) {
0545       DetId id(0);
0546
0547       //Relabel as necessary for HGCAL
0548       if (doHGCAL) {
0549         const uint32_t simId = simHit.id();
0550         if (geometryType_ == 1) {
0551           // no test numbering in new geometry
0552           id = simId;
0553         } else if (isHcal) {
0554           HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd_);
0555           if (hid.subdet() == HcalEndcap)
0556             id = hid;
0557         } else {
0558           int subdet, layer, cell, sec, subsec, zp;
0559           HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
0560           const HGCalDDDConstants *ddd = hgddd_[subdet - 3];
0561           std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo_[subdet - 3]->detectorType());
0562           cell = recoLayerCell.first;
0563           layer = recoLayerCell.second;
0564           // skip simhits with bad barcodes or non-existant layers
0565           if (layer == -1 || simHit.geantTrackId() == 0)
0566             continue;
0567           id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
0568         }
0569       } else {
0570         id = simHit.id();
0571         //Relabel all HCAL hits
0572         if (isHcal) {
0573           HcalDetId hid = HcalHitRelabeller::relabel(simHit.id(), hcddd_);
0574           id = hid;
0575         }
0576       }
0577
0578       if (id == DetId(0)) {
0579         continue;
0580       }
0581       if (simHit.geantTrackId() == 0) {
0582         continue;
0583       }
0584
0585       returnValue.emplace_back(id, &simHit);
0586       simTrackDetIdEnergyMap[simHit.geantTrackId()][id.rawId()] += simHit.energy();
0587       m_detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
0588     }
0589   }  // end of loop over InputTags
0590 }
0591
0592 // Register with the framework
0593 DEFINE_DIGI_ACCUMULATOR(CaloTruthAccumulator);