diff --git a/create_database/BDtaunuGraphWriter.h b/create_database/BDtaunuGraphWriter.h
index 2abe80c..07855d0 100644
--- a/create_database/BDtaunuGraphWriter.h
+++ b/create_database/BDtaunuGraphWriter.h
@@ -7,9 +7,6 @@
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/graphviz.hpp>
 
-#include "BDtaunuMcReader.h"
-#include "utilities/helpers.h"
-
 template <typename GraphT, typename LundMap,
           typename LundPM, typename IdxPM>
 class BDtaunuGraphWriter {
diff --git a/create_database/BDtaunuMcReader.cc b/create_database/BDtaunuMcReader.cc
index be54c00..be7daa0 100644
--- a/create_database/BDtaunuMcReader.cc
+++ b/create_database/BDtaunuMcReader.cc
@@ -1,405 +1,97 @@
-#include "BDtaunuReader.h"
-#include "BDtaunuMcReader.h"
-#include "utilities/helpers.h"
-#include "bdtaunu_definitions.h"
-
 #include <map>
 #include <string>
 #include <cmath>
 #include <cassert>
-#include <algorithm>
+
+#include "utilities/helpers.h"
+#include "bdtaunu_definitions.h"
+#include "BDtaunuReader.h"
+#include "BDtaunuReaderStatus.h"
+#include "BDtaunuMcReader.h"
+#include "McGraphManager.h"
+#include "McBTypeCatalogue.h"
 
 using namespace boost;
 
-const double BDtaunuMcReader::min_photon_energy = 0.02;
 const int BDtaunuMcReader::max_mc_length = 100;
 
-// Static methods to initialize the particle lundId arrays that are
-// needed to determine MC B meson types. 
-// Charged leptons. 
-std::vector<int> BDtaunuMcReader::build_ell() {
-  std::vector<int> ell_temp;
-  ell_temp.push_back(abs(lundIdMap["e-"]));
-  ell_temp.push_back(abs(lundIdMap["mu-"]));
-  ell_temp.push_back(abs(lundIdMap["tau-"]));
-  std::sort(ell_temp.begin(), ell_temp.end());
-  return ell_temp;
-}
-const std::vector<int> BDtaunuMcReader::ell = build_ell();
-
-// Neutrinos. 
-std::vector<int> BDtaunuMcReader::build_nu() {
-  std::vector<int> nu_temp;
-  nu_temp.push_back(abs(lundIdMap["nu_e"]));
-  nu_temp.push_back(abs(lundIdMap["nu_mu"]));
-  nu_temp.push_back(abs(lundIdMap["nu_tau"]));
-  std::sort(nu_temp.begin(), nu_temp.end());
-  return nu_temp;
-}
-const std::vector<int> BDtaunuMcReader::nu = build_nu();
-
-// D mesons. 
-std::vector<int> BDtaunuMcReader::build_dmeson() {
-  std::vector<int> dmeson_temp;
-  dmeson_temp.push_back(abs(lundIdMap["D+"]));
-  dmeson_temp.push_back(abs(lundIdMap["D0"]));
-  std::sort(dmeson_temp.begin(), dmeson_temp.end());
-  return dmeson_temp;
-}
-const std::vector<int> BDtaunuMcReader::dmeson = build_dmeson();
-
-// D* mesons. 
-std::vector<int> BDtaunuMcReader::build_dstar() {
-  std::vector<int> dstar_temp;
-  dstar_temp.push_back(abs(lundIdMap["D*+"]));
-  dstar_temp.push_back(abs(lundIdMap["D*0"]));
-  std::sort(dstar_temp.begin(), dstar_temp.end());
-  return dstar_temp;
-}
-const std::vector<int> BDtaunuMcReader::dstar = build_dstar();
-
-// D** mesons. 
-std::vector<int> BDtaunuMcReader::build_dstarstar() {
-  std::vector<int> dstarstar_temp;
-  dstarstar_temp.push_back(abs(lundIdMap["D_0*+"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D_0*0"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D_1+"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D_10"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D_2*+"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D_2*0"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D'_1+"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D'_10"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D(2S)+"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D(2S)0"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D*(2S)+"]));
-  dstarstar_temp.push_back(abs(lundIdMap["D*(2S)0"]));
-  std::sort(dstarstar_temp.begin(), dstarstar_temp.end());
-  return dstarstar_temp;
-}
-const std::vector<int> BDtaunuMcReader::dstarstar = build_dstarstar();
-
-// Ds(*) mesons.
-std::vector<int> BDtaunuMcReader::build_dstrange() {
-  std::vector<int> dstrange_temp;
-  dstrange_temp.push_back(abs(lundIdMap["D_s+"]));
-  dstrange_temp.push_back(abs(lundIdMap["D_s*+"]));
-  dstrange_temp.push_back(abs(lundIdMap["D_s0*+"]));
-  dstrange_temp.push_back(abs(lundIdMap["D_s1+"]));
-  dstrange_temp.push_back(abs(lundIdMap["D_s2*+"]));
-  dstrange_temp.push_back(abs(lundIdMap["D'_s1+"]));
-  std::sort(dstrange_temp.begin(), dstrange_temp.end());
-  return dstrange_temp;
-}
-const std::vector<int> BDtaunuMcReader::dstrange = build_dstrange();
-
-// pion. 
-std::vector<int> BDtaunuMcReader::build_pion() {
-  std::vector<int> pion_temp;
-  pion_temp.push_back(abs(lundIdMap["pi+"]));
-  pion_temp.push_back(abs(lundIdMap["pi0"]));
-  std::sort(pion_temp.begin(), pion_temp.end());
-  return pion_temp;
-}
-const std::vector<int> BDtaunuMcReader::pion = build_pion();
-
+BDtaunuMcReader::BDtaunuMcReader(
+  const char *root_fname,
+  const char *root_trname) : 
+  BDtaunuReader(root_fname, root_trname) {
 
-BDtaunuMcReader::BDtaunuMcReader() : BDtaunuReader() {
-  Initialize();
-}
+  AllocateBuffer();
+  ClearBuffer();
 
-BDtaunuMcReader::BDtaunuMcReader(const char *root_fname) : 
-  BDtaunuReader(root_fname) {
-  Initialize();
-}
+  mc_graph_manager = McGraphManager(this);
 
-BDtaunuMcReader::BDtaunuMcReader(
-    const char *root_fname,
-    const char *root_trname) : 
-  BDtaunuReader(root_fname, root_trname) {
-  Initialize();
 }
 
 BDtaunuMcReader::~BDtaunuMcReader() {
-  delete[] mcLund;
-  delete[] mothIdx;
-  delete[] dauIdx;
-  delete[] dauLen;
-  delete[] mcenergy;
+  DeleteBuffer();
 }
 
 
-void BDtaunuMcReader::Initialize() {
+void BDtaunuMcReader::AllocateBuffer() {
+
   mcLund = new int[max_mc_length];
   mothIdx = new int[max_mc_length];
   dauIdx = new int[max_mc_length];
   dauLen = new int[max_mc_length];
   mcenergy = new float[max_mc_length];
-  SetBranchAddress();
-  ClearColumnValues();
-}
 
-void BDtaunuMcReader::SetBranchAddress() {
-  BDtaunuReader::SetBranchAddress();
   tr->SetBranchAddress("mcLen", &mcLen);
   tr->SetBranchAddress("mcLund", mcLund);
   tr->SetBranchAddress("mothIdx", mothIdx);
   tr->SetBranchAddress("dauIdx", dauIdx);
   tr->SetBranchAddress("dauLen", dauLen);
   tr->SetBranchAddress("mcenergy", mcenergy);
-}
-
-void BDtaunuMcReader::ClearColumnValues() {
-  BDtaunuReader::ClearColumnValues();
-  mcLen = 0;
-
-  McB1.bflavor = bdtaunu::kUndefinedBFlavor;
-  McB1.mc_idx = -1;
-  McB1.b_mctype = bdtaunu::kUndefinedBMcType;
-  McB1.tau_mctype = bdtaunu::kUndefinedTauMcType;
-  McB1.dtau_max_photon_energy = -1;
-
-  McB2.bflavor = bdtaunu::kUndefinedBFlavor;
-  McB2.mc_idx = -1;
-  McB2.b_mctype = bdtaunu::kUndefinedBMcType;
-  McB2.tau_mctype = bdtaunu::kUndefinedTauMcType;
-  McB2.dtau_max_photon_energy = -1;
-
-  g_mc.clear();
-  mc_vertex_map.clear();
-}
-
-// Find the MC B meson's from the mcLund array of the event. Determine
-// its flavor along the way. 
-void BDtaunuMcReader::FindBMesons() {
-
-  // Scan through the mcLund array. Store B mesons that are found.
-  for (int i = 0; i < mcLen; i++) {
-    if (abs(mcLund[i]) == lundIdMap["B0"] ||
-        abs(mcLund[i]) == lundIdMap["B+"] ) {
-
-      // Error if we find more than two B mesons. 
-      assert(!(McB1.mc_idx > -1 && McB2.mc_idx > -1));
-
-      if (McB1.mc_idx == -1) {
-        McB1.mc_idx = i;
-        if (abs(mcLund[i]) == lundIdMap["B0"]) {
-          McB1.bflavor = bdtaunu::kB0;
-        } else {
-          McB1.bflavor = bdtaunu::kBc;
-        }
-      } else {
-        McB2.mc_idx = i;
-        if (abs(mcLund[i]) == lundIdMap["B0"]) {
-          McB2.bflavor = bdtaunu::kB0;
-        } else {
-          McB2.bflavor = bdtaunu::kBc;
-        }    
-      }
-    }
-  }
-
-  // Error if we find only one B meson. 
-  assert(
-      (McB1.mc_idx == -1 && McB2.mc_idx == -1) ||
-      (McB1.mc_idx > -1 && McB2.mc_idx > -1)
-  );
-
-  // Error if the two B meson's flavor don't agree. Note the condition
-  // is also true for continuum events. 
-  assert(McB1.bflavor == McB2.bflavor);
 
 }
 
-// Determine the B MC types. 
-void BDtaunuMcReader::DetermineBMcType(McBMeson &mcB) {
-
-  // Return if no MC B mesons are in the event. 
-  if (mcB.mc_idx == -1) {
-    mcB.b_mctype = bdtaunu::kCont;
-    return;
-  }
-
-  // Count the number of daughters of the first generation B daughters
-  // and record the lundId of relevant daughters.
-  int n_daughters = 0;
-
-  int n_ell, n_nu;
-  n_ell = n_nu = 0;
-
-  int n_d, n_dstar, n_dstarstar, n_dstrange, n_pi, n_other;
-  n_d = n_dstar = n_dstarstar = n_dstrange = n_pi = n_other = 0;
-
-  int ell_lund, nu_lund;
-  ell_lund = nu_lund = 0;
-  int ell_mcidx = -1;
-
-  // Scan through the entire first generation daughter. 
-  int begin_dauIdx = dauIdx[mcB.mc_idx];
-  int end_dauIdx = begin_dauIdx + dauLen[mcB.mc_idx];
-  for (int i = begin_dauIdx; i < end_dauIdx; i++) {
-
-    // Note down max photon energy, but don't use it to classify B type.
-    if (mcLund[i] == lundIdMap["gamma"]) {
-      if (mcenergy[i] > mcB.dtau_max_photon_energy) {
-        mcB.dtau_max_photon_energy = mcenergy[i];
-      }
-      continue;
-    }
-    n_daughters += 1;
-
-    // Determine the identity of each daughter. 
-    if (std::binary_search(ell.begin(), ell.end(), abs(mcLund[i]))) {
-      n_ell += 1;
-      ell_lund = abs(mcLund[i]);
-      ell_mcidx = i;
-    } else if (std::binary_search(nu.begin(), nu.end(), abs(mcLund[i]))) {
-      n_nu += 1;
-      nu_lund = abs(mcLund[i]);
-    } else if (std::binary_search(dmeson.begin(), dmeson.end(), abs(mcLund[i]))) {
-      n_d += 1;
-    } else if (std::binary_search(dstar.begin(), dstar.end(), abs(mcLund[i]))) {
-      n_dstar += 1;
-    } else if (std::binary_search(dstarstar.begin(), dstarstar.end(), abs(mcLund[i]))) {
-      n_dstarstar += 1;
-    } else if (std::binary_search(dstrange.begin(), dstrange.end(), abs(mcLund[i]))) {
-      n_dstrange += 1;
-    } else if (std::binary_search(pion.begin(), pion.end(), abs(mcLund[i]))) {
-      n_pi += 1;
-    } else {
-      n_other += 1;
-    }
-  }
-
-  // Determine the B MC Type. See bdtaunu_definitions.h for the
-  // definitions. 
-  if (n_ell == 1 && n_nu == 1) {
-    if (n_d + n_dstar + n_dstarstar == 0) {
-      mcB.b_mctype = bdtaunu::kCharmless_SL;
-    } else if (n_dstarstar > 0) {
-      mcB.b_mctype = bdtaunu::kDstarstar_res;
-    } else {
-      assert(n_d + n_dstar == 1);
-      if (n_daughters == 3) {
-        if (ell_lund == abs(lundIdMap["tau-"])) {
-          if (n_d == 1) {
-            mcB.b_mctype = bdtaunu::kDtau;
-          } else {
-            mcB.b_mctype = bdtaunu::kDstartau;
-          }
-        } else {
-          if (n_d == 1) {
-            mcB.b_mctype = bdtaunu::kDl;
-          } else {
-            mcB.b_mctype = bdtaunu::kDstarl;
-          }
-        }
-      } else {
-        mcB.b_mctype = bdtaunu::kDstarstar_nonres;
-      }
-    }
-  } else if (n_ell == 0 && n_nu == 0) {
-    int nD = n_d + n_dstar + n_dstarstar + n_dstrange;
-    if (nD == 0)  {
-      mcB.b_mctype = bdtaunu::k0Charm_Had;
-    } else if (nD == 1) {
-      mcB.b_mctype = bdtaunu::k1Charm_Had;
-    } else if (nD == 2) {
-      mcB.b_mctype = bdtaunu::k2Charm_Had;
-    } else {
-      mcB.b_mctype = bdtaunu::kUndefinedBMcType;
-    }
-  } else {
-    mcB.b_mctype = bdtaunu::kUndefinedBMcType;
-  }
-
-  // Determine tau mc type
-  if (mcB.b_mctype == bdtaunu::kDtau || mcB.b_mctype == bdtaunu::kDstartau) {
-    mcB.tau_mctype = bdtaunu::ktau_h_mc;
-
-    int begin_dauIdx = dauIdx[ell_mcidx];
-    int end_dauIdx = begin_dauIdx + dauLen[ell_mcidx];
-    for (int i = begin_dauIdx; i < end_dauIdx; i++) {
-      if (abs(mcLund[i]) == abs(lundIdMap["e-"])) {
-        mcB.tau_mctype = bdtaunu::ktau_e_mc;
-        break;
-      } else if (abs(mcLund[i]) == abs(lundIdMap["mu-"])) {
-        mcB.tau_mctype = bdtaunu::ktau_mu_mc;
-        break;
-      } else if (abs(mcLund[i]) == abs(lundIdMap["K-"])) {
-        mcB.tau_mctype = bdtaunu::ktau_k_mc;
-        break;
-      } else {
-        continue;
-      }
-    }
-  } else {
-    mcB.tau_mctype = bdtaunu::kUndefinedTauMcType;
-  }
+void BDtaunuMcReader::ClearBuffer() {
+  mcLen = -999;
+  continuum = false;
+  b1_mctype = static_cast<int>(McBTypeCatalogue::BType::null);
+  b2_mctype = static_cast<int>(McBTypeCatalogue::BType::null);
 }
 
-// Main method that fills MC information.
-void BDtaunuMcReader::FillMCInformation() {
-
-  // Find the MC B mesons. 
-  FindBMesons();
-
-  // Determine B MC types. 
-  DetermineBMcType(McB1);
-  DetermineBMcType(McB2);
-
+void BDtaunuMcReader::DeleteBuffer() {
+  delete[] mcLund;
+  delete[] mothIdx;
+  delete[] dauIdx;
+  delete[] dauLen;
+  delete[] mcenergy;
 }
 
 // get the next event
 int BDtaunuMcReader::next_record() {
-  ClearColumnValues();
-
-  int next_record_idx = BDtaunuReader::next_record();
-  if (next_record_idx > -1) {
-    if (!IsMaxCandidateExceeded()) {
-      ConstructMcGraph();
-      FillMCInformation();
-    }
-  }
-
-  return next_record_idx;
-}
 
-void BDtaunuMcReader::ConstructMcGraph() {
+  ClearBuffer();
 
-  property_map<McGraph, vertex_mc_index_t>::type mc_index = get(vertex_mc_index, g_mc);
-  property_map<McGraph, vertex_lund_id_t>::type lund_id = get(vertex_lund_id, g_mc);
+  // read next event from root ntuple into the buffer
+  int reader_status = BDtaunuReader::next_record();
 
-  Vertex u, v;
-  std::map<int, Vertex>::iterator pos;
-  bool inserted;
-
-  for (int i = 0; i < mcLen; i++) {
-    int u_idx = i;
-    boost::tie(pos, inserted) = mc_vertex_map.insert(std::make_pair(u_idx, Vertex()));
-    if (inserted) {
-      u = add_vertex(g_mc);
-      mc_index[u] = u_idx;
-      lund_id[u] = mcLund[i];
-      mc_vertex_map[u_idx] = u;
+  // proceed only when event is not in an error state
+  if (reader_status == bdtaunu::kReadSucceeded) {
+    if (!is_max_mc_exceeded()) {
+      mc_graph_manager.construct_graph();
+      mc_graph_manager.analyze_graph();
+      FillMCInformation();
     } else {
-      u = mc_vertex_map[u_idx];
-    }
-
-    int first_dau_idx = dauIdx[i];
-    for (int j = 0; j < dauLen[i]; j++) {
-
-      int v_idx = first_dau_idx + j;
-      boost::tie(pos, inserted) = mc_vertex_map.insert(std::make_pair(v_idx, Vertex()));
-      if (inserted) {
-        v = add_vertex(g_mc);
-        mc_index[v] = v_idx;
-        lund_id[v] = mcLund[v_idx];
-        mc_vertex_map[v_idx] = v;
-      } else {
-        v = mc_vertex_map[v_idx];
-      }
-
-      add_edge(u, v, g_mc);
+      reader_status = bdtaunu::kMaxMcParticlesExceeded;
     }
   }
+  
+  return reader_status;
+}
+
+void BDtaunuMcReader::FillMCInformation() {
+  if (mc_graph_manager.get_mcY()) 
+    continuum = !(mc_graph_manager.get_mcY()->isBBbar);
+  if (mc_graph_manager.get_mcB1()) 
+    b1_mctype = mc_graph_manager.get_mcB1()->mc_type;
+  if (mc_graph_manager.get_mcB2()) 
+    b2_mctype = mc_graph_manager.get_mcB2()->mc_type;
+  return;
 }
diff --git a/create_database/BDtaunuMcReader.h b/create_database/BDtaunuMcReader.h
index 81b72d3..59ad48f 100644
--- a/create_database/BDtaunuMcReader.h
+++ b/create_database/BDtaunuMcReader.h
@@ -1,106 +1,30 @@
 #ifndef __BDTAUNUMCREADER_H__
 #define __BDTAUNUMCREADER_H__
 
-#include "BDtaunuReader.h"
-#include "bdtaunu_definitions.h"
-
+#include <iostream>
 #include <map>
 #include <string>
 #include <vector>
 
-#include <boost/graph/adjacency_list.hpp>
-
-namespace boost {
-  enum vertex_mc_index_t { vertex_mc_index };
-
-  BOOST_INSTALL_PROPERTY(vertex, mc_index);
-}
-
+#include "bdtaunu_definitions.h"
+#include "BDtaunuReader.h"
+#include "McGraphManager.h"
 
 //! Reads Monte Carlo ntuples and computes truth information. 
 /*! In addition to computing detector response data, this class also
  * computes data related to Monte Carlo truth. */
 class BDtaunuMcReader : public BDtaunuReader {
 
-  typedef boost::adjacency_list<
-  boost::vecS, boost::vecS, boost::directedS,
-  boost::property<boost::vertex_mc_index_t, int,
-  boost::property<boost::vertex_lund_id_t, int>>,
-  boost::property<boost::edge_index_t, int>
-  > McGraph;
-
-  typedef typename boost::graph_traits<McGraph>::vertex_descriptor Vertex;
-
-  McGraph g_mc;
-  std::map<int, Vertex> mc_vertex_map;
-
-  void ConstructMcGraph();
-
-  const static double min_photon_energy;
-  const static int max_mc_length;
-
-  public:
-    McGraph get_mc_graph() const { return g_mc; }
-
-  private:
-    const static std::vector<int> ell;
-    const static std::vector<int> nu;
-    const static std::vector<int> dmeson;
-    const static std::vector<int> dstar;
-    const static std::vector<int> dstarstar;
-    const static std::vector<int> dstrange;
-    const static std::vector<int> pion;
-
-    static std::vector<int> build_ell();
-    static std::vector<int> build_nu();
-    static std::vector<int> build_dmeson();
-    static std::vector<int> build_dstar();
-    static std::vector<int> build_dstarstar();
-    static std::vector<int> build_dstrange();
-    static std::vector<int> build_pion();
-
-	protected:
-    int mcLen;
-    int *mcLund; 
-    int *mothIdx;
-    int *dauIdx, *dauLen;
-    float *mcenergy;
-
-    struct McBMeson {
-      int bflavor;
-      int mc_idx;
-      int b_mctype;
-      int tau_mctype;
-      double dtau_max_photon_energy;
-      McBMeson() :
-        bflavor(bdtaunu::kUndefinedBFlavor), 
-        mc_idx(-1), 
-        b_mctype(bdtaunu::kUndefinedBMcType),
-        tau_mctype(bdtaunu::kUndefinedTauMcType), 
-        dtau_max_photon_energy(-1) {};
-    } McB1, McB2;
-
-  private:
-    void FillMCInformation();
-    void FindBMesons();
-    void DetermineBMcType(McBMeson &mcB);
-
-	protected:
-    virtual void Initialize();
-    virtual void SetBranchAddress();
-    virtual void ClearColumnValues();
+  friend class McGraphManager;
 
   public: 
 
-    //! Default construction undefined. 
-    BDtaunuMcReader();
-
-    //! Construct with specified root file name. TTree name assumed to be "ntp1".
-    BDtaunuMcReader(const char *root_fname);
-
-    //! Construct with specified root file name and TTree name.
-    BDtaunuMcReader(const char *root_fname, const char *root_trname);
-    virtual ~BDtaunuMcReader();
+    // Constructors
+    BDtaunuMcReader() = delete;
+    BDtaunuMcReader(const char *root_fname, const char *root_trname = "ntp1");
+    BDtaunuMcReader(const BDtaunuMcReader&) = delete;
+    BDtaunuMcReader &operator=(const BDtaunuMcReader&) = delete;
+    ~BDtaunuMcReader();
 
     //! Read in the next event. 
     /*! Returns an integer that indexes the event number. Returns -1
@@ -110,29 +34,42 @@ class BDtaunuMcReader : public BDtaunuReader {
      * with the event. */
     virtual int next_record();
 
-    //! B MC type of first truth B. 
-    /*! Returns an int that corresponds to the #BMcType enum in */
-    int get_b1_mctype() const { return McB1.b_mctype; }
+    // Data Accessors
+    bool is_continuum() const { return continuum; }
+    int get_b1_mctype() const { return b1_mctype; }
+    int get_b2_mctype() const { return b2_mctype; }
+
+    // Printer
+    void print_mc_graph(std::ostream &os) const { mc_graph_manager.print(os); }
 
-    //! B MC type of second truth B. 
-    /*! Returns an int that corresponds to the #BMcType enum in */
-    int get_b2_mctype() const { return McB2.b_mctype; }
+  private:
+    const static int max_mc_length;
 
-    //! tau MC type of first truth B 
-    /*! Returns an int corresponding to #TauMcType enum. */
-    int get_b1_tau_mctype() const { return McB1.tau_mctype; }
+	private:
+    int mcLen;
+    int *mcLund; 
+    int *mothIdx;
+    int *dauIdx;
+    int *dauLen;
+    float *mcenergy;
 
-    //! tau MC type of second truth B 
-    /*! Returns an int corresponding to #TauMcType enum. */
-    int get_b2_tau_mctype() const { return McB2.tau_mctype; }
+    bool continuum;
+    int b1_mctype, b2_mctype;
+
+    // Constructor helpers
+    void AllocateBuffer();
+    void DeleteBuffer();
+    void ClearBuffer();
+
+    // Reader status helpers
+    bool is_max_mc_exceeded() { return (mcLen > max_mc_length) ? true : false; }
+
+    // Mutator helpers
+    void FillMCInformation();
 
-    //! Energy of the highest energy photon of first truth B. 
-    /*! Returns -1 if no photons exist. */
-    double get_b1_dtau_max_photon_energy() const { return McB1.dtau_max_photon_energy; }
+    // MC graph helpers
+    McGraphManager mc_graph_manager;
 
-    //! Energy of the highest energy photon of second truth B. 
-    /*! Returns -1 if no photons exist. */
-    double get_b2_dtau_max_photon_energy() const { return McB2.dtau_max_photon_energy; }
 };
 
 #endif
diff --git a/create_database/BDtaunuReader.cc b/create_database/BDtaunuReader.cc
index 3a5f255..5bf5d4f 100644
--- a/create_database/BDtaunuReader.cc
+++ b/create_database/BDtaunuReader.cc
@@ -12,6 +12,7 @@
 #include "utilities/helpers.h"
 #include "RootReader.h"
 #include "BDtaunuReader.h"
+#include "BDtaunuReaderStatus.h"
 #include "UpsilonCandidate.h"
 #include "RecoGraphManager.h"
 
@@ -310,15 +311,20 @@ int BDtaunuReader::next_record() {
   ClearBuffer();
 
   // read next event from root ntuple into the buffer
-  int next_record_idx = RootReader::next_record();
+  int reader_status = RootReader::next_record();
 
   // proceed only when event is not in an error state
-  if ((next_record_idx > -1) && (!is_max_reco_exceeded())) {
-    reco_graph_manager.analyze();
-    FillUpsilonCandidates();
-  }
-
-  return next_record_idx;
+  if (reader_status == bdtaunu::kReadSucceeded) {
+    if (!is_max_reco_exceeded()) {
+      reco_graph_manager.construct_graph();
+      reco_graph_manager.analyze_graph();
+      FillUpsilonCandidates();
+    } else {
+      reader_status = bdtaunu::kMaxRecoCandExceeded;
+    }
+  } 
+  
+  return reader_status;
 }
 
 bool BDtaunuReader::is_max_reco_exceeded() const {
diff --git a/create_database/BDtaunuReader.h b/create_database/BDtaunuReader.h
index f836908..aa61cc7 100644
--- a/create_database/BDtaunuReader.h
+++ b/create_database/BDtaunuReader.h
@@ -21,7 +21,7 @@ class BDtaunuReader : public RootReader {
     BDtaunuReader(const char *root_fname, const char *root_trname = "ntp1");
     BDtaunuReader(const BDtaunuReader&) = delete;
     BDtaunuReader &operator=(const BDtaunuReader&) = delete;
-    ~BDtaunuReader();
+    virtual ~BDtaunuReader();
 
     //! Read in the next event. 
     /*! Returns an integer that indexes the event number. Returns -1
@@ -29,10 +29,7 @@ class BDtaunuReader : public RootReader {
      *
      * Calling this automatically computes all features associated
      * with the event that the analysis is interested in. */
-    int next_record();
-
-    //! Returns true if event exceeds maximum allowable reco particles
-    bool is_max_reco_exceeded() const;
+    virtual int next_record();
 
     // Data Accessors
 
@@ -54,9 +51,10 @@ class BDtaunuReader : public RootReader {
     // Printer
     void print_reco_graph(std::ostream &os) const { reco_graph_manager.print(os); }
 
-  private: 
+  protected:
     static const std::map<int, std::string> lund_to_name;
 
+  private: 
     static const int maximum_Y_candidates;
     static const int maximum_B_candidates;
     static const int maximum_D_candidates;
@@ -107,9 +105,12 @@ class BDtaunuReader : public RootReader {
   private:
 
     // Constructor helpers
-    void AllocateBuffer();
-    void DeleteBuffer();
-    void ClearBuffer();
+    virtual void AllocateBuffer();
+    virtual void DeleteBuffer();
+    virtual void ClearBuffer();
+
+    // Reader status helpers
+    bool is_max_reco_exceeded() const;
 
     // Mutator helpers
     void FillUpsilonCandidates();
diff --git a/create_database/BDtaunuReaderStatus.h b/create_database/BDtaunuReaderStatus.h
new file mode 100644
index 0000000..4fbf757
--- /dev/null
+++ b/create_database/BDtaunuReaderStatus.h
@@ -0,0 +1,15 @@
+#ifndef __BDTAUNUREADERSTATUS_H_
+#define __BDTAUNUREADERSTATUS_H_
+
+namespace bdtaunu {
+
+enum BDtaunuReaderStatus {
+  kReadSucceeded = 0,
+  kEOF = 1,
+  kMaxRecoCandExceeded = 2,
+  kMaxMcParticlesExceeded = 3,
+};
+
+}
+
+#endif
diff --git a/create_database/GraphDef.h b/create_database/GraphDef.h
new file mode 100644
index 0000000..e38c9c7
--- /dev/null
+++ b/create_database/GraphDef.h
@@ -0,0 +1,61 @@
+#ifndef __GRAPHDEF_H_
+#define __GRAPHDEF_H_
+
+#include <boost/graph/adjacency_list.hpp>
+
+namespace boost {
+  enum vertex_lund_id_t { vertex_lund_id };
+  BOOST_INSTALL_PROPERTY(vertex, lund_id);
+}
+
+namespace boost {
+  enum vertex_reco_index_t { vertex_reco_index };
+  enum vertex_block_index_t { vertex_block_index };
+  BOOST_INSTALL_PROPERTY(vertex, reco_index);
+  BOOST_INSTALL_PROPERTY(vertex, block_index);
+}
+
+namespace RecoGraph {
+
+typedef boost::adjacency_list<
+  boost::vecS, boost::vecS, boost::directedS,
+  boost::property<boost::vertex_reco_index_t, int,
+  boost::property<boost::vertex_block_index_t, int,
+  boost::property<boost::vertex_lund_id_t, int>>>,
+  boost::property<boost::edge_index_t, int>
+> Graph;
+
+typedef typename boost::graph_traits<Graph>::vertex_descriptor Vertex;
+
+typedef typename boost::graph_traits<Graph>::adjacency_iterator AdjacencyIterator;
+
+typedef typename boost::property_map<Graph, boost::vertex_lund_id_t>::type LundIdPropertyMap;
+typedef typename boost::property_map<Graph, boost::vertex_reco_index_t>::type RecoIndexPropertyMap;
+typedef typename boost::property_map<Graph, boost::vertex_block_index_t>::type BlockIndexPropertyMap;
+
+}
+
+namespace boost {
+  enum vertex_mc_index_t { vertex_mc_index };
+  BOOST_INSTALL_PROPERTY(vertex, mc_index);
+}
+
+namespace McGraph {
+
+typedef boost::adjacency_list<
+  boost::vecS, boost::vecS, boost::directedS,
+  boost::property<boost::vertex_mc_index_t, int,
+  boost::property<boost::vertex_lund_id_t, int>>,
+  boost::property<boost::edge_index_t, int>
+> Graph;
+
+typedef typename boost::graph_traits<Graph>::vertex_descriptor Vertex;
+
+typedef typename boost::graph_traits<Graph>::adjacency_iterator AdjacencyIterator;
+
+typedef typename boost::property_map<Graph, boost::vertex_lund_id_t>::type LundIdPropertyMap;
+typedef typename boost::property_map<Graph, boost::vertex_mc_index_t>::type McIndexPropertyMap;
+
+}
+
+#endif
diff --git a/create_database/GraphManager.h b/create_database/GraphManager.h
new file mode 100644
index 0000000..359f13d
--- /dev/null
+++ b/create_database/GraphManager.h
@@ -0,0 +1,17 @@
+#ifndef __GRAPHMANAGER_H_
+#define __GRAPHMANAGER_H_
+
+#include <iostream>
+
+class GraphManager {
+
+  public:
+    virtual ~GraphManager() {}
+    virtual void construct_graph() = 0;
+    virtual void analyze_graph() = 0;
+    virtual void print(std::ostream &os) const = 0;
+    virtual void clear() = 0;
+
+};
+
+#endif
diff --git a/create_database/Makefile b/create_database/Makefile
index b9bb579..652f206 100644
--- a/create_database/Makefile
+++ b/create_database/Makefile
@@ -17,15 +17,17 @@ CXXFLAGS += -D__PDT_FILE_PATHNAME='"$(PDT_FILE_PATHNAME)"'
 OS = $(shell uname)
 
 OBJECTS = UpsilonCandidate.o \
-					RootReader.o BDtaunuReader.o \
-					RecoIndexer.o \
-					RecoGraphVisitors.o RecoGraphManager.o
+					RootReader.o BDtaunuReader.o BDtaunuMcReader.o \
+					RecoIndexer.o RecoGraphVisitors.o RecoGraphManager.o \
+					McGraphManager.o McBTypeCatalogue.o McGraphVisitors.o
 
 LIBNAME = libcreatedb.so
 
-HEADER_ONLY = RecoGraph.h RecoParticles.h BDtaunuGraphWriter.h
+HEADER_ONLY = GraphDef.h Particles.h \
+							BDtaunuGraphWriter.h BDtaunuReaderStatus.h \
+							GraphManager.h Trie.h
 
-BINARIES = test2 test #test1 test_bgl
+BINARIES = test_mcreader1 test_mcreader2 #trie_test mcbtrie_test test_reader1 test_reader2 
 
 all : CXXFLAGS += -O3
 all : lib $(BINARIES)
diff --git a/create_database/McBTypeCatalogue.cc b/create_database/McBTypeCatalogue.cc
new file mode 100644
index 0000000..532639f
--- /dev/null
+++ b/create_database/McBTypeCatalogue.cc
@@ -0,0 +1,246 @@
+#include <vector>
+#include <algorithm>
+
+#include "McBTypeCatalogue.h"
+
+McBTypeCatalogue::BType 
+McBTypeCatalogue::search_catalogue(std::vector<int> lund_list) const {
+  std::vector<Alphabet> word;
+  bool hasX = false;
+  for (auto l : lund_list) {
+    Alphabet a = LundToAlphabet(l);
+    if (a == Alphabet::I) {
+      continue;
+    } else if (a == Alphabet::X) {
+      hasX = true;
+    } else {
+      word.push_back(a);
+    }
+  }
+
+  std::sort(word.begin(), word.end());
+  if (hasX) word.push_back(Alphabet::X);
+  word.push_back(Alphabet::null);
+
+  return trie.search_word(word);
+}
+
+void McBTypeCatalogue::RegisterDecays() {
+
+  // nu_tau branch
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, 
+      Alphabet::null
+    }, BType::SL);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, Alphabet::D,
+      Alphabet::null
+    }, BType::Dtau);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, Alphabet::Dstar,
+      Alphabet::null
+    }, BType::Dstartau);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, Alphabet::Dstarstar,
+      Alphabet::null
+    }, BType::Dstarstar_res);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, Alphabet::X,
+      Alphabet::null
+    }, BType::SL);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, 
+      Alphabet::D, Alphabet::X,
+      Alphabet::null
+    }, BType::Dstarstar_nonres);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, 
+      Alphabet::Dstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Dstarstar_nonres);
+
+  trie.add_word({
+      Alphabet::nu_tau, Alphabet::tau, 
+      Alphabet::Dstarstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Dstarstar_res);
+
+  // nu_ell branch
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, 
+      Alphabet::null
+    }, BType::SL);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, Alphabet::D,
+      Alphabet::null
+    }, BType::Dl);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, Alphabet::Dstar,
+      Alphabet::null
+    }, BType::Dstarl);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, Alphabet::Dstarstar,
+      Alphabet::null
+    }, BType::Dstarstar_res);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, Alphabet::X,
+      Alphabet::null
+    }, BType::SL);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, 
+      Alphabet::D, Alphabet::X,
+      Alphabet::null
+    }, BType::Dstarstar_nonres);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, 
+      Alphabet::Dstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Dstarstar_nonres);
+
+  trie.add_word({
+      Alphabet::nu_ell, Alphabet::ell, 
+      Alphabet::Dstarstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Dstarstar_res);
+
+  // hadron branch
+  trie.add_word({
+      Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::X, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstar, Alphabet::X, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstarstar, Alphabet::X, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::D, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::Dstar, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::Dstarstar, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstar, Alphabet::Dstar, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstar, Alphabet::Dstarstar, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstarstar, Alphabet::Dstarstar, 
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::D, Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::Dstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::D, Alphabet::Dstarstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstar, Alphabet::Dstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstar, Alphabet::Dstarstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+  trie.add_word({
+      Alphabet::Dstarstar, Alphabet::Dstarstar, Alphabet::X,
+      Alphabet::null
+    }, BType::Had);
+
+}
+
+McBTypeCatalogue::Alphabet McBTypeCatalogue::LundToAlphabet(int lund) const {
+  switch (abs(lund)) {
+    case 12: // nu_e
+    case 14: // nu_mu
+      return Alphabet::nu_ell;
+      break;
+    case 16: // nu_tau
+      return Alphabet::nu_tau;
+      break;
+    case 11: // e
+    case 13: // mu
+      return Alphabet::ell;
+      break;
+    case 15: // tau
+      return Alphabet::tau;
+      break;
+    case 411: // D+
+    case 421: // D0
+      return Alphabet::D;
+      break;
+    case 413: // D*+
+    case 423: // D*0
+      return Alphabet::Dstar;
+      break;
+    case 10411: // D_0*+
+    case 10421: // D_0*0
+    case 10413: // D_1+
+    case 10423: // D_10
+    case 415: // D_2*+
+    case 425: // D_2*0
+    case 20413: // D'_1+
+    case 20423: // D'_10
+    case 30411: // D(2S)+
+    case 30421: // D(2S)0
+    case 30413: // D*(2S)+
+    case 30423: // D*(2S)0
+      return Alphabet::Dstarstar;
+      break;
+    case 22: 
+      return Alphabet::I;
+      break;
+    default: 
+      return Alphabet::X;
+      break;
+  }
+}
diff --git a/create_database/McBTypeCatalogue.h b/create_database/McBTypeCatalogue.h
new file mode 100644
index 0000000..ee47b76
--- /dev/null
+++ b/create_database/McBTypeCatalogue.h
@@ -0,0 +1,34 @@
+#ifndef _MCBTYPECATALOGUE_H_
+#define _MCBTYPECATALOGUE_H_
+
+#include <vector>
+
+#include "Trie.h"
+
+class McBTypeCatalogue {
+
+  public:
+    enum class Alphabet { 
+      null = 0, nu_ell, nu_tau, ell, tau, 
+      D, Dstar, Dstarstar, X, I,
+    };
+
+    enum class BType { 
+      null = 0, Dtau, Dstartau, Dl, Dstarl, 
+      Dstarstar_res, Dstarstar_nonres, SL, Had, 
+    };
+
+    McBTypeCatalogue() { RegisterDecays(); }
+    ~McBTypeCatalogue() {};
+
+    BType search_catalogue(std::vector<int>) const;
+
+  private:
+    void RegisterDecays();
+    Alphabet LundToAlphabet(int lund) const;
+
+    bdtaunu::Trie<Alphabet, BType, Alphabet::null, BType::null> trie;
+
+};
+
+#endif
diff --git a/create_database/McGraphManager.cc b/create_database/McGraphManager.cc
new file mode 100644
index 0000000..af57ed7
--- /dev/null
+++ b/create_database/McGraphManager.cc
@@ -0,0 +1,134 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include <map>
+
+#include <boost/graph/depth_first_search.hpp>
+#include <boost/graph/graphviz.hpp>
+
+#include "bdtaunu_definitions.h"
+#include "GraphDef.h"
+#include "Particles.h"
+#include "BDtaunuMcReader.h"
+#include "McGraphManager.h"
+#include "BDtaunuGraphWriter.h"
+
+using namespace boost;
+using namespace McGraph;
+
+const std::set<int> McGraphManager::final_state_particles = { 
+  bdtaunu::eLund, bdtaunu::muLund, bdtaunu::piLund, bdtaunu::KLund,
+  bdtaunu::gammaLund, bdtaunu::protonLund, bdtaunu::neutronLund,
+};
+
+bool is_final_state_particle(int lund) {
+  return 
+    (McGraphManager::final_state_particles.find(std::abs(lund)) == 
+     McGraphManager::final_state_particles.end()) ? false : true;
+}
+
+McGraphManager::McGraphManager() : reader(nullptr) { 
+}
+
+McGraphManager::McGraphManager(BDtaunuMcReader *_reader) : reader(_reader) { 
+}
+
+void McGraphManager::clear() {
+  ClearGraph();
+  ClearAnalysis();
+}
+
+void McGraphManager::ClearGraph() {
+  mc_vertex_map.clear();
+  g.clear();
+}
+
+void McGraphManager::construct_graph() {
+
+  ClearGraph();
+
+  assert(reader != nullptr);
+
+  McIndexPropertyMap mc_index = get(vertex_mc_index, g);
+  LundIdPropertyMap lund_id = get(vertex_lund_id, g);
+
+  Vertex u, v;
+  std::map<int, Vertex>::iterator pos;
+  bool inserted;
+
+  for (int i = 0; i < reader->mcLen; i++) {
+
+    int u_idx = i;
+    tie(pos, inserted) = mc_vertex_map.insert(std::make_pair(u_idx, Vertex()));
+    if (inserted) {
+      u = add_vertex(g);
+      mc_index[u] = u_idx;
+      lund_id[u] = (reader->mcLund)[i];
+      mc_vertex_map[u_idx] = u;
+    } else {
+      u = mc_vertex_map[u_idx];
+    }
+
+    int first_dau_idx = (reader->dauIdx)[i];
+    for (int j = 0; j < (reader->dauLen)[i]; j++) {
+
+      int v_idx = first_dau_idx + j;
+      tie(pos, inserted) = mc_vertex_map.insert(std::make_pair(v_idx, Vertex()));
+      if (inserted) {
+        v = add_vertex(g);
+        mc_index[v] = v_idx;
+        lund_id[v] = (reader->mcLund)[v_idx];
+        mc_vertex_map[v_idx] = v;
+      } else {
+        v = mc_vertex_map[v_idx];
+      }
+
+      add_edge(u, v, g);
+    }
+  }
+}
+
+void McGraphManager::ClearAnalysis() {
+  Y_map.clear();
+  B_map.clear();
+}
+
+void McGraphManager::analyze_graph() {
+  ClearAnalysis();
+  depth_first_search(g, visitor(McGraphDfsVisitor(this)));
+  assert(Y_map.size() == 0 || Y_map.size() == 1);
+  assert(B_map.size() == 0 || B_map.size() == 2);
+  return;
+}
+
+const McY* McGraphManager::get_mcY() const {
+  if (Y_map.size()) {
+    return &Y_map.begin()->second;
+  } else {
+    return nullptr;
+  }
+}
+
+const McB* McGraphManager::get_mcB1() const {
+  if (B_map.size()) {
+    return &B_map.begin()->second;
+  } else {
+    return nullptr;
+  }
+}
+
+const McB* McGraphManager::get_mcB2() const {
+  if (B_map.size()) {
+    return &(++B_map.begin())->second;
+  } else {
+    return nullptr;
+  }
+}
+
+void McGraphManager::print(std::ostream &os) const {
+  boost::write_graphviz(
+      os, g, 
+      make_graph_writer(g, BDtaunuMcReader::lund_to_name, 
+                        get(vertex_lund_id, g),
+                        get(vertex_mc_index, g)));
+}
diff --git a/create_database/McGraphManager.h b/create_database/McGraphManager.h
new file mode 100644
index 0000000..dc4ca61
--- /dev/null
+++ b/create_database/McGraphManager.h
@@ -0,0 +1,56 @@
+#ifndef __MCGRAPHMANAGER_H_
+#define __MCGRAPHMANAGER_H_
+
+#include <iostream>
+#include <set>
+#include <map>
+
+#include "GraphManager.h"
+#include "GraphDef.h"
+#include "Particles.h"
+#include "McGraphVisitors.h"
+
+class BDtaunuMcReader;
+
+class McGraphManager : public GraphManager {
+
+  friend bool is_final_state_particle(int lund);
+  friend class McGraphDfsVisitor;
+
+  public:
+    McGraphManager();
+    McGraphManager(BDtaunuMcReader*);
+    McGraphManager(const McGraphManager&) = default;
+    McGraphManager &operator=(const McGraphManager&) = default;
+    ~McGraphManager() {};
+
+    void construct_graph();
+    void analyze_graph();
+    void print(std::ostream &os) const;
+    void clear();
+
+    const McY* get_mcY() const;
+    const McB* get_mcB1() const;
+    const McB* get_mcB2() const;
+
+  private:
+    static const std::set<int> final_state_particles;
+
+  private:
+    BDtaunuMcReader *reader;
+    McGraph::Graph g;
+
+    // graph construction
+    std::map<int, McGraph::Vertex> mc_vertex_map;
+    void ClearGraph();
+
+    // graph analysis
+    std::map<McGraph::Vertex, McY> Y_map;
+    std::map<McGraph::Vertex, McB> B_map;
+    void ClearAnalysis();
+
+};
+
+bool is_final_state_particle(int lund);
+
+#endif
diff --git a/create_database/McGraphVisitors.cc b/create_database/McGraphVisitors.cc
new file mode 100644
index 0000000..e9a1b06
--- /dev/null
+++ b/create_database/McGraphVisitors.cc
@@ -0,0 +1,80 @@
+#include <cmath>
+#include <cassert>
+#include <vector>
+
+#include "GraphDef.h"
+#include "Particles.h"
+#include "McGraphVisitors.h"
+#include "McGraphManager.h"
+#include "McBTypeCatalogue.h"
+
+using namespace boost;
+using namespace McGraph;
+
+McGraphDfsVisitor::McGraphDfsVisitor(McGraphManager *_manager) 
+  : manager(_manager) {
+  lund_map = get(vertex_lund_id, manager->g);
+}
+
+const McBTypeCatalogue McGraphDfsVisitor::mcB_catalogue = McBTypeCatalogue();
+
+void McGraphDfsVisitor::finish_vertex(Vertex u, const Graph &g) {
+  int lund = std::abs(get(lund_map, u));
+  switch (lund) {
+    case bdtaunu::UpsilonLund:
+      AnalyzeY(u, g);
+      break;
+    case bdtaunu::B0Lund:
+    case bdtaunu::BcLund:
+      AnalyzeB(u, g);
+      break;
+    default:
+      return;
+  }
+  return;
+}
+
+
+void McGraphDfsVisitor::AnalyzeY(const Vertex &u, const Graph &g) {
+
+  McY mcY;
+
+  AdjacencyIterator ai, ai_end;
+  for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) {
+
+    int lund = abs(get(lund_map, *ai));
+    switch (lund) {
+      case bdtaunu::B0Lund:
+      case bdtaunu::BcLund:
+        (mcY.B1 == nullptr) ? 
+          (mcY.B1 = &(manager->B_map)[*ai]) : 
+          (mcY.B2 = &(manager->B_map)[*ai]);
+        break;
+      default:
+        mcY.isBBbar = false;
+        return;
+    }
+  }
+
+  (manager->Y_map).insert(std::make_pair(u, mcY));
+}
+
+void McGraphDfsVisitor::AnalyzeB(const Vertex &u, const Graph &g) {
+
+  McB mcB;
+
+  if (abs(get(lund_map, u)) == bdtaunu::B0Lund) {
+    mcB.flavor = bdtaunu::kB0;
+  } else {
+    mcB.flavor = bdtaunu::kBc;
+  }
+
+  std::vector<int> daulund_list;
+  AdjacencyIterator ai, ai_end;
+  for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) {
+    daulund_list.push_back(get(lund_map, *ai));
+  }
+  mcB.mc_type = static_cast<int>(mcB_catalogue.search_catalogue(daulund_list));
+
+  (manager->B_map).insert(std::make_pair(u, mcB));
+}
diff --git a/create_database/McGraphVisitors.h b/create_database/McGraphVisitors.h
new file mode 100644
index 0000000..8162ae8
--- /dev/null
+++ b/create_database/McGraphVisitors.h
@@ -0,0 +1,31 @@
+#ifndef __MCGRAPHVISITOR_H__
+#define __MCGRAPHVISITOR_H__
+
+#include <boost/graph/depth_first_search.hpp>
+
+#include "GraphDef.h"
+#include "McBTypeCatalogue.h"
+
+class McGraphManager;
+
+class McGraphDfsVisitor : public boost::default_dfs_visitor {
+
+  public:
+    McGraphDfsVisitor() = default;
+    McGraphDfsVisitor(McGraphManager*);
+    ~McGraphDfsVisitor() {};
+
+    void finish_vertex(McGraph::Vertex u, const McGraph::Graph &g);
+
+  private:
+    static const McBTypeCatalogue mcB_catalogue;
+
+  private:
+    McGraphManager *manager = nullptr;
+    McGraph::LundIdPropertyMap lund_map;
+
+    void AnalyzeY(const McGraph::Vertex &u, const McGraph::Graph &g);
+    void AnalyzeB(const McGraph::Vertex &u, const McGraph::Graph &g);
+};
+
+#endif
diff --git a/create_database/RecoParticles.h b/create_database/Particles.h
similarity index 60%
rename from create_database/RecoParticles.h
rename to create_database/Particles.h
index f917a14..e85487c 100644
--- a/create_database/RecoParticles.h
+++ b/create_database/Particles.h
@@ -1,8 +1,8 @@
-#ifndef _RECOPARTICLES_H_
-#define _RECOPARTICLES_H_
+#ifndef _PARTICLES_H_
+#define _PARTICLES_H_
 
 #include "bdtaunu_definitions.h"
