TTrackerECALReconModule.cxx

Go to the documentation of this file.

#include <sstream>

#include "TDatum.hxx"
#include "TGeomInfo.hxx"
#include "TPrincipal.h"
#include "TRealDatum.hxx"
#include "TShowerTruthInfo.hxx"
#include "TrackTruthInfo.hxx"

#include "TTrackerECALReconModule.hxx"

ClassImp(ND::TTrackerECALReconModule::TECALReconTrack);
ClassImp(ND::TTrackerECALReconModule::TECALReconShower);
ClassImp(ND::TTrackerECALReconModule::TECALReconObject);
ClassImp(ND::TTrackerECALReconModule::TECALReconCluster);
ClassImp(ND::TTrackerECALReconModule::TECALReconUnmatchedObject);
ClassImp(ND::TTrackerECALReconModule::TECALReconVertexCluster);
ClassImp(ND::TTrackerECALReconModule::TECALReconVertexTrack);
ClassImp(ND::TTrackerECALReconModule::TECALReconVertexCrossing);

#ifdef TTRACKERECALRECONMODULE_USE_EXP
ClassImp(ND::TTrackerECALReconModule::TECALReconHit);
#endif

ClassImp(ND::TTrackerECALReconModule::TECALReconVertexCandidate);

#define CVSTAG \
  "\
  $Name: 7.0$"
#define CVSID \
  "\
  $Id: eventAnalysis TTrackerECALReconModule.cxx,2024/03/20:09:46:11,Alexander_J_Finch,lapw.lancs.ac.uk $"

namespace {
std::ostream& operator<<(std::ostream& os, const TVector3& tv3) {
  return os << "[" << tv3[0] << ", " << tv3[1] << ", " << tv3[2] << "]";
}
std::ostream& operator<<(std::ostream& os, const TLorentzVector& tlv) {
  return os << "[" << tlv[0] << ", " << tlv[1] << ", " << tlv[2] << ", "
            << tlv[3] << "]";
}

#ifdef TTRACKERECALRECONMODULE_USE_EXP

const int MAXHITSTOSAVE = 1000;

/// Pretty prints TECALReconHits
inline std::ostream& operator<<(
    std::ostream& os, const ND::TTrackerECALReconModule::TECALReconHit& hit) {
  return os << "{ Pos: " << hit.Position << ", Charge: " << hit.ChargeDeposit
            << ", GID: " << hit.GeomId << ", IsXYZ: " << hit.IsXYZ
            << " UsedBy: " << hit.UsedByNObjs << " objects }";
}

#endif
}

