From e8d8842e5f62339a3ce8d25e2935660757455d96 Mon Sep 17 00:00:00 2001 From: dchao Date: Tue, 30 Sep 2014 10:40:28 -0700 Subject: [PATCH] Completed truth match algorithm --- .../bdtaunu_tuple_analyzer/BDtaunuDef.h | 4 + .../BDtaunuGraphWriter.h | 137 +++++++++- .../bdtaunu_tuple_analyzer/BDtaunuMcReader.cc | 10 +- .../bdtaunu_tuple_analyzer/BDtaunuMcReader.h | 4 + .../bdtaunu_tuple_analyzer/GraphDef.h | 4 +- .../bdtaunu_tuple_analyzer/McGraphManager.cc | 4 +- .../bdtaunu_tuple_analyzer/McGraphManager.h | 3 + .../TruthMatchManager.cc | 234 ++++++++++++++++-- .../TruthMatchManager.h | 25 +- .../UpsilonCandidate.cc | 4 + .../bdtaunu_tuple_analyzer/UpsilonCandidate.h | 6 + .../bdtaunu_tuple_analyzer/tests/Makefile | 11 +- .../tests/merge_pdf.cmd | 1 - .../tests/truthmatch_test1.cc | 28 +++ .../tests/truthmatch_test2.cc | 56 +++++ .../tests/truthmatch_test3.cc | 215 ++++++++++++++++ 16 files changed, 714 insertions(+), 32 deletions(-) delete mode 100644 create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/merge_pdf.cmd create mode 100644 create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test1.cc create mode 100644 create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test2.cc create mode 100644 create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test3.cc diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuDef.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuDef.h index 233d477..02f17bd 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuDef.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuDef.h @@ -24,13 +24,17 @@ const int DcLund = 411; const int Dstar0Lund = 423; const int DstarcLund = 413; const int KSLund = 310; +const int K0Lund = 311; const int rhoLund = 213; const int pi0Lund = 111; const int KLund = 321; const int piLund = 211; const int eLund = 11; +const int nu_eLund = 12; const int muLund = 13; +const int nu_muLund = 14; const int tauLund = 15; +const int nu_tauLund = 16; const int gammaLund = 22; const int protonLund = 2212; const int neutronLund = 2112; diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuGraphWriter.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuGraphWriter.h index 07855d0..c9f21d9 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuGraphWriter.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuGraphWriter.h @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -41,9 +42,9 @@ void BDtaunuGraphWriter::operator() ( typename LundMap::const_iterator it = lund_map.find(lund_pm[v]); if (it != lund_map.end()) { - out << "[label=\"" << idx_pm[v] << ": " << it->second << "\",color=\"red\"]"; + out << "[label=\"" << idx_pm[v] << ": " << it->second << "\",color=\"blue\"]"; } else { - out << "[label=\"" << idx_pm[v] << ": " << lund_pm[v] << "\",color=\"red\"]"; + out << "[label=\"" << idx_pm[v] << ": " << lund_pm[v] << "\",color=\"blue\"]"; } } @@ -57,4 +58,136 @@ make_graph_writer(const GraphT &g, const LundMap lund_map, return BDtaunuGraphWriter(g, lund_map, lund_pm, idx_pm); } +template +class MyVertexWriter { + + private: + typedef typename boost::graph_traits::vertex_descriptor VertexT; + + public: + MyVertexWriter( + const GraphT &_graph, + const LundMap &_lund_map, + const TruthMatchMap &_tm_map, + const LundPM &_lund_pm, + const IdxPM &_idx_pm) : + g(_graph), + lund_map(_lund_map), + tm_map(_tm_map), + lund_pm(_lund_pm), + idx_pm(_idx_pm) {} + + void operator()(std::ostream &out, const VertexT &v) const; + + private: + GraphT g; + LundMap lund_map; + TruthMatchMap tm_map; + LundPM lund_pm; + IdxPM idx_pm; +}; + +template +void MyVertexWriter::operator() ( + std::ostream &out, const VertexT &v) const { + + typename LundMap::const_iterator it = lund_map.find(lund_pm[v]); + typename TruthMatchMap::const_iterator tm_it = tm_map.find(idx_pm[v]); + if (it != lund_map.end()) { + if (tm_it != tm_map.end() && tm_it->second >= 0) { + //out << "[label=\"" << idx_pm[v] << ": " << it->second << "\",color=\"blue\",penwidth=\"3\",fillcolor=\"cyan\",style=\"filled\"]"; + out << "[label=\"" << idx_pm[v] << ": " << it->second << "\",color=\"red\",penwidth=\"3\",fillcolor=\"lightskyblue\",style=\"filled\"]"; + } else { + out << "[label=\"" << idx_pm[v] << ": " << it->second << "\",color=\"red\"]"; + } + } else { + out << "[label=\"" << idx_pm[v] << ": " << lund_pm[v] << "\",color=\"red\"]"; + } + +} + +template +MyVertexWriter +make_myvertex_writer(const GraphT &g, const LundMap lund_map, const TruthMatchMap truthmatch_map, + const LundPM lund_pm, const IdxPM idx_pm) { + + return MyVertexWriter(g, lund_map, truthmatch_map, lund_pm, idx_pm); +} + +template +class MyEdgeWriter { + + private: + typedef typename boost::graph_traits::edge_descriptor EdgeT; + + public: + MyEdgeWriter( + const GraphT &_graph, + const LundMap &_lund_map, + const TruthMatchMap &_tm_map, + const LundPM &_lund_pm, + const IdxPM &_idx_pm) : + g(_graph), + lund_map(_lund_map), + tm_map(_tm_map), + lund_pm(_lund_pm), + idx_pm(_idx_pm) {} + + void operator()(std::ostream &out, const EdgeT &v) const; + + private: + GraphT g; + LundMap lund_map; + TruthMatchMap tm_map; + LundPM lund_pm; + IdxPM idx_pm; +}; + +template +void MyEdgeWriter::operator() ( + std::ostream &out, const EdgeT &e) const { + + auto u = source(e, g); + auto v = target(e, g); + typename TruthMatchMap::const_iterator u_tm_it, v_tm_it; + u_tm_it = tm_map.find(idx_pm[u]); + v_tm_it = tm_map.find(idx_pm[v]); + if ( + (u_tm_it != tm_map.end() && u_tm_it->second >= 0) && + (v_tm_it != tm_map.end() && v_tm_it->second >= 0) + ) { + out << "[color=\"black\",weight=\"1\", penwidth=\"3\"]"; + } else { + out << "[color=\"grey\"]"; + } + +} + +template +MyEdgeWriter +make_myedge_writer(const GraphT &g, const LundMap lund_map, const TruthMatchMap truthmatch_map, + const LundPM lund_pm, const IdxPM idx_pm) { + + return MyEdgeWriter(g, lund_map, truthmatch_map, lund_pm, idx_pm); +} + +class MyGraphWriter { + + public: + MyGraphWriter(std::string _title) : title(_title) {} + void operator()(std::ostream &out) const { + out << "graph[fontsize=\"32\"]" << std::endl; + out << "labelloc=\"t\"" << std::endl; + out << "label = " << "\"" << title << "\"" << std::endl; + }; + private: + std::string title; +}; + #endif diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.cc index cc2673c..1a371f6 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.cc +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.cc @@ -107,7 +107,7 @@ RootReader::Status BDtaunuMcReader::next_record() { mc_graph_manager.analyze_graph(); // Outsource truth match operations to truth matcher - truth_match_manager.update_graph(reco_graph_manager); + truth_match_manager.update_graph(reco_graph_manager, mc_graph_manager); truth_match_manager.analyze_graph(); // Make derived information ready for access @@ -131,5 +131,13 @@ void BDtaunuMcReader::FillMcInfo() { if (mc_graph_manager.get_mcB2()->tau) b2_tau_mctype = mc_graph_manager.get_mcB2()->tau->mc_type; } + + auto it = upsilon_candidates.begin(); + while (it != upsilon_candidates.end()) { + it->set_truth_match( + truth_match_manager.get_truth_match_status(it->get_reco_index())); + ++it; + } + return; } diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.h index 3d011b6..5ca5400 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/BDtaunuMcReader.h @@ -67,6 +67,10 @@ class BDtaunuMcReader : public BDtaunuReader { //! Printer void print_mc_graph(std::ostream &os) const { mc_graph_manager.print(os); } + void print_contracted_mc_graph(std::ostream &os) const { truth_match_manager.print_mc(os); } + void print_truthmatch_reco_graph(std::ostream &os) const { truth_match_manager.print_reco(os); } + + std::map get_truth_map() const { return truth_match_manager.get_truth_map(); } private: diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/GraphDef.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/GraphDef.h index 213dd39..e2f57d4 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/GraphDef.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/GraphDef.h @@ -37,7 +37,7 @@ namespace RecoGraph { // Boost graph typedefs typedef boost::adjacency_list< - boost::vecS, boost::vecS, boost::directedS, + boost::vecS, boost::vecS, boost::bidirectionalS, boost::property>>, @@ -150,7 +150,7 @@ namespace McGraph { // Boost graph typedefs typedef boost::adjacency_list< - boost::vecS, boost::vecS, boost::directedS, + boost::vecS, boost::vecS, boost::bidirectionalS, boost::property>, boost::property diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.cc index 5169897..fa2ce6d 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.cc +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.cc @@ -153,5 +153,7 @@ void McGraphManager::print(std::ostream &os) const { os, g, make_graph_writer(g, BDtaunuMcReader::lund_to_name, get(vertex_lund_id, g), - get(vertex_mc_index, g))); + get(vertex_mc_index, g)), + boost::default_writer(), + MyGraphWriter("MC Graph")); } diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.h index 6c2553c..5870c35 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/McGraphManager.h @@ -110,6 +110,9 @@ class McGraphManager : public GraphManager { //! Clear cache. void clear(); + //! Get the mc graph + const McGraph::Graph& get_mc_graph() const { return g; } + //! Returns pointer to the MC truth \f$\Upsilon(4S)\f$ if it exists, nullptr otherwise. const McGraph::Y* get_mcY() const; diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.cc index 432ed1c..1e15341 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.cc +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.cc @@ -1,15 +1,23 @@ +#include #include +#include #include +#include #include #include +#include +#include "BDtaunuDef.h" #include "TruthMatchManager.h" #include "BDtaunuMcReader.h" #include "GraphDef.h" #include "RecoGraphManager.h" +#include "McGraphManager.h" +#include "BDtaunuGraphWriter.h" using namespace boost; +using namespace bdtaunu; // TruthMatchMananger // ------------------ @@ -20,18 +28,22 @@ TruthMatchManager::TruthMatchManager() : reader(nullptr) { TruthMatchManager::TruthMatchManager(BDtaunuMcReader *_reader) : reader(_reader) { } -bool TruthMatchManager::is_truth_matched(int reco_idx) const { - std::map::const_iterator it = truth_match.find(reco_idx); +int TruthMatchManager::get_truth_match_status(int reco_idx) const { + std::map::const_iterator it = truth_match.find(reco_idx); assert(it != truth_match.end()); return it->second; } -void TruthMatchManager::update_graph(const RecoGraphManager &reco_graph_manager) { +void TruthMatchManager::update_graph( + const RecoGraphManager &reco_graph_manager, + const McGraphManager &mc_graph_manager) { hMCIdx = reader->hMCIdx; lMCIdx = reader->lMCIdx; gammaMCIdx = reader->gammaMCIdx; reco_graph = reco_graph_manager.get_reco_graph(); reco_indexer = reco_graph_manager.get_reco_indexer(); + mc_graph = mc_graph_manager.get_mc_graph(); + contract_mc_graph(); } void TruthMatchManager::analyze_graph() { @@ -39,6 +51,119 @@ void TruthMatchManager::analyze_graph() { depth_first_search(reco_graph, visitor(TruthMatchDfsVisitor(this))); } +bool TruthMatchManager::is_cleave_vertex( + const McGraph::Vertex &v, + const McGraph::Graph &g, + const McGraph::LundIdPropertyMap &lund_id_pm, + const McGraph::McIndexPropertyMap &mc_idx_pm) const { + + // neutrinos, tau, and K0 + switch (abs(lund_id_pm[v])) { + case nu_eLund: + case nu_muLund: + case nu_tauLund: + case tauLund: + case K0Lund: + return true; + } + + // initial e+e- beam particles + switch (mc_idx_pm[v]) { + case 0: + case 1: + return true; + } + + // elimimate final state particles' daughters + graph_traits::in_edge_iterator ie, ie_end; + tie(ie, ie_end) = in_edges(v, g); + if ((ie != ie_end) && (lund_id_pm[v] != UpsilonLund)) { + assert(ie + 1 == ie_end); + if (is_final_state_particle(lund_id_pm[source(*ie, g)])) { + return true; + } + } + + // eliminate MC added photons + if ((ie != ie_end) && (lund_id_pm[v] == gammaLund)) { + assert(ie + 1 == ie_end); + if (lund_id_pm[source(*ie, g)] != pi0Lund) { + return true; + } + } + + return false; + +} + +void TruthMatchManager::contract_mc_graph() { + + McGraph::Graph &g = mc_graph; + McGraph::McIndexPropertyMap mc_idx_pm = get(vertex_mc_index, g); + McGraph::LundIdPropertyMap lund_id_pm = get(vertex_lund_id, g); + + std::vector to_cleave; + graph_traits::vertex_iterator vi, vi_end; + for (tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + if (is_cleave_vertex(*vi, g, lund_id_pm, mc_idx_pm)) to_cleave.push_back(mc_idx_pm[*vi]); + + for (auto i : to_cleave) { + + graph_traits::vertex_iterator vi_begin, vi_end; + tie(vi_begin, vi_end) = vertices(g); + McGraph::Vertex v = *std::find_if(vi_begin, vi_end, + [i, mc_idx_pm] (const McGraph::Vertex &v) { return (i == mc_idx_pm[v]); }); + + McGraph::Vertex m; + graph_traits::in_edge_iterator ie, ie_end; + tie(ie, ie_end) = in_edges(v, g); + if (ie != ie_end) { + assert(ie + 1 == ie_end); + m = source(*ie, g); + McGraph::Vertex d; + graph_traits::out_edge_iterator oe, oe_end; + for (tie(oe, oe_end) = out_edges(v, g); oe != oe_end; ++oe) { + d = target(*oe, g); + add_edge(m, d, g); + } + } + + clear_vertex(v, g); + remove_vertex(v, g); + + } + +} + +void TruthMatchManager::print_mc(std::ostream &os) const { + boost::write_graphviz( + os, mc_graph, + make_graph_writer(mc_graph, BDtaunuMcReader::lund_to_name, + get(vertex_lund_id, mc_graph), + get(vertex_mc_index, mc_graph)), + boost::default_writer(), + MyGraphWriter("MC Graph with Edge Contraction")); +} + +void TruthMatchManager::print_reco(std::ostream &os) const { + boost::write_graphviz( + os, reco_graph, + make_myvertex_writer(reco_graph, BDtaunuMcReader::lund_to_name, + truth_match, + get(vertex_lund_id, reco_graph), + get(vertex_reco_index, reco_graph)), + make_myedge_writer(reco_graph, BDtaunuMcReader::lund_to_name, + truth_match, + get(vertex_lund_id, reco_graph), + get(vertex_reco_index, reco_graph)), + MyGraphWriter("Reco Graph with Truth Match")); +} + + + + + + // TruthMatchVisitor // ----------------- @@ -46,34 +171,103 @@ TruthMatchDfsVisitor::TruthMatchDfsVisitor() : manager(nullptr) { } TruthMatchDfsVisitor::TruthMatchDfsVisitor(TruthMatchManager *_manager) : manager(_manager) { - lund_pm = get(vertex_lund_id, manager->reco_graph); + reco_lund_pm = get(vertex_lund_id, manager->reco_graph); reco_idx_pm = get(vertex_reco_index, manager->reco_graph); block_idx_pm = get(vertex_block_index, manager->reco_graph); + mc_lund_pm = get(vertex_lund_id, manager->mc_graph); + mc_idx_pm = get(vertex_mc_index, manager->mc_graph); } void TruthMatchDfsVisitor::finish_vertex(RecoGraph::Vertex u, const RecoGraph::Graph &g) { + int lund = std::abs(reco_lund_pm[u]); + switch (lund) { + case bdtaunu::UpsilonLund: + case bdtaunu::B0Lund: + case bdtaunu::BcLund: + case bdtaunu::Dstar0Lund: + case bdtaunu::DstarcLund: + case bdtaunu::D0Lund: + case bdtaunu::DcLund: + case bdtaunu::KSLund: + case bdtaunu::rhoLund: + case bdtaunu::pi0Lund: + MatchCompositeState(u, g); + break; + case bdtaunu::eLund: + case bdtaunu::muLund: + case bdtaunu::piLund: + case bdtaunu::KLund: + case bdtaunu::gammaLund: + MatchFinalState(u); + break; + default: + assert(false); + return; + } +} + +void TruthMatchDfsVisitor::MatchFinalState(const RecoGraph::Vertex &u) { + const int *hitMap = nullptr; + switch (abs(reco_lund_pm[u])) { + case eLund: + case muLund: + hitMap = manager->lMCIdx; + break; + case piLund: + case KLund: + hitMap = manager->hMCIdx; + break; + case gammaLund: + hitMap = manager->gammaMCIdx; + break; + default: + assert(false); + } + + manager->truth_match.insert( + std::make_pair(reco_idx_pm[u], hitMap[block_idx_pm[u]])); + +} + +void TruthMatchDfsVisitor::MatchCompositeState(const RecoGraph::Vertex &u, const RecoGraph::Graph &g) { + + int tm_mc_idx = -1; - bool is_tm = false; + std::vector dau_tm_mc_idx; + RecoGraph::AdjacencyIterator ai, ai_end; + for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) { + dau_tm_mc_idx.push_back((manager->truth_match).find(reco_idx_pm[*ai])->second); + } + std::sort(dau_tm_mc_idx.begin(), dau_tm_mc_idx.end()); + + if (dau_tm_mc_idx[0] < 0) { + manager->truth_match.insert( + std::make_pair(reco_idx_pm[u], tm_mc_idx)); + return; + } - int reco_idx = get(reco_idx_pm, u); + graph_traits::vertex_iterator vi, vi_end; + for (tie(vi, vi_end) = vertices(manager->mc_graph); vi != vi_end; ++vi) { - if ((manager->reco_indexer).is_h_candidate(reco_idx)) { - if (manager->hMCIdx[get(block_idx_pm, u)] >= 0) is_tm = true; - (manager->truth_match).insert(std::make_pair(reco_idx, is_tm)); + if (reco_lund_pm[u] == mc_lund_pm[*vi]) { - } else if ((manager->reco_indexer).is_l_candidate(reco_idx)) { - if (manager->lMCIdx[get(block_idx_pm, u)] >= 0) is_tm = true; - (manager->truth_match).insert(std::make_pair(reco_idx, is_tm)); + std::vector dau_mc_idx; - } else if ((manager->reco_indexer).is_gamma_candidate(reco_idx)) { - if (manager->gammaMCIdx[get(block_idx_pm, u)] >= 0) is_tm = true; - (manager->truth_match).insert(std::make_pair(reco_idx, is_tm)); + McGraph::AdjacencyIterator bi, bi_end; + for (tie(bi, bi_end) = adjacent_vertices(*vi, manager->mc_graph); bi != bi_end; ++bi) { + dau_mc_idx.push_back(mc_idx_pm[*bi]); + } + std::sort(dau_mc_idx.begin(), dau_mc_idx.end()); - } else { - RecoGraph::AdjacencyIterator ai, ai_end; - for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) { - if (!(manager->truth_match).find(get(reco_idx_pm, *ai))->second) is_tm = false; break; + if (dau_tm_mc_idx == dau_mc_idx) { + tm_mc_idx = mc_idx_pm[*vi]; + break; + } } - (manager->truth_match).insert(std::make_pair(reco_idx, is_tm)); + } + + manager->truth_match.insert( + std::make_pair(reco_idx_pm[u], tm_mc_idx)); + } diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.h index 3666339..6f7acd3 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/TruthMatchManager.h @@ -6,6 +6,7 @@ #include "GraphDef.h" #include "RecoGraphManager.h" +#include "McGraphManager.h" class BDtaunuMcReader; @@ -19,14 +20,17 @@ class TruthMatchManager { TruthMatchManager(const TruthMatchManager&) = default; TruthMatchManager &operator=(const TruthMatchManager&) = default; - bool is_truth_matched(int reco_idx) const; + int get_truth_match_status(int reco_idx) const; + std::map get_truth_map() const { return truth_match; } - void update_graph(const RecoGraphManager&); + void update_graph(const RecoGraphManager&, const McGraphManager&); void analyze_graph(); + void print_mc(std::ostream &os) const; + void print_reco(std::ostream &os) const; private: - std::map truth_match; + std::map truth_match; BDtaunuMcReader *reader; const int *hMCIdx; @@ -34,6 +38,14 @@ class TruthMatchManager { const int *gammaMCIdx; RecoGraph::Graph reco_graph; RecoGraph::RecoIndexer reco_indexer; + McGraph::Graph mc_graph; + + void contract_mc_graph(); + bool is_cleave_vertex( + const McGraph::Vertex &v, + const McGraph::Graph &g, + const McGraph::LundIdPropertyMap &lund_id_pm, + const McGraph::McIndexPropertyMap &mc_idx_pm) const; }; class TruthMatchDfsVisitor : public boost::default_dfs_visitor { @@ -47,9 +59,14 @@ class TruthMatchDfsVisitor : public boost::default_dfs_visitor { private: TruthMatchManager *manager; - RecoGraph::LundIdPropertyMap lund_pm; + RecoGraph::LundIdPropertyMap reco_lund_pm; RecoGraph::RecoIndexPropertyMap reco_idx_pm; RecoGraph::BlockIndexPropertyMap block_idx_pm; + McGraph::LundIdPropertyMap mc_lund_pm; + McGraph::McIndexPropertyMap mc_idx_pm; + + void MatchFinalState(const RecoGraph::Vertex &u); + void MatchCompositeState(const RecoGraph::Vertex &u, const RecoGraph::Graph &g); }; #endif diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.cc index fb349cf..e288915 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.cc +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.cc @@ -10,6 +10,7 @@ UpsilonCandidate::UpsilonCandidate() : eventId(""), block_index(-999), reco_index(-1), + truth_match(-1), bflavor(BFlavor::null), eextra50(-999), mmiss_prime2(-999), @@ -45,6 +46,7 @@ UpsilonCandidate::UpsilonCandidate( std::string &_eventId, int _block_index, int _reco_index, + int _truth_match, BFlavor _bflavor, float _eextra50, float _mmiss_prime2, @@ -78,6 +80,7 @@ UpsilonCandidate::UpsilonCandidate( eventId(_eventId), block_index(_block_index), reco_index(_reco_index), + truth_match(_truth_match), bflavor(_bflavor), eextra50(_eextra50), mmiss_prime2(_mmiss_prime2), @@ -124,6 +127,7 @@ void UpsilonCandidate::copy_candidate(const UpsilonCandidate &cand) { eventId = cand.eventId; block_index = cand.block_index; reco_index = cand.reco_index; + truth_match = cand.truth_match; bflavor = cand.bflavor; eextra50 = cand.eextra50; mmiss_prime2 = cand.mmiss_prime2; diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.h b/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.h index 5b29a89..ed61d10 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.h +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/UpsilonCandidate.h @@ -17,6 +17,7 @@ class UpsilonCandidate { std::string &eventId, int block_index, int reco_index, + int truth_match, bdtaunu::BFlavor bflavor, float eextra50, float mmiss_prime2, @@ -64,6 +65,10 @@ class UpsilonCandidate { int get_reco_index() const { return reco_index; } void set_reco_index(int _reco_index) { reco_index = _reco_index; } + //! The index that uniquely identifies the reconstructed candidate within the event. + int get_truth_match() const { return truth_match; } + void set_truth_match(int _truth_match) { truth_match = _truth_match; } + //! \f$ E_{extra} \f$ /*! The energy sum of photons above 50 MeV that are not used in the * reconstruction of this candidate. */ @@ -248,6 +253,7 @@ class UpsilonCandidate { std::string eventId; int block_index; int reco_index; + int truth_match; bdtaunu::BFlavor bflavor; float eextra50; float mmiss_prime2; diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/Makefile b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/Makefile index 444af8a..1d4e51d 100644 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/Makefile +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/Makefile @@ -9,7 +9,7 @@ CXXFLAGS += -fPIC # Contents # -------- -BINARIES = mcreader_test1 mcreader_test2 mcreader_test3 +BINARIES = mcreader_test1 mcreader_test2 mcreader_test3 truthmatch_test1 truthmatch_test2 truthmatch_test3 # Dependencies # ------------ @@ -32,6 +32,8 @@ LDFLAGS += $(shell root-config --libs) # Build Rules # ----------- +.PHONY: all debug pdf pdf_clean clean + $@ : $@.cc $(CXX) $(CXXFLAGS) $< -o $% @@ -44,5 +46,12 @@ debug : $(BINARIES) $(BINARIES) : % : %.cc $(CXX) $(CXXFLAGS) $(INCFLAGS) $(LDFLAGS) -Wl,-rpath,$(TUPLE_READER_LIB_PATH) -o $@ $< +pdf: + printf "%s\n" *.gv | sed s'/^\(.*\).gv/ -Tpdf \1.gv -o \1.pdf/' | xargs -L 1 dot + pdftk *.pdf cat output merged.pdf + +pdf_clean: + rm -f *.gv *.pdf + clean: rm -f *~ *.o *.gv *.pdf `find . -perm +100 -type f` diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/merge_pdf.cmd b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/merge_pdf.cmd deleted file mode 100644 index 7cbf895..0000000 --- a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/merge_pdf.cmd +++ /dev/null @@ -1 +0,0 @@ -ls -l *.gv | sed s'/^.*\(sp[0-9]*.*\.gv\).*$/\1/' | sed s'/\(.*\)\.gv/dot -Tpdf \1.gv -o \1.pdf/' > tmp.cmd; source tmp.cmd; rm tmp.cmd; pdftk *.pdf cat output merged.pdf diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test1.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test1.cc new file mode 100644 index 0000000..ca4bf4f --- /dev/null +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test1.cc @@ -0,0 +1,28 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; + +int main() { + BDtaunuMcReader reader("/Users/dchao/bdtaunu/v4/data/root/signal/aug_12_2014/A/sp11444r1.root"); + int nTM = 0, nevents = 0; + while (reader.next_record() != RootReader::Status::kEOF) { + ++nevents; + vector upsilons = reader.get_upsilon_candidates(); + for (auto ups : upsilons) { + if (ups.get_truth_match() >= 0) ++nTM; + break; + } + } + cout << nTM << " / " << nevents << endl; + return 0; +} diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test2.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test2.cc new file mode 100644 index 0000000..b4df88e --- /dev/null +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test2.cc @@ -0,0 +1,56 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; + +int main() { + ofstream reco_output; + ofstream mc_output; + ofstream mc_contracted_output; + + BDtaunuMcReader reader("/Users/dchao/bdtaunu/v4/data/root/signal/aug_12_2014/A/sp11444r1.root"); + + reader.next_record(); + reco_output.open("1.reco.gv"); + mc_output.open("1.mc.gv"); + mc_contracted_output.open("1.mc_contracted.gv"); + reader.print_mc_graph(mc_output); + reader.print_contracted_mc_graph(mc_contracted_output); + reader.print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader.next_record(); + + reader.next_record(); + reco_output.open("2.reco.gv"); + mc_output.open("2.mc.gv"); + mc_contracted_output.open("2.mc_contracted.gv"); + reader.print_mc_graph(mc_output); + reader.print_contracted_mc_graph(mc_contracted_output); + reader.print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + map truth_map = reader.get_truth_map(); + auto it = truth_map.begin(); + auto it_end = truth_map.end(); + while (it != it_end) { + cout << it->first << " : " << it->second << endl; + ++it; + } + vector upsilons = reader.get_upsilon_candidates(); + cout << upsilons[1].get_reco_index() << endl; + cout << upsilons[1].get_truth_match() << endl; + return 0; +} diff --git a/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test3.cc b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test3.cc new file mode 100644 index 0000000..c292b40 --- /dev/null +++ b/create_database/tuple_reader/bdtaunu_tuple_analyzer/tests/truthmatch_test3.cc @@ -0,0 +1,215 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +using namespace std; + +int main() { + ofstream reco_output; + ofstream mc_output; + ofstream mc_contracted_output; + + BDtaunuMcReader *reader; + + reader = new BDtaunuMcReader("/Users/dchao/bdtaunu/v4/data/root/signal/aug_12_2014/A/sp11444r1.root"); + + reader->next_record(); + reco_output.open("1.reco.gv"); + mc_output.open("1.mc.gv"); + mc_contracted_output.open("1.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("2.reco.gv"); + mc_output.open("2.mc.gv"); + mc_contracted_output.open("2.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("3.reco.gv"); + mc_output.open("3.mc.gv"); + mc_contracted_output.open("3.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("4.reco.gv"); + mc_output.open("4.mc.gv"); + mc_contracted_output.open("4.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + delete reader; + + reader = new BDtaunuMcReader("/Users/dchao/bdtaunu/v4/data/root/signal/aug_12_2014/A/sp11445r1.root"); + + reader->next_record(); + reco_output.open("5.reco.gv"); + mc_output.open("5.mc.gv"); + mc_contracted_output.open("5.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("6.reco.gv"); + mc_output.open("6.mc.gv"); + mc_contracted_output.open("6.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("7.reco.gv"); + mc_output.open("7.mc.gv"); + mc_contracted_output.open("7.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("8.reco.gv"); + mc_output.open("8.mc.gv"); + mc_contracted_output.open("8.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + delete reader; + + reader = new BDtaunuMcReader("/Users/dchao/bdtaunu/v4/data/root/signal/aug_12_2014/A/sp11446r1.root"); + + reader->next_record(); + reco_output.open("9.reco.gv"); + mc_output.open("9.mc.gv"); + mc_contracted_output.open("9.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("10.reco.gv"); + mc_output.open("10.mc.gv"); + mc_contracted_output.open("10.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("11.reco.gv"); + mc_output.open("11.mc.gv"); + mc_contracted_output.open("11.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("12.reco.gv"); + mc_output.open("12.mc.gv"); + mc_contracted_output.open("12.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + delete reader; + + reader = new BDtaunuMcReader("/Users/dchao/bdtaunu/v4/data/root/signal/aug_12_2014/A/sp11447r1.root"); + + reader->next_record(); + reco_output.open("13.reco.gv"); + mc_output.open("13.mc.gv"); + mc_contracted_output.open("13.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("14.reco.gv"); + mc_output.open("14.mc.gv"); + mc_contracted_output.open("14.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("15.reco.gv"); + mc_output.open("15.mc.gv"); + mc_contracted_output.open("15.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + reader->next_record(); + reco_output.open("16.reco.gv"); + mc_output.open("16.mc.gv"); + mc_contracted_output.open("16.mc_contracted.gv"); + reader->print_mc_graph(mc_output); + reader->print_contracted_mc_graph(mc_contracted_output); + reader->print_truthmatch_reco_graph(reco_output); + reco_output.close(); + mc_output.close(); + mc_contracted_output.close(); + + delete reader; + + return 0; +}