Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:04:38

0001 // Author: Aurora Perego, Fabio Cossutti - aurora.perego@cern.ch, fabio.cossutti@ts.infn.it
0002 // Date: 05/2023
0003 
0004 #define DEBUG false
0005 
0006 #if DEBUG
0007 #pragma GCC diagnostic pop
0008 #endif
0009 
0010 #include <algorithm>
0011 #include <iterator>
0012 #include <exception>
0013 #include <memory>
0014 
0015 #include <numeric>  // for std::accumulate
0016 #include <unordered_map>
0017 
0018 #include "FWCore/Framework/interface/ConsumesCollector.h"
0019 #include "FWCore/Framework/interface/ESWatcher.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/ProducesCollector.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/InputTag.h"
0026 
0027 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0028 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0029 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0030 #include "DataFormats/Math/interface/GeantUnits.h"
0031 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0032 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0033 
0034 #include "SimDataFormats/CaloAnalysis/interface/MtdCaloParticle.h"
0035 #include "SimDataFormats/CaloAnalysis/interface/MtdCaloParticleFwd.h"
0036 #include "SimDataFormats/CaloAnalysis/interface/MtdSimCluster.h"
0037 #include "SimDataFormats/CaloAnalysis/interface/MtdSimClusterFwd.h"
0038 #include "SimDataFormats/CaloAnalysis/interface/MtdSimLayerCluster.h"
0039 #include "SimDataFormats/CaloAnalysis/interface/MtdSimLayerClusterFwd.h"
0040 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTrackster.h"
0041 #include "SimDataFormats/CaloAnalysis/interface/MtdSimTracksterFwd.h"
0042 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0043 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0044 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0045 
0046 #include "SimGeneral/MixingModule/interface/DecayGraph.h"
0047 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixMod.h"
0048 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
0049 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0050 #include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
0051 
0052 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0053 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0054 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0055 #include "Geometry/MTDGeometryBuilder/interface/MTDGeomUtil.h"
0056 #include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
0057 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0058 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0059 
0060 namespace {
0061   using Index_t = unsigned;
0062   using Barcode_t = int;
0063   const std::string messageCategoryGraph_("MtdTruthAccumulatorGraphProducer");
0064 }  // namespace
0065 
0066 class MtdTruthAccumulator : public DigiAccumulatorMixMod {
0067 public:
0068   explicit MtdTruthAccumulator(const edm::ParameterSet &config, edm::ProducesCollector, edm::ConsumesCollector &iC);
0069 
0070 private:
0071   void initializeEvent(const edm::Event &event, const edm::EventSetup &setup) override;
0072   void accumulate(const edm::Event &event, const edm::EventSetup &setup) override;
0073   void accumulate(const PileUpEventPrincipal &event, const edm::EventSetup &setup, edm::StreamID const &) override;
0074   void finalizeEvent(edm::Event &event, const edm::EventSetup &setup) override;
0075 
0076   /** @brief Both forms of accumulate() delegate to this templated method. */
0077   template <class T>
0078   void accumulateEvent(const T &event,
0079                        const edm::EventSetup &setup,
0080                        const edm::Handle<edm::HepMCProduct> &hepMCproduct);
0081 
0082   /** @brief Fills the supplied vector with pointers to the SimHits, checking
0083    * for bad modules if required */
0084   template <class T>
0085   void fillSimHits(std::vector<std::pair<uint64_t, const PSimHit *>> &returnValue,
0086                    std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdEnergyMap,
0087                    std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdTimeMap,
0088                    const T &event,
0089                    const edm::EventSetup &setup);
0090 
0091   const std::string messageCategory_;
0092 
0093   std::unordered_map<uint64_t, float> m_detIdToTotalSimEnergy;  // keep track of cell normalizations
0094   std::unordered_multimap<Barcode_t, Index_t> m_simHitBarcodeToIndex;
0095 
0096   /** The maximum bunch crossing BEFORE the signal crossing to create
0097       TrackinParticles for. Use positive values. If set to zero no
0098       previous bunches are added and only in-time, signal and after bunches
0099       (defined by maximumSubsequentBunchCrossing_) are used.
0100   */
0101   const unsigned int maximumPreviousBunchCrossing_;
0102   /** The maximum bunch crossing AFTER the signal crossing to create
0103       TrackinParticles for. E.g. if set to zero only
0104       uses the signal and in time pileup (and previous bunches defined by the
0105       maximumPreviousBunchCrossing_ parameter).
0106   */
0107   const unsigned int maximumSubsequentBunchCrossing_;
0108 
0109   const unsigned int bunchSpacing_;
0110 
0111   const edm::InputTag simTrackLabel_;
0112   const edm::InputTag simVertexLabel_;
0113   edm::Handle<std::vector<SimTrack>> hSimTracks;
0114   edm::Handle<std::vector<SimVertex>> hSimVertices;
0115 
0116   std::vector<edm::InputTag> collectionTags_;
0117   edm::InputTag genParticleLabel_;
0118   /// Needed to add HepMC::GenVertex to SimVertex
0119   edm::InputTag hepMCproductLabel_;
0120   const edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> geomToken_;
0121   const edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0122   // edm::ESWatcher<MTDDigiGeometryRecord> geomWatcher_;
0123 
0124   mtd::MTDGeomUtil geomTools_;
0125 
0126   const double minEnergy_, maxPseudoRapidity_;
0127   const bool premixStage1_;
0128 
0129   bool isEtl_;
0130 
0131 public:
0132   struct OutputCollections {
0133     std::unique_ptr<MtdSimClusterCollection> pSimClusters;
0134     std::unique_ptr<MtdCaloParticleCollection> pCaloParticles;
0135     std::unique_ptr<MtdSimLayerClusterCollection> pMtdSimLayerClusters;
0136     std::unique_ptr<MtdSimTracksterCollection> pMtdSimTracksters;
0137   };
0138 
0139   struct calo_particles {
0140     std::vector<uint32_t> sc_start_;
0141     std::vector<uint32_t> sc_stop_;
0142 
0143     void swap(calo_particles &oth) {
0144       sc_start_.swap(oth.sc_start_);
0145       sc_stop_.swap(oth.sc_stop_);
0146     }
0147 
0148     void clear() {
0149       sc_start_.clear();
0150       sc_stop_.clear();
0151     }
0152   };
0153 
0154 private:
0155   const MTDGeometry *geom = nullptr;
0156   const MTDTopology *topology = nullptr;
0157   OutputCollections output_;
0158   calo_particles m_caloParticles;
0159 };
0160 
0161 /* Graph utility functions */
0162 
0163 namespace {
0164   class CaloParticle_dfs_visitor : public boost::default_dfs_visitor {
0165   public:
0166     CaloParticle_dfs_visitor(MtdTruthAccumulator::OutputCollections &output,
0167                              MtdTruthAccumulator::calo_particles &caloParticles,
0168                              std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex,
0169                              std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdEnergyMap,
0170                              std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdTimeMap,
0171                              std::unordered_map<uint32_t, float> &vertex_time_map,
0172                              Selector selector)
0173         : output_(output),
0174           caloParticles_(caloParticles),
0175           simHitBarcodeToIndex_(simHitBarcodeToIndex),
0176           simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap),
0177           simTrackDetIdTimeMap_(simTrackDetIdTimeMap),
0178           vertex_time_map_(vertex_time_map),
0179           selector_(selector) {}
0180     template <typename Vertex, typename Graph>
0181     void discover_vertex(Vertex u, const Graph &g) {
0182       // If we reach the vertex 0, it means that we are backtracking with respect
0183       // to the first generation of stable particles: simply return;
0184       //    if (u == 0) return;
0185       print_vertex(u, g);
0186       auto const vertex_property = get(vertex_name, g, u);
0187       if (!vertex_property.simTrack)
0188         return;
0189       auto trackIdx = vertex_property.simTrack->trackId();
0190       IfLogDebug(DEBUG, messageCategoryGraph_)
0191           << " Found " << simHitBarcodeToIndex_.count(trackIdx) << " associated simHits" << std::endl;
0192       if (simHitBarcodeToIndex_.count(trackIdx)) {
0193         output_.pSimClusters->emplace_back(*vertex_property.simTrack);
0194         auto &simcluster = output_.pSimClusters->back();
0195         std::unordered_map<uint64_t, float> acc_energy;
0196         for (auto const &hit_and_energy : simTrackDetIdEnergyMap_[trackIdx]) {
0197           acc_energy[hit_and_energy.first] += hit_and_energy.second;
0198         }
0199         for (auto const &hit_and_energy : acc_energy) {
0200           simcluster.addHitAndFraction(hit_and_energy.first, hit_and_energy.second);
0201           simcluster.addHitEnergy(hit_and_energy.second);
0202           simcluster.addHitTime(simTrackDetIdTimeMap_[simcluster.g4Tracks()[0].trackId()][hit_and_energy.first]);
0203         }
0204       }
0205     }
0206     template <typename Edge, typename Graph>
0207     void examine_edge(Edge e, const Graph &g) {
0208       auto src = source(e, g);
0209       auto vertex_property = get(vertex_name, g, src);
0210       if (src == 0 or (vertex_property.simTrack == nullptr)) {
0211         auto edge_property = get(edge_weight, g, e);
0212         IfLogDebug(DEBUG, messageCategoryGraph_) << "Considering CaloParticle: " << edge_property.simTrack->trackId();
0213         if (selector_(edge_property)) {
0214           IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
0215           output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
0216           output_.pCaloParticles->back().addSimTime(vertex_time_map_[(edge_property.simTrack)->vertIndex()]);
0217           caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
0218         }
0219       }
0220     }
0221 
0222     template <typename Edge, typename Graph>
0223     void finish_edge(Edge e, const Graph &g) {
0224       auto src = source(e, g);
0225       auto vertex_property = get(vertex_name, g, src);
0226       if (src == 0 or (vertex_property.simTrack == nullptr)) {
0227         auto edge_property = get(edge_weight, g, e);
0228         if (selector_(edge_property)) {
0229           caloParticles_.sc_stop_.push_back(output_.pSimClusters->size());
0230         }
0231       }
0232     }
0233 
0234   private:
0235     MtdTruthAccumulator::OutputCollections &output_;
0236     MtdTruthAccumulator::calo_particles &caloParticles_;
0237     std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex_;
0238     std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdEnergyMap_;
0239     std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdTimeMap_;
0240     std::unordered_map<uint32_t, float> &vertex_time_map_;
0241     Selector selector_;
0242   };
0243 }  // namespace
0244 
0245 MtdTruthAccumulator::MtdTruthAccumulator(const edm::ParameterSet &config,
0246                                          edm::ProducesCollector producesCollector,
0247                                          edm::ConsumesCollector &iC)
0248     : messageCategory_("MtdTruthAccumulator"),
0249       maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
0250       maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
0251       bunchSpacing_(config.getParameter<unsigned int>("bunchspace")),
0252       simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
0253       simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
0254       collectionTags_(),
0255       genParticleLabel_(config.getParameter<edm::InputTag>("genParticleCollection")),
0256       hepMCproductLabel_(config.getParameter<edm::InputTag>("HepMCProductLabel")),
0257       geomToken_(iC.esConsumes<MTDGeometry, MTDDigiGeometryRecord>()),
0258       mtdtopoToken_(iC.esConsumes<MTDTopology, MTDTopologyRcd>()),
0259       minEnergy_(config.getParameter<double>("MinEnergy")),
0260       maxPseudoRapidity_(config.getParameter<double>("MaxPseudoRapidity")),
0261       premixStage1_(config.getParameter<bool>("premixStage1")) {
0262   iC.consumes<std::vector<SimTrack>>(simTrackLabel_);
0263   iC.consumes<std::vector<SimVertex>>(simVertexLabel_);
0264   iC.consumes<std::vector<reco::GenParticle>>(genParticleLabel_);
0265   iC.consumes<std::vector<int>>(genParticleLabel_);
0266   iC.consumes<std::vector<int>>(hepMCproductLabel_);
0267 
0268   // Fill the collection tags
0269   const edm::ParameterSet &simHitCollectionConfig = config.getParameterSet("simHitCollections");
0270   std::vector<std::string> parameterNames = simHitCollectionConfig.getParameterNames();
0271 
0272   for (auto const &parameterName : parameterNames) {
0273     std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
0274     collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
0275   }
0276 
0277   for (auto const &collectionTag : collectionTags_) {
0278     iC.consumes<std::vector<PSimHit>>(collectionTag);
0279     isEtl_ = (collectionTag.instance().find("FastTimerHitsEndcap") != std::string::npos);
0280   }
0281 
0282   producesCollector.produces<MtdSimClusterCollection>("MergedMtdTruth");
0283   producesCollector.produces<MtdSimLayerClusterCollection>("MergedMtdTruthLC");
0284   producesCollector.produces<MtdSimTracksterCollection>("MergedMtdTruthST");
0285   producesCollector.produces<MtdCaloParticleCollection>("MergedMtdTruth");
0286   if (premixStage1_) {
0287     producesCollector.produces<std::vector<std::pair<uint64_t, float>>>("MergedMtdTruth");
0288   }
0289 }
0290 
0291 void MtdTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
0292   output_.pSimClusters = std::make_unique<MtdSimClusterCollection>();
0293   output_.pCaloParticles = std::make_unique<MtdCaloParticleCollection>();
0294 
0295   output_.pMtdSimLayerClusters = std::make_unique<MtdSimLayerClusterCollection>();
0296   output_.pMtdSimTracksters = std::make_unique<MtdSimTracksterCollection>();
0297 
0298   m_detIdToTotalSimEnergy.clear();
0299 
0300   auto geometryHandle = setup.getTransientHandle(geomToken_);
0301   geom = geometryHandle.product();
0302 
0303   auto topologyHandle = setup.getTransientHandle(mtdtopoToken_);
0304   topology = topologyHandle.product();
0305 
0306   geomTools_.setGeometry(geom);
0307   geomTools_.setTopology(topology);
0308 }
0309 
0310 /** Create handle to edm::HepMCProduct here because event.getByLabel with
0311     edm::HepMCProduct only works for edm::Event but not for
0312     PileUpEventPrincipal; PileUpEventPrincipal::getByLabel tries to call
0313     T::value_type and T::iterator (where T is the type of the object one wants
0314     to get a handle to) which is only implemented for container-like objects
0315     like std::vector but not for edm::HepMCProduct!
0316 */
0317 void MtdTruthAccumulator::accumulate(edm::Event const &event, edm::EventSetup const &setup) {
0318   edm::Handle<edm::HepMCProduct> hepmc;
0319   event.getByLabel(hepMCproductLabel_, hepmc);
0320 
0321   edm::LogInfo(messageCategory_) << " MtdTruthAccumulator::accumulate (signal)";
0322   accumulateEvent(event, setup, hepmc);
0323 }
0324 
0325 void MtdTruthAccumulator::accumulate(PileUpEventPrincipal const &event,
0326                                      edm::EventSetup const &setup,
0327                                      edm::StreamID const &) {
0328   if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
0329       event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
0330     // simply create empty handle as we do not have a HepMCProduct in PU anyway
0331     edm::Handle<edm::HepMCProduct> hepmc;
0332     edm::LogInfo(messageCategory_) << " MtdTruthAccumulator::accumulate (pileup) bunchCrossing="
0333                                    << event.bunchCrossing();
0334     accumulateEvent(event, setup, hepmc);
0335   } else {
0336     edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
0337   }
0338 }
0339 
0340 void MtdTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const &setup) {
0341   using namespace geant_units::operators;
0342 
0343   edm::LogInfo(messageCategory_) << "Adding " << output_.pSimClusters->size() << " SimParticles and "
0344                                  << output_.pCaloParticles->size() << " CaloParticles to the event.";
0345 
0346   // We need to normalize the hits and energies into hits and fractions (since
0347   // we have looped over all pileup events)
0348   // For premixing stage1 we keep the energies, they will be normalized to
0349   // fractions in stage2
0350 
0351   if (premixStage1_) {
0352     auto totalEnergies = std::make_unique<std::vector<std::pair<uint64_t, float>>>();
0353     totalEnergies->reserve(m_detIdToTotalSimEnergy.size());
0354     std::copy(m_detIdToTotalSimEnergy.begin(), m_detIdToTotalSimEnergy.end(), std::back_inserter(*totalEnergies));
0355     std::sort(totalEnergies->begin(), totalEnergies->end());
0356     event.put(std::move(totalEnergies), "MergedMtdTruth");
0357   } else {
0358     for (auto &sc : *(output_.pSimClusters)) {
0359       auto hitsAndEnergies = sc.hits_and_fractions();
0360       sc.clearHitsAndFractions();
0361       sc.clearHitsEnergy();
0362       for (auto &hAndE : hitsAndEnergies) {
0363         const float totalenergy = m_detIdToTotalSimEnergy[hAndE.first];
0364         float fraction = 0.;
0365         if (totalenergy > 0)
0366           fraction = hAndE.second / totalenergy;
0367         else
0368           edm::LogWarning(messageCategory_)
0369               << "TotalSimEnergy for hit " << hAndE.first << " is 0! The fraction for this hit cannot be computed.";
0370         sc.addHitAndFraction(hAndE.first, fraction);
0371         sc.addHitEnergy(hAndE.second);
0372       }
0373     }
0374   }
0375 
0376 #ifdef PRINT_DEBUG
0377   IfLogDebug(DEBUG, messageCategory_) << "SIMCLUSTERS LIST:" << std::endl;
0378   for (const auto &sc : *(output_.pSimClusters)) {
0379     IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "SimCluster from CP with:"
0380                                         << "\n  charge " << sc.charge() << "\n  pdgId  " << sc.pdgId() << "\n  energy "
0381                                         << sc.energy() << " GeV\n  eta    " << sc.eta() << "\n  phi    " << sc.phi()
0382                                         << "\n  number of cells = " << sc.hits_and_fractions().size() << std::endl;
0383     for (unsigned int i = 0; i < sc.hits_and_fractions().size(); ++i) {
0384       DetId id(sc.detIds_and_rows()[i].first);
0385       IfLogDebug(DEBUG, messageCategory_)
0386           << std::fixed << std::setprecision(3) << " hit " << id.rawId() << " on layer " << geomTools_.layer(id)
0387           << " module " << geomTools_.module(id) << " row " << (unsigned int)(sc.detIds_and_rows()[i].second).first
0388           << " col " << (unsigned int)(sc.detIds_and_rows()[i].second).second << " at time "
0389           << sc.hits_and_times()[i].second << " ns" << std::endl;
0390     }
0391     IfLogDebug(DEBUG, messageCategory_) << "--------------\n";
0392   }
0393   IfLogDebug(DEBUG, messageCategory_) << std::endl;
0394 #endif
0395 
0396   // save the SimCluster orphan handle so we can fill the calo particles
0397   auto scHandle = event.put(std::move(output_.pSimClusters), "MergedMtdTruth");
0398 
0399   // reserve for the best case scenario: already splitted
0400   output_.pMtdSimLayerClusters->reserve(scHandle->size());
0401   output_.pMtdSimTracksters->reserve(scHandle->size());
0402 
0403   uint32_t SC_index = 0;
0404   uint32_t LC_index = 0;
0405   for (const auto &sc : *scHandle) {
0406     auto const &hAndF = sc.hits_and_fractions();
0407     auto const &hAndE = sc.hits_and_energies();
0408     auto const &hAndT = sc.hits_and_times();
0409     auto const &hAndR = sc.detIds_and_rows();
0410     // create a vector with the indices of the hits in the simCluster
0411     std::vector<int> indices(hAndF.size());
0412     std::iota(indices.begin(), indices.end(), 0);
0413     // sort the hits indices based on the unique indices created before
0414     std::sort(indices.begin(), indices.end(), [&](int a, int b) { return hAndF[a].first < hAndF[b].first; });
0415 
0416     // now split the sc: loop on the sorted indices and save the first hit in a
0417     // temporary simCluster. If the following hit is in the same module and row (column),
0418     // but next column (row) put it in the temporary simcluster as well, otherwise
0419     // put the temporary simcluster in the collection and start creating a new one
0420     std::vector<uint32_t> LC_indices;
0421     MtdSimLayerCluster tmpLC(sc.g4Tracks()[0]);
0422     int prev = indices[0];
0423     DetId prevId(hAndR[prev].first);
0424 
0425     float SimLCenergy = 0.;
0426     float SimLCx = 0., SimLCy = 0., SimLCz = 0.;
0427 
0428     auto push_back_hit = [&](const int &ind) {
0429       tmpLC.addHitAndFraction(hAndF[ind].first, hAndF[ind].second);
0430       tmpLC.addHitEnergy(hAndE[ind].second);
0431       tmpLC.addHitTime(hAndT[ind].second);
0432     };
0433 
0434     auto update_clu_info = [&](const int &ind) {
0435       double energy = hAndE[ind].second;
0436       auto position =
0437           geomTools_.position((DetId)hAndR[ind].first, (hAndR[ind].second).first, (hAndR[ind].second).second).first;
0438       SimLCenergy += energy;
0439       SimLCx += position.x() * energy;
0440       SimLCy += position.y() * energy;
0441       SimLCz += position.z() * energy;
0442     };
0443 
0444     auto push_back_clu = [&](const uint32_t &SC_index, uint32_t &LC_index) {
0445       tmpLC.addCluEnergy(SimLCenergy);
0446       LocalPoint SimLCpos(SimLCx / SimLCenergy, SimLCy / SimLCenergy, SimLCz / SimLCenergy);
0447       tmpLC.addCluLocalPos(SimLCpos);
0448       SimLCenergy = 0.;
0449       SimLCx = 0.;
0450       SimLCy = 0.;
0451       SimLCz = 0.;
0452       tmpLC.addCluIndex(SC_index);
0453       tmpLC.computeClusterTime();
0454       output_.pMtdSimLayerClusters->push_back(tmpLC);
0455       LC_indices.push_back(LC_index);
0456       LC_index++;
0457       tmpLC.clear();
0458     };
0459 
0460     // fill tmpLC with the first hit
0461     push_back_hit(prev);
0462     update_clu_info(prev);
0463     for (const auto &ind : indices) {
0464       if (ind == indices[0])
0465         continue;
0466       DetId id(hAndR[ind].first);
0467       if (geomTools_.isETL(id) != geomTools_.isETL(prevId) or geomTools_.layer(id) != geomTools_.layer(prevId) or
0468           geomTools_.module(id) != geomTools_.module(prevId) or
0469           ((hAndR[ind].second).first == (hAndR[prev].second).first and
0470            (hAndR[ind].second).second != (hAndR[prev].second).second + 1) or
0471           ((hAndR[ind].second).second == (hAndR[prev].second).second and
0472            (hAndR[ind].second).first != (hAndR[prev].second).first + 1)) {
0473         // the next hit is not adjacent to the previous one, put the current temporary cluster in the collection
0474         // and the hit will be put in an empty temporary cluster
0475         push_back_clu(SC_index, LC_index);
0476       }
0477       // add the hit to the temporary cluster
0478       push_back_hit(ind);
0479       update_clu_info(ind);
0480       prev = ind;
0481       DetId newId(hAndR[prev].first);
0482       prevId = newId;
0483     }
0484     // add the remaining temporary cluster to the collection
0485     push_back_clu(SC_index, LC_index);
0486 
0487     // now the simTrackster: find position and time of the first simHit
0488     // bc right now there is no method to ask the simTrack for pos/time
0489     // at MTD entrance
0490     float timeAtEntrance = 99.;
0491     uint32_t idAtEntrance = 0;
0492     for (uint32_t i = 0; i < (uint32_t)hAndT.size(); i++) {
0493       if (hAndT[i].second < timeAtEntrance) {
0494         timeAtEntrance = hAndT[i].second;
0495         idAtEntrance = i;
0496       }
0497     }
0498 
0499     // sort LCs in the SimTrackster by time
0500     auto &MtdSimLayerClusters = output_.pMtdSimLayerClusters;
0501     std::sort(LC_indices.begin(), LC_indices.end(), [&MtdSimLayerClusters](int i, int j) {
0502       return (*MtdSimLayerClusters)[i].simLCTime() < (*MtdSimLayerClusters)[j].simLCTime();
0503     });
0504 
0505     GlobalPoint posAtEntrance = geomTools_
0506                                     .position((DetId)hAndR[idAtEntrance].first,
0507                                               (hAndR[idAtEntrance].second).first,
0508                                               (hAndR[idAtEntrance].second).second)
0509                                     .second;
0510     output_.pMtdSimTracksters->emplace_back(sc, LC_indices, timeAtEntrance, posAtEntrance);
0511     SC_index++;
0512   }
0513 
0514 #ifdef PRINT_DEBUG
0515   IfLogDebug(DEBUG, messageCategory_) << "SIMLAYERCLUSTERS LIST: \n";
0516   for (auto &sc : *output_.pMtdSimLayerClusters) {
0517     IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "SimLayerCluster with:"
0518                                         << "\n  CP charge " << sc.charge() << "\n  CP pdgId  " << sc.pdgId()
0519                                         << "\n  CP energy " << sc.energy() << " GeV\n  CP eta    " << sc.eta()
0520                                         << "\n  CP phi    " << sc.phi()
0521                                         << "\n  number of cells = " << sc.hits_and_fractions().size() << std::endl;
0522     for (unsigned int i = 0; i < sc.hits_and_fractions().size(); ++i) {
0523       DetId id(sc.detIds_and_rows()[i].first);
0524       IfLogDebug(DEBUG, messageCategory_)
0525           << std::fixed << std::setprecision(3) << " hit " << sc.detIds_and_rows()[i].first << " on layer "
0526           << geomTools_.layer(id) << " at time " << sc.hits_and_times()[i].second << " ns" << std::endl;
0527     }
0528     IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "  Cluster time " << sc.simLCTime()
0529                                         << " ns \n Cluster pos" << sc.simLCPos() << " cm\n"
0530                                         << std::fixed << std::setprecision(6) << " Cluster energy "
0531                                         << convertUnitsTo(0.001_MeV, sc.simLCEnergy()) << " MeV" << std::endl;
0532     IfLogDebug(DEBUG, messageCategory_) << "--------------\n";
0533   }
0534   IfLogDebug(DEBUG, messageCategory_) << std::endl;
0535 
0536   IfLogDebug(DEBUG, messageCategory_) << "SIMTRACKSTERS LIST: \n";
0537   for (auto &sc : *output_.pMtdSimTracksters) {
0538     IfLogDebug(DEBUG, messageCategory_) << std::fixed << std::setprecision(3) << "SimTrackster with:"
0539                                         << "\n  CP charge " << sc.charge() << "\n  CP pdgId  " << sc.pdgId()
0540                                         << "\n  CP energy " << sc.energy() << " GeV\n  CP eta    " << sc.eta()
0541                                         << "\n  CP phi    " << sc.phi()
0542                                         << "\n number of layer clusters = " << sc.numberOfClusters()
0543                                         << "\n time of first simhit " << sc.time() << " ns\n position of first simhit"
0544                                         << sc.position() << "cm" << std::endl;
0545     IfLogDebug(DEBUG, messageCategory_) << "  LCs indices: ";
0546     for (const auto &lc : sc.clusters())
0547       IfLogDebug(DEBUG, messageCategory_) << lc << ", ";
0548     IfLogDebug(DEBUG, messageCategory_) << "\n--------------\n";
0549   }
0550   IfLogDebug(DEBUG, messageCategory_) << std::endl;
0551 #endif
0552 
0553   event.put(std::move(output_.pMtdSimLayerClusters), "MergedMtdTruthLC");
0554   event.put(std::move(output_.pMtdSimTracksters), "MergedMtdTruthST");
0555 
0556   // now fill the calo particles
0557   for (unsigned i = 0; i < output_.pCaloParticles->size(); ++i) {
0558     auto &cp = (*output_.pCaloParticles)[i];
0559     for (unsigned j = m_caloParticles.sc_start_[i]; j < m_caloParticles.sc_stop_[i]; ++j) {
0560       edm::Ref<MtdSimClusterCollection> ref(scHandle, j);
0561       cp.addSimCluster(ref);
0562     }
0563   }
0564   event.put(std::move(output_.pCaloParticles), "MergedMtdTruth");
0565 
0566   calo_particles().swap(m_caloParticles);
0567 
0568   std::unordered_map<uint64_t, float>().swap(m_detIdToTotalSimEnergy);
0569   std::unordered_multimap<Barcode_t, Index_t>().swap(m_simHitBarcodeToIndex);
0570 }
0571 
0572 template <class T>
0573 void MtdTruthAccumulator::accumulateEvent(const T &event,
0574                                           const edm::EventSetup &setup,
0575                                           const edm::Handle<edm::HepMCProduct> &hepMCproduct) {
0576   edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
0577   edm::Handle<std::vector<int>> hGenParticleIndices;
0578 
0579   event.getByLabel(simTrackLabel_, hSimTracks);
0580   event.getByLabel(simVertexLabel_, hSimVertices);
0581 
0582   event.getByLabel(genParticleLabel_, hGenParticles);
0583   event.getByLabel(genParticleLabel_, hGenParticleIndices);
0584 
0585   std::vector<std::pair<uint64_t, const PSimHit *>> simHitPointers;
0586   std::unordered_map<int, std::map<uint64_t, float>> simTrackDetIdEnergyMap;
0587   std::unordered_map<int, std::map<uint64_t, float>> simTrackDetIdTimeMap;
0588   fillSimHits(simHitPointers, simTrackDetIdEnergyMap, simTrackDetIdTimeMap, event, setup);
0589 
0590   // Clear maps from previous event fill them for this one
0591   m_simHitBarcodeToIndex.clear();
0592   for (unsigned int i = 0; i < simHitPointers.size(); ++i) {
0593     m_simHitBarcodeToIndex.emplace(simHitPointers[i].second->trackId(), i);
0594   }
0595 
0596   auto const &tracks = *hSimTracks;
0597   auto const &vertices = *hSimVertices;
0598   std::unordered_map<int, int> trackid_to_track_index;
0599   std::unordered_map<uint32_t, float> vertex_time_map;
0600   DecayChain decay;
0601 
0602   for (uint32_t i = 0; i < vertices.size(); i++) {
0603     vertex_time_map[i] = vertices[i].position().t() * 1e9 + event.bunchCrossing() * bunchSpacing_;
0604   }
0605 
0606   IfLogDebug(DEBUG, messageCategory_) << " TRACKS" << std::endl;
0607   int idx = 0;
0608   for (auto const &t : tracks) {
0609     IfLogDebug(DEBUG, messageCategory_) << " " << idx << "\t" << t.trackId() << "\t" << t << std::endl;
0610     trackid_to_track_index[t.trackId()] = idx;
0611     idx++;
0612   }
0613 
0614   /**
0615   Build the main decay graph and assign the SimTrack to each edge. The graph
0616   built here will only contain the particles that have a decay vertex
0617   associated to them. In order to recover also the particles that will not
0618   decay, we need to keep track of the SimTrack used here and add, a-posteriori,
0619   the ones not used, associating a ghost vertex (starting from the highest
0620   simulated vertex number), in order to build the edge and identify them
0621   immediately as stable (i.e. not decayed).
0622 
0623   To take into account the multi-bremsstrahlung effects in which a single
0624   particle is emitting photons in different vertices **keeping the same
0625   track index**, we also collapsed those vertices into 1 unique vertex. The
0626   other approach of fully representing the decay chain keeping the same
0627   track index would have the problem of over-counting the contributions of
0628   that track, especially in terms of hits.
0629 
0630   The 2 auxiliary vectors are structured as follow:
0631 
0632   1. used_sim_tracks is a vector that has the same size as the overall
0633      number of simulated tracks. The associated integer is the vertexId of
0634      the **decaying vertex for that track**.
0635   2. collapsed_vertices is a vector that has the same size as the overall
0636      number of simulated vertices. The vector's index is the vertexId
0637      itself, the associated value is the vertexId of the vertex on which
0638      this should collapse.
0639   */
0640   idx = 0;
0641   std::vector<int> used_sim_tracks(tracks.size(), 0);
0642   std::vector<int> collapsed_vertices(vertices.size(), 0);
0643   IfLogDebug(DEBUG, messageCategory_) << " VERTICES" << std::endl;
0644   for (auto const &v : vertices) {
0645     IfLogDebug(DEBUG, messageCategory_) << " " << idx++ << "\t" << v << std::endl;
0646     if (v.parentIndex() != -1) {
0647       auto const trk_idx = trackid_to_track_index[v.parentIndex()];
0648       auto origin_vtx = tracks[trk_idx].vertIndex();
0649       if (used_sim_tracks[trk_idx]) {
0650         // collapse the vertex into the original first vertex we saw associated
0651         // to this track. Omit adding the edge in order to avoid double
0652         // counting of the very same particles  and its associated hits.
0653         collapsed_vertices[v.vertexId()] = used_sim_tracks[trk_idx];
0654         continue;
0655       }
0656       // Perform the actual vertex collapsing, if needed.
0657       if (collapsed_vertices[origin_vtx])
0658         origin_vtx = collapsed_vertices[origin_vtx];
0659       add_edge(origin_vtx,
0660                v.vertexId(),
0661                EdgeProperty(&tracks[trk_idx], simTrackDetIdEnergyMap[v.parentIndex()].size(), 0),
0662                decay);
0663       used_sim_tracks[trk_idx] = v.vertexId();
0664     }
0665   }
0666   // Build the motherParticle property to each vertex
0667   auto const &vertexMothersProp = get(vertex_name, decay);
0668   // Now recover the particles that did not decay. Append them with an index
0669   // bigger than the size of the generated vertices.
0670   int offset = vertices.size();
0671   for (size_t i = 0; i < tracks.size(); ++i) {
0672     if (!used_sim_tracks[i]) {
0673       auto origin_vtx = tracks[i].vertIndex();
0674       // Perform the actual vertex collapsing, if needed.
0675       if (collapsed_vertices[origin_vtx])
0676         origin_vtx = collapsed_vertices[origin_vtx];
0677       add_edge(
0678           origin_vtx, offset, EdgeProperty(&tracks[i], simTrackDetIdEnergyMap[tracks[i].trackId()].size(), 0), decay);
0679       // The properties for "fake" vertices associated to stable particles have
0680       // to be set inside this loop, since they do not belong to the vertices
0681       // collection and would be skipped by that loop (coming next)
0682       put(vertexMothersProp, offset, VertexProperty(&tracks[i], 0));
0683       offset++;
0684     }
0685   }
0686   for (auto const &v : vertices) {
0687     if (v.parentIndex() != -1) {
0688       // Skip collapsed_vertices
0689       if (collapsed_vertices[v.vertexId()])
0690         continue;
0691       put(vertexMothersProp, v.vertexId(), VertexProperty(&tracks[trackid_to_track_index[v.parentIndex()]], 0));
0692     }
0693   }
0694   SimHitsAccumulator_dfs_visitor vis;
0695   depth_first_search(decay, visitor(vis));
0696   CaloParticle_dfs_visitor caloParticleCreator(
0697       output_,
0698       m_caloParticles,
0699       m_simHitBarcodeToIndex,
0700       simTrackDetIdEnergyMap,
0701       simTrackDetIdTimeMap,
0702       vertex_time_map,
0703       [&](EdgeProperty &edge_property) -> bool {
0704         // Apply selection on SimTracks in order to promote them to be
0705         // CaloParticles. The function returns TRUE if the particle satisfies
0706         // the selection, FALSE otherwise. Therefore the correct logic to select
0707         // the particle is to ask for TRUE as return value.
0708         return (edge_property.cumulative_simHits != 0 and !edge_property.simTrack->noGenpart() and
0709                 edge_property.simTrack->momentum().E() > minEnergy_ and
0710                 std::abs(edge_property.simTrack->momentum().Eta()) < maxPseudoRapidity_);
0711       });
0712   depth_first_search(decay, visitor(caloParticleCreator));
0713 
0714 #if DEBUG
0715   boost::write_graphviz(std::cout,
0716                         decay,
0717                         make_label_writer(make_transform_value_property_map(&graphviz_vertex, get(vertex_name, decay))),
0718                         make_label_writer(make_transform_value_property_map(&graphviz_edge, get(edge_weight, decay))));
0719 #endif
0720 }
0721 
0722 template <class T>
0723 void MtdTruthAccumulator::fillSimHits(std::vector<std::pair<uint64_t, const PSimHit *>> &returnValue,
0724                                       std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdEnergyMap,
0725                                       std::unordered_map<int, std::map<uint64_t, float>> &simTrackDetIdTimeMap,
0726                                       const T &event,
0727                                       const edm::EventSetup &setup) {
0728   using namespace geant_units::operators;
0729   using namespace angle_units::operators;
0730   for (auto const &collectionTag : collectionTags_) {
0731     edm::Handle<std::vector<PSimHit>> hSimHits;
0732     event.getByLabel(collectionTag, hSimHits);
0733 
0734     for (auto const &simHit : *hSimHits) {
0735       DetId id(0);
0736 
0737       // --- Use only hits compatible with the in-time bunch-crossing
0738       if (simHit.tof() < 0 || simHit.tof() > 25.)
0739         continue;
0740 
0741       id = simHit.detUnitId();
0742 
0743       if (id == DetId(0)) {
0744         edm::LogWarning(messageCategory_) << "Invalid DetId for the current simHit!";
0745         continue;
0746       }
0747 
0748       if (simHit.trackId() == 0) {
0749         continue;
0750       }
0751 
0752       returnValue.emplace_back(id, &simHit);
0753 
0754       // get an unique id: for BTL the detId is unique (one for each crystal), for ETL the detId is not enough
0755       // also row and column are needed. An unique number is created from detId, row, col
0756       // Get row and column
0757       const auto &position = simHit.localPosition();
0758       LocalPoint simscaled(convertMmToCm(position.x()), convertMmToCm(position.y()), convertMmToCm(position.z()));
0759       std::pair<uint8_t, uint8_t> pixel = geomTools_.pixelInModule(id, simscaled);
0760       // create the unique id
0761       uint64_t uniqueId = static_cast<uint64_t>(id.rawId()) << 32;
0762       uniqueId |= pixel.first << 16;
0763       uniqueId |= pixel.second;
0764 
0765       simTrackDetIdEnergyMap[simHit.trackId()][uniqueId] += simHit.energyLoss();
0766       m_detIdToTotalSimEnergy[uniqueId] += simHit.energyLoss();
0767       // --- Get the time of the first SIM hit in the cell
0768       if (simTrackDetIdTimeMap[simHit.trackId()][uniqueId] == 0. ||
0769           simHit.tof() < simTrackDetIdTimeMap[simHit.trackId()][uniqueId]) {
0770         simTrackDetIdTimeMap[simHit.trackId()][uniqueId] = simHit.tof();
0771       }
0772 
0773 #ifdef PRINT_DEBUG
0774       IfLogDebug(DEBUG, messageCategory_)
0775           << "hitId " << id.rawId() << " from track " << simHit.trackId() << " in layer " << geomTools_.layer(id)
0776           << ", module " << geomTools_.module(id) << ", pixel ( " << (int)geomTools_.pixelInModule(id, simscaled).first
0777           << ", " << (int)geomTools_.pixelInModule(id, simscaled).second << " )\n global pos(cm) "
0778           << geomTools_.globalPosition(id, simscaled) << ", time(ns) " << simHit.tof() << ", energy(MeV) "
0779           << convertUnitsTo(0.001_MeV, simHit.energyLoss()) << std::endl;
0780 #endif
0781     }  // end of loop over simHits
0782   }    // end of loop over InputTags
0783 }
0784 
0785 // Register with the framework
0786 DEFINE_DIGI_ACCUMULATOR(MtdTruthAccumulator);