namespace ND {
namespace TTrackerECALReconModule {

OutputManager::OutputManager(const char* name, const char* title) {
  SetNameTitle(name, title);
  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDoSaveHoughResults = false;
  fDoVertexing = false;
  fDescription = "Tracker ECAL recon output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;
  fIsMC = false;
  fTotal = 0;
  fVertexRecon = NULL;
  // Tree variables
  fReconObjects =
      new TClonesArray("ND::TTrackerECALReconModule::TECALReconObject", 200);
  fUnmatchedObjects = new TClonesArray(
      "ND::TTrackerECALReconModule::TECALReconUnmatchedObject", 200);
  fReconVertexClusters = new TClonesArray(
      "ND::TTrackerECALReconModule::TECALReconVertexCluster", 200);
#ifdef TTRACKERECALRECONMODULE_USE_EXP
  fDoSaveHitInfo = false;
  fUsedHits = new TClonesArray("ND::TTrackerECALReconModule::TECALReconHit",
                               MAXHITSTOSAVE);
#endif
  fDoSaveVertexCandidates = false;
  fVertexCandidates = new TClonesArray(
      "ND::TTrackerECALReconModule::TECALReconVertexCandidate", 100);
  fVertexCandidateTracks = new TClonesArray(
      "ND::TTrackerECALReconModule::TECALReconVertexTrack", 100);
}

OutputManager::~OutputManager() {}

Bool_t OutputManager::Configure(std::string& option) {
#ifdef TTRACKERECALRECONMODULE_USE_EXP
  if (option == "SaveHitInfo") {
    fDoSaveHitInfo = true;
    ND280NamedLog("TTrackerECALReconModule", "Will save hit info");
    return true;
  }
#endif
  if (option == "RunVertexing") {
    fDoVertexing = true;
    ND280NamedLog("TTrackerECALReconModule",
                  "Will re-run vertexing code reconstruction.");
    return true;
  }
  if (option == "SaveTheHoff") {
    fDoSaveHoughResults = true;
    ND280NamedLog("TTrackerECALReconModule",
                  "Will save raw Hough transform results for this run.");
    return true;
  }
  if (option == "SaveVertexCandy") {
    fDoSaveVertexCandidates = true;
    ND280NamedLog("TTrackerECALReconModule",
                  "Will save reconstructed vertex "
                  "candidate info");
    return true;
  }
  return false;
}

Bool_t OutputManager::ProcessFirstEvent(ND::TND280Event& event) {
#ifdef TTRACKERECALRECONMODULE_USE_EXP
  if (fDoSaveHitInfo) {
    ND280NamedLog("TTrackerECALReconModule", "Saving Hit Info for this run.");
  }
#endif
  if (fDoSaveVertexCandidates) {
    ND280NamedLog("TTrackerECALReconModule",
                  "Saving reconstructed vertex "
                  "candidates for this run.");
  }
  if (fDoSaveHoughResults) {
    ND280NamedLog("TTrackerECALReconModule",
                  "Saving raw Hough transform results for this run.");
  }
  if (fDoVertexing) {
    ND280NamedLog("TTrackerECALReconModule", "Running ECal-iso vertexing");

    // Initialise the vertex recon
    fVertexRecon = new ND::TECalVertexing();
  }

  return true;
}

bool OutputManager::IsFullEventWorthSaving(ND::TND280Event& event) {
  return fNReconObjects;
}

void OutputManager::InitializeModule() {}
void OutputManager::InitializeBranches() {
  fOutputTree->Branch("NReconObject", &fNReconObjects, "NReconObjects/I",
      fBufferSize)
  ->SetTitle(" Number of  TECALRecon Objects");
  
  fOutputTree->Branch("NReconTrackLike", &fNReconTrackLike, "NReconTrackLike/I",
      fBufferSize)
  ->SetTitle("Number of TrackLike objects ");
  
  fOutputTree->Branch("NReconShowerLike", &fNReconShowerLike,
      "NReconShowerLike/I", fBufferSize)
  ->SetTitle("Number of ShowerLike objects ");
  
  fOutputTree->Branch("NUnmatchedObject", &fNUnmatchedObjects,
      "NUnmatchedObjects/I", fBufferSize)
  ->SetTitle("Number of  Unmatched Objects ");
  
  // A TClonesArray of TECALReconObject objects
  fOutputTree->Branch("ReconObject", &fReconObjects, fBufferSize, fSplitLevel)
  ->SetTitle(" The TECALRecon Objects");
  
  fOutputTree->Branch("UnmatchedObject", &fUnmatchedObjects, fBufferSize,
      fSplitLevel)
  ->SetTitle("The Unmatched Objects ");
  
  if (fDoSaveHoughResults) {
    fOutputTree->Branch("NVertexCluster", &fNReconVertexClusters,
        "NVertexCluster/I", fBufferSize)
    ->SetTitle(" ");
    
    fOutputTree->Branch("VertexCluster", &fReconVertexClusters, fBufferSize,
        fSplitLevel)
    ->SetTitle("The Reconstructed Vertices ");
    
  }

#ifdef TTRACKERECALRECONMODULE_USE_EXP
  fOutputTree->Branch("NUsedHits", &fNUsedHits, "NUsedHits/I", fBufferSize);
  fOutputTree->Branch("UsedHits", &fUsedHits, fBufferSize, fSplitLevel);
#endif

  if (fDoVertexing || fDoSaveVertexCandidates) {
    fOutputTree->Branch("NVertexCandidates", &fNVertexCandidates,
                        "NVertexCandidates/I", fBufferSize);
    fOutputTree->Branch("VertexCandidates", &fVertexCandidates, fBufferSize,
                        fSplitLevel);
    fOutputTree->Branch("NVertexCandidateTracks", &fNVertexCandidateTracks,
                        "NVertexCandidateTracks/I", fBufferSize);
    fOutputTree->Branch("VertexCandidateTracks", &fVertexCandidateTracks,
                        fBufferSize, fSplitLevel);
  }
}

void OutputManager::FillHoughTransformInformation(
    ND::THandle<ND::TAlgorithmResult> ecalFinalResult,
    ND::THandle<ND::TG4TrajectoryContainer> trajectories) 
{
  ND::THandle<ND::TReconObjectContainer> vertices_container =
      ecalFinalResult->GetResultsContainer("ECalHoughResults");
  if (!vertices_container) {
    ND280NamedInfo("TTrackerECALReconModule", "No vertex recon container.");
    return;
  }
  for (ND::TReconObjectContainer::iterator vertIt = vertices_container->begin();
       vertIt != vertices_container->end(); vertIt++) {
    ND::THandle<ND::TReconVertex> top_vertex = (*vertIt);
    if (!top_vertex) {
      continue;
    }

    ND280NamedInfo("TTrackerECALReconModule",
                   "============================================");
    ND280NamedInfo("TTrackerECALReconModule",
                   "======Processing Hough Vertex Cluster("
                       << top_vertex->GetUniqueID() << ")=======");
    ND280NamedInfo("TTrackerECALReconModule",
                   "============================================");

    TECALReconVertex vd;
    if (fDoVertexing) {  // If we are saving the vertex candidates
      vd = fVertexRecon->ProcessHoughResult(top_vertex);

      if (vd.IsValid()) {  // If this is a valid reconstructed vertex candidate
        (void)TECALVertexCandidateFactory(
            vd, top_vertex, fIsMC, (*fVertexCandidates)[fNVertexCandidates++]);

        for (UInt_t trki = 0; trki < vd.NTracks(); ++trki) {
          ND::THandle<ND::TReconPID> trk = vd.Tracks_at(trki);
          if (!trk) {
            continue;
          }
          if (!fSavedVertexTrackIds.count(trk->GetUniqueID())) {
            fSavedVertexTrackIds.insert(
                std::make_pair(trk->GetUniqueID(), fNVertexCandidateTracks));
            TECALReconVertexTrackFactory(
                trk, fIsMC,
                (*fVertexCandidateTracks)[fNVertexCandidateTracks++],
                trajectories);
            ND280NamedVerbose(
                "TTrackerECALReconModule",
                "Saved track constituent of vertex candidate with ID: "
                    << trk->GetUniqueID());
          }
        }

        // Flatten and add any secondary vertex candidates in this cluster.
        for (std::vector<TECALReconVertex>::iterator rv_it =
                 vd.SecondaryVertices_begin();
             rv_it != vd.SecondaryVertices_end(); ++rv_it) {
          (void)TECALVertexCandidateFactory(
              (*rv_it), top_vertex, fIsMC,
              (*fVertexCandidates)[fNVertexCandidates++]);

          for (UInt_t trki = 0; trki < (*rv_it).NTracks(); ++trki) {
            ND::THandle<ND::TReconPID> trk = (*rv_it).Tracks_at(trki);
            if (!trk) {
              continue;
            }
            if (!fSavedVertexTrackIds.count(trk->GetUniqueID())) {
              fSavedVertexTrackIds.insert(
                  std::make_pair(trk->GetUniqueID(), fNVertexCandidateTracks));
              TECALReconVertexTrackFactory(
                  trk, fIsMC,
                  (*fVertexCandidateTracks)[fNVertexCandidateTracks++],
                  trajectories);
              ND280NamedVerbose("TTrackerECALReconModule",
                                "Saved track constituent of secondary vertex "
                                "candidate with ID: "
                                    << trk->GetUniqueID());
            }
          }
        }
      }
    }  // end Do vertexing

    //**************************************************************************
    // Do Hough
    //**************************************************************************
    // Now organise the constituents (tracks and crossings) into
    // TReconObjectContainers
    ND::THandle<ND::TReconObjectContainer> constituents =
        top_vertex->GetConstituents();
    double ClusterTotalCharge = 0;
    for (ND::THitSelection::iterator hitIt = top_vertex->GetHits()->begin();
         hitIt != top_vertex->GetHits()->end(); hitIt++) {
      ClusterTotalCharge += (*hitIt)->GetCharge();
    }
    if (!constituents) {
      // There was not enough information in the 2 views to reconstruct
      // anything useful. Record NTracks and NCrossings as 0 and move on
      continue;
    }
    ND::THandle<ND::TReconObjectContainer> tracks(
        new ND::TReconObjectContainer("vertexRecon_tracks"));
    ND::THandle<ND::TReconObjectContainer> crossings(
        new ND::TReconObjectContainer("vertexRecon_crossings"));
    for (ND::TReconObjectContainer::iterator rocIt = constituents->begin();
         rocIt != constituents->end(); rocIt++) {
      ND::THandle<ND::TReconPID> track = (*rocIt);
      ND::THandle<ND::TReconVertex> vertex = (*rocIt);
      if (track) {
        tracks->push_back(track);
      } else if (vertex) {
        crossings->push_back(vertex);
      } else {
        ND280NamedError("TTrackerECALReconModule",
                        "Found a VertexCluster "
                        "constituent that is not a track or vertex :(");
      }
    }

    TECALReconVertexCluster* vertex_cluster = NULL;
    if (fDoSaveHoughResults) {
      // Create the TClonesArray member
      vertex_cluster = new ((*fReconVertexClusters)[fNReconVertexClusters++])
          TECALReconVertexCluster;
      // Now fill it from the top level 'vertex'
      vertex_cluster->UniqueID = top_vertex->GetUniqueID();
      vertex_cluster->NHits = top_vertex->GetHits()->size();
      // Calculate the total charge
      vertex_cluster->TotalCharge = ClusterTotalCharge;
      for (ND::THitSelection::iterator hitIt = top_vertex->GetHits()->begin();
           hitIt != top_vertex->GetHits()->end(); hitIt++) {
#ifdef TTRACKERECALRECONMODULE_USE_EXP
        SaveUsedHit(*hitIt, top_vertex->GetUniqueID());
#endif
      }
      vertex_cluster->Position = top_vertex->GetPosition();
      vertex_cluster->Module = ND::TGeomInfo::Get()
                                   .ECAL()
                                   .GetModule(top_vertex->GetHits()->front())
                                   .GetName();
      vertex_cluster->N2DTracksMax =
          top_vertex->Get<ND::TIntegerDatum>("N2DProngsMax")->GetValue();
      vertex_cluster->N2DTracksMin =
          top_vertex->Get<ND::TIntegerDatum>("N2DProngsMin")->GetValue();
      // Calculate the hit level information
      // Use sets to easily organise the hit information
      std::set<Int_t> layer_numbers, para_bar_numbers, perp_bar_numbers;
      // Grab the TGeomInfo ECal geometry now as we are going to make repeated
      // calls to it.
      ND::TECALGeom* ecal_geom = &(ND::TGeomInfo::Get().ECAL());
      for (ND::THitSelection::iterator hitIt = top_vertex->GetHits()->begin();
           hitIt != top_vertex->GetHits()->end(); hitIt++) {
        ND::THandle<ND::THit> hit = (*hitIt);
        Int_t layer_number = ecal_geom->GetLayerNumber(hit);
        Int_t bar_number = ecal_geom->GetBarNumber(hit);

        // insert the layer number into the set
        layer_numbers.insert(layer_number);
        // Do the same for the bars
        if (hit->IsZHit())
          perp_bar_numbers.insert(bar_number);
        else
          para_bar_numbers.insert(bar_number);
      }
      // The bar and layer numbers should be automatically organised
      // Remember to check that the sets have entries
      if (layer_numbers.size() == 0) {
        vertex_cluster->MinLayerHit = -9999;
        vertex_cluster->MaxLayerHit = -9999;
      } else {
        std::set<Int_t>::iterator setIt = layer_numbers.begin();
        vertex_cluster->MinLayerHit = (*setIt);
        setIt = layer_numbers.end();
        setIt--;
        vertex_cluster->MaxLayerHit = (*setIt);
      }
      if (perp_bar_numbers.size() == 0) {
        vertex_cluster->MinPerpBarHit = -9999;
        vertex_cluster->MaxPerpBarHit = -9999;
      } else {
        std::set<Int_t>::iterator setIt = perp_bar_numbers.begin();
        vertex_cluster->MinPerpBarHit = (*setIt);
        setIt = perp_bar_numbers.end();
        setIt--;
        vertex_cluster->MaxPerpBarHit = (*setIt);
      }
      if (para_bar_numbers.size() == 0) {
        vertex_cluster->MinParaBarHit = -9999;
        vertex_cluster->MaxParaBarHit = -9999;
      } else {
        std::set<Int_t>::iterator setIt = para_bar_numbers.begin();
        vertex_cluster->MinParaBarHit = (*setIt);
        setIt = para_bar_numbers.end();
        setIt--;
        vertex_cluster->MaxParaBarHit = (*setIt);
      }

      // Get the truth information for the cluster
      ND::TShowerTruthInfo shower_truth(top_vertex->GetHits(),
                                        ND::TShowerTruthInfo::kXYZ, 1, 1, 1);
      // Start with single trajectories
      std::vector<ND::TShowerTruthInfo::TruthInfo> results =
          shower_truth.GetSingleTrajectories();
      if (results.size() > 0) {
        vertex_cluster->G4IDSingle = results[0].ParticleID;
        vertex_cluster->CompletenessSingle = results[0].Completeness_NumHits;
        vertex_cluster->CleanlinessSingle = results[0].Cleanliness_NumHits;
      } else {
        vertex_cluster->G4IDSingle = -1;
        vertex_cluster->CompletenessSingle = -1.;
        vertex_cluster->CleanlinessSingle = -1.;
      }
      // Now the primary trajectories
      results = shower_truth.GetPrimaryTrajectories();
      if (results.size() > 0) {
        vertex_cluster->G4IDPrimary = results[0].ParticleID;
        vertex_cluster->CompletenessPrimary = results[0].Completeness_NumHits;
        vertex_cluster->CleanlinessPrimary = results[0].Cleanliness_NumHits;
      } else {
        vertex_cluster->G4IDPrimary = -1;
        vertex_cluster->CompletenessPrimary = -1.;
        vertex_cluster->CleanlinessPrimary = -1.;
      }
      // Finally the recursive trajectories
      results = shower_truth.GetRecursiveTrajectories();
      if (results.size() > 0) {
        vertex_cluster->G4IDRecursive = results[0].ParticleID;
        vertex_cluster->CompletenessRecursive = results[0].Completeness_NumHits;
        vertex_cluster->CleanlinessRecursive = results[0].Cleanliness_NumHits;
      } else {
        vertex_cluster->G4IDRecursive = -1;
        vertex_cluster->CompletenessRecursive = -1.;
        vertex_cluster->CleanlinessRecursive = -1.;
      }

      // Initialise the number of tracks and crossings to 0
      vertex_cluster->NTracks = 0;
      vertex_cluster->NCrossings = 0;

      // Now process the crossings
      vertex_cluster->NCrossings = 0;

      for (ND::TReconObjectContainer::iterator crossIt = crossings->begin();
           crossIt != crossings->end(); crossIt++) {
        ND::THandle<ND::TReconVertex> input_crossing = (*crossIt);
        TECALReconVertexCrossing* vertex_crossing =
            new (((*vertex_cluster->Crossings))[vertex_cluster->NCrossings++])
                TECALReconVertexCrossing;

        vertex_crossing->UniqueID = input_crossing->GetUniqueID();
        Int_t NHits = input_crossing->GetHits()->size();
        vertex_crossing->NHits = NHits;
        vertex_crossing->HitCharges = new Float_t[NHits];
        vertex_crossing->HitPositionsX = new Float_t[NHits];
        vertex_crossing->HitPositionsY = new Float_t[NHits];
        vertex_crossing->HitPositionsZ = new Float_t[NHits];
        vertex_crossing->HitTimes = new Float_t[NHits];
        vertex_crossing->HitLayerNumbers = new Int_t[NHits];
        vertex_crossing->HitBarNumbers = new Int_t[NHits];

        for (ND::THitSelection::iterator hitIt =
                 input_crossing->GetHits()->begin();
             hitIt != input_crossing->GetHits()->end(); hitIt++) {
          int index = std::distance(input_crossing->GetHits()->begin(), hitIt);
          vertex_crossing->HitCharges[index] = (*hitIt)->GetCharge();
          TVector3 hit_position = (*hitIt)->GetPosition();

          vertex_crossing->HitPositionsX[index] = hit_position.X();
          vertex_crossing->HitPositionsY[index] = hit_position.Y();
          vertex_crossing->HitPositionsZ[index] = hit_position.Z();
          vertex_crossing->HitTimes[index] = (*hitIt)->GetTime();
          vertex_crossing->HitLayerNumbers[index] =
              ecal_geom->GetLayerNumber((*hitIt));
          vertex_crossing->HitBarNumbers[index] =
              ecal_geom->GetBarNumber((*hitIt));

          // Now get the track IDs which created this crossing
          vertex_crossing->NTracks = input_crossing->GetConstituents()->size();
          vertex_crossing->TrackIDs = new UInt_t[vertex_crossing->NTracks];

          for (ND::TReconObjectContainer::iterator trackIt =
                   input_crossing->GetConstituents()->begin();
               trackIt != input_crossing->GetConstituents()->end(); trackIt++) {
            int index = std::distance(
                input_crossing->GetConstituents()->begin(), trackIt);
            vertex_crossing->TrackIDs[index] = (*trackIt)->GetUniqueID();
          }
        }
        vertex_crossing->Position = input_crossing->GetPosition();
        vertex_crossing->DOCA =
            input_crossing->Get<ND::TRealDatum>("distance_of_closest_approach")
                ->GetValue();
      }
      // Now start processing the tracks
      vertex_cluster->NTracks = 0;
    }

    double sumAllTracksCharge = 0;

    for (ND::TReconObjectContainer::iterator trackIt = tracks->begin();
         trackIt != tracks->end(); trackIt++) {
      ND::THandle<ND::TReconPID> input_track = (*trackIt);
      // Create the TReconVertexTrack member on the TClonesArray

      if (fDoVertexing) {
        if (!fSavedVertexTrackIds.count(input_track->GetUniqueID())) {
          fSavedVertexTrackIds.insert(std::make_pair(input_track->GetUniqueID(),
                                                     fNVertexCandidateTracks));

          TECALReconVertexTrackFactory(
              input_track, fIsMC,
              (*fVertexCandidateTracks)[fNVertexCandidateTracks++],
              trajectories);
          ND280NamedVerbose("TTrackerECALReconModule",
                            "Saved track constituent of vertex "
                            "cluster with ID: "
                                << input_track->GetUniqueID());
        }
      }  // end save vertexing tracks

      if (fDoSaveHoughResults) {
        TECALReconVertexTrack* vertex_track = NULL;
        if (fDoVertexing) {
          vertex_track = new (
              ((*vertex_cluster->Tracks))[vertex_cluster->NTracks++])
              TECALReconVertexTrack(
                  static_cast<ND::TTrackerECALReconModule::
                                  TECALReconVertexTrack const&>(
                      *fVertexCandidateTracks->At(
                          fSavedVertexTrackIds[input_track->GetUniqueID()])));
        } else {
          TECALReconVertexTrackFactory(
              input_track, fIsMC,
              (*vertex_cluster->Tracks)[vertex_cluster->NTracks++],
              trajectories);

          vertex_track =
              static_cast<ND::TTrackerECALReconModule::TECALReconVertexTrack*>(
                  vertex_cluster->Tracks->At(vertex_cluster->NTracks - 1));
        }

        sumAllTracksCharge += vertex_track->TotalCharge;

// Now we need to get the crossing IDs to associate to the track
// This is a little more complicated than the reverse situation we
// handled in the vertex_crossing section as crossings are NOT stored as
// constituents of the track
// The chosen method is to loop through all of the crossing and see if
// their constituent tracks match the track we are looking at :(
// You can work out how many crossings there are based on how many
// tracks
// are reconstructed but the below method gives you this number
// automatically

#ifdef USE_FULL_VERTEX_TRACK_INFO
        // Store the IDs in a temporary vector
        std::vector<int> tempIDs;
        // Loop over the crossings
        for (ND::TReconObjectContainer::iterator crossIt = crossings->begin();
             crossIt != crossings->end(); crossIt++) {
          ND::THandle<ND::TReconVertex> crossing = (*crossIt);
          // Now loop over the crossing constituents
          for (ND::TReconObjectContainer::iterator trackIt =
                   crossing->GetConstituents()->begin();
               trackIt != crossing->GetConstituents()->end(); trackIt++) {
            ND::THandle<ND::TReconPID> track = (*trackIt);
            // If the constituent unique id matches the track we are currently
            // processing, the crossing was created using the track
            if (track->GetUniqueID() == vertex_track->UniqueID) {
              tempIDs.push_back(crossing->GetUniqueID());  // Store the crossing
                                                           // ID
              break;  // Track won't be more than one constituent of the
                      // crossing
            }
          }
        }
        // Now that we have stored the IDs in a vector, we can easily store the
        // IDs in the TECALReconVertexTrack
        vertex_track->NCrossings = tempIDs.size();
        vertex_track->CrossingIDs = new UInt_t[vertex_track->NCrossings];
        for (unsigned int i = 0; i < tempIDs.size(); i++) {
          vertex_track->CrossingIDs[i] = tempIDs[i];
        }
#endif

        if (sumAllTracksCharge >= vertex_cluster->TotalCharge) {
          ND280NamedWarn(
              "TTrackerECALReconModule",
              "Vertex cluster charge: " << vertex_cluster->TotalCharge
                                        << ", is less than the sum of the "
                                           "constituent track charge: "
                                        << sumAllTracksCharge);
          vertex_cluster->TotalCharge = sumAllTracksCharge;
        }
      }  // end if save hough result

      if (fDoVertexing) {
        for (std::map<UInt_t, int>::iterator ts_it = vd.tracksStatus.begin();
             ts_it != vd.tracksStatus.end(); ++ts_it) {
          std::pair<UInt_t, int> const& trkSt = (*ts_it);

          /// This track was an intermediate object created during prong
          /// merging.
          if (trkSt.second == TECalVertexing::kReMergedTr) {
            continue;
          }

          if (fSavedVertexTrackIds.count(trkSt.first)) {
            ND::TTrackerECALReconModule::TECALReconVertexTrack& VT =
                static_cast<
                    ND::TTrackerECALReconModule::TECALReconVertexTrack&>(
                    *fVertexCandidateTracks->At(
                        fSavedVertexTrackIds[trkSt.first]));
            VT.Status = trkSt.second;
            ND280NamedVerbose("TTrackerECALReconModule",
                              "Assigning saved track status: " << trkSt.second);

          } else {
            ND280NamedSevere(
                "TTrackerECALReconModule",
                "Failed to find vertex track with ID: " << trkSt.first);
            throw;
          }
        }
      }  // end fDoVertexing
    }    // end vertex tracks loop
  }      // end vertex loop
}

int GetTIntegerDatumValue(ND::THandle<ND::TReconBase> trb, char const* name,
                          int const& defaultValue) {
  if (!trb->Get<ND::TIntegerDatum>(name)) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Could not retrieve TIntegerDatum: \"" << name << "\".");
    return defaultValue;
  }
  return trb->Get<ND::TIntegerDatum>(name)->GetValue();
}