-#include "RecoGraph.h"
+#include "GraphDef.h"
 
 struct RecoB;
 struct RecoD;
@@ -16,7 +16,7 @@ struct RecoY {
 
 struct RecoB {
   RecoB() = default;
-  int flavor;
+  int flavor = bdtaunu::kUndefinedBFlavor;
   RecoD *D = nullptr;
   RecoLepton *Lepton = nullptr;
 };
@@ -33,4 +33,19 @@ struct RecoLepton {
   int tau_mode = bdtaunu::kUndefinedTauMode;
 };
 
+struct McB;
+
+struct McY {
+  McY() = default;
+  bool isBBbar = true;
+  McB *B1 = nullptr;
+  McB *B2 = nullptr;
+};
+
+struct McB {
+  McB() = default;
+  int flavor = bdtaunu::kUndefinedBFlavor;
+  int mc_type = bdtaunu::kUndefinedBMcType;
+};
+
 #endif
diff --git a/create_database/RecoGraphManager.cc b/create_database/RecoGraphManager.cc
index 22668ae..0951d5a 100644
--- a/create_database/RecoGraphManager.cc
+++ b/create_database/RecoGraphManager.cc
@@ -6,8 +6,9 @@
 #include <boost/graph/depth_first_search.hpp>
 #include <boost/graph/graphviz.hpp>
 
+#include "GraphDef.h"
+#include "Particles.h"
 #include "BDtaunuReader.h"
-#include "RecoGraph.h"
 #include "RecoGraphManager.h"
 #include "BDtaunuGraphWriter.h"
 
@@ -25,41 +26,13 @@ void RecoGraphManager::clear() {
   ClearAnalysis();
 }
 
-void RecoGraphManager::analyze() {
-  ConstructGraph();
-  AnalyzeGraph();
-  return;
-}
-
-const RecoY* RecoGraphManager::get_recoY(int i) const { 
-
-  std::map<int, RecoGraph::Vertex>::const_iterator reco_vertex_it;
-  reco_vertex_it = reco_vertex_map.find(reco_indexer(bdtaunu::UpsilonLund, i));
-  assert(reco_vertex_it != reco_vertex_map.end());
-
-  std::map<Vertex, RecoY>::const_iterator y_it;
-  y_it = Y_map.find(reco_vertex_it->second);
-  assert(y_it != Y_map.end());
-
-  return &y_it->second;
-}
-
-void RecoGraphManager::print(std::ostream &os) const {
-  boost::write_graphviz(
-      os, g, 
-      make_graph_writer(g, BDtaunuReader::lund_to_name, 
-                        get(vertex_lund_id, g),
-                        get(vertex_reco_index, g)));
-}
-
-
 void RecoGraphManager::ClearGraph() {
   reco_vertex_map.clear();
   reco_indexer.clear();
   g.clear();
 }
 
-void RecoGraphManager::ConstructGraph() {
+void RecoGraphManager::construct_graph() {
 
   ClearGraph();
 
@@ -87,8 +60,10 @@ void RecoGraphManager::ConstructGraph() {
   AddCandidates(reader->nh, reader->hLund, hdauIdx, hdauLund);
   AddCandidates(reader->nl, reader->lLund, ldauIdx, ldauLund);
 
+  return;
 }
 
+
 void RecoGraphManager::ClearAnalysis() {
   Y_map.clear();
   B_map.clear();
@@ -96,14 +71,34 @@ void RecoGraphManager::ClearAnalysis() {
   Lepton_map.clear();
 }
 
-void RecoGraphManager::AnalyzeGraph() {
-
+void RecoGraphManager::analyze_graph() {
   ClearAnalysis();
   depth_first_search(g, visitor(RecoGraphDfsVisitor(this)));
-
   return;
 }
 
+const RecoY* RecoGraphManager::get_recoY(int i) const { 
+
+  std::map<int, RecoGraph::Vertex>::const_iterator reco_vertex_it;
+  reco_vertex_it = reco_vertex_map.find(reco_indexer(bdtaunu::UpsilonLund, i));
+  assert(reco_vertex_it != reco_vertex_map.end());
+
+  std::map<Vertex, RecoY>::const_iterator y_it;
+  y_it = Y_map.find(reco_vertex_it->second);
+  assert(y_it != Y_map.end());
+
+  return &y_it->second;
+}
+
+void RecoGraphManager::print(std::ostream &os) const {
+  boost::write_graphviz(
+      os, g, 
+      make_graph_writer(g, BDtaunuReader::lund_to_name, 
+                        get(vertex_lund_id, g),
+                        get(vertex_reco_index, g)));
+}
+
+
 void RecoGraphManager::AddCandidates(
     int nCand,
     int *CandLund,
diff --git a/create_database/RecoGraphManager.h b/create_database/RecoGraphManager.h
index b843668..2ec4a86 100644
--- a/create_database/RecoGraphManager.h
+++ b/create_database/RecoGraphManager.h
@@ -5,14 +5,15 @@
 #include <vector>
 #include <map>
 
-#include "RecoGraph.h"
+#include "GraphManager.h"
+#include "GraphDef.h"
+#include "Particles.h"
 #include "RecoIndexer.h"
 #include "RecoGraphVisitors.h"
-#include "RecoParticles.h"
 
 class BDtaunuReader;
 
-class RecoGraphManager {
+class RecoGraphManager : public GraphManager {
 
   friend class RecoGraphDfsVisitor;
 
@@ -23,14 +24,14 @@ class RecoGraphManager {
     RecoGraphManager &operator=(const RecoGraphManager&) = default;
     ~RecoGraphManager() {};
 
-    void analyze();
+    void construct_graph();
+    void analyze_graph();
+    void print(std::ostream &os) const;
     void clear();
 
     int get_reco_index(int lund, int i) const { return reco_indexer(lund, i); }
     const RecoY* get_recoY(int i) const;
 
-    void print(std::ostream &os) const;
-
   private:
     BDtaunuReader *reader;
     RecoGraph::Graph g;
@@ -38,19 +39,17 @@ class RecoGraphManager {
     // graph construction
     RecoIndexer reco_indexer;
     std::map<int, RecoGraph::Vertex> reco_vertex_map;
+    void ClearGraph();
     void AddCandidates(
         int nCand, int *CandLund,
         std::vector<int*> &CandDauIdx, 
         std::vector<int*> &CandDauLund);
-    void ConstructGraph();
-    void ClearGraph();
 
     // graph analysis
     std::map<RecoGraph::Vertex, RecoY> Y_map;
     std::map<RecoGraph::Vertex, RecoB> B_map;
     std::map<RecoGraph::Vertex, RecoD> D_map;
     std::map<RecoGraph::Vertex, RecoLepton> Lepton_map;
-    void AnalyzeGraph();
     void ClearAnalysis();
 };
 
diff --git a/create_database/RecoGraphVisitors.cc b/create_database/RecoGraphVisitors.cc
index 29bfc18..2586cef 100644
--- a/create_database/RecoGraphVisitors.cc
+++ b/create_database/RecoGraphVisitors.cc
@@ -1,10 +1,10 @@
 #include <cmath>
 #include <cassert>
 
-#include "RecoGraph.h"
+#include "GraphDef.h"
+#include "Particles.h"
 #include "RecoGraphVisitors.h"
 #include "RecoGraphManager.h"
-#include "RecoParticles.h"
 
 using namespace boost;
 using namespace RecoGraph;
diff --git a/create_database/RecoGraphVisitors.h b/create_database/RecoGraphVisitors.h
index 4f0dbed..f14d615 100644
--- a/create_database/RecoGraphVisitors.h
+++ b/create_database/RecoGraphVisitors.h
@@ -3,7 +3,7 @@
 
 #include <boost/graph/depth_first_search.hpp>
 
-#include "RecoGraph.h"
+#include "GraphDef.h"
 
 class RecoGraphManager;
 
diff --git a/create_database/RootReader.cc b/create_database/RootReader.cc
index 2db7aab..adbd773 100644
--- a/create_database/RootReader.cc
+++ b/create_database/RootReader.cc
@@ -3,6 +3,7 @@
 #include <iostream>
 
 #include "RootReader.h"
+#include "BDtaunuReaderStatus.h"
 
 RootReader::RootReader() : 
   tfile(0), tr(0), 
@@ -53,8 +54,8 @@ void RootReader::PrepareTreeFile(const char *root_fname,
 int RootReader::next_record() {
   if (record_index < total_records) {
     tr->GetEntry(record_index++);
-    return record_index;
+    return bdtaunu::kReadSucceeded;
   } else {
-    return -1;
+    return bdtaunu::kEOF;
   }
 }
diff --git a/create_database/Trie.h b/create_database/Trie.h
new file mode 100644
index 0000000..f3700ff
--- /dev/null
+++ b/create_database/Trie.h
@@ -0,0 +1,97 @@
+#ifndef _BDTAUNU_TRIE_H_
+#define _BDTAUNU_TRIE_H_
+
+#include <iostream>
+#include <vector>
+#include <algorithm>
+
+namespace bdtaunu {
+
+template<typename alphaT, typename valueT,
+         alphaT alphaNull, valueT valueNull>
+class Trie {
+
+  private:
+    struct node {
+      alphaT alpha;
+      valueT value;
+      std::vector<node*> links;
+
+      node() : alpha(alphaT()), value(valueT()) {};
+      node(alphaT a, valueT v = valueNull) : alpha(a), value(v) {};
+      ~node() { for (auto l : links) delete l; }
+    };
+
+    node *root;
+    node *CopyNode(const node*);
+
+  public:
+    Trie() { root = new node; }
+    ~Trie() { delete root; }
+    Trie(const Trie &t) { root = CopyNode(t.root); }
+    Trie &operator=(const Trie&);
+
+    void add_word(std::vector<alphaT> word, valueT value);
+    valueT search_word(std::vector<alphaT> word) const;
+};
+
+
+template<typename alphaT, typename valueT,
+         alphaT alphaNull, valueT valueNull>
+typename Trie<alphaT, valueT, alphaNull, valueNull>::node*
+ Trie<alphaT, valueT, alphaNull, valueNull>::CopyNode(const node *t) {
+  node *n = new node(t->alpha, t->value);
+  for (auto l : t->links) (n->links).push_back(CopyNode(l));
+  return n;
+}
+
+template<typename alphaT, typename valueT,
+         alphaT alphaNull, valueT valueNull>
+Trie<alphaT, valueT, alphaNull, valueNull>&
+Trie<alphaT, valueT, alphaNull, valueNull>::operator=(const Trie &target) {
+  if (this != &target) {
+    delete root;
+    root = CopyNode(target.root);
+  }
+  return *this;
+}
+
+
+template<typename alphaT, typename valueT,
+         alphaT alphaNull, valueT valueNull>
+void Trie<alphaT, valueT, alphaNull, valueNull>::add_word(
+    std::vector<alphaT> word, valueT value) {
+  node *p = root;
+  for (auto alpha : word) {
+    auto e_it = find_if((p->links).begin(), (p->links).end(), 
+              [alpha] (const node *n) { return (n->alpha == alpha); });
+
+    if (e_it == (p->links).end()) {
+      node *n = new node(alpha);
+      if (alpha == alphaNull) n->value = value;
+      (p->links).push_back(n);
+      p = n;
+    } else {
+      p = *e_it;
+    }
+  }
+}
+
+
+template<typename alphaT, typename valueT,
+         alphaT alphaNull, valueT valueNull>
+valueT Trie<alphaT, valueT, alphaNull, valueNull>::search_word(
+    std::vector<alphaT> word) const { 
+  node *p = root;
+  for (auto alpha : word) {
+    auto e_it = find_if((p->links).begin(), (p->links).end(), 
+              [alpha] (const node *n) { return (n->alpha == alpha); });
+    if (e_it == (p->links.end())) return valueNull; 
+    p = *e_it;
+  }
+  return p->value;
+}
+
+}
+
+#endif
diff --git a/include/bdtaunu_definitions.h b/include/bdtaunu_definitions.h
index efb4764..4bae042 100644
--- a/include/bdtaunu_definitions.h
+++ b/include/bdtaunu_definitions.h
@@ -25,6 +25,8 @@ const int piLund = 211;
 const int eLund = 11;
 const int muLund = 13;
 const int gammaLund = 22;
+const int protonLund = 2212;
+const int neutronLund = 2112;
 
 //! B meson flavors.
 enum BFlavor {