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