void OutputManager::FillECalIsoVertexingInformation(
    ND::THandle<ND::TAlgorithmResult> ecalFinalResult,
    ND::THandle<ND::TG4TrajectoryContainer> trajectories) {
  ND::THandle<ND::TReconObjectContainer> ECalVertices =
      ecalFinalResult->GetResultsContainer("ECalVertexingResults");

  if ((!ECalVertices) || (!ECalVertices->size())) {
    return;
  }

  for (size_t vd_it = 0; vd_it < ECalVertices->size(); ++vd_it) {
    ND::THandle<ND::TReconVertex> ECalVtx = ECalVertices->at(vd_it);

    if (!ECalVtx->GetConstituents() || !ECalVtx->GetConstituents()->size()) {
      return;
    }

    UInt_t VCUID = GetTIntegerDatumValue(ECalVtx, "OVCUID", -1);

    (void)TECALVertexCandidateFactory(
        ECalVtx, VCUID, fIsMC, (*fVertexCandidates)[fNVertexCandidates++]);

    for (UInt_t const_it = 0; const_it < ECalVtx->GetConstituents()->size();
         ++const_it) {
      ND::THandle<ND::TReconPID> trk = ECalVtx->GetConstituents()->at(const_it);
      if (trk) {
        if (!fSavedVertexTrackIds.count(trk->GetUniqueID())) {
          fSavedVertexTrackIds.insert(
              std::make_pair(trk->GetUniqueID(), fNVertexCandidateTracks));
          TECALReconVertexTrackFactory(
              trk, fIsMC, (*fVertexCandidateTracks)[fNVertexCandidateTracks++],
              trajectories);
          ND280NamedVerbose(
              "TTrackerECALReconModule",
              "Saved track constituent of vertex candidate with ID: "
                  << trk->GetUniqueID());
        }
        continue;
      }
      ND::THandle<ND::TReconVertex> vtx =
          ECalVtx->GetConstituents()->at(const_it);
      if (vtx) {
        (void)TECALVertexCandidateFactory(
            ECalVtx, VCUID, fIsMC, (*fVertexCandidates)[fNVertexCandidates++]);

        for (UInt_t const_it = 0; const_it < ECalVtx->GetConstituents()->size();
             ++const_it) {
          ND::THandle<ND::TReconPID> trk =
              ECalVtx->GetConstituents()->at(const_it);
          if (trk) {
            if (!fSavedVertexTrackIds.count(trk->GetUniqueID())) {
              fSavedVertexTrackIds.insert(
                  std::make_pair(trk->GetUniqueID(), fNVertexCandidateTracks));
              TECALReconVertexTrackFactory(
                  trk, fIsMC,
                  (*fVertexCandidateTracks)[fNVertexCandidateTracks++],
                  trajectories);
              ND280NamedVerbose(
                  "TTrackerECALReconModule",
                  "Saved track constituent of vertex candidate with ID: "
                      << trk->GetUniqueID());
            }
          }
        }
      }  // end if secondary vertex
    }
  }
}

void OutputManager::FillVertexingInformation(
    ND::THandle<ND::TAlgorithmResult> ecalFinalResult,
    ND::THandle<ND::TG4TrajectoryContainer> trajectories) {
  if (fDoSaveHoughResults || fDoVertexing) {
    FillHoughTransformInformation(ecalFinalResult, trajectories);
  }
  if (fDoSaveVertexCandidates) {
    FillECalIsoVertexingInformation(ecalFinalResult, trajectories);
  }
}

TECALReconUnmatchedObject* TECALReconUnmatchedObjectFactory(
    ND::THandle<ND::TReconCluster> cluster, bool IsMC, void* pos) {
  if (!cluster) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Attempting to fill "
                   "TECALReconUnmatchedObject from bad TReconCluster.");
    return NULL;
  }
  TECALReconUnmatchedObject* unmatchedObject = NULL;
  if (pos != NULL) {
    unmatchedObject = new (pos) TECALReconUnmatchedObject;
  } else {
    unmatchedObject = new TECALReconUnmatchedObject;
  }

  FillHitInfo(unmatchedObject, cluster->GetHits());

  if (IsMC) {  // Fill the truth-matching info.
    ND::TShowerTruthInfo showerTruth(cluster);

    std::vector<ND::TShowerTruthInfo::TruthInfo> results =
        showerTruth.GetPrimaryTrajectories();
    if (results.size()) {
      unmatchedObject->G4ID_Primary = results[0].ParticleID;
    }

    results = showerTruth.GetRecursiveTrajectories();
    if (results.size()) {
      unmatchedObject->G4ID_Recursive = results[0].ParticleID;
    }

    results = showerTruth.GetSingleTrajectories();
    if (results.size()) {
      unmatchedObject->G4ID_Single = results[0].ParticleID;
    }
  }
  return unmatchedObject;
}

void OutputManager::ProcessOtherContainer(
    ND::THandle<ND::TReconObjectContainer> other) {
  if (!other) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Recieved bad \"other\" container.");
    return;
  }

  for (ND::TReconObjectContainer::iterator cl_it = other->begin();
       cl_it != other->end(); ++cl_it) {
    ND::THandle<ND::TReconCluster> clust = (*cl_it);

    if (!clust) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "Found an invalid "
                     "TReconCluster in the \"other\" container.");
      continue;
    }

    if (clust->GetHits() && clust->GetHits()->size()) {
      TECALReconUnmatchedObjectFactory(
          clust, fIsMC, (*fUnmatchedObjects)[fNUnmatchedObjects++]);
    }
  }
}

TECALReconObject* TECALReconObjectFactory(ND::THandle<ND::TReconTrack> track,
                                          ND::THandle<ND::TReconPID> pid,
                                          bool IsMC, void* pos, bool Allocate) {
  ND280NamedInfo("TTrackerECALReconModule",
                 "Found " << track->GetConstituents()->size()
                          << " PID-track-constituent cluster constituent(s).");

  TECALReconObject* obj = NULL;
  if (Allocate) {
    if (pos != NULL) {
      obj = new (pos) TECALReconObject;
    } else {
      obj = new TECALReconObject;
    }
  } else {
    obj = static_cast<TECALReconObject*>(pos);
  }

  obj->FilledAsTrack = true;
  obj->FilledAsShower = false;

  /// This is to identically emulate previous functionality, but it seems
  /// a bit odd, if this is ever more than one then info has been thrown away
  /// without any regard to which set is better.
  int NReconClusters = 0;
  for (ND::TReconObjectContainer::iterator constit_clus_tit =
           track->GetConstituents()->begin();
       constit_clus_tit != track->GetConstituents()->end();
       ++constit_clus_tit) {
    ND::THandle<ND::TReconCluster> cluster = (*constit_clus_tit);
    if (!cluster) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "PID-track-constituent cluster constituent was bad.");
      continue;
    }
    obj->Cluster.NDOF = cluster->GetNDOF();
    obj->Cluster.Position = cluster->GetPosition();
    obj->Cluster.PositionVar = cluster->GetPositionVariance();
    obj->Cluster.EDeposit = cluster->GetEDeposit();

    NReconClusters++;
  }

  // If tracks have hits grab some hit lvl information.  It may
  // make sense in the future to remove some of this info when
  // we have a working hit lvl information module.

  // Save the unique object ID to allow matching to global recon objects
  obj->UniqueID = pid->GetUniqueID();

  // Calculate the pointing vector.
  TVector3 Pointing(-9999, -9999, -9999);
  if (track->GetHits()) {
    ND::THandle<ND::THitSelection> hitsel = track->GetHits();

    FillHitInfo(obj, hitsel);
    Pointing = CalculatePointing(hitsel);
    if (IsMC) {
      obj->G4ID = TrackTruthInfo::GetG4TrajIDHits(*hitsel);

      // Fill the truth-matching information from the new truth-matching
      // algorithms.
      ND::TShowerTruthInfo showerTruth(pid);

      std::vector<ND::TShowerTruthInfo::TruthInfo> results =
          showerTruth.GetPrimaryTrajectories();
      if (results.size() > 0) {
        obj->G4ID_Primary = results[0].ParticleID;
        obj->Completeness_Primary = results[0].Completeness_NumHits;
        obj->Cleanliness_Primary = results[0].Cleanliness_NumHits;
      }

      results = showerTruth.GetRecursiveTrajectories();
      if (results.size() > 0) {
        obj->G4ID_Recursive = results[0].ParticleID;
        obj->Completeness_Recursive = results[0].Completeness_NumHits;
        obj->Cleanliness_Recursive = results[0].Cleanliness_NumHits;
      }
      if (results.size() > 1) {
        obj->G4ID_Recursive2 = results[1].ParticleID;
      }
      if (results.size() > 2) {
        obj->G4ID_Recursive3 = results[2].ParticleID;
      }
      if (results.size() > 3) {
        obj->G4ID_Recursive4 = results[3].ParticleID;
      }

      results = showerTruth.GetSingleTrajectories();
      if (results.size() > 0) {
        obj->G4ID_Single = results[0].ParticleID;
        obj->Completeness_Single = results[0].Completeness_NumHits;
        obj->Cleanliness_Single = results[0].Cleanliness_NumHits;
      }
    }
  }
  obj->Pointing = Pointing;

  // Collecting information stored in a TReconTrack..
  obj->Track.Direction = track->GetDirection();
  obj->Track.EDeposit = track->GetEDeposit();
  obj->Track.Position = track->GetPosition();
  obj->Track.PositionVar = track->GetPositionVariance();
  obj->Track.BackPosition =
      FindBackPosition(track->GetHits(), obj->MostDownStreamLayerHit);
  obj->Track.Width = track->GetWidth();

  // Adding the Information stored in Datums in the Track
  // This for example includes parametres used in the Pid and
  // Energy Recon.
  FillDatumInfo(obj, track);

  // Adding properties of the Alternative EMShower Hypothesis
  ND280NamedInfo("TTrackerECALReconModule", "Track PID had "
                                                << pid->GetAlternates().size()
                                                << " alternates.");
  for (ND::TReconObjectContainer::const_iterator altIter =
           pid->GetAlternates().begin();
       altIter != pid->GetAlternates().end(); ++altIter) {
    ND::THandle<ND::TReconPID> altPID = (*altIter);
    if (!altPID) {
      ND280NamedWarn("TTrackerECALReconModule", "Alternate PID is invalid.");
      continue;
    }

    if (!altPID->GetConstituents() || !altPID->GetConstituents()->size()) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "Alternate Shower PID has no constituents.");
      continue;
    }

    for (ND::TReconObjectContainer::iterator alt_constit_it =
             altPID->GetConstituents()->begin();
         alt_constit_it != altPID->GetConstituents()->end(); ++alt_constit_it) {
      ND::THandle<ND::TReconShower> shower_alt = (*alt_constit_it);
      if (!shower_alt || !shower_alt->GetEDeposit()) {
        continue;
      }

      obj->Shower.ConeAngle = shower_alt->GetConeAngle();
      obj->Shower.Direction = shower_alt->GetDirection();
      obj->Shower.EDeposit = shower_alt->GetEDeposit();
      obj->Shower.Position = shower_alt->GetPosition();
      obj->Shower.PositionVar = shower_alt->GetPositionVariance();
      obj->Shower.BackPosition =
          FindBackPosition(shower_alt->GetHits(), obj->MostDownStreamLayerHit);
    }
  }
  return obj;
}

