Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-05-01 06:57:28 

0001 #define DEBUG false
0002 #if DEBUG
0003 // boost optional (used by boost graph) results in some false positives with
0004 // -Wmaybe-uninitialized
0005 #pragma GCC diagnostic push
0006 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
0007 #endif
0008 
0009 // BOOST GRAPH LIBRARY
0010 #include <boost/graph/adjacency_list.hpp>
0011 #include <boost/graph/breadth_first_search.hpp>
0012 #include <boost/graph/depth_first_search.hpp>
0013 #include <boost/graph/graphviz.hpp>
0014 
0015 #if DEBUG
0016 #pragma GCC diagnostic pop
0017 #endif
0018 
0019 #include <iterator>
0020 #include <memory>
0021 
0022 #include <numeric>  // for std::accumulate
0023 #include <unordered_map>
0024 
0025 #include "FWCore/Framework/interface/ConsumesCollector.h"
0026 #include "FWCore/Framework/interface/ESWatcher.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/EventSetup.h"
0029 #include "FWCore/Framework/interface/ProducesCollector.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033 
0034 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0035 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0036 #include "DataFormats/HcalDetId/interface/HcalTestNumbering.h"
0037 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0038 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0039 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0040 
0041 #include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
0042 #include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
0043 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
0044 #include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
0045 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0046 #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h"
0047 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0048 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0049 
0050 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixMod.h"
0051 #include "SimGeneral/MixingModule/interface/DigiAccumulatorMixModFactory.h"
0052 #include "SimGeneral/MixingModule/interface/PileUpEventPrincipal.h"
0053 #include "SimGeneral/TrackingAnalysis/interface/EncodedTruthId.h"
0054 
0055 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0056 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0057 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0058 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0059 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0060 
0061 namespace {
0062   using Index_t = unsigned;
0063   using Barcode_t = int;
0064   const std::string messageCategoryGraph_("CaloTruthAccumulatorGraphProducer");
0065 }  // namespace
0066 
0067 using boost::add_edge;
0068 using boost::adjacency_list;
0069 using boost::directedS;
0070 using boost::edge;
0071 using boost::edge_weight;
0072 using boost::edge_weight_t;
0073 using boost::listS;
0074 using boost::property;
0075 using boost::vecS;
0076 using boost::vertex;
0077 using boost::vertex_name;
0078 using boost::vertex_name_t;
0079 
0080 /* GRAPH DEFINITIONS
0081 
0082    The graphs represent the full decay chain.
0083 
0084    The parent-child relationship is the natural one, following "time".
0085 
0086    Each edge has a property (edge_weight_t) that holds a const pointer to the
0087    SimTrack that connects the 2 vertices of the edge, the number of simHits
0088    associated to that simTrack and the cumulative number of simHits of itself
0089    and of all its children. Only simHits within the selected detectors are
0090    taken into account. The cumulative property is filled during the dfs
0091    exploration of the graph: if not explored the number is 0.
0092 
0093    Each vertex has a property (vertex_name_t) that holds a const pointer to the
0094    SimTrack that originated that vertex and the cumulative number of simHits of
0095    all its outgoing edges. The cumulative property is filled during the dfs
0096    exploration of the graph: if not explored the number is 0.
0097 
0098    Stable particles are recovered/added in a second iterations and are linked
0099    to ghost vertices with an offset starting from the highest generated vertex.
0100 
0101    Multiple decays of a single particle that retains its original trackId are
0102    merged into one unique vertex (the first encountered) in order to avoid
0103    multiple counting of its associated simHits (if any).
0104 
0105 */
0106 struct EdgeProperty {
0107   EdgeProperty(const SimTrack *t, int h, int c) : simTrack(t), simHits(h), cumulative_simHits(c) {}
0108   const SimTrack *simTrack;
0109   int simHits;
0110   int cumulative_simHits;
0111 };
0112 
0113 struct VertexProperty {
0114   VertexProperty() : simTrack(nullptr), cumulative_simHits(0) {}
0115   VertexProperty(const SimTrack *t, int c) : simTrack(t), cumulative_simHits(c) {}
0116   VertexProperty(const VertexProperty &other)
0117       : simTrack(other.simTrack), cumulative_simHits(other.cumulative_simHits) {}
0118   const SimTrack *simTrack;
0119   int cumulative_simHits;
0120 };
0121 
0122 using EdgeParticleClustersProperty = property<edge_weight_t, EdgeProperty>;
0123 using VertexMotherParticleProperty = property<vertex_name_t, VertexProperty>;
0124 using DecayChain = adjacency_list<listS, vecS, directedS, VertexMotherParticleProperty, EdgeParticleClustersProperty>;
0125 
0126 class CaloTruthAccumulator : public DigiAccumulatorMixMod {
0127 public:
0128   explicit CaloTruthAccumulator(const edm::ParameterSet &config, edm::ProducesCollector, edm::ConsumesCollector &iC);
0129 
0130 private:
0131   void initializeEvent(const edm::Event &event, const edm::EventSetup &setup) override;
0132   void accumulate(const edm::Event &event, const edm::EventSetup &setup) override;
0133   void accumulate(const PileUpEventPrincipal &event, const edm::EventSetup &setup, edm::StreamID const &) override;
0134   void finalizeEvent(edm::Event &event, const edm::EventSetup &setup) override;
0135 
0136   /** @brief Both forms of accumulate() delegate to this templated method. */
0137   template <class T>
0138   void accumulateEvent(const T &event,
0139                        const edm::EventSetup &setup,
0140                        const edm::Handle<edm::HepMCProduct> &hepMCproduct);
0141 
0142   /** @brief Fills the supplied vector with pointers to the SimHits, checking
0143    * for bad modules if required */
0144   template <class T>
0145   void fillSimHits(std::vector<std::pair<DetId, const PCaloHit *>> &returnValue,
0146                    std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0147                    const T &event,
0148                    const edm::EventSetup &setup);
0149 
0150   const std::string messageCategory_;
0151 
0152   std::unordered_map<Index_t, float> m_detIdToTotalSimEnergy;  // keep track of cell normalizations
0153   std::unordered_multimap<Barcode_t, Index_t> m_simHitBarcodeToIndex;
0154 
0155   /** The maximum bunch crossing BEFORE the signal crossing to create
0156       TrackinParticles for. Use positive values. If set to zero no
0157       previous bunches are added and only in-time, signal and after bunches
0158       (defined by maximumSubsequentBunchCrossing_) are used.
0159   */
0160   const unsigned int maximumPreviousBunchCrossing_;
0161   /** The maximum bunch crossing AFTER the signal crossing to create
0162       TrackinParticles for. E.g. if set to zero only
0163       uses the signal and in time pileup (and previous bunches defined by the
0164       maximumPreviousBunchCrossing_ parameter).
0165   */
0166   const unsigned int maximumSubsequentBunchCrossing_;
0167 
0168   const edm::InputTag simTrackLabel_;
0169   const edm::InputTag simVertexLabel_;
0170   edm::Handle<std::vector<SimTrack>> hSimTracks;
0171   edm::Handle<std::vector<SimVertex>> hSimVertices;
0172 
0173   std::vector<edm::InputTag> collectionTags_;
0174   edm::InputTag genParticleLabel_;
0175   /// Needed to add HepMC::GenVertex to SimVertex
0176   edm::InputTag hepMCproductLabel_;
0177   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0178   edm::ESWatcher<CaloGeometryRecord> geomWatcher_;
0179 
0180   const double minEnergy_, maxPseudoRapidity_;
0181   const bool premixStage1_;
0182 
0183 public:
0184   struct OutputCollections {
0185     std::unique_ptr<SimClusterCollection> pSimClusters;
0186     std::unique_ptr<CaloParticleCollection> pCaloParticles;
0187   };
0188 
0189   struct calo_particles {
0190     std::vector<uint32_t> sc_start_;
0191     std::vector<uint32_t> sc_stop_;
0192 
0193     void swap(calo_particles &oth) {
0194       sc_start_.swap(oth.sc_start_);
0195       sc_stop_.swap(oth.sc_stop_);
0196     }
0197 
0198     void clear() {
0199       sc_start_.clear();
0200       sc_stop_.clear();
0201     }
0202   };
0203 
0204 private:
0205   const HGCalTopology *hgtopo_[3] = {nullptr, nullptr, nullptr};
0206   const HGCalDDDConstants *hgddd_[3] = {nullptr, nullptr, nullptr};
0207   const HcalDDDRecConstants *hcddd_ = nullptr;
0208   OutputCollections output_;
0209   calo_particles m_caloParticles;
0210   // geometry type (0 pre-TDR; 1 TDR)
0211   int geometryType_;
0212   bool doHGCAL;
0213 };
0214 
0215 /* Graph utility functions */
0216 
0217 namespace {
0218   template <typename Edge, typename Graph, typename Visitor>
0219   void accumulateSimHits_edge(Edge &e, const Graph &g, Visitor *v) {
0220     auto const edge_property = get(edge_weight, g, e);
0221     v->total_simHits += edge_property.simHits;
0222     IfLogDebug(DEBUG, messageCategoryGraph_)
0223         << " Examining edges " << e << " --> particle " << edge_property.simTrack->type() << "("
0224         << edge_property.simTrack->trackId() << ")"
0225         << " with SimClusters: " << edge_property.simHits << " Accumulated SimClusters: " << v->total_simHits
0226         << std::endl;
0227   }
0228   template <typename Vertex, typename Graph>
0229   void print_vertex(Vertex &u, const Graph &g) {
0230     auto const vertex_property = get(vertex_name, g, u);
0231     IfLogDebug(DEBUG, messageCategoryGraph_) << " At " << u;
0232     // The Mother of all vertices has **no** SimTrack associated.
0233     if (vertex_property.simTrack)
0234       IfLogDebug(DEBUG, messageCategoryGraph_) << " [" << vertex_property.simTrack->type() << "]"
0235                                                << "(" << vertex_property.simTrack->trackId() << ")";
0236     IfLogDebug(DEBUG, messageCategoryGraph_) << std::endl;
0237   }
0238 
0239 // Graphviz output functions will only be generated in DEBUG mode
0240 #if DEBUG
0241   std::string graphviz_vertex(const VertexProperty &v) {
0242     std::ostringstream oss;
0243     oss << "{id: " << (v.simTrack ? v.simTrack->trackId() : 0) << ",\\ntype: " << (v.simTrack ? v.simTrack->type() : 0)
0244         << ",\\nchits: " << v.cumulative_simHits << "}";
0245     return oss.str();
0246   }
0247 
0248   std::string graphviz_edge(const EdgeProperty &e) {
0249     std::ostringstream oss;
0250     oss << "[" << (e.simTrack ? e.simTrack->trackId() : 0) << "," << (e.simTrack ? e.simTrack->type() : 0) << ","
0251         << e.simHits << "," << e.cumulative_simHits << "]";
0252     return oss.str();
0253   }
0254 #endif
0255 
0256   class SimHitsAccumulator_dfs_visitor : public boost::default_dfs_visitor {
0257   public:
0258     int total_simHits = 0;
0259     template <typename Edge, typename Graph>
0260     void examine_edge(Edge e, const Graph &g) {
0261       accumulateSimHits_edge(e, g, this);
0262     }
0263     template <typename Edge, typename Graph>
0264     void finish_edge(Edge e, const Graph &g) {
0265       auto const edge_property = get(edge_weight, g, e);
0266       auto src = source(e, g);
0267       auto trg = target(e, g);
0268       auto cumulative = edge_property.simHits + get(vertex_name, g, trg).cumulative_simHits +
0269                         (get(vertex_name, g, src).simTrack ? get(vertex_name, g, src).cumulative_simHits
0270                                                            : 0);  // when we hit the root vertex we have to stop
0271                                                                   // adding back its contribution.
0272       auto const src_vertex_property = get(vertex_name, g, src);
0273       put(get(vertex_name, const_cast<Graph &>(g)), src, VertexProperty(src_vertex_property.simTrack, cumulative));
0274       put(get(edge_weight, const_cast<Graph &>(g)),
0275           e,
0276           EdgeProperty(edge_property.simTrack, edge_property.simHits, cumulative));
0277       IfLogDebug(DEBUG, messageCategoryGraph_)
0278           << " Finished edge: " << e << " Track id: " << get(edge_weight, g, e).simTrack->trackId()
0279           << " has accumulated " << cumulative << " hits" << std::endl;
0280       IfLogDebug(DEBUG, messageCategoryGraph_) << " SrcVtx: " << src << "\t" << get(vertex_name, g, src).simTrack
0281                                                << "\t" << get(vertex_name, g, src).cumulative_simHits << std::endl;
0282       IfLogDebug(DEBUG, messageCategoryGraph_) << " TrgVtx: " << trg << "\t" << get(vertex_name, g, trg).simTrack
0283                                                << "\t" << get(vertex_name, g, trg).cumulative_simHits << std::endl;
0284     }
0285   };
0286 
0287   using Selector = std::function<bool(EdgeProperty &)>;
0288 
0289   class CaloParticle_dfs_visitor : public boost::default_dfs_visitor {
0290   public:
0291     CaloParticle_dfs_visitor(CaloTruthAccumulator::OutputCollections &output,
0292                              CaloTruthAccumulator::calo_particles &caloParticles,
0293                              std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex,
0294                              std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0295                              Selector selector)
0296         : output_(output),
0297           caloParticles_(caloParticles),
0298           simHitBarcodeToIndex_(simHitBarcodeToIndex),
0299           simTrackDetIdEnergyMap_(simTrackDetIdEnergyMap),
0300           selector_(selector) {}
0301     template <typename Vertex, typename Graph>
0302     void discover_vertex(Vertex u, const Graph &g) {
0303       // If we reach the vertex 0, it means that we are backtracking with respect
0304       // to the first generation of stable particles: simply return;
0305       //    if (u == 0) return;
0306       print_vertex(u, g);
0307       auto const vertex_property = get(vertex_name, g, u);
0308       if (!vertex_property.simTrack)
0309         return;
0310       auto trackIdx = vertex_property.simTrack->trackId();
0311       IfLogDebug(DEBUG, messageCategoryGraph_)
0312           << " Found " << simHitBarcodeToIndex_.count(trackIdx) << " associated simHits" << std::endl;
0313       if (simHitBarcodeToIndex_.count(trackIdx)) {
0314         output_.pSimClusters->emplace_back(*vertex_property.simTrack);
0315         auto &simcluster = output_.pSimClusters->back();
0316         std::unordered_map<uint32_t, float> acc_energy;
0317         for (auto const &hit_and_energy : simTrackDetIdEnergyMap_[trackIdx]) {
0318           acc_energy[hit_and_energy.first] += hit_and_energy.second;
0319         }
0320         for (auto const &hit_and_energy : acc_energy) {
0321           simcluster.addRecHitAndFraction(hit_and_energy.first, hit_and_energy.second);
0322         }
0323       }
0324     }
0325     template <typename Edge, typename Graph>
0326     void examine_edge(Edge e, const Graph &g) {
0327       auto src = source(e, g);
0328       auto vertex_property = get(vertex_name, g, src);
0329       if (src == 0 or (vertex_property.simTrack == nullptr)) {
0330         auto edge_property = get(edge_weight, g, e);
0331         IfLogDebug(DEBUG, messageCategoryGraph_) << "Considering CaloParticle: " << edge_property.simTrack->trackId();
0332         if (selector_(edge_property)) {
0333           IfLogDebug(DEBUG, messageCategoryGraph_) << "Adding CaloParticle: " << edge_property.simTrack->trackId();
0334           output_.pCaloParticles->emplace_back(*(edge_property.simTrack));
0335           caloParticles_.sc_start_.push_back(output_.pSimClusters->size());
0336         }
0337       }
0338     }
0339 
0340     template <typename Edge, typename Graph>
0341     void finish_edge(Edge e, const Graph &g) {
0342       auto src = source(e, g);
0343       auto vertex_property = get(vertex_name, g, src);
0344       if (src == 0 or (vertex_property.simTrack == nullptr)) {
0345         auto edge_property = get(edge_weight, g, e);
0346         if (selector_(edge_property)) {
0347           caloParticles_.sc_stop_.push_back(output_.pSimClusters->size());
0348         }
0349       }
0350     }
0351 
0352   private:
0353     CaloTruthAccumulator::OutputCollections &output_;
0354     CaloTruthAccumulator::calo_particles &caloParticles_;
0355     std::unordered_multimap<Barcode_t, Index_t> &simHitBarcodeToIndex_;
0356     std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap_;
0357     Selector selector_;
0358   };
0359 }  // namespace
0360 
0361 CaloTruthAccumulator::CaloTruthAccumulator(const edm::ParameterSet &config,
0362                                            edm::ProducesCollector producesCollector,
0363                                            edm::ConsumesCollector &iC)
0364     : messageCategory_("CaloTruthAccumulator"),
0365       maximumPreviousBunchCrossing_(config.getParameter<unsigned int>("maximumPreviousBunchCrossing")),
0366       maximumSubsequentBunchCrossing_(config.getParameter<unsigned int>("maximumSubsequentBunchCrossing")),
0367       simTrackLabel_(config.getParameter<edm::InputTag>("simTrackCollection")),
0368       simVertexLabel_(config.getParameter<edm::InputTag>("simVertexCollection")),
0369       collectionTags_(),
0370       genParticleLabel_(config.getParameter<edm::InputTag>("genParticleCollection")),
0371       hepMCproductLabel_(config.getParameter<edm::InputTag>("HepMCProductLabel")),
0372       geomToken_(iC.esConsumes()),
0373       minEnergy_(config.getParameter<double>("MinEnergy")),
0374       maxPseudoRapidity_(config.getParameter<double>("MaxPseudoRapidity")),
0375       premixStage1_(config.getParameter<bool>("premixStage1")),
0376       geometryType_(-1),
0377       doHGCAL(config.getParameter<bool>("doHGCAL")) {
0378   producesCollector.produces<SimClusterCollection>("MergedCaloTruth");
0379   producesCollector.produces<CaloParticleCollection>("MergedCaloTruth");
0380   if (premixStage1_) {
0381     producesCollector.produces<std::vector<std::pair<unsigned int, float>>>("MergedCaloTruth");
0382   }
0383 
0384   iC.consumes<std::vector<SimTrack>>(simTrackLabel_);
0385   iC.consumes<std::vector<SimVertex>>(simVertexLabel_);
0386   iC.consumes<std::vector<reco::GenParticle>>(genParticleLabel_);
0387   iC.consumes<std::vector<int>>(genParticleLabel_);
0388   iC.consumes<std::vector<int>>(hepMCproductLabel_);
0389 
0390   // Fill the collection tags
0391   const edm::ParameterSet &simHitCollectionConfig = config.getParameterSet("simHitCollections");
0392   std::vector<std::string> parameterNames = simHitCollectionConfig.getParameterNames();
0393 
0394   for (auto const &parameterName : parameterNames) {
0395     std::vector<edm::InputTag> tags = simHitCollectionConfig.getParameter<std::vector<edm::InputTag>>(parameterName);
0396     collectionTags_.insert(collectionTags_.end(), tags.begin(), tags.end());
0397   }
0398 
0399   for (auto const &collectionTag : collectionTags_) {
0400     iC.consumes<std::vector<PCaloHit>>(collectionTag);
0401   }
0402 }
0403 
0404 void CaloTruthAccumulator::initializeEvent(edm::Event const &event, edm::EventSetup const &setup) {
0405   output_.pSimClusters = std::make_unique<SimClusterCollection>();
0406   output_.pCaloParticles = std::make_unique<CaloParticleCollection>();
0407 
0408   m_detIdToTotalSimEnergy.clear();
0409 
0410   if (geomWatcher_.check(setup)) {
0411     auto const &geom = setup.getData(geomToken_);
0412     const HGCalGeometry *eegeom = nullptr, *fhgeom = nullptr, *bhgeomnew = nullptr;
0413     const HcalGeometry *bhgeom = nullptr;
0414     bhgeom = static_cast<const HcalGeometry *>(geom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0415 
0416     if (doHGCAL) {
0417       eegeom = static_cast<const HGCalGeometry *>(
0418           geom.getSubdetectorGeometry(DetId::HGCalEE, ForwardSubdetector::ForwardEmpty));
0419       // check if it's the new geometry
0420       if (eegeom) {
0421         geometryType_ = 1;
0422         fhgeom = static_cast<const HGCalGeometry *>(
0423             geom.getSubdetectorGeometry(DetId::HGCalHSi, ForwardSubdetector::ForwardEmpty));
0424         bhgeomnew = static_cast<const HGCalGeometry *>(
0425             geom.getSubdetectorGeometry(DetId::HGCalHSc, ForwardSubdetector::ForwardEmpty));
0426       } else {
0427         geometryType_ = 0;
0428         eegeom = static_cast<const HGCalGeometry *>(geom.getSubdetectorGeometry(DetId::Forward, HGCEE));
0429         fhgeom = static_cast<const HGCalGeometry *>(geom.getSubdetectorGeometry(DetId::Forward, HGCHEF));
0430         bhgeom = static_cast<const HcalGeometry *>(geom.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0431       }
0432       hgtopo_[0] = &(eegeom->topology());
0433       hgtopo_[1] = &(fhgeom->topology());
0434       if (bhgeomnew)
0435         hgtopo_[2] = &(bhgeomnew->topology());
0436 
0437       for (unsigned i = 0; i < 3; ++i) {
0438         if (hgtopo_[i])
0439           hgddd_[i] = &(hgtopo_[i]->dddConstants());
0440       }
0441     }
0442 
0443     if (bhgeom) {
0444       hcddd_ = bhgeom->topology().dddConstants();
0445     }
0446   }
0447 }
0448 
0449 /** Create handle to edm::HepMCProduct here because event.getByLabel with
0450     edm::HepMCProduct only works for edm::Event but not for
0451     PileUpEventPrincipal; PileUpEventPrincipal::getByLabel tries to call
0452     T::value_type and T::iterator (where T is the type of the object one wants
0453     to get a handle to) which is only implemented for container-like objects
0454     like std::vector but not for edm::HepMCProduct!
0455 */
0456 void CaloTruthAccumulator::accumulate(edm::Event const &event, edm::EventSetup const &setup) {
0457   edm::Handle<edm::HepMCProduct> hepmc;
0458   event.getByLabel(hepMCproductLabel_, hepmc);
0459 
0460   edm::LogInfo(messageCategory_) << " CaloTruthAccumulator::accumulate (signal)";
0461   accumulateEvent(event, setup, hepmc);
0462 }
0463 
0464 void CaloTruthAccumulator::accumulate(PileUpEventPrincipal const &event,
0465                                       edm::EventSetup const &setup,
0466                                       edm::StreamID const &) {
0467   if (event.bunchCrossing() >= -static_cast<int>(maximumPreviousBunchCrossing_) &&
0468       event.bunchCrossing() <= static_cast<int>(maximumSubsequentBunchCrossing_)) {
0469     // simply create empty handle as we do not have a HepMCProduct in PU anyway
0470     edm::Handle<edm::HepMCProduct> hepmc;
0471     edm::LogInfo(messageCategory_) << " CaloTruthAccumulator::accumulate (pileup) bunchCrossing="
0472                                    << event.bunchCrossing();
0473     accumulateEvent(event, setup, hepmc);
0474   } else {
0475     edm::LogInfo(messageCategory_) << "Skipping pileup event for bunch crossing " << event.bunchCrossing();
0476   }
0477 }
0478 
0479 void CaloTruthAccumulator::finalizeEvent(edm::Event &event, edm::EventSetup const &setup) {
0480   edm::LogInfo(messageCategory_) << "Adding " << output_.pSimClusters->size() << " SimParticles and "
0481                                  << output_.pCaloParticles->size() << " CaloParticles to the event.";
0482 
0483   // We need to normalize the hits and energies into hits and fractions (since
0484   // we have looped over all pileup events)
0485   // For premixing stage1 we keep the energies, they will be normalized to
0486   // fractions in stage2
0487 
0488   if (premixStage1_) {
0489     auto totalEnergies = std::make_unique<std::vector<std::pair<unsigned int, float>>>();
0490     totalEnergies->reserve(m_detIdToTotalSimEnergy.size());
0491     std::copy(m_detIdToTotalSimEnergy.begin(), m_detIdToTotalSimEnergy.end(), std::back_inserter(*totalEnergies));
0492     std::sort(totalEnergies->begin(), totalEnergies->end());
0493     event.put(std::move(totalEnergies), "MergedCaloTruth");
0494   } else {
0495     for (auto &sc : *(output_.pSimClusters)) {
0496       auto hitsAndEnergies = sc.hits_and_fractions();
0497       sc.clearHitsAndFractions();
0498       sc.clearHitsEnergy();
0499       for (auto &hAndE : hitsAndEnergies) {
0500         const float totalenergy = m_detIdToTotalSimEnergy[hAndE.first];
0501         float fraction = 0.;
0502         if (totalenergy > 0)
0503           fraction = hAndE.second / totalenergy;
0504         else
0505           edm::LogWarning(messageCategory_)
0506               << "TotalSimEnergy for hit " << hAndE.first << " is 0! The fraction for this hit cannot be computed.";
0507         sc.addRecHitAndFraction(hAndE.first, fraction);
0508         sc.addHitEnergy(hAndE.second);
0509       }
0510     }
0511   }
0512 
0513   // save the SimCluster orphan handle so we can fill the calo particles
0514   auto scHandle = event.put(std::move(output_.pSimClusters), "MergedCaloTruth");
0515 
0516   // now fill the calo particles
0517   for (unsigned i = 0; i < output_.pCaloParticles->size(); ++i) {
0518     auto &cp = (*output_.pCaloParticles)[i];
0519     for (unsigned j = m_caloParticles.sc_start_[i]; j < m_caloParticles.sc_stop_[i]; ++j) {
0520       edm::Ref<SimClusterCollection> ref(scHandle, j);
0521       cp.addSimCluster(ref);
0522     }
0523   }
0524 
0525   event.put(std::move(output_.pCaloParticles), "MergedCaloTruth");
0526 
0527   calo_particles().swap(m_caloParticles);
0528 
0529   std::unordered_map<Index_t, float>().swap(m_detIdToTotalSimEnergy);
0530   std::unordered_multimap<Barcode_t, Index_t>().swap(m_simHitBarcodeToIndex);
0531 }
0532 
0533 template <class T>
0534 void CaloTruthAccumulator::accumulateEvent(const T &event,
0535                                            const edm::EventSetup &setup,
0536                                            const edm::Handle<edm::HepMCProduct> &hepMCproduct) {
0537   edm::Handle<std::vector<reco::GenParticle>> hGenParticles;
0538   edm::Handle<std::vector<int>> hGenParticleIndices;
0539 
0540   event.getByLabel(simTrackLabel_, hSimTracks);
0541   event.getByLabel(simVertexLabel_, hSimVertices);
0542 
0543   event.getByLabel(genParticleLabel_, hGenParticles);
0544   event.getByLabel(genParticleLabel_, hGenParticleIndices);
0545 
0546   std::vector<std::pair<DetId, const PCaloHit *>> simHitPointers;
0547   std::unordered_map<int, std::map<int, float>> simTrackDetIdEnergyMap;
0548   fillSimHits(simHitPointers, simTrackDetIdEnergyMap, event, setup);
0549 
0550   // Clear maps from previous event fill them for this one
0551   m_simHitBarcodeToIndex.clear();
0552   for (unsigned int i = 0; i < simHitPointers.size(); ++i) {
0553     m_simHitBarcodeToIndex.emplace(simHitPointers[i].second->geantTrackId(), i);
0554   }
0555 
0556   auto const &tracks = *hSimTracks;
0557   auto const &vertices = *hSimVertices;
0558   std::unordered_map<int, int> trackid_to_track_index;
0559   DecayChain decay;
0560   int idx = 0;
0561 
0562   IfLogDebug(DEBUG, messageCategory_) << " TRACKS" << std::endl;
0563   for (auto const &t : tracks) {
0564     IfLogDebug(DEBUG, messageCategory_) << " " << idx << "\t" << t.trackId() << "\t" << t << std::endl;
0565     trackid_to_track_index[t.trackId()] = idx;
0566     idx++;
0567   }
0568 
0569   /**
0570   Build the main decay graph and assign the SimTrack to each edge. The graph
0571   built here will only contain the particles that have a decay vertex
0572   associated to them. In order to recover also the particles that will not
0573   decay, we need to keep track of the SimTrack used here and add, a-posteriori,
0574   the ones not used, associating a ghost vertex (starting from the highest
0575   simulated vertex number), in order to build the edge and identify them
0576   immediately as stable (i.e. not decayed).
0577 
0578   To take into account the multi-bremsstrahlung effects in which a single
0579   particle is emitting photons in different vertices **keeping the same
0580   track index**, we also collapsed those vertices into 1 unique vertex. The
0581   other approach of fully representing the decay chain keeping the same
0582   track index would have the problem of over-counting the contributions of
0583   that track, especially in terms of hits.
0584 
0585   The 2 auxiliary vectors are structured as follow:
0586 
0587   1. used_sim_tracks is a vector that has the same size as the overall
0588      number of simulated tracks. The associated integer is the vertexId of
0589      the **decaying vertex for that track**.
0590   2. collapsed_vertices is a vector that has the same size as the overall
0591      number of simulated vertices. The vector's index is the vertexId
0592      itself, the associated value is the vertexId of the vertex on which
0593      this should collapse.
0594   */
0595   idx = 0;
0596   std::vector<int> used_sim_tracks(tracks.size(), 0);
0597   std::vector<int> collapsed_vertices(vertices.size(), 0);
0598   IfLogDebug(DEBUG, messageCategory_) << " VERTICES" << std::endl;
0599   for (auto const &v : vertices) {
0600     IfLogDebug(DEBUG, messageCategory_) << " " << idx++ << "\t" << v << std::endl;
0601     if (v.parentIndex() != -1) {
0602       auto trk_idx = trackid_to_track_index[v.parentIndex()];
0603       auto origin_vtx = tracks[trk_idx].vertIndex();
0604       if (used_sim_tracks[trk_idx]) {
0605         // collapse the vertex into the original first vertex we saw associated
0606         // to this track. Omit adding the edge in order to avoid double
0607         // counting of the very same particles  and its associated hits.
0608         collapsed_vertices[v.vertexId()] = used_sim_tracks[trk_idx];
0609         continue;
0610       }
0611       // Perform the actual vertex collapsing, if needed.
0612       if (collapsed_vertices[origin_vtx])
0613         origin_vtx = collapsed_vertices[origin_vtx];
0614       add_edge(origin_vtx,
0615                v.vertexId(),
0616                EdgeProperty(&tracks[trk_idx], simTrackDetIdEnergyMap[v.parentIndex()].size(), 0),
0617                decay);
0618       used_sim_tracks[trk_idx] = v.vertexId();
0619     }
0620   }
0621   // Build the motherParticle property to each vertex
0622   auto const &vertexMothersProp = get(vertex_name, decay);
0623   // Now recover the particles that did not decay. Append them with an index
0624   // bigger than the size of the generated vertices.
0625   int offset = vertices.size();
0626   for (size_t i = 0; i < tracks.size(); ++i) {
0627     if (!used_sim_tracks[i]) {
0628       auto origin_vtx = tracks[i].vertIndex();
0629       // Perform the actual vertex collapsing, if needed.
0630       if (collapsed_vertices[origin_vtx])
0631         origin_vtx = collapsed_vertices[origin_vtx];
0632       add_edge(
0633           origin_vtx, offset, EdgeProperty(&tracks[i], simTrackDetIdEnergyMap[tracks[i].trackId()].size(), 0), decay);
0634       // The properties for "fake" vertices associated to stable particles have
0635       // to be set inside this loop, since they do not belong to the vertices
0636       // collection and would be skipped by that loop (coming next)
0637       put(vertexMothersProp, offset, VertexProperty(&tracks[i], 0));
0638       offset++;
0639     }
0640   }
0641   for (auto const &v : vertices) {
0642     if (v.parentIndex() != -1) {
0643       // Skip collapsed_vertices
0644       if (collapsed_vertices[v.vertexId()])
0645         continue;
0646       put(vertexMothersProp, v.vertexId(), VertexProperty(&tracks[trackid_to_track_index[v.parentIndex()]], 0));
0647     }
0648   }
0649   SimHitsAccumulator_dfs_visitor vis;
0650   depth_first_search(decay, visitor(vis));
0651   CaloParticle_dfs_visitor caloParticleCreator(
0652       output_,
0653       m_caloParticles,
0654       m_simHitBarcodeToIndex,
0655       simTrackDetIdEnergyMap,
0656       [&](EdgeProperty &edge_property) -> bool {
0657         // Apply selection on SimTracks in order to promote them to be
0658         // CaloParticles. The function returns TRUE if the particle satisfies
0659         // the selection, FALSE otherwise. Therefore the correct logic to select
0660         // the particle is to ask for TRUE as return value.
0661         return (edge_property.cumulative_simHits != 0 and !edge_property.simTrack->noGenpart() and
0662                 edge_property.simTrack->momentum().E() > minEnergy_ and
0663                 std::abs(edge_property.simTrack->momentum().Eta()) < maxPseudoRapidity_);
0664       });
0665   depth_first_search(decay, visitor(caloParticleCreator));
0666 
0667 #if DEBUG
0668   boost::write_graphviz(std::cout,
0669                         decay,
0670                         make_label_writer(make_transform_value_property_map(&graphviz_vertex, get(vertex_name, decay))),
0671                         make_label_writer(make_transform_value_property_map(&graphviz_edge, get(edge_weight, decay))));
0672 #endif
0673 }
0674 
0675 template <class T>
0676 void CaloTruthAccumulator::fillSimHits(std::vector<std::pair<DetId, const PCaloHit *>> &returnValue,
0677                                        std::unordered_map<int, std::map<int, float>> &simTrackDetIdEnergyMap,
0678                                        const T &event,
0679                                        const edm::EventSetup &setup) {
0680   for (auto const &collectionTag : collectionTags_) {
0681     edm::Handle<std::vector<PCaloHit>> hSimHits;
0682     const bool isHcal = (collectionTag.instance().find("HcalHits") != std::string::npos);
0683     event.getByLabel(collectionTag, hSimHits);
0684 
0685     for (auto const &simHit : *hSimHits) {
0686       DetId id(0);
0687 
0688       //Relabel as necessary for HGCAL
0689       if (doHGCAL) {
0690         const uint32_t simId = simHit.id();
0691         if (geometryType_ == 1) {
0692           // no test numbering in new geometry
0693           id = simId;
0694         } else if (isHcal) {
0695           HcalDetId hid = HcalHitRelabeller::relabel(simId, hcddd_);
0696           if (hid.subdet() == HcalEndcap)
0697             id = hid;
0698         } else {
0699           int subdet, layer, cell, sec, subsec, zp;
0700           HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell);
0701           const HGCalDDDConstants *ddd = hgddd_[subdet - 3];
0702           std::pair<int, int> recoLayerCell = ddd->simToReco(cell, layer, sec, hgtopo_[subdet - 3]->detectorType());
0703           cell = recoLayerCell.first;
0704           layer = recoLayerCell.second;
0705           // skip simhits with bad barcodes or non-existant layers
0706           if (layer == -1 || simHit.geantTrackId() == 0)
0707             continue;
0708           id = HGCalDetId((ForwardSubdetector)subdet, zp, layer, subsec, sec, cell);
0709         }
0710       } else {
0711         id = simHit.id();
0712         //Relabel all HCAL hits
0713         if (isHcal) {
0714           HcalDetId hid = HcalHitRelabeller::relabel(simHit.id(), hcddd_);
0715           id = hid;
0716         }
0717       }
0718 
0719       if (id == DetId(0)) {
0720         continue;
0721       }
0722       if (simHit.geantTrackId() == 0) {
0723         continue;
0724       }
0725 
0726       returnValue.emplace_back(id, &simHit);
0727       simTrackDetIdEnergyMap[simHit.geantTrackId()][id.rawId()] += simHit.energy();
0728       m_detIdToTotalSimEnergy[id.rawId()] += simHit.energy();
0729     }
0730   }  // end of loop over InputTags
0731 }
0732 
0733 // Register with the framework
0734 DEFINE_DIGI_ACCUMULATOR(CaloTruthAccumulator);
[ 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!