TP0DReconModule.cxx

Go to the documentation of this file.

#include <functional>
#include <algorithm>
#include <cassert>
#include <string>

#include "TString.h"
#include "TVector3.h"

#include "TAlgorithmResult.hxx"
#include "TDatum.hxx"
#include "TG4PrimaryVertex.hxx"
#include "TGeomInfo.hxx"
#include "TGeometryId.hxx"
#include "THandle.hxx"
#include "THitSelection.hxx"
#include "TND280Event.hxx"
#include "TND280Log.hxx"
#include "TPIDState.hxx"
#include "TP0DGeom.hxx"
#include "TRealDatum.hxx"
#include "TIntegerDatum.hxx"
#include "TReconCluster.hxx"
#include "TReconNode.hxx"
#include "TReconPID.hxx"
#include "TReconShower.hxx"
#include "TReconTrack.hxx"
#include "TReconVertex.hxx"
#include "TReconState.hxx"
#include "TShowerState.hxx"
#include "TTrackState.hxx"
#include "TVertexState.hxx"
#include "TG4HitSegment.hxx"
#include "TTFBChannelId.hxx"
#include "retrieveHitTruthInfo.hxx"

#include "TP0DReconModule.hxx"
#include "oaEvent_version.h"

#if (oaEvent_MAJOR_VERSION>8)
#include "TReconStateUtilities.hxx"
#include "TShowerState.hxx"
#include "TPIDState.hxx"
using ND::TReconStateUtilities::GetPosition;
using ND::TReconStateUtilities::GetPositionVariance;
using ND::TReconStateUtilities::GetPositionIndex;
#endif
using namespace ND;

ClassImp(ND::TP0DReconModule::TP0DAlgoRes);
ND::TP0DReconModule::TP0DAlgoRes::~TP0DAlgoRes(){}

ClassImp(ND::TP0DReconModule::TP0DVertex);
ND::TP0DReconModule::TP0DVertex::~TP0DVertex(){}

ClassImp(ND::TP0DReconModule::TP0DParticle);
ND::TP0DReconModule::TP0DParticle::~TP0DParticle() {}

ClassImp(ND::TP0DReconModule::TP0DTrack);
ND::TP0DReconModule::TP0DTrack::~TP0DTrack(){}

ClassImp(ND::TP0DReconModule::TP0DShower);
ND::TP0DReconModule::TP0DShower::~TP0DShower(){}

ClassImp(ND::TP0DReconModule::TP0DNode);
ND::TP0DReconModule::TP0DNode::~TP0DNode(){}

ClassImp(ND::TP0DReconModule::TP0DCluster);
ND::TP0DReconModule::TP0DCluster::~TP0DCluster() {}

ClassImp(ND::TP0DReconModule::TP0DHit);

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

ND::TP0DReconModule::TP0DReconModule(const char *name,
                                     const char *title) {
  SetNameTitle(name, title);

  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "P0D Recon Output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;

  // Tree variables
  fNAlgoResults = 0;
  fNVertices = 0;
  fNParticles = 0;
  fNShowers = 0;
  //  Commented out this line to remove extraneous eventAnalysis output
  //  std::cout << "fNShowers:\tSta\t" << fNShowers << std::endl;
  fNTracks = 0;
  fNNodes = 0;
  fNClusters = 0;
  fNHits = 0;

  fAlgoResults = new TClonesArray("ND::TP0DReconModule::TP0DAlgoRes", 200);
  fVertices    = new TClonesArray("ND::TP0DReconModule::TP0DVertex", 200);
  fParticles   = new TClonesArray("ND::TP0DReconModule::TP0DParticle", 200);
  fShowers     = new TClonesArray("ND::TP0DReconModule::TP0DShower", 200);
  fTracks      = new TClonesArray("ND::TP0DReconModule::TP0DTrack", 200);
  fNodes       = new TClonesArray("ND::TP0DReconModule::TP0DNode", 200);
  fHits        = new TClonesArray("ND::TP0DReconModule::TP0DHit", 200);
  fClusters    = new TClonesArray("ND::TP0DReconModule::TP0DCluster", 200);

  fRejectAlgoResultList.clear();
  fRejectAlgoResultList.push_back(TPRegexp("TP0DTrack$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DExpandSeeds$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0D3DTrackMatching$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DPairwiseVertex$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0D3DShowerMatching$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DShareCharge$"));
  fRejectAlgoResultList.push_back(TPRegexp("-Expand$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DPairwiseVertexSingle$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DTrackVertex$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DTrackPID$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DShowerInputs$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DShowerVertex$"));
  fRejectAlgoResultList.push_back(TPRegexp("TP0DShowerPID$"));
}


ND::TP0DReconModule::~TP0DReconModule() {}

Bool_t ND::TP0DReconModule::ProcessFirstEvent(ND::TND280Event&) {
  return true;
}


void ND::TP0DReconModule::InitializeModule() {}