TECALReconObject* TECALReconObjectFactory(ND::THandle<ND::TReconShower> shower,
                                          ND::THandle<ND::TReconPID> pid,
                                          bool IsMC, void* pos, bool Allocate) {
  ND280NamedInfo("TTrackerECALReconModule",
                 "Found " << shower->GetConstituents()->size()
                          << " PID-shower-constituent cluster constituent(s).");

  TECALReconObject* obj = NULL;
  if (Allocate) {
    if (pos != NULL) {
      obj = new (pos) TECALReconObject;
    } else {
      obj = new TECALReconObject;
    }
  } else {
    obj = static_cast<TECALReconObject*>(pos);
  }

  obj->FilledAsTrack = false;
  obj->FilledAsShower = true;

  /// This is to identically emulate previous functionality, but it seems
  /// a bit odd, if this is ever more than one then info has been thrown away
  /// without any regard to which set is better.
  int NReconClusters = 0;
  for (ND::TReconObjectContainer::iterator constit_clus_tit =
           shower->GetConstituents()->begin();
       constit_clus_tit != shower->GetConstituents()->end();
       ++constit_clus_tit) {
    ND::THandle<ND::TReconCluster> cluster = (*constit_clus_tit);
    if (!cluster) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "PID-shower-constituent cluster constituent was bad.");
      continue;
    }
    obj->Cluster.NDOF = cluster->GetNDOF();
    obj->Cluster.Position = cluster->GetPosition();
    obj->Cluster.PositionVar = cluster->GetPositionVariance();
    obj->Cluster.EDeposit = cluster->GetEDeposit();

    NReconClusters++;
  }

  // If showers has hits grab some hit lvl information.  It may
  // make sense in the future to remove some of this info when
  // we have a working hit lvl information module.

  // Save the unique object ID to allow matching to global recon objects
  obj->UniqueID = pid->GetUniqueID();

  // Calculate the pointing vector.
  TVector3 Pointing(-9999, -9999, -9999);

  if (shower->GetHits()) {
    ND::THandle<ND::THitSelection> hitsel = shower->GetHits();

    FillHitInfo(obj, hitsel);
    Pointing = CalculatePointing(hitsel);
    if (IsMC) {
      obj->G4ID = TrackTruthInfo::GetG4TrajIDHits(*hitsel);

      // Fill the truth-matching information from the new truth-matching
      // algorithms.
      ND::TShowerTruthInfo showerTruth(pid);

      std::vector<ND::TShowerTruthInfo::TruthInfo> results =
          showerTruth.GetPrimaryTrajectories();
      if (results.size() > 0) {
        obj->G4ID_Primary = results[0].ParticleID;
        obj->Completeness_Primary = results[0].Completeness_NumHits;
        obj->Cleanliness_Primary = results[0].Cleanliness_NumHits;
      }
      results = showerTruth.GetRecursiveTrajectories();
      if (results.size() > 0) {
        obj->G4ID_Recursive = results[0].ParticleID;
        obj->Completeness_Recursive = results[0].Completeness_NumHits;
        obj->Cleanliness_Recursive = results[0].Cleanliness_NumHits;
      }
      if (results.size() > 1) {
        obj->G4ID_Recursive2 = results[1].ParticleID;
      }
      if (results.size() > 2) {
        obj->G4ID_Recursive3 = results[2].ParticleID;
      }
      if (results.size() > 3) {
        obj->G4ID_Recursive4 = results[3].ParticleID;
      }
      results = showerTruth.GetSingleTrajectories();
      if (results.size() > 0) {
        obj->G4ID_Single = results[0].ParticleID;
        obj->Completeness_Single = results[0].Completeness_NumHits;
        obj->Cleanliness_Single = results[0].Cleanliness_NumHits;
      }
    }
  }

  obj->Pointing = Pointing;
  obj->Shower.ConeAngle = shower->GetConeAngle();
  obj->Shower.Direction = shower->GetDirection();
  obj->Shower.EDeposit = shower->GetEDeposit();
  obj->Shower.Position = shower->GetPosition();
  obj->Shower.PositionVar = shower->GetPositionVariance();
  obj->Shower.BackPosition =
      FindBackPosition(shower->GetHits(), obj->MostDownStreamLayerHit);

  // Adding the Information stored in Datums in the Shower
  // This for example includes parametres using in TMVA and Energy
  // fitting.
  FillDatumInfo(obj, shower);

  ND280NamedInfo("TTrackerECALReconModule", "Shower PID had "
                                                << pid->GetAlternates().size()
                                                << " alternates.");
  for (ND::TReconObjectContainer::const_iterator altIter =
           pid->GetAlternates().begin();
       altIter != pid->GetAlternates().end(); ++altIter) {
    ND::THandle<ND::TReconPID> altPID = (*altIter);
    if (!altPID) {
      ND280NamedWarn("TTrackerECALReconModule", "Alternate PID is invalid.");
      continue;
    }

    if (!altPID->GetConstituents() || !altPID->GetConstituents()->size()) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "Alternate Track PID has no constituents.");
      continue;
    }

    for (ND::TReconObjectContainer::iterator alt_constit_it =
             altPID->GetConstituents()->begin();
         alt_constit_it != altPID->GetConstituents()->end(); ++alt_constit_it) {
      ND::THandle<ND::TReconTrack> track_alt = (*alt_constit_it);
      if (!track_alt || !track_alt->GetEDeposit()) {
        continue;
      }

      obj->Track.Direction = track_alt->GetDirection();
      obj->Track.EDeposit = track_alt->GetEDeposit();
      obj->Track.Position = track_alt->GetPosition();
      obj->Track.PositionVar = track_alt->GetPositionVariance();
      obj->Track.BackPosition =
          FindBackPosition(track_alt->GetHits(), obj->MostDownStreamLayerHit);
      obj->Track.Width = track_alt->GetWidth();
    }
  }
  return obj;
}

void OutputManager::ProcessFinalContainer(
    ND::THandle<ND::TReconObjectContainer> final) {
  if (!final) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Recieved bad \"final\" container.");
    return;
  }

  // Loop looking for TReconPID's in the top recon conatiner
  for (ND::TReconObjectContainer::iterator pid_it = final->begin();
       pid_it != final->end(); ++pid_it) {
    ND::THandle<ND::TReconPID> pid = (*pid_it);

    if (!pid) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "Found bad PID in the "
                     "\"final\" container.");
      continue;
    }
    if (!pid->GetConstituents() || !pid->GetConstituents()->size()) {
      ND280NamedWarn("TTrackerECALReconModule", "PID has no constituents.");
      continue;
    }
    ND280NamedInfo("TTrackerECALReconModule",
                   "PID has " << pid->GetConstituents()->size()
                              << " constituents.");

    for (ND::TReconObjectContainer::iterator constit_it =
             pid->GetConstituents()->begin();
         constit_it != pid->GetConstituents()->end(); ++constit_it) {
      ND::THandle<ND::TReconTrack> track_constit = (*constit_it);
      ND::THandle<ND::TReconShower> shower_constit = (*constit_it);

      // Only the implicit dynamic casts that are valid will result in a new
      // object.
      if (track_constit) {
        if (!track_constit->GetEDeposit()) {
          ND280NamedWarn("TTrackerECALReconModule",
                         "PID-constituent track deposited no energy.");
          continue;
        }
        // Making a TECALReconObject ClonesArray entry...
        TECALReconObjectFactory(track_constit, pid, fIsMC,
                                (*fReconObjects)[fNReconObjects++]);
        fTotal++;
        fNReconTrackLike++;

        ND280NamedInfo("TTrackerECALReconModule",
                       "Added information about a Track PID constituent.");
      } else if (shower_constit) {
        if (!shower_constit->GetEDeposit()) {
          ND280NamedWarn("TTrackerECALReconModule",
                         "PID-constituent shower deposited no energy.");
          continue;
        }

        TECALReconObjectFactory(shower_constit, pid, fIsMC,
                                (*fReconObjects)[fNReconObjects++]);
        fTotal++;
        fNReconShowerLike++;

        ND280NamedInfo("TTrackerECALReconModule",
                       "Added information about a Shower PID constituent.");
      }
    }
  }
}

void OutputManager::PerEvReset() {
  fNReconObjects = 0;
  fNUnmatchedObjects = 0;
  fNReconTrackLike = 0;
  fNReconShowerLike = 0;
  fNReconVertexClusters = 0;

#ifdef TTRACKERECALRECONMODULE_USE_EXP
  fNUsedHits = 0;
  fSavedHitMap.clear();
  fUsedHits->Clear();
#endif
  fNVertexCandidates = 0;
  fNVertexCandidateTracks = 0;
  fSavedVertexTrackIds.clear();
  fVertexCandidates->Clear();
  fVertexCandidateTracks->Clear();
  fReconObjects->Clear();
  fUnmatchedObjects->Clear();
  fReconVertexClusters->Clear();
  ND280NamedInfo("TTrackerECALReconModule", "===============================");
  ND280NamedInfo("TTrackerECALReconModule", "=========Clearing Event========");
  ND280NamedInfo("TTrackerECALReconModule", "===============================");
}

bool OutputManager::FillTree(ND::TND280Event& event) {
  ND280NamedInfo("TTrackerECALReconModule",
                 "Processing event: " << event.GetContext().GetEvent());
  PerEvReset();

  // Getting the final Algorithm Result from ecalRecon.
  ND::THandle<ND::TAlgorithmResult> ecalReconResult = event.GetFit("ecalRecon");

  if (!ecalReconResult) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Couldn't find the ecalRecon algorithm result.");
    return true;
  }
  ND::THandle<ND::TReconObjectContainer> ecalFinal =
      ecalReconResult->GetResultsContainer("final");

  if (!ecalFinal) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Couldn't find a valid ecalRecon \"final\" container");
  } else {
    ProcessFinalContainer(ecalFinal);
  }
  ND::THandle<ND::TReconObjectContainer> ecalOther =
      ecalReconResult->GetResultsContainer("other");
  if (!ecalOther) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Couldn't find a valid ecalRecon \"other\" container");
  } else {
    ProcessOtherContainer(ecalOther);
  }

  ND::THandle<ND::TG4TrajectoryContainer> trajectories(NULL);
  if (fIsMC) {
    trajectories =
        event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
  }

  FillVertexingInformation(ecalReconResult, trajectories);
  return true;
}

TECALReconObject::TECALReconObject() {
  G4ID_Primary = -1;
  Completeness_Primary = -1;
  Cleanliness_Primary = -1;
  G4ID_Recursive = -1;
  Completeness_Recursive = -1;
  Cleanliness_Recursive = -1;
  G4ID_Single = -1;
  Completeness_Single = -1;
  Cleanliness_Single = -1;
}

