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