void ND::TP0DReconModule::InitializeBranches() {

  // The number of objects that have been saved.
  fOutputTree->Branch("NAlgoResults", &fNAlgoResults,
      "NAlgoResults/S", fBufferSize)
  ->SetTitle("The number of Algorithm Results ");
  fOutputTree->Branch("NVertices",    &fNVertices,
      "NVertices/S",    fBufferSize)
  ->SetTitle(" The number of added vertices");
  fOutputTree->Branch("NParticles",   &fNParticles,
      "NParticles/S",   fBufferSize)
  ->SetTitle(" The number of particles ");
  fOutputTree->Branch("NShowers",     &fNShowers,
      "NShowers/S",     fBufferSize)
  ->SetTitle(" The number of showers");
  fOutputTree->Branch("NTracks",      &fNTracks,
      "NTracks/S",      fBufferSize)
  ->SetTitle("  The number of tracks");
  fOutputTree->Branch("NNodes",       &fNNodes,
      "NNodes/S",       fBufferSize)
  ->SetTitle(" The number of track nodes");
  fOutputTree->Branch("NHits",        &fNHits,
      "NHits/S",        fBufferSize)
  ->SetTitle("The number of hits that are saved ");
  fOutputTree->Branch("NClusters",    &fNClusters,
      "NClusters/S",    fBufferSize)
  ->SetTitle("The number of clusters ");

  // A branch for the clone array of P0D object information.
  fOutputTree->Branch("AlgoResults", &fAlgoResults,
      fBufferSize, fSplitLevel)
  ->SetTitle("The TP0DAlgoRes vector of Algorithm Results. ");
  fOutputTree->Branch("Vertices",    &fVertices,
      fBufferSize, fSplitLevel)
  ->SetTitle("  The TP0DVertex vector of vertices.");
  fOutputTree->Branch("Particles",   &fParticles,
      fBufferSize, fSplitLevel)
  ->SetTitle("The TP0DParticle vector of particles. ");
  fOutputTree->Branch("Showers",     &fShowers,
      fBufferSize, fSplitLevel)
  ->SetTitle("The TP0DShower vector of showers. ");
  fOutputTree->Branch("Tracks",      &fTracks,
      fBufferSize, fSplitLevel)
  ->SetTitle(" The TP0DTrack vector of tracks.");
  fOutputTree->Branch("Nodes",       &fNodes,
      fBufferSize, fSplitLevel)
  ->SetTitle("  The TP0DNode vector of nodes.");
  fOutputTree->Branch("Hits",        &fHits,
      fBufferSize, fSplitLevel)
  ->SetTitle("The TP0DHit vector of hits. ");
  fOutputTree->Branch("Clusters",    &fClusters,
      fBufferSize, fSplitLevel)
  ->SetTitle("The TP0DCluster vector of clusters ");
}


bool ND::TP0DReconModule::FillTree(ND::TND280Event& event) {
  // Clear variables
  fNAlgoResults = 0;
  fNVertices    = 0;
  fNParticles   = 0;
  fNShowers     = 0;
  fNTracks      = 0;
  fNNodes       = 0;
  fNClusters    = 0;
  fNHits        = 0;

  fAlgoResults->Clear();
  fVertices->Clear();
  fParticles->Clear();
  fShowers->Clear();
  fTracks->Clear();
  fNodes->Clear();
  fClusters->Clear();
  fHits->Clear();

  fTempHitMap.clear();

  //###############################################################

  // Get the main result to save.
  ND::THandle<ND::TAlgorithmResult> p0dRecon = event.GetFit("p0dRecon");

  if (!p0dRecon) {
    ND280Info("No p0dRecon in Event");
    return true;
  }

  FillAlgorithmResult(p0dRecon, -1);

  return true;
}

template <class T>
void TP0DReconModule::FillBaseObject(T basePtr,
                                     ND::THandle<ND::TReconBase> baseObject,
                                     bool saveHits){
  // Algorithm Name
  basePtr->AlgorithmName = baseObject->GetAlgorithmName();

  // Hit Based Information
  ND::THandle<ND::THitSelection> baseHits = baseObject->GetHits();

  if(baseHits){
    basePtr->Cycle = TP0DReconModule::GetCycle(baseHits);
    basePtr->NHits = (short)baseHits->size();

    if(saveHits){
      ND::THitSelection::const_iterator hit;
      for(hit = baseHits->begin();
          hit != baseHits->end(); ++hit){
        short hitId = FillHit(*hit);
        if(hitId >= 0)
          basePtr->Hits.push_back(hitId);
      }
    }
  } else {
    basePtr->Cycle = -1;
    basePtr->NHits = -1;
  }
  // Unique ID
  basePtr->UniqueID      = baseObject->GetUniqueID();

  // Daughter objects
  ND::THandle<ND::TReconObjectContainer> objCont =
    baseObject->GetConstituents();
  if(objCont) {
    ND::TReconObjectContainer::const_iterator constituent;
    for (constituent = objCont->begin();
         constituent != objCont->end();
         constituent++){

      ObjectID conID = FillReconObject(*constituent, saveHits);
      if(conID.type == ObjectID::kVertex){
        basePtr->Vertices.push_back(conID.id);
      }
      else if(conID.type == ObjectID::kParticle){
        basePtr->Particles.push_back(conID.id);
      }
      else if(conID.type == ObjectID::kTrack){
        basePtr->Tracks.push_back(conID.id);
      }
      else if(conID.type == ObjectID::kShower){
        basePtr->Showers.push_back(conID.id);
      }
      else if(conID.type == ObjectID::kCluster){
        basePtr->Clusters.push_back(conID.id);
      }
    }
  }

  // Nodes (probably only used by tracks, showers and particles, but
  // technically valid for all TReconBase objects.
  const ND::TReconNodeContainer& nodeCont = baseObject->GetNodes();
  ND::TReconNodeContainer::const_iterator node;
  for (node = nodeCont.begin();
       node != nodeCont.end();
       node++){

    int nodeID = FillNode(*node, saveHits);
    basePtr->Nodes.push_back(nodeID);
  }

  // Muon Decay Clusters
  for (TDataVector::iterator d = baseObject->begin();
       d != baseObject->end(); ++d) {

    if ( !(*d) ){continue;}
    std::string sub_class_name = (*d)->ClassName();
    std::string fullname = (*d)->GetFullName().Data();

    if ( (sub_class_name.find("ND::TReconObjectContainer") !=
                                             std::string::npos) &&
        (fullname.find("//unnamed/decays") != std::string::npos) ) {

      ND::THandle<ND::TReconObjectContainer> clusters
        = (*d)->Get<ND::TReconObjectContainer>(".");
      if(!clusters){continue;}

      ND::TReconObjectContainer::const_iterator clus;
      for (clus = clusters->begin();
           clus != clusters->end();
           clus++){

        ObjectID conID = FillReconObject(*clus, saveHits);
        if(conID.type == ObjectID::kCluster){
          basePtr->Clusters.push_back(conID.id);
        }
        else{
          ND280Error("Muon Decay object not a cluster");
        }
      }
    }
  }
};

