CaloAnalysis/plugins/MtdTruthAccumulator.cc

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