void FillHitInfo(TECALReconObject* obj, ND::THandle<ND::THitSelection> hitsel) {
  // Setting hit lvl quantities to appropriate start values
  obj->MostUpStreamLayerHit = 36;
  obj->MostDownStreamLayerHit = -1;
  obj->MaxPerpBarHit = -1;
  obj->MaxParaBarHit = -1;
  obj->MinBarHit = 100;
  obj->NLayersHit = 0;
  obj->NHits = 0;
  obj->TotalHitCharge = 0;
  obj->AverageHitTime = 0;
  obj->NIsXHits = 0;
  obj->NIsYHits = 0;
  obj->NIsZHits = 0;
  obj->ObjectLength = 0.0;
  obj->AverageZPosition = 0.0;

  std::string module;

  bool layer[36] = {false};
  std::vector<TVector3> layerHitPos(36);
  std::vector<double> layerCharge(36, 0.0);

  // Vector to store the layer for each hit
  std::vector<int> layerhits;
  layerhits.reserve(hitsel->size());
  double HitCharge = -1;

  // Looping through hits and filling the Hit info
  for (ND::THitSelection::iterator hit_it = hitsel->begin();
       hit_it != hitsel->end(); ++hit_it) {
    ND::THandle<ND::TReconHit> reconHit = (*hit_it);

    if (!reconHit) {
      ND280NamedWarn("TTrackerECALReconModule", "TReconHit is invalid.");
      continue;
    }

    TVector3 HitPos = reconHit->GetContributor(0)->GetPosition();
    ND::TGeometryId HitGID = reconHit->GetContributor(0)->GetGeomId();

    Int_t layerHit = ND::TGeomInfo::Get().ECAL().GetLayerNumber(HitGID);
    Int_t barHit = ND::TGeomInfo::Get().ECAL().GetBarNumber(HitGID);

    layerhits.push_back(layerHit);

    if (hit_it == hitsel->begin()) {
      module = ND::TGeomInfo::Get().ECAL().GetModule(HitGID).GetName();
    }

    // Average Z Position
    obj->AverageZPosition += HitPos.Z();

    // Charge weighted layer positon
    TVector3 chargePos = reconHit->GetPosition() * reconHit->GetCharge();
    layerHitPos[layerHit] += chargePos;
    layerCharge[layerHit] += reconHit->GetCharge();

    // Using the array to collect the number of layers hit.
    if (!layer[layerHit]) {
      layer[layerHit] = true;
      obj->NLayersHit++;
    }
    // Filling numbers on hits in different types of bars.
    if (reconHit->GetContributor(0)->IsXHit()) {
      obj->NIsXHits++;
    }
    if (reconHit->GetContributor(0)->IsYHit()) {
      obj->NIsYHits++;
    }
    if (reconHit->GetContributor(0)->IsZHit()) {
      obj->NIsZHits++;
    }
    // Obtaining the max and min bars hit
    if (barHit < obj->MinBarHit) {
      obj->MinBarHit = barHit;
    }
    if ((barHit > obj->MaxParaBarHit) &&
        (!reconHit->GetContributor(0)->IsZHit())) {
      obj->MaxParaBarHit = barHit;
    }
    if ((barHit > obj->MaxPerpBarHit) &&
        reconHit->GetContributor(0)->IsZHit()) {
      obj->MaxPerpBarHit = barHit;
    }
    // Obtaining the first and last layers hit.
    if (layerHit < obj->MostUpStreamLayerHit) {
      obj->MostUpStreamLayerHit = layerHit;
    }
    if (layerHit > obj->MostDownStreamLayerHit) {
      obj->MostDownStreamLayerHit = layerHit;
    }
    // Collecting other useful hit information
    obj->NHits++;
    obj->TotalHitCharge += reconHit->GetCharge();
    obj->AverageHitTime += (reconHit->GetTime() * reconHit->GetCharge());
    // Save the layer with the highest hit charge
    if (reconHit->GetCharge() > HitCharge) {
      HitCharge = reconHit->GetCharge();
      obj->MaxHitChargeLayer = (layerHit + 1);
    }
  }  // End hit loop

  obj->Module = module;
  obj->AverageHitTime /= obj->TotalHitCharge;
  obj->AverageZPosition /= obj->NHits;

  // Charge weighted layer position
  for (int i = 0; i < 36; ++i) {
    if (layerCharge[i] != 0.0) {
      layerHitPos[i] *= (1 / layerCharge[i]);
    }
  }

  // Now calculate object length
  TVector3 objLength = (layerHitPos[obj->MostDownStreamLayerHit] -
                        layerHitPos[obj->MostUpStreamLayerHit]);
  obj->ObjectLength = objLength.Mag();

  // Sort the layer for each hit by increasing order
  sort(layerhits.begin(), layerhits.end());
  // Temporary variables to be used
  // used to calculate the number of hits in the same layer.
  // This is used to count the number of layers with 1,2,3,4,5,6 or 7 hits
  int seclayer = -1;
  // used to calculate the exact number of hits in the same layer.
  // This is used to know the number of hits in a layer with many hits
  int doublelayer = -1;
  // used as a flag to count the first,second,third and last layer with many
  // hits
  int doublelayersec = -1;
  // counts the first layer with two or more hits
  int layertwohits = -1;
  // counts the last layer with two or more hits
  int lastlayermanyhits = -1;
  // used as a flag to find when there are more than one hit in a layer
  int layerA = -1;
  // used as a flag to stop multiple counting of first,second third and last
  // layer with many hits
  int layerB = -1;
  // counts the number of layers with at least two hits
  int NLayers2Hits = -1;
  // counts the number of layers with at least three hits
  int NLayers3Hits = -1;
  // the number of hits in the layer that received the most hits
  int maxhits = -1;

  // loop over all the ordered hits
  for (unsigned int i = 0; i < layerhits.size(); i++) {
    int layer = layerhits.at(i);

    // this hit is in a different layer to the previous hit
    if (layer > layerA) {
      layerA = layer;
      seclayer = -1;
    } else {  // found a hit in the same layer as the previous hit
      // count the number of hits in this layer
      seclayer++;
      // do the same again (this time so we know exactly how many hits there
      // were in this layer
      doublelayer++;
      lastlayermanyhits = layer;
      // Now we just want to know if this is the first,second or third layer
      // which had two or more hits
      if (layer > layerB) {
        layerB = layer;
        doublelayersec++;
        if (doublelayersec == 0) {
          layertwohits = layer;
        }
      }

      // see if this layer has had the most hits
      if (seclayer > maxhits) {
        maxhits = seclayer;
      }

      // see if this layer has at least 2,3,4,5,6 or 7 hits
      if (seclayer == 0) {
        NLayers2Hits++;
      } else if (seclayer == 1) {
        NLayers3Hits++;
      }
    }
  }

  // Save the variables - The first layer is always starting from 1
  obj->NHitsAtLayersWithManyHits = (doublelayer + 1);
  obj->NLayersTwoHits = (NLayers2Hits + 1);
  obj->NLayersThreeHits = (NLayers3Hits + 1);
  obj->FirstLayerManyHits = (layertwohits + 1);
  obj->LastLayerManyHits = (lastlayermanyhits + 1);
  obj->MaxHitsInALayer = (maxhits + 1);
}

double GetTRealDatumValue(ND::THandle<ND::TReconBase> trb, char const* name,
                          double const& defaultValue) {
  if (!trb->Get<ND::TRealDatum>(name)) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Could not retrieve TRealDatum: \"" << name << "\".");
    return defaultValue;
  }
  return trb->Get<ND::TRealDatum>(name)->GetValue();
}

void FillDatumInfo(TECALReconObject* obj, ND::THandle<ND::TReconBase> trb) {
  if (!trb) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "TReconBase was invalid in FillDatumInfo.");
    return;
  }

  // =========================Pid Vars==========================
  //                         ----------

  // Circularity
  // Combined and in each view
  obj->PID_Circularity = GetTRealDatumValue(trb, "Circularity", -1);

  // Containment
  // This isn't a PID variable, but it's easy to save it the same way.
  obj->Containment = GetTRealDatumValue(trb, "Containment", -1);

  obj->PID_Angle = GetTRealDatumValue(trb, "Angle", -1);

  obj->PID_TruncatedMaxRatio = GetTRealDatumValue(trb, "TruncatedMaxRatio", -1);

  obj->PID_TransverseChargeRatio =
      GetTRealDatumValue(trb, "TransverseChargeRatio", -1);

  // FrontBackRatio. Note - Assumes particle is coming from the
  // tracker into the ECal
  if (trb->Get<ND::TRealDatum>("FrontBackRatio")) {
    obj->PID_FrontBackRatio =
        trb->Get<ND::TRealDatum>("FrontBackRatio")->GetValue();
  } else {
    ND::THandle<ND::THitSelection> hits = trb->Get<ND::THitSelection>("hits");
    double FrontBackRatio = -1;
    if (hits) {
      FrontBackRatio = CalculateFrontBackRatio(hits);
    }
    obj->PID_FrontBackRatio = FrontBackRatio;
  }

  obj->PID_ShowerAngle = GetTRealDatumValue(trb, "ShowerAngle", -1);

  obj->PID_Asymmetry = GetTRealDatumValue(trb, "Asymmetry", -1);

  obj->PID_LLR_Quality = GetTRealDatumValue(trb, "LLR_Quality", -9999.0);  
  
  obj->PID_LLR_MIP_EM = GetTRealDatumValue(trb, "LLR_MIP_EM", -9999.0);

  obj->PID_LLR_MIP_Pion = GetTRealDatumValue(trb, "LLR_MIP_Pion", -9999.0);

  obj->PID_LLR_EM_HIP = GetTRealDatumValue(trb, "LLR_EM_HIP", -9999.0);

  obj->PID_LLR_MIP_EM_LowMomentum =
      GetTRealDatumValue(trb, "LLR_MIP_EM_LowMomentum", -9999.0);

  ND280NamedInfo("TTrackerECALReconModule",
                 "TReconBase has " << trb->GetConstituents()->size()
                                   << " constituents.");
  for (ND::TReconObjectContainer::iterator trb_it =
           trb->GetConstituents()->begin();
       trb_it != trb->GetConstituents()->end(); ++trb_it) {
    ND::THandle<ND::TReconCluster> cluster_constit = (*trb_it);
    if (!cluster_constit) {
      ND280NamedWarn("TTrackerECALReconModule",
                     "TReconBase cluster constituent"
                     " is invalid.");
      continue;
    }

    // ====================3D Clustering Variables========================
    //                    ------------------

    // 3D Clustering Seed Type
    if (cluster_constit->Get<ND::TIntegerDatum>("3DSeedType")) {
      obj->Clustering_3DSeedType =
          cluster_constit->Get<ND::TIntegerDatum>("3DSeedType")->GetValue();
    } else {
      obj->Clustering_3DSeedType = -1;  // not a 3D cluster, sed to minus 1.
    }

    // ====================Energy Variables========================
    //                    ------------------

    // EM Energy Fit Result
    if (cluster_constit->Get<ND::TRealDatum>("EMEnergyFitResult_Corrected")) {
      obj->EMEnergyFit_Result =
          (*cluster_constit->Get<ND::TRealDatum>("EMEnergyFitResult_Corrected")
                ->begin());
    } else {
      obj->EMEnergyFit_Result = -1;
    }

    obj->EMEnergyFit_Uncertainty =
        GetTRealDatumValue(trb, "EMEnergyFitUncertainty_Corrected", -1);

    // EM Energy Fit Likelihoods
    ND::THandle<ND::TRealDatum> emefl =
        cluster_constit->Get<ND::TRealDatum>("EMEnergyFitLikelihood");
    if (emefl && (emefl->size() >= 3)) {
      obj->EMEnergyFit_Likelihood_energyGrad = emefl->at(0);
      obj->EMEnergyFit_Likelihood_qsum_like = emefl->at(1);
      obj->EMEnergyFit_Likelihood_energy_qsumGrad = emefl->at(2);
    } else {
      obj->EMEnergyFit_Likelihood_energyGrad = -1;
      obj->EMEnergyFit_Likelihood_energy_qsumGrad = -1;
      obj->EMEnergyFit_Likelihood_qsum_like = -1;
    }

    // EM Energy Fit Parameters
    ND::THandle<ND::TRealDatum> emefp =
        cluster_constit->Get<ND::TRealDatum>("EMEnergyFitParameters");
    if (emefp && (emefp->size() >= 5)) {
      obj->EMEnergyFit_Para_QSum = emefp->at(0);
      obj->EMEnergyFit_Para_QMean = emefp->at(1);
      obj->EMEnergyFit_Para_QRMS = emefp->at(2);
      obj->EMEnergyFit_Para_QSkew = emefp->at(3);
      obj->EMEnergyFit_Para_QMax = emefp->at(4);
    } else {
      obj->EMEnergyFit_Para_QSum = -1;
      obj->EMEnergyFit_Para_QMean = -1;
      obj->EMEnergyFit_Para_QRMS = -1;
      obj->EMEnergyFit_Para_QSkew = -1;
      obj->EMEnergyFit_Para_QMax = -1;
    }

    // ====================Thrust Variables========================
    //                    ------------------
    ND::THandle<ND::TRealDatum> thrustVect =
        cluster_constit->Get<ND::TRealDatum>("Thrust");
    if (thrustVect && (thrustVect->size() >= 7)) {
      obj->Thrust = thrustVect->at(0);
      obj->ThrustOrigin.SetXYZ(thrustVect->at(1), thrustVect->at(2),
                               thrustVect->at(3));
      obj->ThrustAxis.SetXYZ(thrustVect->at(4), thrustVect->at(5),
                             thrustVect->at(6));
    } else {
      obj->Thrust = -1.;
    }

    // ====================Matching Likelihood========================
    //                    ------------------

    obj->MatchingLikelihood = GetTRealDatumValue(trb, "Likelihood", -1);
    break;  // Only care about the first (good) cluster constituent.
  }
}

TLorentzVector FindBackPosition(ND::THandle<ND::THitSelection> hits,
                                int MostDownStreamLayerHit) {
  TLorentzVector pos;
  double totalCharge(0);

  // Iterate over hits to get the total charge in the back layer.
  for (ND::THitSelection::iterator hit_it = hits->begin();
       hit_it != hits->end(); ++hit_it) {
    int layer =
        ND::TGeomInfo::Get().ECAL().GetLayerNumber((*hit_it)->GetGeomId());

    if (layer == MostDownStreamLayerHit) {
      totalCharge += (*hit_it)->GetCharge();
    }
  }

  // Iterate over hits again to find a charge weighted position.
  for (ND::THitSelection::iterator hit_it = hits->begin();
       hit_it != hits->end(); ++hit_it) {
    int layer =
        ND::TGeomInfo::Get().ECAL().GetLayerNumber((*hit_it)->GetGeomId());

    if (layer == MostDownStreamLayerHit) {
      TLorentzVector temp((*hit_it)->GetPosition(), (*hit_it)->GetTime());
      temp *= ((*hit_it)->GetCharge() / totalCharge);
      pos += temp;
    }
  }
  return pos;
}