short
ND::TP0DReconModule::FillAlgorithmResult(
    const ND::THandle<ND::TAlgorithmResult> algoRes,
    short parent){

  TString fullName = algoRes->GetFullName();
  TString baseName = fullName;
  baseName.Remove(0,baseName.Last('/') + 1);

  // To save space, we're not going to save a number of the lower
  // algorithm results.  It is important to remember, that as soon as
  // a result is not saved, none of it's daughter results will be
  // saved either.
  for(std::vector<TPRegexp>::iterator cut = fRejectAlgoResultList.begin();
      cut != fRejectAlgoResultList.end(); ++cut){

    if(cut->MatchB(fullName, "", 0, 1)){
       ND280Verbose("Rejecting:\t" << baseName
                    << " (" << cut->GetPattern() << ")");
       return -1;
    }
  }

  ND280Verbose("Fill Algorithm Result:\t" << baseName);

  int algoResID = fNAlgoResults;

  TP0DAlgoRes* p0dAlgoResult =
    new ((*fAlgoResults)[fNAlgoResults++]) TP0DAlgoRes;

  p0dAlgoResult->FullName      = fullName;
  p0dAlgoResult->AlgorithmName = baseName;
  p0dAlgoResult->Parent        = parent;

  // Get the used hits from the algo result
  ND::THandle<ND::THitSelection> usedHits = algoRes->GetHitSelection("used");

  // Save the used hits to a cluster
  if(usedHits && (!usedHits->empty())){
    ND::THitSelection::const_iterator hit;

    for(hit = usedHits->begin();
        hit != usedHits->end(); ++hit){

      short hitId = FillHit(*hit);
      if(hitId >= 0){
        p0dAlgoResult->Hits.push_back(hitId);
      }
    }
    ND::THandle<ND::TReconCluster> usedCluster(new ND::TReconCluster);
    usedCluster->FillFromHits((baseName+"Used").Data(), *usedHits);

    ObjectID clusterID = FillReconObject(usedCluster, true);
    p0dAlgoResult->UsedHitCluster = clusterID.id;
  }
  else{ p0dAlgoResult->UsedHitCluster     = -1;}

  // Get the unused hits from the algo result
  ND::THandle<ND::THitSelection> unusedHits
      = algoRes->GetHitSelection("unused");

  // Save the unused hits to a cluster
  if(unusedHits && (!unusedHits->empty())){
    ND::THandle<ND::TReconCluster> unusedCluster(new ND::TReconCluster);
    unusedCluster->FillFromHits((baseName+"Unused").Data(), *unusedHits);

    // We don't want to save the unused hits. The 'false' prevents the
    ObjectID clusterID = FillReconObject(unusedCluster, false);
    p0dAlgoResult->UnusedHitCluster = clusterID.id;
  }
  else{ p0dAlgoResult->UnusedHitCluster     = -1;}

  // Get the final recon object container
  ND::THandle<ND::TReconObjectContainer> final
    = algoRes->GetResultsContainer("final");
  if (final) {
    ND::TReconObjectContainer::const_iterator reconObject;
    for (reconObject = final->begin();
         reconObject != final->end(); reconObject++){

      ObjectID objID = FillReconObject(*reconObject, true);
      if(objID.type == ObjectID::kVertex){
        p0dAlgoResult->Vertices.push_back(objID.id);
      }
      else if(objID.type == ObjectID::kParticle){
        p0dAlgoResult->Particles.push_back(objID.id);
      }
      else if(objID.type == ObjectID::kTrack){
        p0dAlgoResult->Tracks.push_back(objID.id);
      }
      else if(objID.type == ObjectID::kShower){
        p0dAlgoResult->Showers.push_back(objID.id);
      }
      else if(objID.type == ObjectID::kCluster){
        p0dAlgoResult->Clusters.push_back(objID.id);
      }
    }
  }

  // Get the clusters of non-reconstructed hits.  This only happens in the
  // cycle reconstruction.
  ND::THandle<ND::TReconObjectContainer> clusters
    = algoRes->GetResultsContainer("clusters");
  if (clusters) {
    ND::TReconObjectContainer::const_iterator reconObject;
    for (reconObject = clusters->begin();
         reconObject != clusters->end(); reconObject++){

      ObjectID objID = FillReconObject(*reconObject, true);
      if (objID.type == ObjectID::kVertex) {
        p0dAlgoResult->Vertices.push_back(objID.id);
      }
      else if (objID.type == ObjectID::kParticle) {
        p0dAlgoResult->Particles.push_back(objID.id);
      }
      else if (objID.type == ObjectID::kTrack) {
        p0dAlgoResult->Tracks.push_back(objID.id);
      }
      else if (objID.type == ObjectID::kShower) {
        p0dAlgoResult->Showers.push_back(objID.id);
      }
      else if (objID.type == ObjectID::kCluster) {
        p0dAlgoResult->Clusters.push_back(objID.id);
      }
    }
  }

  for (TDataVector::iterator d = algoRes->begin();
       d != algoRes->end(); ++d) {
    if (!(*d)){continue;}
    std::string class_name = (*d)->ClassName();

    if (class_name.find("ND::TAlgorithmResult") != std::string::npos) {
      ND::THandle<ND::TAlgorithmResult> subResult =
        (*d)->Get<ND::TAlgorithmResult>(".");
      if(!subResult){continue;}
      int daughterID = FillAlgorithmResult(subResult, algoResID);
      if(daughterID != -1){
        p0dAlgoResult->AlgoResults.push_back(daughterID);
      }
    }
  }
  return algoResID;
}

