Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-12-20 03:14:16 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!