TVector3 CalculatePointing(ND::THandle<ND::THitSelection> hits) {
  TPrincipal pca(3, "ND");

  // Add the charge weighted hits to the pca.
  for (ND::THitSelection::iterator hit_it = hits->begin();
       hit_it != hits->end(); ++hit_it) {
    double row[3] = {(*hit_it)->GetPosition().X(), (*hit_it)->GetPosition().Y(),
                     (*hit_it)->GetPosition().Z()};
    for (double q = 0; q < (*hit_it)->GetCharge(); q += 1.0) {
      pca.AddRow(row);
    }
  }

  // Find the principal components of the shower.
  pca.MakePrincipals();

  // Find the centre of the shower.
  //   The P2X function does this:
  //   Get the mean position of the shower.
  //    for (Int_t i = 0; i < fNumberOfVariables; i++){
  //      x[i] = fMeanValues(i);
  //      for (Int_t j = 0; j < nTest; j++)
  //        x[i] += p[j] * (fIsNormalised ? fSigmas(i) : 1)
  //      * fEigenVectors(i,j);
  //    }
  //   Give p[j] = 0, so actually just return x[i], which is the
  //   mean value over all data points for x, y and z.
  double posP[] = {0., 0., 0.};
  double posX[3];
  pca.P2X(posP, posX, 3);
  TVector3 pcaCentre;
  pcaCentre.SetXYZ(posX[0], posX[1], posX[2]);

  // Set the primary pca axis.
  TVector3 pcaAxis;
  pcaAxis.SetXYZ((*pca.GetEigenVectors())[0][0], (*pca.GetEigenVectors())[1][0],
                 (*pca.GetEigenVectors())[2][0]);

  // Use the PCA info to calcualte pointing.
  double sumcharge = 0;
  double numerator = 0;
  for (ND::THitSelection::iterator hit = hits->begin(); hit != hits->end();
       ++hit) {
    sumcharge += (*hit)->GetCharge();
    numerator +=
        ((*hit)->GetCharge() * pcaAxis.Dot((*hit)->GetPosition() - pcaCentre));
  }

  TVector3 pointing;
  pointing = ((numerator / sumcharge) * pcaAxis);

  return pointing;
}

TECALReconUnmatchedObject::TECALReconUnmatchedObject() {
  G4ID_Primary = -1;
  G4ID_Recursive = -1;
  G4ID_Single = -1;
}

void FillHitInfo(TECALReconUnmatchedObject* obj,
                 ND::THandle<ND::THitSelection> hits) {
  // Initialise the variables.
  obj->NHits = 0;
  obj->View = -1;
  obj->TotalHitCharge = 0.0;
  obj->AverageHitTime = 0.0;
  obj->MostUpStreamLayerHit = 0;
  obj->MostDownStreamLayerHit = 0;

  obj->NHits = hits->size();

  int MinLayer(34);
  int MaxLayer(-1);

  // Loop through the hits.
  for (ND::THitSelection::iterator hit_it = hits->begin();
       hit_it != hits->end(); ++hit_it) {
    int layer =
        ND::TGeomInfo::Get().ECAL().GetLayerNumber((*hit_it)->GetGeomId());

    if (layer < MinLayer) {
      MinLayer = layer;
    }
    if (layer > MaxLayer) {
      MaxLayer = layer;
    }

    obj->TotalHitCharge += (*hit_it)->GetCharge();
    obj->AverageHitTime += ((*hit_it)->GetTime() * (*hit_it)->GetCharge());
  }

  // Finish charge weighting the time.
  obj->AverageHitTime /= obj->TotalHitCharge;

  // Layer Info
  obj->MostUpStreamLayerHit = MinLayer;
  obj->MostDownStreamLayerHit = MaxLayer;

  double MinLayerCharge(0.0);
  double MaxLayerCharge(0.0);

  // Find Charge Weighted Front and Back Positions
  for (ND::THitSelection::iterator hit_it = hits->begin();
       hit_it != hits->end(); ++hit_it) {
    TVector3 hitPos = (*hit_it)->GetPosition();
    int layer =
        ND::TGeomInfo::Get().ECAL().GetLayerNumber((*hit_it)->GetGeomId());

    if (layer == MinLayer) {
      obj->FrontPosition += (hitPos * (*hit_it)->GetCharge());
      MinLayerCharge += (*hit_it)->GetCharge();
    }
    if (layer == MaxLayer) {
      obj->BackPosition += (hitPos * (*hit_it)->GetCharge());
      MaxLayerCharge += (*hit_it)->GetCharge();
    }
  }
  obj->FrontPosition *= (1.0 / MinLayerCharge);
  obj->BackPosition *= (1.0 / MaxLayerCharge);

  // Find the view from one of the layers used.
  if (!(MinLayer % 2)) {
    obj->View = 1;
  } else {
    obj->View = 0;
  }

  obj->G4ID = TrackTruthInfo::GetG4TrajIDHits(*hits);
}

double CalculateFrontBackRatio(ND::THandle<ND::THitSelection> hits) {
  // create a vector to store the hits in
  std::vector<std::pair<TVector3, double> > positionAndCharge;
  double FrontBackRatio = 0;
  if (hits && hits->size()) {
    // determine principal PCA axis
    TPrincipal pca(3, "");

    for (ND::THitSelection::iterator hit_it = hits->begin();
         hit_it != hits->end(); ++hit_it) {
      const TVector3& pos = (*hit_it)->GetPosition();
      double charge = (*hit_it)->GetCharge();
      int repeatEntries = static_cast<int>(charge / 0.1);
      for (int i = 0; i < repeatEntries; ++i) {
        double row[3] = {pos.X(), pos.Y(), pos.Z()};
        pca.AddRow(row);
      }
    }
    pca.MakePrincipals();

    // Rotate hits into PCA space, and fill vector with the new positions
    for (ND::THitSelection::iterator hit_it = hits->begin();
         hit_it != hits->end(); ++hit_it) {
      TVector3 position = (*hit_it)->GetPosition();
      ND::THandle<ND::THit> h = *hit_it;
      double charge = h->GetCharge();

      double x[3] = {position.X(), position.Y(), position.Z()};
      double p[3];
      pca.X2P(x, p);
      // define Z-axis of new co-ordinate system to be the direction along the
      // principal PCA axis
      position.SetZ(p[0]);
      position.SetY(p[1]);
      position.SetX(p[2]);

      std::pair<TVector3, double> p_c(position, charge);
      positionAndCharge.push_back(p_c);
    }

    // define four blocks, first to fourth
    double first(0);
    double second(0);
    double third(0);
    double fourth(0);

    int size = positionAndCharge.size();
    // set initial values for the ends
    double zmin = positionAndCharge.at(0).first.Z();
    double zmax = positionAndCharge.at(0).first.Z();

    // find the actual maximum and minimum "z" values
    //(i.e. the ends of the track)
    for (int i(0); i < size; i++) {
      if (positionAndCharge.at(i).first.Z() > zmax) {
        zmax = positionAndCharge.at(i).first.Z();
      }
      if (positionAndCharge.at(i).first.Z() < zmin) {
        zmin = positionAndCharge.at(i).first.Z();
      }
    }

    // calculate the lengh of the track along the PCA axis
    double length = zmax - zmin;

    // loop through the hits, calculating which block it falls in, and adding
    // its charge to that block
    for (int i = 0; i < size; i++) {
      if ((positionAndCharge.at(i).first.Z() - zmin) > (3 * length / 4)) {
        fourth += positionAndCharge.at(i).second;
      } else if ((positionAndCharge.at(i).first.z() - zmin) > (length / 2)) {
        third += positionAndCharge.at(i).second;
      } else if ((positionAndCharge.at(i).first.z() - zmin) > (length / 4)) {
        second += positionAndCharge.at(i).second;
      } else {
        first += positionAndCharge.at(i).second;
      }
    }

    // if the filling has succeeded, calculate the ratio
    if (first != 0 && fourth != 0) {
      FrontBackRatio = fourth / first;
    }
  }
  return FrontBackRatio;
}

TECALReconVertexCluster::TECALReconVertexCluster()
    : UniqueID(0),
      NHits(0),
      TotalCharge(0),
      Position(0, 0, 0, 0),
      Module(""),
      MinLayerHit(0),
      MaxLayerHit(0),
      MinParaBarHit(0),
      MaxParaBarHit(0),
      MinPerpBarHit(0),
      MaxPerpBarHit(0),
      N2DTracksMax(0),
      N2DTracksMin(0),
      G4IDSingle(0),
      CompletenessSingle(0),
      CleanlinessSingle(0),
      G4IDPrimary(0),
      CompletenessPrimary(0),
      CleanlinessPrimary(0),
      G4IDRecursive(0),
      CompletenessRecursive(0),
      CleanlinessRecursive(0) {
  NTracks = 0;
  NCrossings = 0;
  Tracks = new TClonesArray(
      "ND::TTrackerECALReconModule::TECALReconVertexTrack", 50);
  Crossings = new TClonesArray(
      "ND::TTrackerECALReconModule::TECALReconVertexCrossing", 200);
}

TECALReconVertexTrack::TECALReconVertexTrack()
    : Status(-1),
      UniqueID(0),
      NHits(0),
      MatchingLikelihood(0),
      FrontPosition(0, 0, 0, 0),
      FrontDirection(0, 0, 0),
      FrontLayerNumber(0),
      BackPosition(0, 0, 0, 0),
      BackDirection(0, 0, 0),
      BackLayerNumber(0),
      TotalCharge(0),
      AverageHitTime(0),
      G4IDSingle(0),
#ifdef USE_FULL_VERTEX_TRACK_INFO
      CompletenessSingle(0),
      CleanlinessSingle(0),
      Single_FrontPosition(0, 0, 0),
      Single_FrontDirection(0, 0, 0),
      Single_FrontG4PointDist(0),
      Single_BackPosition(0, 0, 0),
      Single_BackDirection(0, 0, 0),
      Single_BackG4PointDist(0),
      G4IDPrimary(0),
      CompletenessPrimary(0),
      CleanlinessPrimary(0),
      G4IDRecursive(0),
      CompletenessRecursive(0),
      CleanlinessRecursive(0),
#endif
      PID_LLR_MIP_HIP(-9999),
      PID_LLR_MIP_HIP_VACut(-9999) {
#ifdef USE_FULL_VERTEX_TRACK_INFO
  SubTrack_UniqueID[0] = 0;
  SubTrack_UniqueID[1] = 0;
  SubTrack_NHits[0] = 0;
  SubTrack_NHits[1] = 0;
  SubTrack_Angle[0] = 0;
  SubTrack_Angle[1] = 0;
  SubTrack_FrontLayerNumber[0] = 0;
  SubTrack_FrontLayerNumber[1] = 0;
  SubTrack_BackLayerNumber[0] = 0;
  SubTrack_BackLayerNumber[1] = 0;
  SubTrack_DOCA[0] = 0;
  SubTrack_DOCA[1] = 0;
  SubTrack_TotalCharge[0] = 0;
  SubTrack_TotalCharge[1] = 0;
  SubTrack_G4IDSingle[0] = 0;
  SubTrack_G4IDSingle[1] = 0;
  SubTrack_CompletenessSingle[0] = 0;
  SubTrack_CompletenessSingle[1] = 0;
  SubTrack_CleanlinessSingle[0] = 0;
  SubTrack_CleanlinessSingle[1] = 0;
  SubTrack_G4IDPrimary[0] = 0;
  SubTrack_G4IDPrimary[1] = 0;
  SubTrack_CompletenessPrimary[0] = 0;
  SubTrack_CompletenessPrimary[1] = 0;
  SubTrack_CleanlinessPrimary[0] = 0;
  SubTrack_CleanlinessPrimary[1] = 0;
  SubTrack_G4IDRecursive[0] = 0;
  SubTrack_G4IDRecursive[1] = 0;
  SubTrack_CompletenessRecursive[0] = 0;
  SubTrack_CompletenessRecursive[1] = 0;
  SubTrack_CleanlinessRecursive[0] = 0;
  SubTrack_CleanlinessRecursive[1] = 0;
  NCrossings = 0;
  CrossingIDs = NULL;
#endif
}

TECALReconVertexTrack::TECALReconVertexTrack(TECALReconVertexTrack const& other)
    : TObject(other),
      Status(other.Status),
      UniqueID(other.UniqueID),
      NHits(other.NHits),
      MatchingLikelihood(other.MatchingLikelihood),
      FrontPosition(other.FrontPosition),
      FrontDirection(other.FrontDirection),
      FrontLayerNumber(other.FrontLayerNumber),
      BackPosition(other.BackPosition),
      BackDirection(other.BackDirection),
      BackLayerNumber(other.BackLayerNumber),
      TotalCharge(other.TotalCharge),
      AverageHitTime(other.AverageHitTime),
      G4IDSingle(other.G4IDSingle),
#ifdef USE_FULL_VERTEX_TRACK_INFO
      CompletenessSingle(other.CompletenessSingle),
      CleanlinessSingle(other.CleanlinessSingle),
      Single_FrontPosition(other.Single_FrontPosition),
      Single_FrontDirection(other.Single_FrontDirection),
      Single_FrontG4PointDist(other.Single_FrontG4PointDist),
      Single_BackPosition(other.Single_BackPosition),
      Single_BackDirection(other.Single_BackDirection),
      Single_BackG4PointDist(other.Single_BackG4PointDist),
      G4IDPrimary(other.G4IDPrimary),
      CompletenessPrimary(other.CompletenessPrimary),
      CleanlinessPrimary(other.CleanlinessPrimary),
      G4IDRecursive(other.G4IDRecursive),
      CompletenessRecursive(other.CompletenessRecursive),
      CleanlinessRecursive(other.CleanlinessRecursive),