ND::TP0DReconModule::ObjectID
ND::TP0DReconModule::FillReconObject(
  const ND::THandle<ND::TReconBase> reconObject,
  bool saveHits = true){

  std::string reco_class_name = reconObject->ClassName();

  if(reco_class_name.find("ND::TReconVertex") != std::string::npos){
    ND280Verbose("Fill Vertex");
    return FillVertexObject(reconObject, saveHits);
  }
  else if(reco_class_name.find("ND::TReconPID") != std::string::npos){
    ND280Verbose("Fill PID");
    return FillParticleObject(reconObject, saveHits);
  }
  else if(reco_class_name.find("ND::TReconTrack") != std::string::npos){
    ND280Verbose("Fill Track");
    return FillTrackObject(reconObject, saveHits);
  }
  else if(reco_class_name.find("ND::TReconShower") != std::string::npos){
    ND280Verbose("Fill Shower");
    return FillShowerObject(reconObject, saveHits);
  }
  else if(reco_class_name.find("ND::TReconCluster") != std::string::npos){
    ND280Verbose("Fill Cluster");
    return FillClusterObject(reconObject, saveHits);
  }

  ND280Error("Unknown object passed to TP0DReconModule");
  ObjectID blank;
  return blank;
}

ND::TP0DReconModule::ObjectID
ND::TP0DReconModule::FillVertexObject(
  const ND::THandle<ND::TReconVertex> vertex,
  bool saveHits = true){

  ObjectID vertexID;
  vertexID.id   = fNVertices;
  vertexID.type = ObjectID::kVertex;

  TP0DVertex* p0dVertex =
    new ((*fVertices)[fNVertices++]) TP0DVertex;

  FillBaseObject(p0dVertex, vertex, saveHits);

  ND::THandle<ND::TVertexState> vertexState = vertex->GetState();

  p0dVertex->Status        = vertex->GetStatus();
  p0dVertex->Quality       = vertex->GetQuality();
  p0dVertex->NDOF          = vertex->GetNDOF();
  p0dVertex->Position      = vertexState->GetPosition();
  p0dVertex->PosVariance   = vertexState->GetPositionVariance();
  p0dVertex->Fiducial =
    ND::TGeomInfo::P0D().WaterFiducialDistance(p0dVertex->Position.X(),
                                               p0dVertex->Position.Y(),
                                               p0dVertex->Position.Z());

  p0dVertex->ValidDimensions =
    CountValidDimensions(p0dVertex->PosVariance.Vect());

  p0dVertex->Truth_TrajIDs.clear();
  p0dVertex->Truth_HitCount.clear();
  p0dVertex->Truth_ChargeShare.clear();

  ND::THandle<ND::THitSelection> vertexHits = vertex->GetHits();
  if (vertexHits){

    std::map< int, std::pair<int, float> > hitInfo = HitTruthInfo(vertexHits);
    std::map< int, std::pair< int, float > >::iterator parentID;

    for(parentID = hitInfo.begin(); parentID != hitInfo.end(); ++parentID){
      p0dVertex->Truth_TrajIDs.push_back((*parentID).first);
      p0dVertex->Truth_HitCount.push_back((*parentID).second.first);
      p0dVertex->Truth_ChargeShare.push_back((*parentID).second.second);
    }

    p0dVertex->Truth_PrimaryTrajIDs = HitTruthPrimaryInfo(vertexHits);
  }

  return vertexID;
}

