diff --git a/create_database/BDtaunuGraphWriter.h b/create_database/BDtaunuGraphWriter.h
new file mode 100644
index 0000000..2abe80c
--- /dev/null
+++ b/create_database/BDtaunuGraphWriter.h
@@ -0,0 +1,63 @@
+#ifndef __BDTAUNUGRAPHWRITER_H_
+#define __BDTAUNUGRAPHWRITER_H_
+
+#include <iostream> 
+#include <fstream> 
+
+#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 {
+
+  private:
+    typedef typename boost::graph_traits<GraphT>::vertex_descriptor VertexT;
+
+  public:
+    BDtaunuGraphWriter(
+        const GraphT &_graph, 
+        const LundMap &_lund_map,
+        const LundPM &_lund_pm, 
+        const IdxPM &_idx_pm) : 
+        g(_graph), 
+        lund_map(_lund_map), 
+        lund_pm(_lund_pm), 
+        idx_pm(_idx_pm) {}
+
+    void operator()(std::ostream &out, const VertexT &v) const;
+
+  private:
+    GraphT g;
+    LundMap lund_map;
+    LundPM lund_pm;
+    IdxPM idx_pm;
+};
+
+template <typename GraphT, typename LundMap,
+          typename LundPM, typename IdxPM>
+void BDtaunuGraphWriter<GraphT, LundMap, LundPM, IdxPM>::operator() (
+    std::ostream &out, const VertexT &v) const {
+
+  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\"]";
+  } else {
+    out << "[label=\"" << idx_pm[v] << ": " << lund_pm[v] << "\",color=\"red\"]";
+  }
+
+}
+
+template <typename GraphT, typename LundMap,
+          typename LundPM, typename IdxPM>
+BDtaunuGraphWriter<GraphT, LundMap, LundPM, IdxPM>
+make_graph_writer(const GraphT &g, const LundMap lund_map, 
+                  const LundPM lund_pm, const IdxPM idx_pm) {
+
+  return BDtaunuGraphWriter<GraphT, LundMap, LundPM, IdxPM>(g, lund_map, lund_pm, idx_pm);
+}
+
+#endif
diff --git a/create_database/BDtaunuReader.cc b/create_database/BDtaunuReader.cc
index 3cc9ce4..3a5f255 100644
--- a/create_database/BDtaunuReader.cc
+++ b/create_database/BDtaunuReader.cc
@@ -1,5 +1,3 @@
-#include <TTree.h>
-
 #include <iostream>
 #include <fstream>
 #include <string>
@@ -11,17 +9,14 @@
 #include <cassert>
 
 #include "bdtaunu_definitions.h"
+#include "utilities/helpers.h"
 #include "RootReader.h"
 #include "BDtaunuReader.h"
-#include "RecoIndexer.h"
 #include "UpsilonCandidate.h"
-#include "RecoGraphVisitors.h"
-#include "utilities/helpers.h"
-
-#include <boost/graph/depth_first_search.hpp>
-#include <boost/graph/breadth_first_search.hpp>
+#include "RecoGraphManager.h"
 
-using namespace boost;
+// Lund to particle name map needed for printing.
+const std::map<int, std::string> BDtaunuReader::lund_to_name = bdtaunu::LundToNameMap();
 
 // The maximum number of candidates allowed in an event. This should
 // be consistent with the number set in BtaTupleMaker. 
@@ -33,25 +28,23 @@ const int BDtaunuReader::maximum_h_candidates = 100;
 const int BDtaunuReader::maximum_l_candidates = 100;
 const int BDtaunuReader::maximum_gamma_candidates = 100;
 