#endif
      PID_LLR_MIP_HIP(other.PID_LLR_MIP_HIP),
      PID_LLR_MIP_HIP_VACut(other.PID_LLR_MIP_HIP_VACut) {
#ifdef USE_FULL_VERTEX_TRACK_INFO
  SubTrack_UniqueID[0] = other.SubTrack_UniqueID[0];
  SubTrack_UniqueID[1] = other.SubTrack_UniqueID[1];
  SubTrack_NHits[0] = other.SubTrack_NHits[0];
  SubTrack_NHits[1] = other.SubTrack_NHits[1];
  SubTrack_Angle[0] = other.SubTrack_Angle[0];
  SubTrack_Angle[1] = other.SubTrack_Angle[1];
  SubTrack_FrontLayerNumber[0] = other.SubTrack_FrontLayerNumber[0];
  SubTrack_FrontLayerNumber[1] = other.SubTrack_FrontLayerNumber[1];
  SubTrack_BackLayerNumber[0] = other.SubTrack_BackLayerNumber[0];
  SubTrack_BackLayerNumber[1] = other.SubTrack_BackLayerNumber[1];
  SubTrack_DOCA[0] = other.SubTrack_DOCA[0];
  SubTrack_DOCA[1] = other.SubTrack_DOCA[1];
  SubTrack_TotalCharge[0] = other.SubTrack_TotalCharge[0];
  SubTrack_TotalCharge[1] = other.SubTrack_TotalCharge[1];
  SubTrack_G4IDSingle[0] = other.SubTrack_G4IDSingle[0];
  SubTrack_G4IDSingle[1] = other.SubTrack_G4IDSingle[1];
  SubTrack_CompletenessSingle[0] = other.SubTrack_CompletenessSingle[0];
  SubTrack_CompletenessSingle[1] = other.SubTrack_CompletenessSingle[1];
  SubTrack_CleanlinessSingle[0] = other.SubTrack_CleanlinessSingle[0];
  SubTrack_CleanlinessSingle[1] = other.SubTrack_CleanlinessSingle[1];
  SubTrack_G4IDPrimary[0] = other.SubTrack_G4IDPrimary[0];
  SubTrack_G4IDPrimary[1] = other.SubTrack_G4IDPrimary[1];
  SubTrack_CompletenessPrimary[0] = other.SubTrack_CompletenessPrimary[0];
  SubTrack_CompletenessPrimary[1] = other.SubTrack_CompletenessPrimary[1];
  SubTrack_CleanlinessPrimary[0] = other.SubTrack_CleanlinessPrimary[0];
  SubTrack_CleanlinessPrimary[1] = other.SubTrack_CleanlinessPrimary[1];
  SubTrack_G4IDRecursive[0] = other.SubTrack_G4IDRecursive[0];
  SubTrack_G4IDRecursive[1] = other.SubTrack_G4IDRecursive[1];
  SubTrack_CompletenessRecursive[0] = other.SubTrack_CompletenessRecursive[0];
  SubTrack_CompletenessRecursive[1] = other.SubTrack_CompletenessRecursive[1];
  SubTrack_CleanlinessRecursive[0] = other.SubTrack_CleanlinessRecursive[0];
  SubTrack_CleanlinessRecursive[1] = other.SubTrack_CleanlinessRecursive[1];
  NCrossings = other.NCrossings;
  if (NCrossings) {
    CrossingIDs = new UInt_t[NCrossings];
    for (Int_t c_it = 0; c_it < NCrossings; ++c_it) {
      CrossingIDs[c_it] = other.CrossingIDs[c_it];
    }
  } else {
    CrossingIDs = NULL;
  }
#endif
}

TECALReconVertexTrack* TECALReconVertexTrackFactory(
    ND::THandle<ND::TReconPID> trk, bool IsMC, void* pos,
    ND::THandle<ND::TG4TrajectoryContainer> trajectories, bool Allocate) {
  TECALReconVertexTrack* obj = NULL;
  if (Allocate) {
    if (pos != NULL) {
      obj = new (pos) TECALReconVertexTrack;
    } else {
      obj = new TECALReconVertexTrack;
    }
  } else {
    obj = static_cast<TECALReconVertexTrack*>(pos);
  }

  obj->UniqueID = trk->GetUniqueID();
  obj->NHits = trk->GetHits()->size();
  obj->MatchingLikelihood =
      GetTRealDatumValue(trk, "matching_likelihood", -9999);

  // Get the positions and directions of the track
  // Get the start first
  ND::THandle<ND::TPIDState> front_state = trk->GetNodes().front()->GetState();
  obj->FrontPosition = front_state->GetPosition();
  obj->FrontDirection = front_state->GetDirection();
  obj->FrontLayerNumber = ND::TGeomInfo::Get().ECAL().GetLayerNumber(
      trk->GetNodes().front()->GetObject()->GetHits()->front());

  // Now the end
  ND::THandle<ND::TPIDState> back_state = trk->GetNodes().back()->GetState();
  obj->BackPosition = back_state->GetPosition();
  obj->BackDirection = back_state->GetDirection();
  obj->BackLayerNumber = ND::TGeomInfo::Get().ECAL().GetLayerNumber(
      trk->GetNodes().back()->GetObject()->GetHits()->front());

  // Calculate the total charge
  obj->TotalCharge = 0;
  obj->AverageHitTime = 0;
  for (ND::THitSelection::iterator hitIt = trk->GetHits()->begin();
       hitIt != trk->GetHits()->end(); hitIt++) {
    obj->TotalCharge += (*hitIt)->GetCharge();
    obj->AverageHitTime += ((*hitIt)->GetTime() * (*hitIt)->GetCharge());
  }
  obj->AverageHitTime /= obj->TotalCharge;

  obj->StraightTrackLength =
      (obj->FrontPosition - obj->BackPosition).Vect().Mag();
  obj->NodeTrackLength = 0;

  TLorentzVector LastPos;
  for (ND::TReconNodeContainer::iterator n_it = trk->GetNodes().begin();
       n_it != trk->GetNodes().end(); ++n_it) {
    if (n_it != trk->GetNodes().begin()) {
      obj->NodeTrackLength +=
          (LastPos -
           ND::THandle<ND::TPIDState>((*n_it)->GetState())->GetPosition())
              .Vect()
              .Mag();
    }
    LastPos = ND::THandle<ND::TPIDState>((*n_it)->GetState())->GetPosition();
  }

// Now lets get the 2D constituents of the track and store those bits of
// information.  The current set up is that each 3D track has TWO 2D
// constituents.  To keep the file size increase to a minimum we are
// going to store all 2D information as direct data members of the
// TECALReconVertexTracks, rather than creating entirely new classes
// to store the info
// Lets get the constituents and fill

#ifdef USE_FULL_VERTEX_TRACK_INFO
  for (ND::TReconObjectContainer::iterator rocIt =
           trk->GetConstituents()->begin();
       rocIt != trk->GetConstituents()->end(); rocIt++) {
    // We are going to need an index to fill the arrays, so get one now
    int index = std::distance(trk->GetConstituents()->begin(), rocIt);

    // This is meant to be a loop over the two 1D hough tracks.
    if (index > 1) {
      break;
    }

    // All 2D constituents are stored as TReconClusters
    ND::THandle<ND::TReconCluster> input_constituent = (*rocIt);

    // This will happen if the track is a merged track as its constituents are
    // TReconTracks
    if (!input_constituent) {
      continue;
    }

    obj->SubTrack_UniqueID[index] = input_constituent->GetUniqueID();
    obj->SubTrack_NHits[index] = input_constituent->GetHits()->size();
    obj->SubTrack_Angle[index] =
        input_constituent->Get<ND::TIntegerDatum>("hough_angle")->GetValue();
    obj->SubTrack_DOCA[index] =
        input_constituent->Get<ND::TIntegerDatum>("hough_distance")->GetValue();
    obj->SubTrack_TotalCharge[index] = 0;
    obj->SubTrack_FrontLayerNumber[index] = 9999;
    obj->SubTrack_BackLayerNumber[index] = -9999;
    for (ND::THitSelection::iterator hitIt =
             input_constituent->GetHits()->begin();
         hitIt != input_constituent->GetHits()->end(); hitIt++) {
      // Calculate the total charge
      obj->SubTrack_TotalCharge[index] += (*hitIt)->GetCharge();
      // Calculate the front and back layer numbers
      int Hit_LayerNumber = ND::TGeomInfo::Get().ECAL().GetLayerNumber(*hitIt);
      if (Hit_LayerNumber < obj->SubTrack_FrontLayerNumber[index]) {
        obj->SubTrack_FrontLayerNumber[index] = Hit_LayerNumber;
      }
      if (Hit_LayerNumber > obj->SubTrack_BackLayerNumber[index]) {
        obj->SubTrack_BackLayerNumber[index] = Hit_LayerNumber;
      }
    }

    if (!IsMC) {
      continue;
    }

    ND::TShowerTruthInfo showerTruth(input_constituent);
    // Start with single trajectories
    std::vector<ND::TShowerTruthInfo::TruthInfo> results =
        showerTruth.GetSingleTrajectories();
    if (results.size() > 0) {
      obj->SubTrack_G4IDSingle[index] = results[0].ParticleID;
      obj->SubTrack_CompletenessSingle[index] = results[0].Completeness_NumHits;
      obj->SubTrack_CleanlinessSingle[index] = results[0].Cleanliness_NumHits;

    } else {
      obj->SubTrack_G4IDSingle[index] = -1;
      obj->SubTrack_CompletenessSingle[index] = -1.;
      obj->SubTrack_CleanlinessSingle[index] = -1.;
    }
    // Now the primary trajectories
    results = showerTruth.GetPrimaryTrajectories();
    if (results.size() > 0) {
      obj->SubTrack_G4IDPrimary[index] = results[0].ParticleID;
      obj->SubTrack_CompletenessPrimary[index] =
          results[0].Completeness_NumHits;
      obj->SubTrack_CleanlinessPrimary[index] = results[0].Cleanliness_NumHits;
    } else {
      obj->SubTrack_G4IDPrimary[index] = -1;
      obj->SubTrack_CompletenessPrimary[index] = -1.;
      obj->SubTrack_CleanlinessPrimary[index] = -1.;
    }
    // Now the recursive trajectories
    results = showerTruth.GetRecursiveTrajectories();
    if (results.size() > 0) {
      obj->SubTrack_G4IDRecursive[index] = results[0].ParticleID;
      obj->SubTrack_CompletenessRecursive[index] =
          results[0].Completeness_NumHits;
      obj->SubTrack_CleanlinessRecursive[index] =
          results[0].Cleanliness_NumHits;
    } else {
      obj->SubTrack_G4IDRecursive[index] = -1;
      obj->SubTrack_CompletenessRecursive[index] = -1.;
      obj->SubTrack_CleanlinessRecursive[index] = -1.;
    }
  }  // End loop ovr 2D constituents

  // For what this is used for, we don't care
  obj->NCrossings = 0;
  obj->CrossingIDs = 0;
#endif

  obj->PID_LLR_MIP_EM = GetTRealDatumValue(trk, "LLR_MIP_EM", -9999);
  obj->PID_LLR_EM_HIP = GetTRealDatumValue(trk, "LLR_EM_HIP", -9999);

  obj->PID_LLR_MIP_HIP = GetTRealDatumValue(trk, "LLR_MIP_HIP", -9999);
  obj->PID_LLR_MIP_HIP_VACut =
      GetTRealDatumValue(trk, "LLR_MIP_HIP_VACut", -9999);

  // Containment
  // This isn't a PID variable, but it's easy to save it the same way.
  obj->Containment = GetTRealDatumValue(trk, "Containment", -1);

  if (!IsMC) {
    return obj;
  }

  // Get the truth information for the track
  ND::TShowerTruthInfo showerTruth(trk);
  // Start with single trajectories
  std::vector<ND::TShowerTruthInfo::TruthInfo> results =
      showerTruth.GetSingleTrajectories();
  if (results.size() > 0) {
    obj->G4IDSingle = results[0].ParticleID;

#ifdef USE_FULL_VERTEX_TRACK_INFO
    obj->CompletenessSingle = results[0].Completeness_NumHits;
    obj->CleanlinessSingle = results[0].Cleanliness_NumHits;

    if (trajectories) {
      ND::THandle<ND::TG4Trajectory> traj =
          trajectories->GetTrajectory(results[0].ParticleID);
      if (traj) {
        TVector3 Scratch = obj->FrontPosition.Vect();
        ND::TG4TrajectoryPoint const& g4pointF =
            TrackTruthInfo::GetCloserG4Point(traj, &Scratch,
                                             obj->Single_FrontG4PointDist);
        obj->Single_FrontPosition = g4pointF.GetPosition().Vect();
        obj->Single_FrontDirection = g4pointF.GetMomentum().Unit();

        Scratch = obj->BackPosition.Vect();
        ND::TG4TrajectoryPoint const& g4pointB =
            TrackTruthInfo::GetCloserG4Point(traj, &Scratch,
                                             obj->Single_BackG4PointDist);
        obj->Single_BackPosition = g4pointB.GetPosition().Vect();
        obj->Single_BackDirection = g4pointB.GetMomentum().Unit();
      }
    }
#endif
  } else {
    obj->G4IDSingle = -1;
#ifdef USE_FULL_VERTEX_TRACK_INFO
    obj->CompletenessSingle = -1.;
    obj->CleanlinessSingle = -1.;
#endif
  }
#ifdef USE_FULL_VERTEX_TRACK_INFO
  // Now the primary trajectories
  results = showerTruth.GetPrimaryTrajectories();
  if (results.size() > 0) {
    obj->G4IDPrimary = results[0].ParticleID;
    obj->CompletenessPrimary = results[0].Completeness_NumHits;
    obj->CleanlinessPrimary = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDPrimary = -1;
    obj->CompletenessPrimary = -1.;
    obj->CleanlinessPrimary = -1.;
  }
  // Finally the recursive trajectories
  results = showerTruth.GetRecursiveTrajectories();
  if (results.size() > 0) {
    obj->G4IDRecursive = results[0].ParticleID;
    obj->CompletenessRecursive = results[0].Completeness_NumHits;
    obj->CleanlinessRecursive = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDRecursive = -1;
    obj->CompletenessRecursive = -1.;
    obj->CleanlinessRecursive = -1.;
  }
#endif
  return obj;
}