ND::TP0DReconModule::ObjectID
ND::TP0DReconModule::FillParticleObject(
  const ND::THandle<ND::TReconPID> particle,
  bool saveHits = true){

  ObjectID particleID;
  particleID.id   = fNParticles;
  particleID.type = ObjectID::kParticle;

  TP0DParticle* p0dParticle =
    new ((*fParticles)[fNParticles++]) TP0DParticle;

  FillBaseObject(p0dParticle, particle, saveHits);

  ND::THandle<ND::TPIDState> particleState = particle->GetState();

  p0dParticle->Status        = particle->GetStatus();
  p0dParticle->Quality       = particle->GetQuality();
  p0dParticle->NDOF          = particle->GetNDOF();
  p0dParticle->Position      = particleState->GetPosition();
  p0dParticle->PosVariance   = particleState->GetPositionVariance();
  p0dParticle->Direction     = particleState->GetDirection();
  p0dParticle->DirVariance   = particleState->GetDirectionVariance();
  p0dParticle->Momentum      = particle->GetMomentum();
  p0dParticle->Charge        = particle->GetCharge();

  p0dParticle->ValidDimensions =
    CountValidDimensions(p0dParticle->PosVariance.Vect());

  const ND::THandle<ND::THitSelection> particleHits = particle->GetHits();
  if (particleHits){

    std::map< int, std::pair<int, float> > hitInfo = HitTruthInfo(particleHits);
    std::map< int, std::pair< int, float > >::iterator parentID;

    for(parentID = hitInfo.begin(); parentID != hitInfo.end(); ++parentID){
      p0dParticle->Truth_TrajIDs.push_back((*parentID).first);
      p0dParticle->Truth_HitCount.push_back((*parentID).second.first);
      p0dParticle->Truth_ChargeShare.push_back((*parentID).second.second);
    }

    p0dParticle->Truth_PrimaryTrajIDs = HitTruthPrimaryInfo(particleHits);
  }

  p0dParticle->SideDeposit = 0.0;
  p0dParticle->EndDeposit  = 0.0;

  ND::THitSelection::const_iterator hit;
  for (hit = particleHits->begin();
       hit != particleHits->end();
       ++hit) {

    ND::TGeometryId geomid = (*hit)->GetGeomId();
    if (!ND::GeomId::P0D::IsP0D(geomid)) continue;
    int p0dule = ND::TGeomInfo::P0D().GetBar(geomid).GetP0Dule();
    int layer  = ND::TGeomInfo::P0D().GetBar(geomid).GetLayer();
    int bar    = ND::TGeomInfo::P0D().GetBar(geomid).GetNumber();

    if(p0dule > 38){
      p0dParticle->EndDeposit += (*hit)->GetCharge();
    }

    if(bar < 4 || (layer == 0 && bar > 121) || (layer == 1 && bar > 129)){
      p0dParticle->SideDeposit += (*hit)->GetCharge();
    }
  }

  // The  TReconPID  may  have  extra  PID  variables  stored  as  TReal-  and
  // TIntegerData.  Check to  see if they exist, and store  the values if they
  // do.

  // std::copy is required when the input vectors are doubles, but we
  // only store floats.  It is not required for integers, but it is
  // easier to keep the code the same.

  // Currently we only have TRealDatums as inputs (although some could
  // be stored as TIntegerDatums), but the code is included for
  // completeness.

  for(ND::TDataVector::iterator dvIter = particle->begin();
      dvIter != particle->end(); ++dvIter){
    ND::TRealDatum *rDatum = dynamic_cast<ND::TRealDatum*>(*dvIter);
    if(rDatum){
      // Save the datum name, but remove the //unnamed/ at the beginning.
      TString datumName = rDatum->GetFullName();
      Ssiz_t lastSlash = datumName.Last('/');
      TSubString pidAlgoName = datumName(lastSlash + 1, datumName.Length());

      std::vector<float> tempVector;
      std::copy(rDatum->begin(),
                rDatum->end(),
                std::back_inserter(tempVector));

      p0dParticle->realPIDNames.push_back(pidAlgoName.Data());
      p0dParticle->realPIDValues.push_back(tempVector);
    }


    ND::TIntegerDatum *iDatum = dynamic_cast<ND::TIntegerDatum*>(*dvIter);
    if(iDatum){
      // Save the datum name, but remove the //unnamed/ at the beginning.
      TString datumName = iDatum->GetFullName();
      Ssiz_t lastSlash = datumName.Last('/');
      TSubString pidAlgoName = datumName(lastSlash + 1, datumName.Length());

      std::vector<short> tempVector;
      std::copy(iDatum->begin(),
                iDatum->end(),
                std::back_inserter(tempVector));

      p0dParticle->integerPIDNames.push_back(pidAlgoName.Data());
      p0dParticle->integerPIDValues.push_back(tempVector);
    }
  }

  // The most  likely PID is  available through particle->GetParticle,
  // but it is also stored as one of the alternates.  These are sorted
  // in order of PID, so copy them directly into the TTree.
  ND::TReconObjectContainer::const_iterator alternate;

  for(alternate = particle->GetAlternates().begin();
      alternate != particle->GetAlternates().end();
      alternate++){

    ND::THandle<ND::TReconPID> alt = (*alternate);
    if(alt){
      p0dParticle->PID.push_back(alt->GetParticleId());
      p0dParticle->PID_weight.push_back(alt->GetPIDWeight());
    }
    else {ND280Warn("TReconPID alternate is not a TReconPID - Ignored");}
  }

  return particleID;
}

ND::TP0DReconModule::ObjectID
ND::TP0DReconModule::FillTrackObject(const ND::THandle<ND::TReconTrack> track,
                                     bool saveHits = true){
  ObjectID trackID;
  trackID.id   = fNTracks;
  trackID.type = ObjectID::kTrack;

  TP0DTrack* p0dTrack =
    new ((*fTracks)[fNTracks++]) TP0DTrack;

  FillBaseObject(p0dTrack, track, saveHits);

  ND::THandle<ND::TTrackState> trackState = track->GetState();

  p0dTrack->Status        = track->GetStatus();
  p0dTrack->Quality       = track->GetQuality();
  p0dTrack->NDOF          = track->GetNDOF();
  p0dTrack->Position      = trackState->GetPosition();
  p0dTrack->PosVariance   = trackState->GetPositionVariance();
  p0dTrack->Direction     = trackState->GetDirection();
  p0dTrack->DirVariance   = trackState->GetDirectionVariance();
  p0dTrack->EDeposit      = trackState->GetEDeposit();

  p0dTrack->ValidDimensions =
    CountValidDimensions(p0dTrack->PosVariance.Vect());

  const ND::THandle<ND::THitSelection> trackHits = track->GetHits();
  if (trackHits) {

    std::map< int, std::pair<int, float> > hitInfo = HitTruthInfo(trackHits);
    std::map< int, std::pair< int, float > >::iterator parentID;

    for(parentID = hitInfo.begin(); parentID != hitInfo.end(); ++parentID){
      p0dTrack->Truth_TrajIDs.push_back((*parentID).first);
      p0dTrack->Truth_HitCount.push_back((*parentID).second.first);
      p0dTrack->Truth_ChargeShare.push_back((*parentID).second.second);
    }

    p0dTrack->Truth_PrimaryTrajIDs = HitTruthPrimaryInfo(trackHits);
  }

  p0dTrack->SideDeposit = 0.0;
  p0dTrack->EndDeposit  = 0.0;

  ND::THitSelection::const_iterator hit;
  for (hit = trackHits->begin();
       hit != trackHits->end();
       ++hit) {

    ND::TGeometryId geomid = (*hit)->GetGeomId();
    if (!ND::GeomId::P0D::IsP0D(geomid)) continue;
    int p0dule = ND::TGeomInfo::P0D().GetBar(geomid).GetP0Dule();
    int layer  = ND::TGeomInfo::P0D().GetBar(geomid).GetLayer();
    int bar    = ND::TGeomInfo::P0D().GetBar(geomid).GetNumber();
    if(p0dule > 38){
      p0dTrack->EndDeposit += (*hit)->GetCharge();
    }
    if(bar < 4 || (layer == 0 && bar > 121) || (layer == 1 && bar > 129)){
      p0dTrack->SideDeposit += (*hit)->GetCharge();
    }

  }

  const ND::TReconNodeContainer& nodeCont = track->GetNodes();

  // Get the track length from the nodes
  p0dTrack->Length = 0.0;
  if (nodeCont.size() < 1) return trackID;

  ND::TReconNodeContainer::const_iterator node = nodeCont.begin();
  ND::THandle<ND::TTrackState> nodeState = (*node)->GetState();
  TVector3 frontPos = nodeState->GetPosition().Vect();
  ++node;

  while(node != nodeCont.end()){

    // Check the node is valid
    if (!(*node)) {
      ++node;
      continue;
    }
    nodeState = (*node)->GetState();

    p0dTrack->Length += (frontPos - nodeState->GetPosition().Vect()).Mag();
    frontPos = nodeState->GetPosition().Vect();

    ++node;
  }

  return trackID;
}

