MixingModule/interface/DecayGraph.h

0001 #ifndef SimGeneral_MixingModule_DecayGraph_h
0002 #define SimGeneral_MixingModule_DecayGraph_h
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "SimDataFormats/Track/interface/SimTrack.h"
0006
0007 #if DEBUG
0008 // boost optional (used by boost graph) results in some false positives with
0009 // // -Wmaybe-uninitialized
0010 #pragma GCC diagnostic push
0011 #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
0012 #endif
0013
0014 // BOOST GRAPH LIBRARY
0015 #include <boost/graph/adjacency_list.hpp>
0016 #include <boost/graph/breadth_first_search.hpp>
0017 #include <boost/graph/depth_first_search.hpp>
0018 #include <boost/graph/graphviz.hpp>
0019
0020 using boost::add_edge;
0021 using boost::adjacency_list;
0022 using boost::directedS;
0023 using boost::edge;
0024 using boost::edge_weight;
0025 using boost::edge_weight_t;
0026 using boost::listS;
0027 using boost::property;
0028 using boost::vecS;
0029 using boost::vertex;
0030 using boost::vertex_name;
0031 using boost::vertex_name_t;
0032
0033 /* GRAPH DEFINITIONS
0034
0035    The graphs represent the full decay chain.
0036
0037    The parent-child relationship is the natural one, following "time".
0038
0039    Each edge has a property (edge_weight_t) that holds a const pointer to the
0040    SimTrack that connects the 2 vertices of the edge, the number of simHits
0041    associated to that simTrack and the cumulative number of simHits of itself
0042    and of all its children. Only simHits within the selected detectors are
0043    taken into account. The cumulative property is filled during the dfs
0044    exploration of the graph: if not explored the number is 0.
0045
0046    Each vertex has a property (vertex_name_t) that holds a const pointer to the
0047    SimTrack that originated that vertex and the cumulative number of simHits of
0048    all its outgoing edges. The cumulative property is filled during the dfs
0049    exploration of the graph: if not explored the number is 0.
0050
0051    Stable particles are recovered/added in a second iterations and are linked
0052    to ghost vertices with an offset starting from the highest generated vertex.
0053
0054    Multiple decays of a single particle that retains its original trackId are
0055    merged into one unique vertex (the first encountered) in order to avoid
0056    multiple counting of its associated simHits (if any).
0057
0058 */
0059 struct EdgeProperty {
0060   EdgeProperty(const SimTrack *t, int h, int c) : simTrack(t), simHits(h), cumulative_simHits(c) {}
0061   const SimTrack *simTrack;
0062   int simHits;
0063   int cumulative_simHits;
0064 };
0065
0066 struct VertexProperty {
0067   VertexProperty() : simTrack(nullptr), cumulative_simHits(0) {}
0068   VertexProperty(const SimTrack *t, int c) : simTrack(t), cumulative_simHits(c) {}
0069   const SimTrack *simTrack;
0070   int cumulative_simHits;
0071 };
0072
0073 using EdgeParticleClustersProperty = property<edge_weight_t, EdgeProperty>;
0074 using VertexMotherParticleProperty = property<vertex_name_t, VertexProperty>;
0075 using DecayChain = adjacency_list<listS, vecS, directedS, VertexMotherParticleProperty, EdgeParticleClustersProperty>;
0076
0077 namespace {
0078   extern const std::string messageCategoryGraph_;
0079
0080   template <typename Edge, typename Graph, typename Visitor>
0081   void accumulateSimHits_edge(Edge &e, const Graph &g, Visitor *v) {
0082     auto const edge_property = get(edge_weight, g, e);
0083     v->total_simHits += edge_property.simHits;
0084     IfLogDebug(DEBUG, messageCategoryGraph_)
0085         << " Examining edges " << e << " --> particle " << edge_property.simTrack->type() << "("
0086         << edge_property.simTrack->trackId() << ")"
0087         << " with SimClusters: " << edge_property.simHits << " Accumulated SimClusters: " << v->total_simHits
0088         << std::endl;
0089   }
0090   template <typename Vertex, typename Graph>
0091   void print_vertex(Vertex &u, const Graph &g) {
0092     auto const vertex_property = get(vertex_name, g, u);
0093     IfLogDebug(DEBUG, messageCategoryGraph_) << " At " << u;
0094     // The Mother of all vertices has **no** SimTrack associated.
0095     if (vertex_property.simTrack)
0096       IfLogDebug(DEBUG, messageCategoryGraph_) << " [" << vertex_property.simTrack->type() << "]"
0097                                                << "(" << vertex_property.simTrack->trackId() << ")";
0098     IfLogDebug(DEBUG, messageCategoryGraph_) << std::endl;
0099   }
0100
0101 // Graphviz output functions will only be generated in DEBUG mode
0102 #if DEBUG
0103   std::string graphviz_vertex(const VertexProperty &v) {
0104     std::ostringstream oss;
0105     oss << "{id: " << (v.simTrack ? v.simTrack->trackId() : 0) << ",\\ntype: " << (v.simTrack ? v.simTrack->type() : 0)
0106         << ",\\nchits: " << v.cumulative_simHits << "}";
0107     return oss.str();
0108   }
0109
0110   std::string graphviz_edge(const EdgeProperty &e) {
0111     std::ostringstream oss;
0112     oss << "[" << (e.simTrack ? e.simTrack->trackId() : 0) << "," << (e.simTrack ? e.simTrack->type() : 0) << ","
0113         << e.simHits << "," << e.cumulative_simHits << "]";
0114     return oss.str();
0115   }
0116 #endif
0117
0118   class SimHitsAccumulator_dfs_visitor : public boost::default_dfs_visitor {
0119   public:
0120     int total_simHits = 0;
0121     template <typename Edge, typename Graph>
0122     void examine_edge(Edge e, const Graph &g) {
0123       accumulateSimHits_edge(e, g, this);
0124     }
0125     template <typename Edge, typename Graph>
0126     void finish_edge(Edge e, const Graph &g) {
0127       auto const edge_property = get(edge_weight, g, e);
0128       auto src = source(e, g);
0129       auto trg = target(e, g);
0130       auto cumulative = edge_property.simHits + get(vertex_name, g, trg).cumulative_simHits +
0131                         (get(vertex_name, g, src).simTrack ? get(vertex_name, g, src).cumulative_simHits
0132                                                            : 0);  // when we hit the root vertex we have to stop
0133                                                                   // adding back its contribution.
0134       auto const src_vertex_property = get(vertex_name, g, src);
0135       put(get(vertex_name, const_cast<Graph &>(g)), src, VertexProperty(src_vertex_property.simTrack, cumulative));
0136       put(get(edge_weight, const_cast<Graph &>(g)),
0137           e,
0138           EdgeProperty(edge_property.simTrack, edge_property.simHits, cumulative));
0139       IfLogDebug(DEBUG, messageCategoryGraph_)
0140           << " Finished edge: " << e << " Track id: " << get(edge_weight, g, e).simTrack->trackId()
0141           << " has accumulated " << cumulative << " hits" << std::endl;
0142       IfLogDebug(DEBUG, messageCategoryGraph_) << " SrcVtx: " << src << "\t" << get(vertex_name, g, src).simTrack
0143                                                << "\t" << get(vertex_name, g, src).cumulative_simHits << std::endl;
0144       IfLogDebug(DEBUG, messageCategoryGraph_) << " TrgVtx: " << trg << "\t" << get(vertex_name, g, trg).simTrack
0145                                                << "\t" << get(vertex_name, g, trg).cumulative_simHits << std::endl;
0146     }
0147   };
0148
0149   using Selector = std::function<bool(EdgeProperty &)>;
0150 }  // namespace
0151
0152 #endif