TECALReconVertexCrossing::TECALReconVertexCrossing()
    : UniqueID(0),
      NHits(0),
      HitCharges(NULL),
      HitPositionsX(NULL),
      HitPositionsY(NULL),
      HitPositionsZ(NULL),
      HitTimes(NULL),
      HitLayerNumbers(NULL),
      HitBarNumbers(NULL),
      Position(0, 0, 0, 0),
      DOCA(0),
      NTracks(0),
      TrackIDs(NULL) {}

#ifdef TTRACKERECALRECONMODULE_USE_EXP

void OutputManager::SaveUsedHit(ND::THandle<ND::THit> hit, unsigned int UID,
                                bool overwritePrevious) {
  if (!fDoSaveHitInfo) {
    return;
  }
  if (fNUsedHits == MAXHITSTOSAVE) {
    ND280NamedWarn("TTrackerECALReconModule", "Recieved too many ecalHits...");
    return;
  }
  unsigned long hitHash = TECALHitHasher(hit);
  ND280NamedInfo("TTrackerECALReconModule", "HH: " << hitHash);

  TECALReconHit* trHit = 0;
  if (fSavedHitMap.count(hitHash)) {  // already have the hit, then get it.
    trHit = static_cast<TECALReconHit*>(fUsedHits->At(fSavedHitMap[hitHash]));
  }

  // If we do not want to overwrite previous versions of this hit, then bail if
  // we see it again
  if (trHit && (!overwritePrevious)) {
    ND280NamedDebug("TTrackerECALReconModule",
                    "Already have Hit info (Hash:"
                        << hitHash << " ), Adding this association: " << UID);
    trHit->AddAssocUID(UID);
    ND280NamedInfo("TTrackerECALReconModule", "\t" << trHit->Position << ", "
                                                   << trHit->ChargeDeposit);
    return;
  }

  if (!trHit) {  // We don't already have it, make a new one
    trHit = new ((*fUsedHits)[fNUsedHits++]) TECALReconHit();
    fSavedHitMap.insert(std::make_pair(hitHash, fNUsedHits - 1));
    trHit->GeomId = hit->GetGeomId().AsInt();
    trHit->IsXYZ =
        ((hit->IsXHit()) | (hit->IsYHit() << 1) | (hit->IsZHit() << 2));
  } else {
    ND280NamedInfo("TTrackerECALReconModule", "overwriting hit: " << (*trHit));
  }

  // Overwrite old info or write to a new hit
  trHit->AddAssocUID(UID);
  trHit->Position = TLorentzVector(hit->GetPosition(), hit->GetTime());
  trHit->PosUncertainty = hit->GetUncertainty();
  trHit->ChargeDeposit = hit->GetCharge();
  trHit->IsXYZ =
      ((hit->IsXHit()) | (hit->IsYHit() << 1) | (hit->IsZHit() << 2));

  ND280NamedInfo("TTrackerECALReconModule",
                 "Adding hit info Hit(#" << (fNUsedHits - 1) << "), "
                                         << (*trHit) << " to UID: " << UID);
  ND280NamedDebug("TTrackerECALReconModule", "\tHitHash (" << hitHash << ")");
}

TECALReconHit::TECALReconHit()
    : GeomId(-1),
      UsedByNObjs(0),
      ObjUIDs(NULL),
      Position(0, 0, 0, 0),
      PosUncertainty(0, 0, 0),
      ChargeDeposit(0),
      IsXYZ(0) {
  ObjUIDs = new UInt_t[TECALReconHit::MaxAssocs];
}

TECALReconHit::~TECALReconHit() {
  if (ObjUIDs) {
    delete ObjUIDs;
  }
}

void TECALReconHit::AddAssocUID(unsigned int tuid) {
  if (UsedByNObjs == MaxAssocs) {
    ND280NamedWarn("TTrackerECALReconModule",
                   "Attempting to add more object"
                   " associations than we have space for: "
                       << MaxAssocs << ".");
    return;
  }
  if (UsedByNObjs == 2) {
    ND280NamedWarn(
        "TTrackerECALReconModule",
        "This hit is already in use"
        " by two objects, likely one VertexCluster and one VertexTrack.");
    ND280NamedWarn("TTrackerECALReconModule",
                   "+Attempting to add associat"
                   "ion: "
                       << tuid);
    for (int it = 0; it < UsedByNObjs; ++it) {
      ND280NamedWarn("TTrackerECALReconModule", "Used By: " << ObjUIDs[it]);
    }
  }
  ObjUIDs[UsedByNObjs++] = tuid;
}

#endif

TECALReconVertexCandidate::TECALReconVertexCandidate() {
  VertexClusterUID = 0;
  Position = TLorentzVector(0, 0, 0, 0);
  ECalModule = "";
  NTracks = 0;
  NHits = 0;
  G4IDSingle = -1;
  CompletenessSingle = -1;
  CleanlinessSingle = -1;
  G4IDPrimary = -1;
  CompletenessPrimary = -1;
  CleanlinessPrimary = -1;
  G4IDRecursive = -1;
  CompletenessRecursive = -1;
  CleanlinessRecursive = -1;
}

TECALReconVertexCandidate* TECALVertexCandidateFactory(
    const TECALReconVertex& rv, ND::THandle<ND::TReconVertex> top_vertex,
    bool IsMC, void* pos) {
  UInt_t VCUID = top_vertex->GetUniqueID();
  ND280NamedInfo("TTrackerECALReconModule",
                 "Adding vertex candidate info:"
                 " Association: "
                     << VCUID << ", Pos: " << rv.GetPosition() << ", uses "
                     << rv.NTracks() << " tracks.");

  TECALReconVertexCandidate* obj = NULL;
  if (pos != NULL) {
    obj = new (pos) TECALReconVertexCandidate;
  } else {
    ND280NamedWarn("TTrackerECALReconModule",
                   "TECALVertexCandidateFactory making new "
                   "TECALReconVertexCandidate instance with non-positional new "
                   "operator.");
    obj = new TECALReconVertexCandidate;
  }

  obj->NHits = rv.VertCandyHits->size();
  obj->ECalModule = ND::TGeomInfo::Get()
                        .ECAL()
                        .GetModule(rv.VertCandyHits->front())
                        .GetName();

  obj->NTracks = rv.NTracks();
  obj->TrackUIDs = new UInt_t[obj->NTracks];

  int tctr = 0;
  for (std::vector<ND::THandle<ND::TReconPID> >::const_iterator
           trk_it = rv.Tracks_cbegin();
       trk_it != rv.Tracks_cend(); ++trk_it, ++tctr) {
    obj->TrackUIDs[tctr] = (*trk_it)->GetUniqueID();
  }

  obj->Position = rv.GetPosition();
  obj->VertexClusterUID = VCUID;

  // Here be truth matching routines.
  if (!IsMC) {
    return obj;
  }

  // Now we have semi-sensibly separate vertex candidate hits
  // Get the truth information for the cluster
  ND::TShowerTruthInfo shower_truth(rv.VertCandyHits,
                                    ND::TShowerTruthInfo::kXYZ, 1, 1, 1);
  // Start with single trajectories
  std::vector<ND::TShowerTruthInfo::TruthInfo> results =
      shower_truth.GetSingleTrajectories();
  if (results.size()) {
    obj->G4IDSingle = results[0].ParticleID;
    obj->CompletenessSingle = results[0].Completeness_NumHits;
    obj->CleanlinessSingle = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDSingle = -1;
    obj->CompletenessSingle = -1;
    obj->CleanlinessSingle = -1;
  }
  // Now the primary trajectories
  results = shower_truth.GetPrimaryTrajectories();
  if (results.size()) {
    obj->G4IDPrimary = results[0].ParticleID;
    obj->CompletenessPrimary = results[0].Completeness_NumHits;
    obj->CleanlinessPrimary = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDPrimary = -1;
    obj->CompletenessPrimary = -1;
    obj->CleanlinessPrimary = -1;
  }
  // Now the primary trajectories
  results = shower_truth.GetRecursiveTrajectories();
  if (results.size()) {
    obj->G4IDRecursive = results[0].ParticleID;
    obj->CompletenessRecursive = results[0].Completeness_NumHits;
    obj->CleanlinessRecursive = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDRecursive = -1;
    obj->CompletenessRecursive = -1;
    obj->CleanlinessRecursive = -1;
  }
  return obj;
}

TECALReconVertexCandidate* TECALVertexCandidateFactory(
    ND::THandle<ND::TReconVertex> ECalVtx, UInt_t VCUID, bool IsMC, void* pos) {
  ND280NamedInfo("TTrackerECALReconModule",
                 "Adding vertex candidate info:"
                 " Association: "
                     << VCUID << ", Pos: " << ECalVtx->GetPosition());

  TECALReconVertexCandidate* obj = NULL;
  if (pos != NULL) {
    obj = new (pos) TECALReconVertexCandidate;
  } else {
    ND280NamedWarn("TTrackerECALReconModule",
                   "TECALVertexCandidateFactory making new "
                   "TECALReconVertexCandidate instance with non-positional new "
                   "operator.");
    obj = new TECALReconVertexCandidate;
  }

  obj->NHits = ECalVtx->GetHits()->size();
  obj->ECalModule = ND::TGeomInfo::Get()
                        .ECAL()
                        .GetModule(ECalVtx->GetHits()->front())
                        .GetName();

  obj->NTracks = 0;
  for (UInt_t const_it = 0; const_it < ECalVtx->GetConstituents()->size();
       ++const_it) {
    ND::THandle<ND::TReconPID> trk = ECalVtx->GetConstituents()->at(const_it);
    if (trk) {
      obj->NTracks++;
    }
  }
  obj->TrackUIDs = new UInt_t[obj->NTracks];
  size_t tit = 0;
  for (UInt_t const_it = 0; const_it < ECalVtx->GetConstituents()->size();
       ++const_it) {
    ND::THandle<ND::TReconPID> trk = ECalVtx->GetConstituents()->at(const_it);
    if (trk) {
      obj->TrackUIDs[tit++] = trk->GetUniqueID();
    }
  }

  obj->Position = ECalVtx->GetPosition();
  obj->VertexClusterUID = VCUID;

  // Here be truth matching routines.
  if (!IsMC) {
    return obj;
  }

  // Get the truth information for the cluster
  ND::TShowerTruthInfo shower_truth(ECalVtx->GetHits(),
                                    ND::TShowerTruthInfo::kXYZ, 1, 1, 1);
  // Start with single trajectories
  std::vector<ND::TShowerTruthInfo::TruthInfo> results =
      shower_truth.GetSingleTrajectories();
  if (results.size()) {
    obj->G4IDSingle = results[0].ParticleID;
    obj->CompletenessSingle = results[0].Completeness_NumHits;
    obj->CleanlinessSingle = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDSingle = -1;
    obj->CompletenessSingle = -1;
    obj->CleanlinessSingle = -1;
  }
  // Now the primary trajectories
  results = shower_truth.GetPrimaryTrajectories();
  if (results.size()) {
    obj->G4IDPrimary = results[0].ParticleID;
    obj->CompletenessPrimary = results[0].Completeness_NumHits;
    obj->CleanlinessPrimary = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDPrimary = -1;
    obj->CompletenessPrimary = -1;
    obj->CleanlinessPrimary = -1;
  }
  // Now the primary trajectories
  results = shower_truth.GetRecursiveTrajectories();
  if (results.size()) {
    obj->G4IDRecursive = results[0].ParticleID;
    obj->CompletenessRecursive = results[0].Completeness_NumHits;
    obj->CleanlinessRecursive = results[0].Cleanliness_NumHits;
  } else {
    obj->G4IDRecursive = -1;
    obj->CompletenessRecursive = -1;
    obj->CleanlinessRecursive = -1;
  }
  return obj;
}

TECALReconTrack::~TECALReconTrack() {}
TECALReconShower::~TECALReconShower() {}
TECALReconCluster::~TECALReconCluster() {}
TECALReconObject::~TECALReconObject() {}
TECALReconUnmatchedObject::~TECALReconUnmatchedObject() {}
}
}