ND::TP0DReconModule::ObjectID
ND::TP0DReconModule::FillShowerObject(
  const ND::THandle<ND::TReconShower> shower,
  bool saveHits = true){

  ObjectID showerID;
  showerID.id   = fNShowers;
  showerID.type = ObjectID::kShower;

  TP0DShower* p0dShower =
    new ((*fShowers)[fNShowers++]) TP0DShower;

  FillBaseObject(p0dShower, shower, saveHits);

  // For Shower objects, the cluster information of the nodes is vital for the
  // pi0 analysis.  However, the cluster information for most nodes is
  // useless.  Therefore, rather than adding a "FillClusterObject" to each
  // "FillNode", we will do a parallel fill, within the
  // FillShowerObject. N.B. FillNode has already been run within
  // FillBaseObject.

  const ND::TReconNodeContainer& nodeCont = shower->GetNodes();
  ND::TReconNodeContainer::const_iterator node;
  for (node = nodeCont.begin();
       node != nodeCont.end();
       node++){
    ND::THandle<ND::TReconCluster> cluster = (*node)->GetObject();
    ND280Verbose("Fill Cluster - Shower");
    ObjectID clusterID = FillClusterObject(cluster, saveHits);
    p0dShower->Clusters.push_back(clusterID.id);
  }

  ND::THandle<ND::TShowerState> showerState = shower->GetState();

  p0dShower->Status        = shower->GetStatus();
  p0dShower->Quality       = shower->GetQuality();
  p0dShower->NDOF          = shower->GetNDOF();
  p0dShower->Position      = showerState->GetPosition();
  p0dShower->PosVariance   = showerState->GetPositionVariance();
  p0dShower->Direction     = showerState->GetDirection();
  p0dShower->DirVariance   = showerState->GetDirectionVariance();
  p0dShower->EDeposit      = showerState->GetEDeposit();
  p0dShower->Cone          = showerState->GetCone()[0];

  p0dShower->ValidDimensions =
    CountValidDimensions(p0dShower->PosVariance.Vect());

  p0dShower->Width         = 0.0;
  p0dShower->Length        = 0.0;
  p0dShower->SideDeposit   = 0.0;
  p0dShower->EndDeposit    = 0.0;

  const ND::THandle<ND::THitSelection> showerHits = shower->GetHits();

  double totalCharge = 0.0;
  float avgLen = 0.0;
  float avgLen2 = 0.0;

  if (showerHits) {
    std::map< int, std::pair<int, float> > hitInfo = HitTruthInfo(showerHits);

    std::map< int, std::pair< int, float > >::iterator parentID;
    for(parentID = hitInfo.begin(); parentID != hitInfo.end(); ++parentID){
      p0dShower->Truth_TrajIDs.push_back((*parentID).first);
      p0dShower->Truth_HitCount.push_back((*parentID).second.first);
      p0dShower->Truth_ChargeShare.push_back((*parentID).second.second);
    }

    p0dShower->Truth_PrimaryTrajIDs = HitTruthPrimaryInfo(showerHits);

    ND::THitSelection::const_iterator hit;

    for (hit = showerHits->begin();
         hit != showerHits->end();
         ++hit) {

      // The distance vector between the hit and the shower position
      TVector3 diff = (*hit)->GetPosition() - p0dShower->Position.Vect();

      // The distance vector along the shower direction
      TVector3 para = (diff * p0dShower->Direction) * p0dShower->Direction;

      // The distance vector normal to the shower direction
      TVector3 perp = diff - para;

      if ((*hit)->IsXHit()) {
        diff.SetY(0.0);
        perp.SetY(0.0);
      }
      else {
        diff.SetX(0.0);
        perp.SetX(0.0);
      }

      float q =  (*hit)->GetCharge();
      totalCharge += q;

      p0dShower->Width += q * perp.Mag();

      avgLen  += q * para.Mag();
      avgLen2 += q * (para * para);

      // Get the side and end charge
      ND::TGeometryId geomid = (*hit)->GetGeomId();
      if (!ND::GeomId::P0D::IsP0D(geomid)) {continue;}
      int p0dule = ND::TGeomInfo::P0D().GetBar(geomid).GetP0Dule();
      int layer  = ND::TGeomInfo::P0D().GetBar(geomid).GetLayer();
      int bar    = ND::TGeomInfo::P0D().GetBar(geomid).GetNumber();

      if(p0dule > 38){
        p0dShower->EndDeposit += (*hit)->GetCharge();
      }
      if(bar < 4 || (layer == 0 && bar > 121) || (layer == 1 && bar > 129)){
        p0dShower->SideDeposit += (*hit)->GetCharge();
      }

    }
  }

  p0dShower->Width /= std::max(totalCharge,1.0);

  avgLen  /= std::max(totalCharge,1.0);
  avgLen2 /= std::max(totalCharge,1.0);
  double len2 = avgLen2 - avgLen*avgLen;

  p0dShower->Length = std::sqrt(std::max(0.0, len2));

  return showerID;
}