-std::map<std::string, int> BDtaunuReader::lundIdMap = bdtaunu::NameToLundMap();
+BDtaunuReader::BDtaunuReader(
+    const char *root_fname, 
+    const char *root_trname) : RootReader(root_fname, root_trname) {
 
-BDtaunuReader::BDtaunuReader() : RootReader() {
-  Initialize();
-}
+  AllocateBuffer();
+  ClearBuffer();
 
-BDtaunuReader::BDtaunuReader(const char *root_fname) : RootReader(root_fname) {
-  Initialize();
+  reco_graph_manager = RecoGraphManager(this);
 }
 
-BDtaunuReader::BDtaunuReader(
-    const char *root_fname, 
-    const char *root_trname) : RootReader(root_fname, root_trname) {
-  Initialize();
+BDtaunuReader::~BDtaunuReader() {
+  DeleteBuffer();
 }
 
 
 // All constructors call this function to initialize the class members.
-void BDtaunuReader::Initialize() {
+void BDtaunuReader::AllocateBuffer() {
 
   // Allocate space to read in arrays from ntuples. 
   YBPairMmissPrime2 = new float[maximum_Y_candidates];
@@ -127,28 +120,7 @@ void BDtaunuReader::Initialize() {
   ld2Lund = new int[maximum_l_candidates];
   ld3Lund = new int[maximum_l_candidates];
 
-  YdauIdx = { Yd1Idx, Yd2Idx };
-  YdauLund = { Yd1Lund, Yd2Lund };
-  BdauIdx = { Bd1Idx, Bd2Idx, Bd3Idx, Bd4Idx };
-  BdauLund = { Bd1Lund, Bd2Lund, Bd3Lund, Bd4Lund };
-  DdauIdx = { Dd1Idx, Dd2Idx, Dd3Idx, Dd4Idx, Dd5Idx };
-  DdauLund = { Dd1Lund, Dd2Lund, Dd3Lund, Dd4Lund, Dd5Lund };
-  CdauIdx = { Cd1Idx, Cd2Idx };
-  CdauLund = { Cd1Lund, Cd2Lund };
-  hdauIdx = { hd1Idx, hd2Idx };
-  hdauLund = { hd1Lund, hd2Lund };
-  ldauIdx = { ld1Idx, ld2Idx, ld3Idx };
-  ldauLund = { ld1Lund, ld2Lund, ld3Lund };
-
   // Specify the variables where each ntuple branch should be read into. 
-  SetBranchAddress();
-  
-  // Zero out variables where ntuple data are read into. 
-  ClearColumnValues();
-}
-
-void BDtaunuReader::SetBranchAddress() {
-
   tr->SetBranchAddress("platform", &platform);
   tr->SetBranchAddress("partition", &partition);
   tr->SetBranchAddress("upperID", &upperID);
@@ -237,7 +209,7 @@ void BDtaunuReader::SetBranchAddress() {
 
 }
 
-void BDtaunuReader::ClearColumnValues() {
+void BDtaunuReader::ClearBuffer() {
   platform = -999;
   partition = -999;
   upperID = -999;
@@ -251,15 +223,10 @@ void BDtaunuReader::ClearColumnValues() {
   nh = -999;
   nl = -999;
   ngamma = -999;
-
-  eventId.clear();
   upsilon_candidates.clear();
-  reco_indexer.clear();
-  reco_vertex_map.clear();
-  g.clear();
 }
 
-BDtaunuReader::~BDtaunuReader() {
+void BDtaunuReader::DeleteBuffer() {
 
   delete[] YBPairMmissPrime2;
   delete[] YBPairEextra50;
@@ -336,162 +303,87 @@ BDtaunuReader::~BDtaunuReader() {
 
 }
 
-// Determine whether the maximum number of allowed candidates 
-// have been exceeded. Events that exceed the maximum are not 
-// well behaved and should be skipped.
-bool BDtaunuReader::IsMaxCandidateExceeded() const {
-    if ( 
-        (nY < maximum_Y_candidates) &&
-        (nB < maximum_B_candidates) &&
-        (nD < maximum_D_candidates) &&
-        (nC < maximum_C_candidates) &&
-        (nh < maximum_h_candidates) &&
-        (nl < maximum_l_candidates) &&
-        (ngamma < maximum_gamma_candidates) 
-       ) {
-      return false;
-    } else {
-      return true;
-    }
-}
-
 // Read in the next event in the ntuple. It computes all relevant
 // information as a side effect. 
 int BDtaunuReader::next_record() {
 
-  // Zero out garbage information from the previous event. 
-  ClearColumnValues();
+  ClearBuffer();
 
-  // Read in the next event in the ntuple. The variables now contain
-  // information from the new event. 
+  // read next event from root ntuple into the buffer
   int next_record_idx = RootReader::next_record();
 
-  // Compute higher level information that is not directly avaible. 
-  if ((next_record_idx > -1) && (!IsMaxCandidateExceeded())) {
-
-    eventId = std::to_string(platform) + ":" 
-              + std::to_string(partition) + ":" 
-              + std::to_string(upperID) + "/" 
-              + std::to_string(lowerID);
-
-    reco_indexer.set({nY, nB, nD, nC, nh, nl, ngamma});
-
-    ConstructRecoGraph();
-
-    // Construct the Y(4S) candidate list for this event. 
-    // Maximum number of candidates must not be exceeded.
-    FillUpsilonList();
+  // 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;
 }
 
-void BDtaunuReader::ConstructRecoGraph() {
-  AddCandidatesToGraph(nY, YLund, YdauIdx, YdauLund);
-  AddCandidatesToGraph(nB, BLund, BdauIdx, BdauLund);
-  AddCandidatesToGraph(nD, DLund, DdauIdx, DdauLund);
-  AddCandidatesToGraph(nC, CLund, CdauIdx, CdauLund);
-  AddCandidatesToGraph(nh, hLund, hdauIdx, hdauLund);
-  AddCandidatesToGraph(nl, lLund, ldauIdx, ldauLund);
-}
-
-void BDtaunuReader::AddCandidatesToGraph(
-    int nCand,
-    int *CandLund,
-    std::vector<int*> &CandDauIdx,
-    std::vector<int*> &CandDauLund) {
-
-  property_map<RecoGraph, vertex_reco_index_t>::type reco_index = get(vertex_reco_index, g);
-  property_map<RecoGraph, vertex_block_index_t>::type block_index = get(vertex_block_index, g);
-  property_map<RecoGraph, vertex_lund_id_t>::type lund_id = get(vertex_lund_id, g);
-
-  Vertex u, v;
-  std::map<int, Vertex>::iterator pos;
-  bool inserted;
-
-  for (int i = 0; i < nCand; i++) {
-
-    int u_idx = reco_indexer(CandLund[i], i);
-    boost::tie(pos, inserted) = reco_vertex_map.insert(std::make_pair(u_idx, Vertex()));
-    if (inserted) {
-      u = add_vertex(g);
-      reco_index[u] = u_idx;
-      block_index[u] = i;
-      lund_id[u] = CandLund[i];
-      reco_vertex_map[u_idx] = u;
+bool BDtaunuReader::is_max_reco_exceeded() const {
+    if ( 
+        (nY < maximum_Y_candidates) &&
+        (nB < maximum_B_candidates) &&
+        (nD < maximum_D_candidates) &&
+        (nC < maximum_C_candidates) &&
+        (nh < maximum_h_candidates) &&
+        (nl < maximum_l_candidates) &&
+        (ngamma < maximum_gamma_candidates) 
+       ) {
+      return false;
     } else {
-      u = reco_vertex_map[u_idx];
+      return true;
     }
+}
 
-    for (std::vector<int*>::size_type j = 0; j < CandDauIdx.size(); j++) {
-      if (CandDauIdx[j][i] == -1) {
-        break;
-      } else {
-
-        int v_idx = reco_indexer(CandDauLund[j][i], CandDauIdx[j][i]);
-        boost::tie(pos, inserted) = reco_vertex_map.insert(std::make_pair(v_idx, Vertex()));
-        if (inserted) {
-          v = add_vertex(g);
-          reco_index[v] = v_idx;
-          block_index[v] = CandDauIdx[j][i];
-          lund_id[v] = CandDauLund[j][i];
-          reco_vertex_map[v_idx] = v;
-        } else {
-          v = reco_vertex_map[v_idx];
-        }
-
-        add_edge(u, v, g);
-      }
-    }
-  }
+std::string BDtaunuReader::get_eventId() const { 
+  return std::to_string(platform) 
+         + ":" + std::to_string(partition) 
+         + ":" + std::to_string(upperID) 
+         + "/" + std::to_string(lowerID);
 }
 
+void BDtaunuReader::FillUpsilonCandidates() {
+
+  for (int i = 0; i < nY; i++) {
+
+    UpsilonCandidate ups;
+
+    ups.set_eventId(get_eventId());
+    ups.set_block_index(i);
+    ups.set_reco_index(reco_graph_manager.get_reco_index(bdtaunu::UpsilonLund, i));
+    ups.set_bflavor(reco_graph_manager.get_recoY(i)->tagB->flavor);
+    ups.set_eextra50(YBPairEextra50[i]);
+    ups.set_mmiss_prime2(YBPairMmissPrime2[i]);
+    ups.set_cosThetaT(YBPairCosThetaT[i]);
+    ups.set_tag_lp3(YTagBlP3MagCM[i]);
+    ups.set_tag_cosBY(YTagBCosBY[i]);
+    ups.set_tag_cosThetaDl(YTagBCosThetaDlCM[i]);
+    ups.set_tag_Dmass(YTagBDMass[i]);
+    ups.set_tag_deltaM(YTagBDstarDeltaM[i]);
+    ups.set_tag_cosThetaDSoft(YTagBCosThetaDSoftCM[i]);
+    ups.set_tag_softP3MagCM(YTagBsoftP3MagCM[i]);
+    ups.set_tag_d_mode(reco_graph_manager.get_recoY(i)->tagB->D->D_mode);
+    ups.set_tag_dstar_mode(reco_graph_manager.get_recoY(i)->tagB->D->Dstar_mode);
+    ups.set_l_ePidMap(eSelectorsMap[lTrkIdx[reco_graph_manager.get_recoY(i)->tagB->Lepton->l_block_idx]]);
+    ups.set_l_muPidMap(muSelectorsMap[lTrkIdx[reco_graph_manager.get_recoY(i)->tagB->Lepton->l_block_idx]]);
+    ups.set_sig_hp3(YSigBhP3MagCM[i]);
+    ups.set_sig_cosBY(YSigBCosBY[i]);
+    ups.set_sig_cosThetaDtau(YSigBCosThetaDtauCM[i]);
+    ups.set_sig_vtxB(YSigBVtxProbB[i]);
+    ups.set_sig_Dmass(YSigBDMass[i]);
+    ups.set_sig_deltaM(YSigBDstarDeltaM[i]);
+    ups.set_sig_cosThetaDSoft(YSigBCosThetaDSoftCM[i]);
+    ups.set_sig_softP3MagCM(YSigBsoftP3MagCM[i]);
+    ups.set_sig_hmass(YSigBhMass[i]);
+    ups.set_sig_vtxh(YSigBVtxProbh[i]);
+    ups.set_sig_d_mode(reco_graph_manager.get_recoY(i)->sigB->D->D_mode);
+    ups.set_sig_dstar_mode(reco_graph_manager.get_recoY(i)->sigB->D->Dstar_mode);
+    ups.set_sig_tau_mode(reco_graph_manager.get_recoY(i)->sigB->Lepton->tau_mode);
+    ups.set_h_ePidMap(eSelectorsMap[hTrkIdx[reco_graph_manager.get_recoY(i)->sigB->Lepton->pi_block_idx]]);
+    ups.set_h_muPidMap(muSelectorsMap[hTrkIdx[reco_graph_manager.get_recoY(i)->sigB->Lepton->pi_block_idx]]);
 
-// This function takes the decay tree information read from the ntuple
-// (which are read into the arrays new'd in Initialize()), and creates
-// the UpsilonList of the event. It constructs an UpsilonCandidate
-// object for each candidate and fills it with correctly computed
-// higher level features.
-void BDtaunuReader::FillUpsilonList() {
-
-  auto lund_map = get(vertex_lund_id, g);
-  std::vector<YDecayProperties> upsilon_decays;
-  depth_first_search(g, visitor(reco_graph_dfs_visitor(lund_map, upsilon_decays)));
-
-  auto block_idx_map = get(vertex_block_index, g);
-  auto reco_idx_map = get(vertex_reco_index, g);
-  for (auto y : upsilon_decays) {
-    int cand_idx = get(block_idx_map, y.Y);
-    int reco_idx = get(reco_idx_map, y.Y);
-    int l_ePidMap = eSelectorsMap[lTrkIdx[get(block_idx_map, y.l)]];
-    int l_muPidMap = muSelectorsMap[lTrkIdx[get(block_idx_map, y.l)]];
-    int h_ePidMap = eSelectorsMap[hTrkIdx[get(block_idx_map, y.tau_pi)]];
-    int h_muPidMap = muSelectorsMap[hTrkIdx[get(block_idx_map, y.tau_pi)]];
-
-    UpsilonCandidate ups(eventId, cand_idx, reco_idx, 
-        YBPairEextra50[cand_idx], YBPairMmissPrime2[cand_idx], 
-        YTagBlP3MagCM[cand_idx], YSigBhP3MagCM[cand_idx], 
-        YTagBCosBY[cand_idx], YSigBCosBY[cand_idx], 
-        YTagBCosThetaDlCM[cand_idx], YSigBCosThetaDtauCM[cand_idx], 
-        YSigBVtxProbB[cand_idx], 
-        YBPairCosThetaT[cand_idx], 
-        YTagBDMass[cand_idx], YTagBDstarDeltaM[cand_idx], 
-        YTagBCosThetaDSoftCM[cand_idx], YTagBsoftP3MagCM[cand_idx],
-        YSigBDMass[cand_idx], YSigBDstarDeltaM[cand_idx], 
-        YSigBCosThetaDSoftCM[cand_idx], YSigBsoftP3MagCM[cand_idx],
-        YSigBhMass[cand_idx], YSigBVtxProbh[cand_idx], 
-        y.bflavor,
-        y.tag_dstar_mode, y.tag_d_mode, 
-        y.sig_dstar_mode, y.sig_d_mode, 
-        y.tau_mode, 
-        l_ePidMap, l_muPidMap, 
-        h_ePidMap, h_muPidMap); 
     upsilon_candidates.push_back(ups);
   }
 }
-
-BDtaunuReader::RecoGraph BDtaunuReader::get_reco_subgraph(int reco_idx) {
-  RecoGraph subg;
-  breadth_first_search(g, reco_vertex_map[reco_idx], visitor(reco_graph_bfs_visitor(g, subg)));
-  return subg;
-}
diff --git a/create_database/BDtaunuReader.h b/create_database/BDtaunuReader.h
index b05bce8..f836908 100644
--- a/create_database/BDtaunuReader.h
+++ b/create_database/BDtaunuReader.h
@@ -1,67 +1,61 @@
 #ifndef __BDTAUNUREADER_H__
 #define __BDTAUNUREADER_H__
 
-#include <TTree.h>
-
 #include <string>
 #include <vector>
 #include <map>
+#include <iostream>
 
 #include "RootReader.h"
 #include "UpsilonCandidate.h"
-#include "RecoIndexer.h"
-
-#include <boost/graph/adjacency_list.hpp>
-
-namespace boost {
-  enum vertex_reco_index_t { vertex_reco_index };
-  enum vertex_block_index_t { vertex_block_index };
-  enum vertex_lund_id_t { vertex_lund_id };
-
-  BOOST_INSTALL_PROPERTY(vertex, reco_index);
-  BOOST_INSTALL_PROPERTY(vertex, block_index);
-  BOOST_INSTALL_PROPERTY(vertex, lund_id);
-}
-
-template <typename T> class RecoGraphDfsVisitor;
-class RecoGraphBfsVisitor;
-
-//! Reads ntuples produced by BtaTupleMaker and computes data related to detector response. 
-/*! This class is responsible for assembling all data that can be
- * derived from the detector's response to an event. More
- * specifically, it contains the following types of information: 
- * - Macro level event infromation. Examples:
- *   - Fox-Wolfram moment. 
- *   - nTrks. 
- * - \f$\Upsilon(4S)\f$ candidate information. Examples:
- *   - \f$E_{extra}\f$.
- *   - \f$m_D\f$ of \f$B_{sig}\f$.
- *   - Reconstructed type. See #CandType enum. 
- *
- * All other event information that does not belong in the above
- * categories must be computed in a class that derives from this one.
- * For example, Monte Carlo information should be computed in a
- * subclass. */
+#include "RecoGraphManager.h"
+
 class BDtaunuReader : public RootReader {
 
-  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>
-  > RecoGraph;
+  friend class RecoGraphManager;
+
+  public: 
+
+    // Constructors
+    BDtaunuReader() = delete;
+    BDtaunuReader(const char *root_fname, const char *root_trname = "ntp1");
+    BDtaunuReader(const BDtaunuReader&) = delete;
+    BDtaunuReader &operator=(const BDtaunuReader&) = delete;
+    ~BDtaunuReader();
+
+    //! Read in the next event. 
+    /*! Returns an integer that indexes the event number. Returns -1
+     * when all events have been read. 
+     *
+     * 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;
+
+    // Data Accessors
 
-  typedef typename boost::graph_traits<RecoGraph>::vertex_descriptor Vertex;
+    //! Babar event Id. 
+    std::string get_eventId() const;
+
+    //! nTRK defined in BtaTupleMaker. 
+    int get_nTrk() const { return nTrk; }
+
+    //! Number of \f$\Upsilon(4S)\f$ candidates associated with this event. 
+    int get_nY() const { return nY; }
 
-  struct YDecayProperties;
-  template <typename T> friend class RecoGraphDfsVisitor;
-  template <typename T> friend RecoGraphDfsVisitor<T> reco_graph_dfs_visitor(T, std::vector<YDecayProperties>&);
-  friend class RecoGraphBfsVisitor;
-  friend RecoGraphBfsVisitor reco_graph_bfs_visitor(const RecoGraph&, RecoGraph&);
+    //! Second Fox-Wolfram moment. 
+    float get_R2All() const { return R2All; }
 
-  protected: 
-    static std::map<std::string, int> lundIdMap;
+    //! Return list of \f$\Upsilon(4S)\f$ candidates in this event.
+    const std::vector<UpsilonCandidate> &get_upsilon_candidates() const { return upsilon_candidates; }
+
+    // Printer
+    void print_reco_graph(std::ostream &os) const { reco_graph_manager.print(os); }
+
+  private: 
+    static const std::map<int, std::string> lund_to_name;
 
     static const int maximum_Y_candidates;
     static const int maximum_B_candidates;
@@ -72,6 +66,7 @@ class BDtaunuReader : public RootReader {
     static const int maximum_gamma_candidates;
 
   private: 
+    // read from root tuples
     int platform, partition, upperID, lowerID;
     int nTrk;
     float R2All;
@@ -91,15 +86,6 @@ class BDtaunuReader : public RootReader {
     int *lTrkIdx, *hTrkIdx;
     int *eSelectorsMap, *muSelectorsMap, *KSelectorsMap, *piSelectorsMap;
 
-    struct YDecayProperties {
-      Vertex Y, l, tau_pi, tau_pi0;
-      int bflavor;
-      int tag_d_mode, sig_d_mode;
-      int tag_dstar_mode, sig_dstar_mode;
-      int tau_mode;
-    };
-
-  protected: 
     int nY, nB, nD, nC, nh, nl, ngamma; 
     int *YLund, *BLund, *DLund, *CLund, *hLund, *lLund, *gammaLund;
     int *Yd1Idx, *Yd2Idx;
@@ -115,73 +101,21 @@ class BDtaunuReader : public RootReader {
     int *hd1Lund, *hd2Lund;
     int *ld1Lund, *ld2Lund, *ld3Lund;
 
-    std::vector<int*> YdauIdx, YdauLund;
-    std::vector<int*> BdauIdx, BdauLund;
-    std::vector<int*> DdauIdx, DdauLund;
-    std::vector<int*> CdauIdx, CdauLund;
-    std::vector<int*> hdauIdx, hdauLund;
-    std::vector<int*> ldauIdx, ldauLund;
-
-  private: 
-    std::string eventId;
+    // computed from root tuples
     std::vector<UpsilonCandidate> upsilon_candidates;
-    RecoIndexer reco_indexer;
-    std::map<int, Vertex> reco_vertex_map;
-    RecoGraph g;
-
-    void ConstructRecoGraph();
-    void AddCandidatesToGraph(
-        int nCand, int *CandLund,
-        std::vector<int*> &CandDauIdx,
-        std::vector<int*> &CandDauLund);
-
-    void FillUpsilonList();
-
-  protected: 
-
-    virtual void Initialize();
-    virtual void SetBranchAddress();
-    virtual void ClearColumnValues();
-
-  public: 
 
-    bool IsMaxCandidateExceeded() const;
+  private:
 
-    //! Default construction undefined. 
-    BDtaunuReader();
-
-    //! Construct with specified root file name. TTree name assumed to be "ntp1".
-    BDtaunuReader(const char *root_fname);
-
-    //! Construct with specified root file name and TTree name.
-    BDtaunuReader(const char *root_fname, const char *root_trname);
-    virtual ~BDtaunuReader();
-
-    //! Read in the next event. 
-    /*! Returns an integer that indexes the event number. Returns -1
-     * when all events have been read. 
-     *
-     * Calling this automatically computes all features associated
-     * with the event that the analysis is interested in. */
-    virtual int next_record();
-
-    //! Babar event Id. 
-    std::string get_eventId() const { return eventId; }
-
-    //! nTRK defined in BtaTupleMaker. 
-    int get_nTrk() const { return nTrk; }
-
-    //! Number of \f$\Upsilon(4S)\f$ candidates associated with this event. 
-    int get_nY() const { return nY; }
-
-    //! Second Fox-Wolfram moment. 
-    float get_R2All() const { return R2All; }
+    // Constructor helpers
+    void AllocateBuffer();
+    void DeleteBuffer();
+    void ClearBuffer();
 
-    //! Return list of \f$\Upsilon(4S)\f$ candidates in this event. 
-    const std::vector<UpsilonCandidate> &get_candidate_list() const { return upsilon_candidates; }
+    // Mutator helpers
+    void FillUpsilonCandidates();
 
-    RecoGraph get_reco_graph() const { return g; }
-    RecoGraph get_reco_subgraph(int reco_index);
+    // Reco graph helpers
+    RecoGraphManager reco_graph_manager;
 
 };
 
diff --git a/create_database/Makefile b/create_database/Makefile
index 4d71fb6..b9bb579 100644
--- a/create_database/Makefile
+++ b/create_database/Makefile
@@ -16,19 +16,16 @@ CXXFLAGS += -D__PDT_FILE_PATHNAME='"$(PDT_FILE_PATHNAME)"'
 
 OS = $(shell uname)
 
-OBJECTS = RecoIndexer.o \
-          RootReader.o BDtaunuReader.o \
-          BDtaunuMcReader.o \
-          UpsilonCandidate.o \
-          SQLiteTableBuilder.o EventSQLiteTableBuilder.o McEventSQLiteTableBuilder.o \
-					CandidateSQLiteTableBuilder.o McCandidateSQLiteTableBuilder.o \
-					EventStatusSQLiteTableBuilder.o OptimalCandidateSQLiteTableBuilder.o
+OBJECTS = UpsilonCandidate.o \
+					RootReader.o BDtaunuReader.o \
+					RecoIndexer.o \
+					RecoGraphVisitors.o RecoGraphManager.o
 
 LIBNAME = libcreatedb.so
 
-HEADER_ONLY = RecoGraphVisitors.h
+HEADER_ONLY = RecoGraph.h RecoParticles.h BDtaunuGraphWriter.h
 
-BINARIES = test2 #test test1 test_bgl
+BINARIES = test2 test #test1 test_bgl
 
 all : CXXFLAGS += -O3
 all : lib $(BINARIES)
diff --git a/create_database/RecoGraph.h b/create_database/RecoGraph.h
new file mode 100644
index 0000000..a7bf280
--- /dev/null
+++ b/create_database/RecoGraph.h
@@ -0,0 +1,36 @@
+#ifndef __RECOGRAPH_H_
+#define __RECOGRAPH_H_
+
+#include <boost/graph/adjacency_list.hpp>
+
+namespace boost {
+  enum vertex_reco_index_t { vertex_reco_index };
+  enum vertex_block_index_t { vertex_block_index };
+  enum vertex_lund_id_t { vertex_lund_id };
+
+  BOOST_INSTALL_PROPERTY(vertex, reco_index);
+  BOOST_INSTALL_PROPERTY(vertex, block_index);
+  BOOST_INSTALL_PROPERTY(vertex, lund_id);
+}
+
+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;
+
+}
+
+#endif
diff --git a/create_database/RecoGraphManager.cc b/create_database/RecoGraphManager.cc
new file mode 100644
index 0000000..22668ae
--- /dev/null
+++ b/create_database/RecoGraphManager.cc
@@ -0,0 +1,156 @@
+#include <cassert>
+#include <vector>
+#include <map>
+#include <iostream>
+
+#include <boost/graph/depth_first_search.hpp>
+#include <boost/graph/graphviz.hpp>
+
+#include "BDtaunuReader.h"
+#include "RecoGraph.h"
+#include "RecoGraphManager.h"
+#include "BDtaunuGraphWriter.h"
+
+using namespace boost;
+using namespace RecoGraph;
+
+RecoGraphManager::RecoGraphManager() : reader(nullptr) { 
+}
+
+RecoGraphManager::RecoGraphManager(BDtaunuReader *_reader) : reader(_reader) { 
+}
+
+void RecoGraphManager::clear() {
+  ClearGraph();
+  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() {
+
+  ClearGraph();
+
+  assert(reader != nullptr);
+  reco_indexer.set({reader->nY, reader->nB, reader->nD, reader->nC, 
+                    reader->nh, reader->nl, reader->ngamma});
+
+  std::vector<int*> YdauIdx{ reader->Yd1Idx, reader->Yd2Idx };
+  std::vector<int*> YdauLund{ reader->Yd1Lund, reader->Yd2Lund };
+  std::vector<int*> BdauIdx{ reader->Bd1Idx, reader->Bd2Idx, reader->Bd3Idx, reader->Bd4Idx };
+  std::vector<int*> BdauLund{ reader->Bd1Lund, reader->Bd2Lund, reader->Bd3Lund, reader->Bd4Lund };
+  std::vector<int*> DdauIdx{ reader->Dd1Idx, reader->Dd2Idx, reader->Dd3Idx, reader->Dd4Idx, reader->Dd5Idx };
+  std::vector<int*> DdauLund{ reader->Dd1Lund, reader->Dd2Lund, reader->Dd3Lund, reader->Dd4Lund, reader->Dd5Lund };
+  std::vector<int*> CdauIdx{ reader->Cd1Idx, reader->Cd2Idx };
+  std::vector<int*> CdauLund{ reader->Cd1Lund, reader->Cd2Lund };
+  std::vector<int*> hdauIdx{ reader->hd1Idx, reader->hd2Idx };
+  std::vector<int*> hdauLund{ reader->hd1Lund, reader->hd2Lund };
+  std::vector<int*> ldauIdx{ reader->ld1Idx, reader->ld2Idx, reader->ld3Idx };
+  std::vector<int*> ldauLund{ reader->ld1Lund, reader->ld2Lund, reader->ld3Lund };
+
+  AddCandidates(reader->nY, reader->YLund, YdauIdx, YdauLund);
+  AddCandidates(reader->nB, reader->BLund, BdauIdx, BdauLund);
+  AddCandidates(reader->nD, reader->DLund, DdauIdx, DdauLund);
+  AddCandidates(reader->nC, reader->CLund, CdauIdx, CdauLund);
+  AddCandidates(reader->nh, reader->hLund, hdauIdx, hdauLund);
+  AddCandidates(reader->nl, reader->lLund, ldauIdx, ldauLund);
+
+}
+
+void RecoGraphManager::ClearAnalysis() {
+  Y_map.clear();
+  B_map.clear();
+  D_map.clear();
+  Lepton_map.clear();
+}
+
+void RecoGraphManager::AnalyzeGraph() {
+
+  ClearAnalysis();
+  depth_first_search(g, visitor(RecoGraphDfsVisitor(this)));
+
+  return;
+}
+
+void RecoGraphManager::AddCandidates(
+    int nCand,
+    int *CandLund,
+    std::vector<int*> &CandDauIdx,
+    std::vector<int*> &CandDauLund) {
+
+  RecoIndexPropertyMap reco_index = get(vertex_reco_index, g);
+  BlockIndexPropertyMap block_index = get(vertex_block_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 < nCand; i++) {
+
+    int u_idx = reco_indexer(CandLund[i], i);
+    tie(pos, inserted) = reco_vertex_map.insert(std::make_pair(u_idx, Vertex()));
+    if (inserted) {
+      u = add_vertex(g);
+      reco_index[u] = u_idx;
+      block_index[u] = i;
+      lund_id[u] = CandLund[i];
+      reco_vertex_map[u_idx] = u;
+    } else {
+      u = reco_vertex_map[u_idx];
+    }
+
+    for (std::vector<int*>::size_type j = 0; j < CandDauIdx.size(); j++) {
+      if (CandDauIdx[j][i] == -1) {
+        break;
+      } else {
+
+        int v_idx = reco_indexer(CandDauLund[j][i], CandDauIdx[j][i]);
+        tie(pos, inserted) = reco_vertex_map.insert(std::make_pair(v_idx, Vertex()));
+        if (inserted) {
+          v = add_vertex(g);
+          reco_index[v] = v_idx;
+          block_index[v] = CandDauIdx[j][i];
+          lund_id[v] = CandDauLund[j][i];
+          reco_vertex_map[v_idx] = v;
+        } else {
+          v = reco_vertex_map[v_idx];
+        }
+
+        add_edge(u, v, g);
+      }
+    }
+  }
+}
diff --git a/create_database/RecoGraphManager.h b/create_database/RecoGraphManager.h
new file mode 100644
index 0000000..b843668
--- /dev/null
+++ b/create_database/RecoGraphManager.h
@@ -0,0 +1,57 @@
+#ifndef __RECOGRAPHMANAGER_H_
+#define __RECOGRAPHMANAGER_H_
+
+#include <iostream>
+#include <vector>
+#include <map>
+
+#include "RecoGraph.h"
+#include "RecoIndexer.h"
+#include "RecoGraphVisitors.h"
+#include "RecoParticles.h"
+
+class BDtaunuReader;
+
+class RecoGraphManager {
+
+  friend class RecoGraphDfsVisitor;
+
+  public:
+    RecoGraphManager();
+    RecoGraphManager(BDtaunuReader*);
+    RecoGraphManager(const RecoGraphManager&) = default;
+    RecoGraphManager &operator=(const RecoGraphManager&) = default;
+    ~RecoGraphManager() {};
+
+    void analyze();
+    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;
+
+    // graph construction
+    RecoIndexer reco_indexer;
+    std::map<int, RecoGraph::Vertex> reco_vertex_map;
+    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();
+};
+
+#endif
diff --git a/create_database/RecoGraphVisitors.cc b/create_database/RecoGraphVisitors.cc
new file mode 100644
index 0000000..29bfc18
--- /dev/null
+++ b/create_database/RecoGraphVisitors.cc
@@ -0,0 +1,288 @@
+#include <cmath>
+#include <cassert>
+
+#include "RecoGraph.h"
+#include "RecoGraphVisitors.h"
+#include "RecoGraphManager.h"
+#include "RecoParticles.h"
+
+using namespace boost;
+using namespace RecoGraph;
+
+RecoGraphDfsVisitor::RecoGraphDfsVisitor(RecoGraphManager *_manager) 
+  : manager(_manager) {
+  lund_map = get(vertex_lund_id, manager->g);
+  block_idx_map = get(vertex_block_index, manager->g);
+}
+
+void RecoGraphDfsVisitor::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;
+    case bdtaunu::Dstar0Lund:
+    case bdtaunu::DstarcLund:
+      AnalyzeDstar(u, g);
+      break;
+    case bdtaunu::D0Lund:
+    case bdtaunu::DcLund:
+      AnalyzeD(u, g);
+      break;
+    case bdtaunu::piLund:
+    case bdtaunu::rhoLund:
+    case bdtaunu::eLund:
+    case bdtaunu::muLund:
+      AnalyzeLepton(u, g);
+      break;
+    default:
+      return;
+  }
+  return;
+}
+
+void RecoGraphDfsVisitor::AnalyzeD(const Vertex &u, const Graph &g) {
+
+  RecoD recoD;
+
+  int n_daughters, n_K, n_Ks, n_pi, n_pi0;
+  n_daughters = n_K = n_Ks = n_pi = n_pi0 = 0;
+
+  AdjacencyIterator ai, ai_end;
+  for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) {
+    n_daughters += 1;
+    int lund = abs(get(lund_map, *ai));
+    switch (lund) {
+      case bdtaunu::KLund:
+        n_K += 1;
+        break;
+      case bdtaunu::KSLund:
+        n_Ks += 1;
+        break;
+      case bdtaunu::piLund:
+        n_pi += 1;
+        break;
+      case bdtaunu::pi0Lund:
+        n_pi0 += 1;
+        break;
+      default:
+        assert(false);
+        return;
+    }
+  }
+
+  int mode;
+  if (std::abs(get(lund_map, u)) == bdtaunu::DcLund) {
+    if (n_daughters == 3 && n_K == 1 && n_pi == 2) {
+      mode = bdtaunu::kDc_Kpipi;
+    } else if (n_daughters == 4 && n_K == 1 && n_pi == 2 && n_pi0 == 1) {
+      mode = bdtaunu::kDc_Kpipipi0;
+    } else if (n_daughters == 2 && n_K == 1 && n_Ks == 1) {
+      mode = bdtaunu::kDc_KsK;
+    } else if (n_daughters == 2 && n_Ks == 1 && n_pi == 1) {
+      mode = bdtaunu::kDc_Kspi;
+    } else if (n_daughters == 3 && n_Ks == 1 && n_pi == 1 && n_pi0 == 1) {
+      mode = bdtaunu::kDc_Kspipi0;
+    } else if (n_daughters == 4 && n_Ks == 1 && n_pi == 3) {
+      mode = bdtaunu::kDc_Kspipipi;
+    } else if (n_daughters == 3 && n_K == 2 && n_pi == 1) {
+      mode = bdtaunu::kDc_KKpi;
+    } else {
+      mode = bdtaunu::kUndefinedDMode;
+    }
+  } else {
+    if (n_daughters == 2 && n_K == 1 && n_pi == 1) {
+      mode = bdtaunu::kD0_Kpi;
+    } else if (n_daughters == 3 && n_K == 1 && n_pi == 1 && n_pi0 == 1) {
+      mode = bdtaunu::kD0_Kpipi0;
+    } else if (n_daughters == 4 && n_K == 1 && n_pi == 3) {
+      mode = bdtaunu::kD0_Kpipipi;
+    } else if (n_daughters == 5 && n_K == 1 && n_pi == 3 && n_pi0 == 1) {
+      mode = bdtaunu::kD0_Kpipipipi0;
+    } else if (n_daughters == 3 && n_Ks == 1 && n_pi == 2) {
+      mode = bdtaunu::kD0_Kspipi;
+    } else if (n_daughters == 4 && n_Ks == 1 && n_pi == 2 && n_pi0 == 1) {
+      mode = bdtaunu::kD0_Kspipipi0;
+    } else if (n_daughters == 2 && n_Ks == 1 && n_pi0 == 1) {
+      mode = bdtaunu::kD0_Kspi0;
+    } else if (n_daughters == 2 && n_K == 2) {
+      mode = bdtaunu::kD0_KK;
+    } else {
+      mode = bdtaunu::kUndefinedDMode;
+    }
+  }
+
+  recoD.D_mode = mode;
+
+  (manager->D_map).insert(std::make_pair(u, recoD));
+}
+
+void RecoGraphDfsVisitor::AnalyzeDstar(const Vertex &u, const Graph &g) {
+
+  RecoD recoD;
+
+  int n_daughters, n_D0, n_Dc, n_pi, n_pi0, n_gamma;
+  n_daughters = n_D0 = n_Dc = n_pi = n_pi0 = n_gamma = 0;
+
+  AdjacencyIterator ai, ai_end;
+  for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) {
+    n_daughters += 1;
+    int lund = abs(get(lund_map, *ai));
+    switch (lund) {
+      case bdtaunu::D0Lund:
+        n_D0 += 1;
+        recoD.D_mode = (manager->D_map)[*ai].D_mode;
+        break;
+      case bdtaunu::DcLund:
+        n_Dc += 1;
+        recoD.D_mode = (manager->D_map)[*ai].D_mode;
+        break;
+      case bdtaunu::piLund:
+        n_pi += 1;
+        break;
+      case bdtaunu::pi0Lund:
+        n_pi0 += 1;
+        break;
+      case bdtaunu::gammaLund:
+        n_gamma += 1;
+        break;
+      default:
+        assert(false);
+        return;
+    }
+  }
+
+  int mode;
+  if (n_daughters == 2 && n_D0 == 1 && n_pi0 == 1) {
+    mode = bdtaunu::kDstar0_D0pi0;
+  } else if (n_daughters == 2 && n_D0 == 1 && n_gamma == 1) {
+    mode = bdtaunu::kDstar0_D0gamma;
+  } else if (n_daughters == 2 && n_D0 == 1 && n_pi == 1) {
+    mode = bdtaunu::kDstarc_D0pi;
+  } else if (n_daughters == 2 && n_Dc == 1 && n_pi0 == 1) {
+    mode = bdtaunu::kDstarc_Dcpi0;
+  } else if (n_daughters == 2 && n_Dc == 1 && n_gamma == 1) {
+    mode = bdtaunu::kDstarc_Dcgamma;
+  } else {
+    mode = bdtaunu::kUndefinedDstarMode;
+  }
+
+  recoD.Dstar_mode = mode;
+
+  (manager->D_map).insert(std::make_pair(u, recoD));
+
+}
+
+void RecoGraphDfsVisitor::AnalyzeLepton(const Vertex &u, const Graph &g) {
+
+  RecoLepton recoLepton;
+
+  AdjacencyIterator ai, ai_end;
+  int lund = abs(get(lund_map, u));
+  switch (lund) {
+
+    case bdtaunu::eLund:
+      recoLepton.l_block_idx = block_idx_map[u];
+      recoLepton.pi_block_idx = -1;
+      recoLepton.tau_mode = bdtaunu::ktau_e;
+
+    case bdtaunu::muLund:
+      recoLepton.l_block_idx = block_idx_map[u];
+      recoLepton.pi_block_idx = -1;
+      recoLepton.tau_mode = bdtaunu::ktau_mu;
+      break;
+
+    case bdtaunu::piLund:
+      recoLepton.l_block_idx = -1;
+      recoLepton.pi_block_idx = block_idx_map[u];
+      recoLepton.tau_mode = bdtaunu::ktau_pi;
+      break;
+
+    case bdtaunu::rhoLund:
+      recoLepton.l_block_idx = -1;
+      for (tie(ai, ai_end) = adjacent_vertices(u, g); ai != ai_end; ++ai) {
+        if (abs(get(lund_map, *ai)) == bdtaunu::piLund) {
+          recoLepton.pi_block_idx = (manager->Lepton_map)[*ai].pi_block_idx;
+          break;
+        }
+      }
+      recoLepton.tau_mode = bdtaunu::ktau_rho;
+      break;
+
+    default:
+      assert(false);
+      return;
+  }
+
+  (manager->Lepton_map).insert(std::make_pair(u, recoLepton));
+}
+
+
+void RecoGraphDfsVisitor::AnalyzeB(const Vertex &u, const Graph &g) {
+
+  RecoB recoB;
+
+  if (abs(get(lund_map, u)) == bdtaunu::B0Lund) {
+    recoB.flavor = bdtaunu::kB0;
+  } else {
+    recoB.flavor = bdtaunu::kBc;
+  }
+
+  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::D0Lund:
+      case bdtaunu::DcLund:
+      case bdtaunu::Dstar0Lund:
+      case bdtaunu::DstarcLund:
+        recoB.D = &(manager->D_map)[*ai];
+        break;
+      case bdtaunu::eLund:
+      case bdtaunu::muLund:
+      case bdtaunu::piLund:
+      case bdtaunu::rhoLund:
+        recoB.Lepton = &(manager->Lepton_map)[*ai];
+        break;
+      default:
+        assert(false);
+        return;
+    }
+  }
+
+  (manager->B_map).insert(std::make_pair(u, recoB));
+}
+
+
+
+void RecoGraphDfsVisitor::AnalyzeY(const Vertex &u, const Graph &g) {
+
+  RecoY recoY;
+
+  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:
+        if ((manager->B_map)[*ai].Lepton->l_block_idx >= 0) {
+          recoY.tagB = &(manager->B_map)[*ai];
+        } else {
+          recoY.sigB = &(manager->B_map)[*ai];
+        }
+        break;
+      default:
+        assert(false);
+        return;
+    }
+  }
+
+  (manager->Y_map).insert(std::make_pair(u, recoY));
+}
+
diff --git a/create_database/RecoGraphVisitors.h b/create_database/RecoGraphVisitors.h
index a29b4be..4f0dbed 100644
--- a/create_database/RecoGraphVisitors.h
+++ b/create_database/RecoGraphVisitors.h
@@ -1,477 +1,31 @@
 #ifndef __RECOGRAPHVISITOR_H__
 #define __RECOGRAPHVISITOR_H__
 
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include <vector>
-#include <map>
-#include <utility>
 #include <boost/graph/depth_first_search.hpp>
-#include <boost/graph/breadth_first_search.hpp>
-#include "BDtaunuReader.h"
-#include "bdtaunu_definitions.h"
 
-template <typename LundMap>
-class RecoGraphDfsVisitor : public boost::default_dfs_visitor {
-  private:
-    typedef typename BDtaunuReader::RecoGraph Graph;
-    typedef typename BDtaunuReader::Vertex Vertex;
-    typedef typename BDtaunuReader::YDecayProperties YDecayProperties;
-    typedef typename boost::graph_traits<Graph>::adjacency_iterator AdjacencyIterator;
-
-  public:
-    RecoGraphDfsVisitor(LundMap &l, std::vector<YDecayProperties> &v) : 
-      lund_map(l), upsilon_properties(v) {};
-    void finish_vertex(Vertex u, const Graph &g);
-    std::vector<YDecayProperties> get_upsilons() { return upsilon_properties; }
-
-  private:
-    LundMap lund_map;
-    std::vector<YDecayProperties> &upsilon_properties;
-
-    void ComputeY(const Vertex &u, const Graph &g);
-    void ComputeB(const Vertex &u, const Graph &g);
-    void ComputeDstar(const Vertex &u, const Graph &g);
-    void ComputeD(const Vertex &u, const Graph &g);
-    void ComputeRhoDaughters(Vertex &tau_pi_dau, Vertex &tau_pi0_dau, 
-                             const Vertex &u, const Graph &g);
-
-  private:
-    struct D {
-      int dec_mode;
-      D() : dec_mode(bdtaunu::kUndefinedDMode) {};
-      D(int mode) : dec_mode(mode) {};
-    };
-
-    struct Dstar {
-      int dec_mode;
-      Vertex D_dau;
-      Dstar() : dec_mode(bdtaunu::kUndefinedDMode) {};
-      Dstar(int mode, const Vertex &v) : dec_mode(mode), D_dau(v) {};
-    };
-
-    struct B {
-      int bflavor;
-      bool has_dstardau;
-
-      bool is_tag;
-      Vertex D_dau;
-      Vertex lep_dau;
-
-      int tau_mode;
-      Vertex tau_pi_dau;
-      Vertex tau_pi0_dau;
-
-      B() {};
-      B(int _bflavor, bool _has_dstardau,
-        bool _is_tag, const Vertex &_D, 
-        const Vertex &_lep, int _tau_mode, 
-        const Vertex &_tau_pi, const Vertex &_tau_pi0) :
-        bflavor(_bflavor), has_dstardau(_has_dstardau), 
-        is_tag(_is_tag), D_dau(_D), lep_dau(_lep), 
-        tau_mode(_tau_mode), tau_pi_dau(_tau_pi), 
-        tau_pi0_dau(_tau_pi0) {};
-    };
-
-    std::map<Vertex, D> Dmap;
-    std::map<Vertex, Dstar> Dstarmap;
-    std::map<Vertex, B> Bmap;
-
-};
-
-template <typename LundMap>
-void RecoGraphDfsVisitor<LundMap>::finish_vertex(Vertex u, const Graph &g) {
-  int lund = std::abs(boost::get(lund_map, u));
-  switch (lund) {
-    case bdtaunu::UpsilonLund:
-      ComputeY(u, g);
-      break;
-    case bdtaunu::B0Lund:
-    case bdtaunu::BcLund:
-      ComputeB(u, g);
-      break;
-    case bdtaunu::Dstar0Lund:
-    case bdtaunu::DstarcLund:
-      ComputeDstar(u, g);
-      break;
-    case bdtaunu::D0Lund:
-    case bdtaunu::DcLund:
-      ComputeD(u, g);
-      break;
-    default:
-      return;
-  }
-}
-
-template<typename LundMap>
-void RecoGraphDfsVisitor<LundMap>::ComputeY(const Vertex &u, const Graph &g) {
-
-  YDecayProperties y;
-  y.Y = u;
-
-  AdjacencyIterator ai, ai_end;
-  for (boost::tie(ai, ai_end) = adjacent_vertices(u, g);
-       ai != ai_end; ++ai) {
-
-    auto bi = Bmap.find(*ai);
-    assert(bi != Bmap.end());
-    B b = bi->second;
-
-    if (b.is_tag) {
-      y.bflavor = b.bflavor;
-      y.l = b.lep_dau;
-      if (b.has_dstardau) {
-        auto dstari = Dstarmap.find(b.D_dau);
-        assert(dstari != Dstarmap.end());
-        y.tag_dstar_mode = (dstari->second).dec_mode;
-
-        auto di = Dmap.find((dstari->second).D_dau);
-        assert(di != Dmap.end());
-        y.tag_d_mode = (di->second).dec_mode;
-      } else {
-        auto di = Dmap.find(b.D_dau);
-        assert(di != Dmap.end());
-        y.tag_dstar_mode = bdtaunu::kNoDstar;
-        y.tag_d_mode = (di->second).dec_mode;
-      }
-    } else {
-      y.tau_mode = b.tau_mode;
-      y.tau_pi = b.tau_pi_dau;
-      y.tau_pi0 = b.tau_pi0_dau;
-      if (b.has_dstardau) {
-        auto dstari = Dstarmap.find(b.D_dau);
-        assert(dstari != Dstarmap.end());
-        y.sig_dstar_mode = (dstari->second).dec_mode;
-
-        auto di = Dmap.find((dstari->second).D_dau);
-        assert(di != Dmap.end());
-        y.sig_d_mode = (di->second).dec_mode;
-      } else {
-        auto di = Dmap.find(b.D_dau);
-        assert(di != Dmap.end());
-        y.sig_dstar_mode = bdtaunu::kNoDstar;
-        y.sig_d_mode = (di->second).dec_mode;
-      }
-    }
-  }
-
-  upsilon_properties.push_back(y);
-}
-
-template<typename LundMap>
-void RecoGraphDfsVisitor<LundMap>::ComputeB(const Vertex &u, const Graph &g) {
-
-  assert(Bmap.find(u) == Bmap.end());
-
-  int bflavor;
-  if (abs(boost::get(lund_map, u)) == bdtaunu::B0Lund) {
-    bflavor = bdtaunu::kB0;
-  } else {
-    bflavor = bdtaunu::kBc;
-  }
-
-  bool is_tag, has_dstardau;
-  Vertex D_dau, lep_dau;
-  int tau_mode;
-  Vertex tau_pi_dau, tau_pi0_dau;
-  AdjacencyIterator ai, ai_end;
-  for (boost::tie(ai, ai_end) = adjacent_vertices(u, g);
-       ai != ai_end; 
-       ++ai) {
-
-    int lund = abs(boost::get(lund_map, *ai));
-    switch (lund) {
-      case bdtaunu::D0Lund:
-      case bdtaunu::DcLund:
-        has_dstardau = false;
-        D_dau = *ai;
-        break;
-      case bdtaunu::Dstar0Lund:
-      case bdtaunu::DstarcLund:
-        has_dstardau = true;
-        D_dau = *ai;
-        break;
-      case bdtaunu::eLund:
-      case bdtaunu::muLund:
-        is_tag = true;
-        lep_dau = *ai;
-        tau_mode = bdtaunu::kUndefinedTauMode;
-        break;
-      case bdtaunu::piLund:
-        is_tag = false;
-        lep_dau = *ai;
-        tau_mode = bdtaunu::ktau_pi;
-        tau_pi_dau = *ai;
-        break;
-      case bdtaunu::rhoLund:
-        is_tag = false;
-        lep_dau = *ai;
-        tau_mode = bdtaunu::ktau_rho;
-        ComputeRhoDaughters(tau_pi_dau, tau_pi0_dau, *ai, g);
-        break;
-      default:
-        assert(false);
-        return;
-    }
-  }
-
-  Bmap.insert(std::make_pair(
-        u, B(bflavor, has_dstardau, is_tag, 
-             D_dau, lep_dau, tau_mode, tau_pi_dau, tau_pi0_dau)));
+#include "RecoGraph.h"
 
-}
+class RecoGraphManager;
 
-template<typename LundMap>
-void RecoGraphDfsVisitor<LundMap>::ComputeRhoDaughters(
-    Vertex &tau_pi_dau, Vertex &tau_pi0_dau, 
-    const Vertex &u, const Graph &g) {
-  AdjacencyIterator ai, ai_end;
-  for (boost::tie(ai, ai_end) = adjacent_vertices(u, g);
-       ai != ai_end; 
-       ++ai) {
-    switch(abs(boost::get(lund_map, *ai))) {
-      case bdtaunu::piLund:
-        tau_pi_dau = *ai;
-        break;
-      case bdtaunu::pi0Lund:
-        tau_pi0_dau = *ai;
-        break;
-      default:
-        assert(false);
-    }
-  }
-}
-
-template<typename LundMap>
-void RecoGraphDfsVisitor<LundMap>::ComputeDstar(const Vertex &u, const Graph &g) {
-
-  assert(Dstarmap.find(u) == Dstarmap.end());
-
-  int n_daughters, n_D0, n_Dc, n_pi, n_pi0, n_gamma;
-  n_daughters = n_D0 = n_Dc = n_pi = n_pi0 = n_gamma = 0;
-
-  Vertex v;
-  AdjacencyIterator ai, ai_end;
-  for (boost::tie(ai, ai_end) = adjacent_vertices(u, g);
-       ai != ai_end; 
-       ++ai) {
-    n_daughters += 1;
-    int lund = abs(boost::get(lund_map, *ai));
-    switch (lund) {
-      case bdtaunu::D0Lund:
-        n_D0 += 1;
-        v = *ai;
-        break;
-      case bdtaunu::DcLund:
-        n_Dc += 1;
-        v = *ai;
-        break;
-      case bdtaunu::piLund:
-        n_pi += 1;
-        break;
-      case bdtaunu::pi0Lund:
-        n_pi0 += 1;
-        break;
-      case bdtaunu::gammaLund:
-        n_gamma += 1;
-        break;
-      default:
-        Dstarmap.insert(std::make_pair(u, Dstar()));
-        return;
-    }
-  }
-
-  int mode;
-  if (n_daughters == 2 && n_D0 == 1 && n_pi0 == 1) {
-    mode = bdtaunu::kDstar0_D0pi0;
-  } else if (n_daughters == 2 && n_D0 == 1 && n_gamma == 1) {
-    mode = bdtaunu::kDstar0_D0gamma;
-  } else if (n_daughters == 2 && n_D0 == 1 && n_pi == 1) {
-    mode = bdtaunu::kDstarc_D0pi;
-  } else if (n_daughters == 2 && n_Dc == 1 && n_pi0 == 1) {
-    mode = bdtaunu::kDstarc_Dcpi0;
-  } else if (n_daughters == 2 && n_Dc == 1 && n_gamma == 1) {
-    mode = bdtaunu::kDstarc_Dcgamma;
-  } else {
-    mode = bdtaunu::kUndefinedDstarMode;
-  }
-
-  Dstarmap.insert(std::make_pair(u, Dstar(mode, v)));
-
-}
-
-
-template<typename LundMap>
-void RecoGraphDfsVisitor<LundMap>::ComputeD(const Vertex &u, const Graph &g) {
-
-  assert(Dmap.find(u) == Dmap.end());
-
-  int n_daughters, n_K, n_Ks, n_pi, n_pi0;
-  n_daughters = n_K = n_Ks = n_pi = n_pi0 = 0;
-
-  AdjacencyIterator ai, ai_end;
-  for (boost::tie(ai, ai_end) = adjacent_vertices(u, g);
-       ai != ai_end; 
-       ++ai) {
-    n_daughters += 1;
-    int lund = abs(boost::get(lund_map, *ai));
-    switch (lund) {
-      case bdtaunu::KLund:
-        n_K += 1;
-        break;
-      case bdtaunu::KSLund:
-        n_Ks += 1;
-        break;
-      case bdtaunu::piLund:
-        n_pi += 1;
-        break;
-      case bdtaunu::pi0Lund:
-        n_pi0 += 1;
-        break;
-      default:
-        Dmap.insert(std::make_pair(u, D()));
-        return;
-    }
-  }
-
-  int mode;
-  if (std::abs(boost::get(lund_map, u)) == bdtaunu::DcLund) {
-    if (n_daughters == 3 && n_K == 1 && n_pi == 2) {
-      mode = bdtaunu::kDc_Kpipi;
-    } else if (n_daughters == 4 && n_K == 1 && n_pi == 2 && n_pi0 == 1) {
-      mode = bdtaunu::kDc_Kpipipi0;
-    } else if (n_daughters == 2 && n_K == 1 && n_Ks == 1) {
-      mode = bdtaunu::kDc_KsK;
-    } else if (n_daughters == 2 && n_Ks == 1 && n_pi == 1) {
-      mode = bdtaunu::kDc_Kspi;
-    } else if (n_daughters == 3 && n_Ks == 1 && n_pi == 1 && n_pi0 == 1) {
-      mode = bdtaunu::kDc_Kspipi0;
-    } else if (n_daughters == 4 && n_Ks == 1 && n_pi == 3) {
-      mode = bdtaunu::kDc_Kspipipi;
-    } else if (n_daughters == 3 && n_K == 2 && n_pi == 1) {
-      mode = bdtaunu::kDc_KKpi;
-    } else {
-      mode = bdtaunu::kUndefinedDMode;
-    }
-  } else {
-    if (n_daughters == 2 && n_K == 1 && n_pi == 1) {
-      mode = bdtaunu::kD0_Kpi;
-    } else if (n_daughters == 3 && n_K == 1 && n_pi == 1 && n_pi0 == 1) {
-      mode = bdtaunu::kD0_Kpipi0;
-    } else if (n_daughters == 4 && n_K == 1 && n_pi == 3) {
-      mode = bdtaunu::kD0_Kpipipi;
-    } else if (n_daughters == 5 && n_K == 1 && n_pi == 3 && n_pi0 == 1) {
-      mode = bdtaunu::kD0_Kpipipipi0;
-    } else if (n_daughters == 3 && n_Ks == 1 && n_pi == 2) {
-      mode = bdtaunu::kD0_Kspipi;
-    } else if (n_daughters == 4 && n_Ks == 1 && n_pi == 2 && n_pi0 == 1) {
-      mode = bdtaunu::kD0_Kspipipi0;
-    } else if (n_daughters == 2 && n_Ks == 1 && n_pi0 == 1) {
-      mode = bdtaunu::kD0_Kspi0;
-    } else if (n_daughters == 2 && n_K == 2) {
-      mode = bdtaunu::kD0_KK;
-    } else {
-      mode = bdtaunu::kUndefinedDMode;
-    }
-  }
-
-  Dmap.insert(std::make_pair(u, D(mode)));
-}
-
-template <typename LundMap>
-RecoGraphDfsVisitor<LundMap> reco_graph_dfs_visitor(
-    LundMap m, std::vector<BDtaunuReader::YDecayProperties> &u) {
-  return RecoGraphDfsVisitor<LundMap>(m, u);
-}
-
-/*
-class RecoGraphBfsVisitor : public boost::default_bfs_visitor {
-  private:
-    typedef typename BDtaunuReader::RecoGraph Graph;
-    typedef typename BDtaunuReader::Vertex Vertex;
+class RecoGraphDfsVisitor : public boost::default_dfs_visitor {
 
   public:
-    RecoGraphBfsVisitor(std::vector<Vertex> &v) : 
-      subgraph_vertices(v) {};
-    void discover_vertex(Vertex u, const Graph &g);
-    std::vector<Vertex> get_subgraph_vertices() { return subgraph_vertices; }
-
-  private:
-    std::vector<Vertex> &subgraph_vertices;
-};
-
-void RecoGraphBfsVisitor::discover_vertex(Vertex u, const Graph &g) {
-  subgraph_vertices.push_back(u);
-}
+    RecoGraphDfsVisitor() = default;
+    RecoGraphDfsVisitor(RecoGraphManager*);
+    ~RecoGraphDfsVisitor() {};
 
-RecoGraphBfsVisitor reco_graph_bfs_visitor(
-    std::vector<BDtaunuReader::Vertex> &u) {
-  return RecoGraphBfsVisitor(u);
-}
-*/
+    void finish_vertex(RecoGraph::Vertex u, const RecoGraph::Graph &g);
 
-class RecoGraphBfsVisitor : public boost::default_bfs_visitor {
   private:
-    typedef typename BDtaunuReader::RecoGraph Graph;
-    typedef typename boost::graph_traits<Graph>::vertex_descriptor Vertex;
-    typedef typename boost::graph_traits<Graph>::edge_descriptor Edge;
-    typedef typename boost::graph_traits<Graph>::adjacency_iterator AdjacencyIterator;
-
-  public:
-    RecoGraphBfsVisitor(const Graph &g, Graph &subg) : 
-      graph(g),
-      subgraph(subg) { 
-        reco_index = get(boost::vertex_reco_index, graph);
-        block_index = get(boost::vertex_block_index, graph);
-        lund_id = get(boost::vertex_lund_id, graph); 
-      };
-    void discover_vertex(Vertex u, const Graph &g);
-    void tree_edge(Edge e, const Graph &g);
-
-  private:
-    Graph graph;
-    Graph &subgraph;
-    std::map<int, Vertex> inserted_vertices;
-    typename boost::property_map<Graph, boost::vertex_reco_index_t>::type reco_index;
-    typename boost::property_map<Graph, boost::vertex_block_index_t>::type block_index;
-    typename boost::property_map<Graph, boost::vertex_lund_id_t>::type lund_id; 
+    RecoGraphManager *manager = nullptr;
+    RecoGraph::LundIdPropertyMap lund_map;
+    RecoGraph::BlockIndexPropertyMap block_idx_map;
+
+    void AnalyzeY(const RecoGraph::Vertex &u, const RecoGraph::Graph &g);
+    void AnalyzeB(const RecoGraph::Vertex &u, const RecoGraph::Graph &g);
+    void AnalyzeDstar(const RecoGraph::Vertex &u, const RecoGraph::Graph &g);
+    void AnalyzeD(const RecoGraph::Vertex &u, const RecoGraph::Graph &g);
+    void AnalyzeLepton(const RecoGraph::Vertex &u, const RecoGraph::Graph &g);
 };
 
-void RecoGraphBfsVisitor::discover_vertex(Vertex u, const Graph &g) {
-  if (inserted_vertices.find(reco_index[u]) == inserted_vertices.end()) {
-    Vertex v = add_vertex(subgraph);
-    put(get(boost::vertex_reco_index, subgraph), v, reco_index[u]);
-    put(get(boost::vertex_block_index, subgraph), v, block_index[u]);
-    put(get(boost::vertex_lund_id, subgraph), v, lund_id[u]);
-    inserted_vertices[reco_index[u]] = v;
-  }
-}
-
-void RecoGraphBfsVisitor::tree_edge(Edge e, const Graph &g) {
-
-  Vertex u_bfs = source(e, g);
-  Vertex v_bfs = target(e, g);
-  if (inserted_vertices.find(reco_index[v_bfs]) == inserted_vertices.end()) {
-    Vertex u = inserted_vertices[reco_index[u_bfs]];
-    Vertex v = add_vertex(subgraph);
-    put(get(boost::vertex_reco_index, subgraph), v, reco_index[v_bfs]);
-    put(get(boost::vertex_block_index, subgraph), v, block_index[v_bfs]);
-    put(get(boost::vertex_lund_id, subgraph), v, lund_id[v_bfs]);
-    inserted_vertices[reco_index[v_bfs]] = v;
-    add_edge(u, v, subgraph);
-  } else {
-    Vertex u = inserted_vertices[reco_index[u_bfs]];
-    Vertex v = inserted_vertices[reco_index[v_bfs]];
-    add_edge(u, v, subgraph);
-  }
-}
-
-RecoGraphBfsVisitor reco_graph_bfs_visitor(
-    const BDtaunuReader::RecoGraph &g,
-    BDtaunuReader::RecoGraph &subg) {
-  return RecoGraphBfsVisitor(g, subg);
-}
-
-
 #endif
diff --git a/create_database/RecoIndexer.cc b/create_database/RecoIndexer.cc
index f153464..a337ecb 100644
--- a/create_database/RecoIndexer.cc
+++ b/create_database/RecoIndexer.cc
@@ -16,7 +16,7 @@ RecoIndexer::RecoIndexer(
     nY(_nY), nB(_nB), nD(_nD), nC(_nC), 
     nh(_nh), nl(_nl), ngamma(_ngamma) {};
 
-int RecoIndexer::operator()(int lund, int idx) {
+int RecoIndexer::operator()(int lund, int idx) const {
 
   int abslund = std::abs(lund);
   switch (abslund) {
diff --git a/create_database/RecoIndexer.h b/create_database/RecoIndexer.h
index f4b8ff0..c94c0a9 100644
--- a/create_database/RecoIndexer.h
+++ b/create_database/RecoIndexer.h
@@ -16,7 +16,7 @@ class RecoIndexer {
     RecoIndexer &operator=(const RecoIndexer &r) = default;
     ~RecoIndexer() {};
 
-    int operator()(int lund, int idx);
+    int operator()(int lund, int idx) const;
 
     int total() const { return nY + nB + nD + nC + nh + nl + ngamma; }
     void clear();
diff --git a/create_database/RecoParticles.h b/create_database/RecoParticles.h
new file mode 100644
index 0000000..f917a14
--- /dev/null
+++ b/create_database/RecoParticles.h
@@ -0,0 +1,36 @@
+#ifndef _RECOPARTICLES_H_
+#define _RECOPARTICLES_H_
+
+#include "bdtaunu_definitions.h"
+#include "RecoGraph.h"
+
+struct RecoB;
+struct RecoD;
+struct RecoLepton;
+
+struct RecoY {
+  RecoY() = default;
+  RecoB *tagB = nullptr;
+  RecoB *sigB = nullptr;
+};
+
+struct RecoB {
+  RecoB() = default;
+  int flavor;
+  RecoD *D = nullptr;
+  RecoLepton *Lepton = nullptr;
+};
+
+struct RecoD {
+  RecoD() = default;
+  int D_mode = bdtaunu::kUndefinedDMode;
+  int Dstar_mode = bdtaunu::kUndefinedDstarMode;
+};
+
+struct RecoLepton {
+  int l_block_idx = -1;
+  int pi_block_idx = -1;
+  int tau_mode = bdtaunu::kUndefinedTauMode;
+};
+
+#endif
diff --git a/create_database/RootReader.h b/create_database/RootReader.h
index 64843c6..1ff96de 100644
--- a/create_database/RootReader.h
+++ b/create_database/RootReader.h
@@ -23,7 +23,6 @@ class RootReader {
     int total_records;
 
     void PrepareTreeFile(const char *root_fname, const char *root_trname);
-    virtual void SetBranchAddress() = 0;
 
   public:
 
diff --git a/create_database/UpsilonCandidate.cc b/create_database/UpsilonCandidate.cc
index e17572d..6341ecc 100644
--- a/create_database/UpsilonCandidate.cc
+++ b/create_database/UpsilonCandidate.cc
@@ -4,193 +4,156 @@
 #include "bdtaunu_definitions.h"
 #include "UpsilonCandidate.h"
 
-UpsilonCandidate::UpsilonCandidate() : 
-  eventId(""), event_candidate_index(-999),
-  reco_index(-1), 
-  eextra50(-999), mmiss_prime2(-999), 
-  tag_lp3(-999), sig_hp3(-999),
-  tag_cosBY(-999), sig_cosBY(-999), 
-  tag_cosThetaDl(-999), sig_cosThetaDtau(-999), 
-  sig_vtxB(-999), 
-  cosThetaT(-999), 
-  tag_Dmass(-999), tag_deltaM(-999), 
-  tag_cosThetaDSoft(-999), tag_softP3MagCM(-999),
-  sig_Dmass(-999), sig_deltaM(-999), 
-  sig_cosThetaDSoft(-999), sig_softP3MagCM(-999), 
-  sig_hmass(-999), sig_vtxh(-999), 
+UpsilonCandidate::UpsilonCandidate() :
+  eventId(""),
+  block_index(-999),
+  reco_index(-1),
   bflavor(bdtaunu::kUndefinedBFlavor),
-  tag_dstar_mode(bdtaunu::kUndefinedDstarMode), tag_d_mode(bdtaunu::kUndefinedDMode),
-  sig_dstar_mode(bdtaunu::kUndefinedDstarMode), sig_d_mode(bdtaunu::kUndefinedDMode),
-  sig_tau_mode(bdtaunu::kUndefinedTauMode), 
-  l_ePidMap(0), l_muPidMap(0),
-  h_ePidMap(0), h_muPidMap(0) {
-}
-
-UpsilonCandidate::UpsilonCandidate(
-    std::string &event_id, int candidate_idx, 
-    float eextra50, float mmiss_prime2, 
-    float tag_lp3, float sig_hp3, 
-    float tag_cosBY, float sig_cosBY, 
-    float tag_cosThetaDl, float sig_cosThetaDtau, 
-    float sig_vtxB, 
-    float cosThetaT, 
-    float tag_Dmass, float tag_deltaM, 
-    float tag_cosThetaDSoft, float tag_softP3MagCM,
-    float sig_Dmass, float sig_deltaM, 
-    float sig_cosThetaDSoft, float sig_softP3MagCM, 
-    float sig_hmass, float sig_vtxh, 
-    int bflavor,
-    int tag_dstar_mode, int tag_d_mode,
-    int sig_dstar_mode, int sig_d_mode,
-    int sig_tau_mode) :
-  eventId(event_id), event_candidate_index(candidate_idx),
-  reco_index(-1), 
-  eextra50(eextra50), mmiss_prime2(mmiss_prime2),
-  tag_lp3(tag_lp3), sig_hp3(sig_hp3),
-  tag_cosBY(tag_cosBY), sig_cosBY(sig_cosBY),
-  tag_cosThetaDl(tag_cosThetaDl), sig_cosThetaDtau(sig_cosThetaDtau),
-  sig_vtxB(sig_vtxB),
-  cosThetaT(cosThetaT),
-  tag_Dmass(tag_Dmass), tag_deltaM(tag_deltaM),
-  tag_cosThetaDSoft(tag_cosThetaDSoft), tag_softP3MagCM(tag_softP3MagCM),
-  sig_Dmass(sig_Dmass), sig_deltaM(sig_deltaM),
-  sig_cosThetaDSoft(sig_cosThetaDSoft), sig_softP3MagCM(sig_softP3MagCM),
-  sig_hmass(sig_hmass), sig_vtxh(sig_vtxh),
-  bflavor(bflavor),
-  tag_dstar_mode(tag_dstar_mode), tag_d_mode(tag_d_mode),
-  sig_dstar_mode(sig_dstar_mode), sig_d_mode(sig_d_mode),
-  sig_tau_mode(sig_tau_mode),
-  l_ePidMap(0), l_muPidMap(0),
-  h_ePidMap(0), h_muPidMap(0) {
-}
+  eextra50(-999),
+  mmiss_prime2(-999),
+  cosThetaT(-999),
+  tag_lp3(-999),
+  tag_cosBY(-999),
+  tag_cosThetaDl(-999),
+  tag_Dmass(-999),
+  tag_deltaM(-999),
+  tag_cosThetaDSoft(-999),
+  tag_softP3MagCM(-999),
+  tag_d_mode(bdtaunu::kUndefinedDMode),
+  tag_dstar_mode(bdtaunu::kUndefinedDstarMode),
+  l_ePidMap(0),
+  l_muPidMap(0),
+  sig_hp3(-999),
+  sig_cosBY(-999),
+  sig_cosThetaDtau(-999),
+  sig_vtxB(-999),
+  sig_Dmass(-999),
+  sig_deltaM(-999),
+  sig_cosThetaDSoft(-999),
+  sig_softP3MagCM(-999),
+  sig_hmass(-999),
+  sig_vtxh(-999),
+  sig_d_mode(bdtaunu::kUndefinedDMode),
+  sig_dstar_mode(bdtaunu::kUndefinedDstarMode),
+  sig_tau_mode(bdtaunu::kUndefinedTauMode),
+  h_ePidMap(0),
+  h_muPidMap(0) {}
 
 UpsilonCandidate::UpsilonCandidate(
-    std::string &event_id, int candidate_idx, 
-    float eextra50, float mmiss_prime2, 
-    float tag_lp3, float sig_hp3, 
-    float tag_cosBY, float sig_cosBY, 
-    float tag_cosThetaDl, float sig_cosThetaDtau, 
-    float sig_vtxB, 
-    float cosThetaT, 
-    float tag_Dmass, float tag_deltaM, 
-    float tag_cosThetaDSoft, float tag_softP3MagCM,
-    float sig_Dmass, float sig_deltaM, 
-    float sig_cosThetaDSoft, float sig_softP3MagCM, 
-    float sig_hmass, float sig_vtxh, 
-    int bflavor,
-    int tag_dstar_mode, int tag_d_mode,
-    int sig_dstar_mode, int sig_d_mode,
-    int sig_tau_mode, 
-    int l_ePidMap, int l_muPidMap) :
-  eventId(event_id), event_candidate_index(candidate_idx),
-  reco_index(-1), 
-  eextra50(eextra50), mmiss_prime2(mmiss_prime2),
-  tag_lp3(tag_lp3), sig_hp3(sig_hp3),
-  tag_cosBY(tag_cosBY), sig_cosBY(sig_cosBY),
-  tag_cosThetaDl(tag_cosThetaDl), sig_cosThetaDtau(sig_cosThetaDtau),
-  sig_vtxB(sig_vtxB),
-  cosThetaT(cosThetaT),
-  tag_Dmass(tag_Dmass), tag_deltaM(tag_deltaM),
-  tag_cosThetaDSoft(tag_cosThetaDSoft), tag_softP3MagCM(tag_softP3MagCM),
-  sig_Dmass(sig_Dmass), sig_deltaM(sig_deltaM),
-  sig_cosThetaDSoft(sig_cosThetaDSoft), sig_softP3MagCM(sig_softP3MagCM),
-  sig_hmass(sig_hmass), sig_vtxh(sig_vtxh),
-  bflavor(bflavor),
-  tag_dstar_mode(tag_dstar_mode), tag_d_mode(tag_d_mode),
-  sig_dstar_mode(sig_dstar_mode), sig_d_mode(sig_d_mode),
-  sig_tau_mode(sig_tau_mode),
-  l_ePidMap(l_ePidMap), l_muPidMap(l_muPidMap),
-  h_ePidMap(0), h_muPidMap(0) {
-}
-
-UpsilonCandidate::UpsilonCandidate(
-    std::string &event_id, int candidate_idx, 
-    int reco_idx,
-    float eextra50, float mmiss_prime2, 
-    float tag_lp3, float sig_hp3, 
-    float tag_cosBY, float sig_cosBY, 
-    float tag_cosThetaDl, float sig_cosThetaDtau, 
-    float sig_vtxB, 
-    float cosThetaT, 
-    float tag_Dmass, float tag_deltaM, 
-    float tag_cosThetaDSoft, float tag_softP3MagCM,
-    float sig_Dmass, float sig_deltaM, 
-    float sig_cosThetaDSoft, float sig_softP3MagCM, 
-    float sig_hmass, float sig_vtxh, 
-    int bflavor,
-    int tag_dstar_mode, int tag_d_mode,
-    int sig_dstar_mode, int sig_d_mode,
-    int sig_tau_mode, 
-    int l_ePidMap, int l_muPidMap,
-    int h_ePidMap, int h_muPidMap) :
-  eventId(event_id), event_candidate_index(candidate_idx),
-  reco_index(reco_idx), 
-  eextra50(eextra50), mmiss_prime2(mmiss_prime2),
-  tag_lp3(tag_lp3), sig_hp3(sig_hp3),
-  tag_cosBY(tag_cosBY), sig_cosBY(sig_cosBY),
-  tag_cosThetaDl(tag_cosThetaDl), sig_cosThetaDtau(sig_cosThetaDtau),
-  sig_vtxB(sig_vtxB),
-  cosThetaT(cosThetaT),
-  tag_Dmass(tag_Dmass), tag_deltaM(tag_deltaM),
-  tag_cosThetaDSoft(tag_cosThetaDSoft), tag_softP3MagCM(tag_softP3MagCM),
-  sig_Dmass(sig_Dmass), sig_deltaM(sig_deltaM),
-  sig_cosThetaDSoft(sig_cosThetaDSoft), sig_softP3MagCM(sig_softP3MagCM),
-  sig_hmass(sig_hmass), sig_vtxh(sig_vtxh),
-  bflavor(bflavor),
-  tag_dstar_mode(tag_dstar_mode), tag_d_mode(tag_d_mode),
-  sig_dstar_mode(sig_dstar_mode), sig_d_mode(sig_d_mode),
-  sig_tau_mode(sig_tau_mode),
-  l_ePidMap(l_ePidMap), l_muPidMap(l_muPidMap),
-  h_ePidMap(h_ePidMap), h_muPidMap(h_muPidMap) {
-}
+  std::string &_eventId,
+  int _block_index,
+  int _reco_index,
+  int _bflavor,
+  float _eextra50,
+  float _mmiss_prime2,
+  float _cosThetaT,
+  float _tag_lp3,
+  float _tag_cosBY,
+  float _tag_cosThetaDl,
+  float _tag_Dmass,
+  float _tag_deltaM,
+  float _tag_cosThetaDSoft,
+  float _tag_softP3MagCM,
+  int _tag_d_mode,
+  int _tag_dstar_mode,
+  int _l_ePidMap,
+  int _l_muPidMap,
+  float _sig_hp3,
+  float _sig_cosBY,
+  float _sig_cosThetaDtau,
+  float _sig_vtxB,
+  float _sig_Dmass,
+  float _sig_deltaM,
+  float _sig_cosThetaDSoft,
+  float _sig_softP3MagCM,
+  float _sig_hmass,
+  float _sig_vtxh,
+  int _sig_d_mode,
+  int _sig_dstar_mode,
+  int _sig_tau_mode,
+  int _h_ePidMap,
+  int _h_muPidMap) : 
+  eventId(_eventId),
+  block_index(_block_index),
+  reco_index(_reco_index),
+  bflavor(_bflavor),
+  eextra50(_eextra50),
+  mmiss_prime2(_mmiss_prime2),
+  cosThetaT(_cosThetaT),
+  tag_lp3(_tag_lp3),
+  tag_cosBY(_tag_cosBY),
+  tag_cosThetaDl(_tag_cosThetaDl),
+  tag_Dmass(_tag_Dmass),
+  tag_deltaM(_tag_deltaM),
+  tag_cosThetaDSoft(_tag_cosThetaDSoft),
+  tag_softP3MagCM(_tag_softP3MagCM),
+  tag_d_mode(_tag_d_mode),
+  tag_dstar_mode(_tag_dstar_mode),
+  l_ePidMap(_l_ePidMap),
+  l_muPidMap(_l_muPidMap),
+  sig_hp3(_sig_hp3),
+  sig_cosBY(_sig_cosBY),
+  sig_cosThetaDtau(_sig_cosThetaDtau),
+  sig_vtxB(_sig_vtxB),
+  sig_Dmass(_sig_Dmass),
+  sig_deltaM(_sig_deltaM),
+  sig_cosThetaDSoft(_sig_cosThetaDSoft),
+  sig_softP3MagCM(_sig_softP3MagCM),
+  sig_hmass(_sig_hmass),
+  sig_vtxh(_sig_vtxh),
+  sig_d_mode(_sig_d_mode),
+  sig_dstar_mode(_sig_dstar_mode),
+  sig_tau_mode(_sig_tau_mode),
+  h_ePidMap(_h_ePidMap),
+  h_muPidMap(_h_muPidMap) {}
 
 UpsilonCandidate::UpsilonCandidate(const UpsilonCandidate &cand) {
   copy_candidate(cand);
 }
 
+UpsilonCandidate & UpsilonCandidate::operator=(const UpsilonCandidate &cand) {
+  if (this != &cand) {
+    copy_candidate(cand);
+  }
+  return *this;
+}
+
 void UpsilonCandidate::copy_candidate(const UpsilonCandidate &cand) {
   eventId = cand.eventId;
-  event_candidate_index = cand.event_candidate_index;
+  block_index = cand.block_index;
   reco_index = cand.reco_index;
+  bflavor = cand.bflavor;
   eextra50 = cand.eextra50;
   mmiss_prime2 = cand.mmiss_prime2;
+  cosThetaT = cand.cosThetaT;
   tag_lp3 = cand.tag_lp3;
-  sig_hp3 = cand.sig_hp3;
   tag_cosBY = cand.tag_cosBY;
-  sig_cosBY = cand.sig_cosBY;
   tag_cosThetaDl = cand.tag_cosThetaDl;
-  sig_cosThetaDtau = cand.sig_cosThetaDtau;
-  sig_vtxB = cand.sig_vtxB;
-  cosThetaT = cand.cosThetaT;
   tag_Dmass = cand.tag_Dmass;
   tag_deltaM = cand.tag_deltaM;
   tag_cosThetaDSoft = cand.tag_cosThetaDSoft;
   tag_softP3MagCM = cand.tag_softP3MagCM;
+  tag_d_mode = cand.tag_d_mode;
+  tag_dstar_mode = cand.tag_dstar_mode;
+  l_ePidMap = cand.l_ePidMap;
+  l_muPidMap = cand.l_muPidMap;
+  sig_hp3 = cand.sig_hp3;
+  sig_cosBY = cand.sig_cosBY;
+  sig_cosThetaDtau = cand.sig_cosThetaDtau;
+  sig_vtxB = cand.sig_vtxB;
   sig_Dmass = cand.sig_Dmass;
   sig_deltaM = cand.sig_deltaM;
   sig_cosThetaDSoft = cand.sig_cosThetaDSoft;
   sig_softP3MagCM = cand.sig_softP3MagCM;
   sig_hmass = cand.sig_hmass;
   sig_vtxh = cand.sig_vtxh;
-  bflavor = cand.bflavor;
-  tag_dstar_mode = cand.tag_dstar_mode;
-  tag_d_mode = cand.tag_d_mode;
-  sig_dstar_mode = cand.sig_dstar_mode;
   sig_d_mode = cand.sig_d_mode;
+  sig_dstar_mode = cand.sig_dstar_mode;
   sig_tau_mode = cand.sig_tau_mode;
-  l_ePidMap = cand.l_ePidMap;
-  l_muPidMap = cand.l_muPidMap;
   h_ePidMap = cand.h_ePidMap;
   h_muPidMap = cand.h_muPidMap;
 }
 
-UpsilonCandidate & UpsilonCandidate::operator=(const UpsilonCandidate &cand) {
-  if (this != &cand) {
-    copy_candidate(cand);
-  }
-  return *this;
-}
-
 // Examine the D, D*, and tau modes to determine the candidate type. 
 int UpsilonCandidate::get_cand_type() const {
 
diff --git a/create_database/UpsilonCandidate.h b/create_database/UpsilonCandidate.h
index ea12b58..3d70646 100644
--- a/create_database/UpsilonCandidate.h
+++ b/create_database/UpsilonCandidate.h
@@ -6,120 +6,78 @@
 //! Class representing an \f$\Upsilon(4S)\f$ candidate
 class UpsilonCandidate {
 
-  private:
-    std::string eventId;
-    int event_candidate_index;
-    int reco_index;
-    float eextra50, mmiss_prime2;
-    float tag_lp3, sig_hp3;
-    float tag_cosBY, sig_cosBY;
-    float tag_cosThetaDl, sig_cosThetaDtau;
-    float sig_vtxB;
-    float cosThetaT;
-    float tag_Dmass, tag_deltaM;
-    float tag_cosThetaDSoft, tag_softP3MagCM;
-    float sig_Dmass, sig_deltaM;
-    float sig_cosThetaDSoft, sig_softP3MagCM;
-    float sig_hmass, sig_vtxh;
-    int bflavor;
-    int tag_dstar_mode, tag_d_mode;
-    int sig_dstar_mode, sig_d_mode;
-    int sig_tau_mode;
-    
-    int l_ePidMap, l_muPidMap;
-    int h_ePidMap, h_muPidMap;
-
-    void copy_candidate(const UpsilonCandidate &cand);
-
   public:
     //! Constructs candidate with non-physical attributes.
     UpsilonCandidate();
 
     //! Constructs candidate with specified attributes. 
     UpsilonCandidate(
-        std::string& eventId, int event_candidate_index, 
-        float eextra50, float mmiss_prime2, 
-        float tag_lp3, float sig_hp3, 
-        float tag_cosBY, float sig_cosBY, 
-        float tag_cosThetaDl, float sig_cosThetaDtau, 
-        float sig_vtxB, 
-        float cosThetaT, 
-        float tag_Dmass, float tag_deltaM, 
-        float tag_cosThetaDSoft, float tag_softP3MagCM,
-        float sig_Dmass, float sig_deltaM, 
-        float sig_cosThetaDSoft, float sig_softP3MagCM, 
-        float sig_hmass, float sig_vtxh, 
-        int bflavor, 
-        int tag_dstar_mode, int tag_d_mode,
-        int sig_dstar_mode, int sig_d_mode,
-        int sig_tau_mode);
-
-    //! Constructs candidate with specified attributes. 
-    UpsilonCandidate(
-        std::string& eventId, int event_candidate_index, 
-        float eextra50, float mmiss_prime2, 
-        float tag_lp3, float sig_hp3, 
-        float tag_cosBY, float sig_cosBY, 
-        float tag_cosThetaDl, float sig_cosThetaDtau, 
-        float sig_vtxB, 
-        float cosThetaT, 
-        float tag_Dmass, float tag_deltaM, 
-        float tag_cosThetaDSoft, float tag_softP3MagCM,
-        float sig_Dmass, float sig_deltaM, 
-        float sig_cosThetaDSoft, float sig_softP3MagCM, 
-        float sig_hmass, float sig_vtxh, 
-        int bflavor, 
-        int tag_dstar_mode, int tag_d_mode,
-        int sig_dstar_mode, int sig_d_mode,
-        int sig_tau_mode,
-        int l_ePidMap, int l_muPidMap); 
-
-    //! Constructs candidate with specified attributes. 
-    UpsilonCandidate(
-        std::string& eventId, int event_candidate_index, 
-        int reco_index, 
-        float eextra50, float mmiss_prime2, 
-        float tag_lp3, float sig_hp3, 
-        float tag_cosBY, float sig_cosBY, 
-        float tag_cosThetaDl, float sig_cosThetaDtau, 
-        float sig_vtxB, 
-        float cosThetaT, 
-        float tag_Dmass, float tag_deltaM, 
-        float tag_cosThetaDSoft, float tag_softP3MagCM,
-        float sig_Dmass, float sig_deltaM, 
-        float sig_cosThetaDSoft, float sig_softP3MagCM, 
-        float sig_hmass, float sig_vtxh, 
-        int bflavor, 
-        int tag_dstar_mode, int tag_d_mode,
-        int sig_dstar_mode, int sig_d_mode,
-        int sig_tau_mode,
-        int l_ePidMap, int l_muPidMap,
-        int h_ePidMap, int h_muPidMap); 
-
+      std::string &eventId,
+      int block_index,
+      int reco_index,
+      int bflavor,
+      float eextra50,
+      float mmiss_prime2,
+      float cosThetaT,
+      float tag_lp3,
+      float tag_cosBY,
+      float tag_cosThetaDl,
+      float tag_Dmass,
+      float tag_deltaM,
+      float tag_cosThetaDSoft,
+      float tag_softP3MagCM,
+      int tag_d_mode,
+      int tag_dstar_mode,
+      int l_ePidMap,
+      int l_muPidMap,
+      float sig_hp3,
+      float sig_cosBY,
+      float sig_cosThetaDtau,
+      float sig_vtxB,
+      float sig_Dmass,
+      float sig_deltaM,
+      float sig_cosThetaDSoft,
+      float sig_softP3MagCM,
+      float sig_hmass,
+      float sig_vtxh,
+      int sig_d_mode,
+      int sig_dstar_mode,
+      int sig_tau_mode,
+      int h_ePidMap,
+      int h_muPidMap);
     UpsilonCandidate(const UpsilonCandidate &cand);
+    UpsilonCandidate & operator=(const UpsilonCandidate &cand);
     ~UpsilonCandidate() {};
 
-    UpsilonCandidate & operator=(const UpsilonCandidate &cand);
 
     //! The babar ID of the event this candidate belongs to. 
     std::string get_eventId() const { return eventId; }
+    void set_eventId(std::string _eventId) { eventId = _eventId; }
+
+    //! The index that uniquely identifies the upsilon candidate within the candidate block.
+    int get_block_index() const { return block_index; }
+    void set_block_index(int _block_index) { block_index = _block_index; }
 
-    //! The index that uniquely identifies the candidate within the event.
-    int get_event_candidate_index() const { return event_candidate_index; }
+    //! The index that uniquely identifies the reconstructed candidate within the event.
+    int get_reco_index() const { return reco_index; }
+    void set_reco_index(int _reco_index) { reco_index = _reco_index; }
 
     //! \f$ E_{extra} \f$ 
     /*! The energy sum of photons above 50 MeV that are not used in the
      * reconstruction of this candidate. */
     float get_eextra50() const { return eextra50; }
+    void set_eextra50(float _eextra50) { eextra50 = _eextra50; }
 
     //! \f$ M^2_{miss} \f$
     /*! This is just like the ordinary missing mass squared, except the
      * beam energy information is also incorporated to reduce the
      * uncertainty. */
     float get_mmiss_prime2() const { return mmiss_prime2; }
+    void set_mmiss_prime2(float _mmiss_prime2) { mmiss_prime2 = _mmiss_prime2; }
 
     //! \f$ |\vec{p}_\ell| \f$ of the \f$B_{tag}\f$
     float get_tag_lp3() const { return tag_lp3; }
+    void set_tag_lp3(float _tag_lp3) { tag_lp3 = _tag_lp3; }
 
     //! \f$ |\vec{p}_h| \f$ of the \f$B_{sig}\f$
     /*! \f$h\f$ is either \f$\pi\f$ or \f$\rho\f$. The letter "h"
@@ -127,127 +85,200 @@ class UpsilonCandidate {
      * the reconstructed \f$\tau\f$ since the neutrino is
      * undetectable. */
     float get_sig_hp3() const { return sig_hp3; }
+    void set_sig_hp3(float _sig_hp3) { sig_hp3 = _sig_hp3; }
 
     //! \f$ \cos\theta_{BY} \f$ of the \f$B_{tag}\f$
     float get_tag_cosBY() const { return tag_cosBY; }
+    void set_tag_cosBY(float _tag_cosBY) { tag_cosBY = _tag_cosBY; }
 
     //! \f$ \cos\theta_{BY} \f$ of the \f$B_{sig}\f$
     float get_sig_cosBY() const { return sig_cosBY; }
+    void set_sig_cosBY(float _sig_cosBY) { sig_cosBY = _sig_cosBY; }
 
     //! \f$ \cos\theta_{D\ell} \f$ of the \f$B_{tag}\f$
     float get_tag_cosThetaDl() const { return tag_cosThetaDl; }
+    void set_tag_cosThetaDl(float _tag_cosThetaDl) { tag_cosThetaDl = _tag_cosThetaDl; }
 
     //! \f$ \cos\theta_{D\tau} \f$ of the \f$B_{sig}\f$
     /*! This is actually \f$ \cos\theta_{D\tau} \f$ to be precise,
      * since the \f$\tau\f$ itself is of course unreconstructable. */
     float get_sig_cosThetaDtau() const { return sig_cosThetaDtau; }
+    void set_sig_cosThetaDtau(float _sig_cosThetaDtau) { sig_cosThetaDtau = _sig_cosThetaDtau; }
 
     //! \f$B_{sig}\f$ vertex quality. 
     /*! It is really (1 - pvalue)? of the \f$\chi^2\f$ fit. */
     float get_sig_vtxB() const { return sig_vtxB; }
+    void set_sig_vtxB(float _sig_vtxB) { sig_vtxB = _sig_vtxB; }
 
     //! \f$ \cos\theta_{T} \f$
     /*! This is the cosine of the angle that the thrust vector makes
      * with the beam axis.  */
     float get_cosThetaT() const { return cosThetaT; }
+    void set_cosThetaT(float _cosThetaT) { cosThetaT = _cosThetaT; }
 
     //! Mass of the \f$ D \f$ meson of the \f$B_{tag}\f$.
     float get_tag_Dmass() const { return tag_Dmass; }
+    void set_tag_Dmass(float _tag_Dmass) { tag_Dmass = _tag_Dmass; }
 
     //! \f$\Delta m\f$ of the \f$ D \f$ meson of the \f$B_{tag}\f$.
     float get_tag_deltaM() const { return tag_deltaM; }
+    void set_tag_deltaM(float _tag_deltaM) { tag_deltaM = _tag_deltaM; }
 
     //! \f$ \cos\theta_{D soft} \f$ of the \f$B_{tag}\f$
     /*! This is the cosine of the angle between the daughters of the
      * \f$D^*\f$ decay; namely the angle between the \f$ D \f$ meson
      * and the soft daughter (photon or a \f$\pi^0\f$). */
     float get_tag_cosThetaDSoft() const { return tag_cosThetaDSoft; }
+    void set_tag_cosThetaDSoft(float _tag_cosThetaDSoft) { tag_cosThetaDSoft = _tag_cosThetaDSoft; }
 
     //! \f$ |\vec{p}_{soft}| \f$ of the \f$B_{tag}\f$
     /*! Magnitude of the three momentum of the \f$D^*\f$ meson's soft
      * daughter in the center of mass frame.  */
     float get_tag_softP3MagCM() const { return tag_softP3MagCM; }
+    void set_tag_softP3MagCM(float _tag_softP3MagCM) { tag_softP3MagCM = _tag_softP3MagCM; }
 
     //! Mass of the \f$ D \f$ meson of the \f$B_{sig}\f$.
     float get_sig_Dmass() const { return sig_Dmass; }
+    void set_sig_Dmass(float _sig_Dmass) { sig_Dmass = _sig_Dmass; }
 
     //! \f$\Delta m\f$ of the \f$ D \f$ meson of the \f$B_{sig}\f$.
     float get_sig_deltaM() const { return sig_deltaM; }
+    void set_sig_deltaM(float _sig_deltaM) { sig_deltaM = _sig_deltaM; }
 
     //! \f$ \cos\theta_{D soft} \f$ of the \f$B_{sig}\f$
     /*! This is the cosine of the angle between the daughters of the
      * \f$D^*\f$ decay; namely the angle between the \f$ D \f$ meson
      * and the soft daughter (photon or a \f$\pi^0\f$). */
     float get_sig_cosThetaDSoft() const { return sig_cosThetaDSoft; }
+    void set_sig_cosThetaDSoft(float _sig_cosThetaDSoft) { sig_cosThetaDSoft = _sig_cosThetaDSoft; }
 
     //! \f$ |\vec{p}_{soft}| \f$ of the \f$B_{sig}\f$
     /*! Magnitude of the three momentum of the \f$D^*\f$ meson's soft
      * daughter in the center of mass frame. */
     float get_sig_softP3MagCM() const { return sig_softP3MagCM; }
+    void set_sig_softP3MagCM(float _sig_softP3MagCM) { sig_softP3MagCM = _sig_softP3MagCM; }
 
     //! Mass of the \f$h\f$. 
     float get_sig_hmass() const { return sig_hmass; }
+    void set_sig_hmass(float _sig_hmass) { sig_hmass = _sig_hmass; }
 
     //! \f$h\f$ vertex quality. 
     /*! Vertex quality of the \f$h\f$. In this case, it can only be
      * the \f$\rho^+\f$. */
     float get_sig_vtxh() const { return sig_vtxh; }
+    void set_sig_vtxh(float _sig_vtxh) { sig_vtxh = _sig_vtxh; }
 
     //! \f$B\f$ flavor. 
     /*! Returns an integer that corresponds to the #BFlavor enum in
      * bdtaunu_definitions.h */
     int get_bflavor() const { return bflavor; }
+    void set_bflavor(int _bflavor) { bflavor = _bflavor; }
  
     //! Reconstructed \f$D^*\f$ mode index of the \f$B_{tag}\f$. 
     /*! Returns an integer that corresponds to the #DstarMode enum in
      * bdtaunu_definitions.h */
     int get_tag_dstar_mode() const { return tag_dstar_mode; }
+    void set_tag_dstar_mode(int _tag_dstar_mode) { tag_dstar_mode = _tag_dstar_mode; }
 
     //! Reconstructed \f$D\f$ mode index of the \f$B_{tag}\f$. 
     /*! Returns an integer that corresponds to the #DMode enum in
      * bdtaunu_definitions.h */
     int get_tag_d_mode() const { return tag_d_mode; }
+    void set_tag_d_mode(int _tag_d_mode) { tag_d_mode = _tag_d_mode; }
 
     //! Reconstructed \f$D^*\f$ mode index of the \f$B_{sig}\f$. 
     /*! Returns an integer that corresponds to the #DstarMode enum in
      * bdtaunu_definitions.h */
     int get_sig_dstar_mode() const { return sig_dstar_mode; }
+    void set_sig_dstar_mode(int _sig_dstar_mode) { sig_dstar_mode = _sig_dstar_mode; }
 
     //! Reconstructed \f$D\f$ mode index of the \f$B_{sig}\f$. 
     /*! Returns an integer that corresponds to the #DMode enum in
      * bdtaunu_definitions.h */
     int get_sig_d_mode() const { return sig_d_mode; }
+    void set_sig_d_mode(int _sig_d_mode) { sig_d_mode = _sig_d_mode; }
 
     //! Reconstructed \f$\tau\f$ mode index of the \f$B_{sig}\f$. 
     /*! Returns an integer that corresponds to the #TauMode enum in
      * bdtaunu_definitions.h */
     int get_sig_tau_mode() const { return sig_tau_mode; }
-
-    //! Candidate type. 
-    /*! Returns an int that corresponds to the #CandType enum in
-     * bdtaunu_definitions.h */
-    int get_cand_type() const;
-
-    //! Sample type. 
-    /*! Returns an int that corresponds to the #SampleType enum in
-     * bdtaunu_definitions.h */
-    int get_sample_type() const;
+    void set_sig_tau_mode(int _sig_tau_mode) { sig_tau_mode = _sig_tau_mode; }
 
     //! Electron PID map of tag lepton.
     /*! Bit map is here:
      * http://www.slac.stanford.edu/BFROOT/www/Physics/Tools/Pid/Selectors/r24c/selectors.html
      */
     int get_l_ePidMap() const { return l_ePidMap; }
+    void set_l_ePidMap(int _l_ePidMap) { l_ePidMap = _l_ePidMap; }
+    
+    //! Electron PID map of sig pion.
+    /*! Bit map is here:
+     * http://www.slac.stanford.edu/BFROOT/www/Physics/Tools/Pid/Selectors/r24c/selectors.html
+     */
     int get_h_ePidMap() const { return h_ePidMap; }
+    void set_h_ePidMap(int _h_ePidMap) { h_ePidMap = _h_ePidMap; }
 
     //! Muon PID map of tag lepton.
     /*! Bit map is here:
      * http://www.slac.stanford.edu/BFROOT/www/Physics/Tools/Pid/Selectors/r24c/selectors.html
      */
     int get_l_muPidMap() const { return l_muPidMap; }
+    void set_l_muPidMap(int _l_muPidMap) { l_muPidMap = _l_muPidMap; }
+
+    //! Muon PID map of sig pion.
+    /*! Bit map is here:
+     * http://www.slac.stanford.edu/BFROOT/www/Physics/Tools/Pid/Selectors/r24c/selectors.html
+     */
     int get_h_muPidMap() const { return h_muPidMap; }
+    void set_h_muPidMap(int _h_muPidMap) { h_muPidMap = _h_muPidMap; }
 
-    int get_reco_index() const { return reco_index; }
+    //! Candidate type. 
+    /*! Returns an int that corresponds to the #CandType enum in
+     * bdtaunu_definitions.h */
+    int get_cand_type() const;
+
+    //! Sample type. 
+    /*! Returns an int that corresponds to the #SampleType enum in
+     * bdtaunu_definitions.h */
+    int get_sample_type() const;
+
+  private:
+    std::string eventId;
+    int block_index;
+    int reco_index;
+    int bflavor;
+    float eextra50;
+    float mmiss_prime2;
+    float cosThetaT;
+    float tag_lp3;
+    float tag_cosBY;
+    float tag_cosThetaDl;
+    float tag_Dmass;
+    float tag_deltaM;
+    float tag_cosThetaDSoft;
+    float tag_softP3MagCM;
+    int tag_d_mode;
+    int tag_dstar_mode;
+    int l_ePidMap;
+    int l_muPidMap;
+    float sig_hp3;
+    float sig_cosBY;
+    float sig_cosThetaDtau;
+    float sig_vtxB;
+    float sig_Dmass;
+    float sig_deltaM;
+    float sig_cosThetaDSoft;
+    float sig_softP3MagCM;
+    float sig_hmass;
+    float sig_vtxh;
+    int sig_d_mode;
+    int sig_dstar_mode;
+    int sig_tau_mode;
+    int h_ePidMap;
+    int h_muPidMap;
+
+    // constructor helper
+    void copy_candidate(const UpsilonCandidate &cand);
 };
 
 #endif
diff --git a/create_database/UpsilonList.cc b/create_database/UpsilonList.cc
deleted file mode 100644
index 34edfe6..0000000
--- a/create_database/UpsilonList.cc
+++ /dev/null
@@ -1,82 +0,0 @@
-#include "UpsilonList.h"
-#include "UpsilonCandidate.h"
-
-UpsilonList::UpsilonList() {
-  head = NULL;
-  iter = NULL;
-  tail = NULL;
-}
-
-UpsilonList::UpsilonList(const UpsilonList &upslist) {
-  copy(upslist);
-}
-
-UpsilonList::~UpsilonList() {
-  clear();
-}
-
-// Delete the linked list.
-void UpsilonList::clear() {
-  if (head != NULL) {
-    node *curr = head;
-    while (curr != NULL) {
-      head = curr->next;
-      delete curr;
-      curr = head;
-    }
-    tail = NULL;
-  }
-  iter = NULL;
-}
-
-// Starting from an empty list, create a new list that is a copy of
-// another list. 
-void UpsilonList::copy(const UpsilonList &upslist) {
-  head = NULL;
-  tail = NULL;
-  if (upslist.head != NULL) {
-    node *other_curr = upslist.head;
-    while (other_curr != NULL) {
-      add_candidate(other_curr->candidate);
-      other_curr = other_curr->next;
-    }
-  }
-  iter = NULL;
-}
-
-UpsilonList & UpsilonList::operator=(const UpsilonList &upslist) {
-  if (this != &upslist) {
-    clear();
-    copy(upslist);
-  }
-  return *this;
-}
-
-// Add a new candidate at the end of the list. 
-void UpsilonList::add_candidate(const UpsilonCandidate &cand) {
-  node *new_node = new node(cand);
-  if (head == NULL) {
-    head = new_node;
-    tail = head;
-  } else {
-    tail->next = new_node;
-    tail = new_node;
-  }
-}
-
-// Move the iterator to the next candidate. 
-int UpsilonList::next_candidate() {
-  if (iter == NULL) {
-    if (head == NULL) {
-      return -1;
-    } else {
-      iter = head;
-      return (iter->candidate).get_event_candidate_index();
-    }
-  } else if (iter == tail) {
-    return -1;
-  } else {
-    iter = iter->next;
-    return (iter->candidate).get_event_candidate_index();
-  }
-}
diff --git a/create_database/UpsilonList.h b/create_database/UpsilonList.h
deleted file mode 100644
index 8f2e9e4..0000000
--- a/create_database/UpsilonList.h
+++ /dev/null
@@ -1,49 +0,0 @@
-#ifndef __UPSILONLIST__
-#define __UPSILONLIST__
-
-#include <string>
-
-#include "UpsilonCandidate.h"
-
-//! Linked list that stores UpsilonCandidate objects from a single event. 
-/*! Currently supports only single pass iteration and back
- * insertion. The list has a pointer that tracks the current candidate
- * available for reading/extraction. */
-class UpsilonList {
-  private:
-
-    struct node {
-      UpsilonCandidate candidate;
-      node *next;
-
-      node(const UpsilonCandidate &cand, node *next = NULL) : 
-        candidate(cand), next(next) {}
-    };
-
-    node *head, *tail;
-    node *iter;
-
-    void clear();
-    void copy(const UpsilonList &);
-
-  public:
-    UpsilonList();
-    UpsilonList(const UpsilonList &);
-    ~UpsilonList();
-
-    UpsilonList & operator=(const UpsilonList &);
-
-    //! Insert an UpsilonCandidate at the end of the linked list. 
-    void add_candidate(const UpsilonCandidate &);
-
-    //! Move to the UpsilonCandidate in the list for reading.
-    /*! This simply moves the tracking pointer forward by one. It
-     * returns the new candidate's candidate index in the event;
-     * return -1 if the end of list is reached. */
-    int next_candidate();
-
-    //! Get the current candidate the tracking pointer is pointing to. 
-    const UpsilonCandidate & get_current_candidate() const { return iter->candidate; }
-};
-
-#endif
diff --git a/include/bdtaunu_definitions.h b/include/bdtaunu_definitions.h
index e24a333..efb4764 100644
--- a/include/bdtaunu_definitions.h
+++ b/include/bdtaunu_definitions.h
@@ -68,6 +68,8 @@ enum DstarMode {
 enum TauMode {
   ktau_pi = 1,               /*!< \f$ \tau^+\rightarrow \pi^+ \f$ */
   ktau_rho = 2,              /*!< \f$ \tau^+\rightarrow \rho^+ \f$ */
+  ktau_e = 3,               /*!< \f$ \tau^+\rightarrow \e^+ \f$ */
+  ktau_mu = 4,               /*!< \f$ \tau^+\rightarrow \mu^+ \f$ */
   kUndefinedTauMode = -1,    /*!< Undefined */
 };