ND::TP0DReconModule::ObjectID
ND::TP0DReconModule::FillClusterObject(
  const ND::THandle<ND::TReconCluster> cluster,
  bool saveHits = true){

  ObjectID clusterID;
  clusterID.id   = fNClusters;
  clusterID.type = ObjectID::kCluster;

  TP0DCluster* p0dCluster =
    new ((*fClusters)[fNClusters++]) TP0DCluster;

  FillBaseObject(p0dCluster,cluster, saveHits);

  p0dCluster->NFiducialHits = -1;
  p0dCluster->EDeposit      = cluster->GetEDeposit();
  p0dCluster->Position      = cluster->GetPosition();
  p0dCluster->PosVariance   = cluster->GetPositionVariance();

  p0dCluster->ValidDimensions =
    CountValidDimensions(p0dCluster->PosVariance.Vect());

  TMatrixF moments = cluster->GetMoments();
  if( p0dCluster->arraySize ==
      (moments.GetNrows() * moments.GetNcols()) ){
    std::copy(moments.GetMatrixArray(),
              moments.GetMatrixArray()+p0dCluster->arraySize,
              p0dCluster->Moments);
  }

  const ND::THandle<ND::THitSelection> clusterHits = cluster->GetHits();
  if (clusterHits){

    std::map< int, std::pair<int, float> > hitInfo = HitTruthInfo(clusterHits);
    std::map< int, std::pair< int, float > >::iterator parentID;

    for(parentID = hitInfo.begin(); parentID != hitInfo.end(); ++parentID){
      p0dCluster->Truth_TrajIDs.push_back((*parentID).first);
      p0dCluster->Truth_HitCount.push_back((*parentID).second.first);
      p0dCluster->Truth_ChargeShare.push_back((*parentID).second.second);
    }

    p0dCluster->Truth_PrimaryTrajIDs = HitTruthPrimaryInfo(clusterHits);
    p0dCluster->NFiducialHits = 0;

    for (ND::THitSelection::const_iterator hit = clusterHits->begin();
         hit != clusterHits->end(); ++hit){
      TVector3 hitPosition = (*hit)->GetPosition();
      float fiducial =
      ND::TGeomInfo::P0D().WaterFiducialDistance(hitPosition.X(),
                                                 hitPosition.Y(),
                                                 hitPosition.Z());

      if (fiducial > 0){p0dCluster->NFiducialHits += 1;}
    }
  }
  return clusterID;
}

short ND::TP0DReconModule::FillNode(
  const ND::THandle<ND::TReconNode> node,
  bool saveHits){

  short nodeID = fNNodes;
  ND280Verbose("Fill Node");

  TP0DNode* p0dNode =
    new ((*fNodes)[fNNodes++]) TP0DNode;

  // The node state may be a PIDState, TrackState or ShowerState, but
  // they are all TMPositionDirectionStates, so we access the position
  // and direction regardless.
  ND::THandle<ND::TReconState> s1 = node->GetState();
  if (s1) {
    #if BEFORE_oaEvent(9,0,0)
    const ND::TMPositionDirectionState* posDirState = NULL;
    posDirState =
      dynamic_cast<const ND::TMPositionDirectionState*>(ND::GetPointer(s1));
    if(posDirState){
      p0dNode->Position      = posDirState->GetPosition();
      p0dNode->PosVariance   = posDirState->GetPositionVariance();
      p0dNode->Direction     = posDirState->GetDirection();
      p0dNode->DirVariance   = posDirState->GetDirectionVariance();
    #else
      p0dNode->Position      = GetPosition(s1);
      p0dNode->PosVariance   = GetPositionVariance(s1);
      p0dNode->Direction     = ND::TPIDState::GetStateDirection(s1);
      p0dNode->DirVariance   = ND::TPIDState::GetStateDirectionVariance(s1);
    #endif
      p0dNode->ValidDimensions =
        CountValidDimensions(p0dNode->PosVariance.Vect());
    #if BEFORE_oaEvent(9,0,0)
    }
    #endif
  }

  // Every node should have an object

  ND::THandle<ND::TReconBase> nodeObject = node->GetObject();
  if(nodeObject){
    ND::THandle<ND::TReconState> objState = nodeObject->GetState();
    if(objState){
    #if BEFORE_oaEvent(9,0,0)
      const ND::TMEDepositState* eDepState = NULL;
      eDepState =
      dynamic_cast<const ND::TMEDepositState*>(ND::GetPointer(s1));
      if(eDepState){
        p0dNode->EDeposit = eDepState->GetEDeposit();
      }
     #else 
        p0dNode->EDeposit = ND::TPIDState::GetStateEDeposit(objState);
     #endif
    }

    const ND::THandle<ND::THitSelection> nodeHits = nodeObject->GetHits();
    if (nodeHits){
      ND::THitSelection::const_iterator hit;

      for(hit = nodeHits->begin();
          hit != nodeHits->end(); ++hit){

        short hitId = FillHit(*hit);
        if(hitId >= 0){
          p0dNode->Hits.push_back(hitId);
        }
      }

      std::map< int, std::pair<int, float> > hitInfo = HitTruthInfo(nodeHits);
      std::map< int, std::pair< int, float > >::iterator parentID;

      for(parentID = hitInfo.begin(); parentID != hitInfo.end(); ++parentID){
        p0dNode->Truth_TrajIDs.push_back((*parentID).first);
        p0dNode->Truth_HitCount.push_back((*parentID).second.first);
        p0dNode->Truth_ChargeShare.push_back((*parentID).second.second);
      }

      p0dNode->Truth_PrimaryTrajIDs = HitTruthPrimaryInfo(nodeHits);
    }
  }
  return nodeID;
}

short
ND::TP0DReconModule::FillHit(const ND::THandle<ND::THit> hit){
  if (!hit) {return -1;}
  ND::THandle<ND::TSingleHit> sh = hit;
  if (!sh) {return -1;}

  // Check if this Hit has already been filled, if so pass back the
  // reference, instead of creating a new one.
  UInt_t chanIdUInt = sh->GetChannelId().AsUInt();
  std::map<UInt_t, short>::const_iterator mapIter;
  mapIter = fTempHitMap.find(chanIdUInt);
  if( mapIter != fTempHitMap.end() ){
    return mapIter->second;
  }

  // We only want to save P0D hits at this time.  If a hit is from
  // outside the P0D, skip and return -1.
  if(sh->GetChannelId().GetSubDetector() != ND::TChannelId::kP0D){
    return -1;
  }

  // Otherwise, make a new entry and save it in the TempHitMap
  short hitID = fNHits;
  fTempHitMap[chanIdUInt] = hitID;

  TP0DHit* p0dHit =
    new ((*fHits)[fNHits++]) TP0DHit;

  p0dHit->ChanID = sh->GetChannelId().AsUInt();
  p0dHit->GeomID = sh->GetGeomId().AsInt();
  p0dHit->Charge = sh->GetCharge();
  p0dHit->Time   = sh->GetTime();

  return hitID;
}

std::map< int, std::pair<int, float> >
ND::TP0DReconModule::HitTruthInfo(const ND::THandle<ND::THitSelection> hits){

  std::map< int, std::pair<int, float> > result;

  // then loop over everything and log the trajectory ids.
  for (ND::THitSelection::const_iterator hit = hits->begin();
       hit != hits->end(); hit++) {

    float hitCharge = (*hit)->GetCharge();
    std::map< int, int > hitCount;
    std::map< int, float > chargeShare;

    std::vector<ND::TG4VHit*> hitContribs = HitTruthInfo::GetHitTruthInfo(*hit);

    for (std::vector<ND::TG4VHit*>::const_iterator g4Hit = hitContribs.begin();
         g4Hit != hitContribs.end(); g4Hit++) {

      ND::TG4HitSegment* g4HitCast = dynamic_cast<ND::TG4HitSegment*>(*g4Hit);

      // For each contributing trajectory, increase the hit count and
      // give it an equal share of the energy deposit (the share may
      // not be equal, but that information is not currently
      // available).
      std::vector<int>::const_iterator trajIDIter;
      for (trajIDIter  = g4HitCast->GetContributors().begin();
           trajIDIter != g4HitCast->GetContributors().end();
           ++trajIDIter){

        hitCount[(*trajIDIter)]++;
        chargeShare[(*trajIDIter)] +=
          g4HitCast->GetEnergyDeposit() / g4HitCast->GetContributors().size();
      }
    }

    for (std::map< int, int >::const_iterator hC = hitCount.begin();
         hC != hitCount.end(); ++hC){
      // If this hit  had a contribution from a given  ID, then increment that
      // ID's hitcount.
      if( (*hC).second > 0) {result[(*hC).first].first++;}
    }

    // Count up the  contributed charge (this won't necessarily  match the hit
    // Charge).
    float totalChargeShare = 0;
    for (std::map< int, float >::const_iterator cS = chargeShare.begin();
         cS != chargeShare.end(); ++cS){
      totalChargeShare += (*cS).second;
    }

    // For each entry, determine the contributed share.
    for (std::map< int, float >::const_iterator cS = chargeShare.begin();
         cS != chargeShare.end(); ++cS){
      result[(*cS).first].second += hitCharge * (*cS).second / totalChargeShare;
    }
  }
  return result;
}

std::vector<int>
ND::TP0DReconModule::HitTruthPrimaryInfo(
  const ND::THandle<ND::THitSelection> hits){

  std::vector<int> result;
  // Loop over all hits, storing each primary ID which contributed
  for (ND::THitSelection::const_iterator hit = hits->begin();
       hit != hits->end(); hit++) {
    std::vector<ND::TG4VHit*> hitContribs = HitTruthInfo::GetHitTruthInfo(*hit);

    for (std::vector<ND::TG4VHit*>::const_iterator g4Hit = hitContribs.begin();
         g4Hit != hitContribs.end(); g4Hit++) {
      ND::TG4HitSegment* g4HitCast = dynamic_cast<ND::TG4HitSegment*>(*g4Hit);

      result.push_back(g4HitCast->GetPrimaryId());
    }
  }

  // Sort the list of primary IDs
  std::sort(result.begin(), result.end());

  // Remove duplicate IDs
  std::vector<int>::iterator newEnd =
    std::unique(result.begin(), result.end());
  result.erase(newEnd, result.end());

  return result;
}

short
ND::TP0DReconModule::GetCycle(const ND::THandle<ND::THitSelection> hits){
  // As of July 2012, the p0dRecon results may contain objects from
  // outside the P0D.  Currently, this is only due to the Incremental
  // matching with the TPC.  Therefore, we consider all potential
  // hits, looking for a TFB hit.

  // Access the channel ID for each hit
  ND::TTFBChannelId chanID;//(hits->front()->GetChannelId());
  ND::THitSelection::const_iterator hit;

  for(hit = hits->begin(); hit != hits->end(); ++hit){
    ND::THandle<ND::TSingleHit> sh = (*hit);
    if (!sh){continue;}
    chanID = ND::TTFBChannelId(sh->GetChannelId());
    if( chanID.IsValid() && chanID.IsTFBChannel() &&
        (chanID.GetSubDetector() == ND::TChannelId::kP0D) ){
      break;
    }
  }

  // If the ID is invalid, return -1 as the cycle
  if(!chanID.IsValid() || !chanID.IsTFBChannel()){
    //ND280Error("Invalid P0D Event containing no cycles!");
    return -1;
  }

  // Otherwise, return the cycle of the channel ID
  short cycle = (short) chanID.GetCapacitor();
  return cycle;
}

short
ND::TP0DReconModule::CountValidDimensions(TVector3 posVar){
  short valid = 0;
  if (std::abs(posVar.X()) < 1E+6) {valid += 1;}
  if (std::abs(posVar.Y()) < 1E+6) {valid += 2;}
  if (std::abs(posVar.Z()) < 1E+6) {valid += 4;}
  return valid;
}