TSFGReconModule.cxx

Go to the documentation of this file.

#include "TSFGReconModule.hxx"

#include <TAlgorithmResult.hxx>
#include <TG4PrimaryVertex.hxx>
#include <TGeomInfo.hxx>
#include <TGeometryId.hxx>
#include <THandle.hxx>
#include <THitSelection.hxx>
#include <TND280Event.hxx>
#include <TND280Log.hxx>
#include <TSFGGeom.hxx>
#include <TReconNode.hxx>
#include <TReconCluster.hxx>
#include <TReconShower.hxx>
#include <TReconTrack.hxx>
#include <TReconVertex.hxx>
#include <TReconPID.hxx>
#include <TReconState.hxx>
#include <TShowerState.hxx>
#include <TTrackState.hxx>
#include <TVertexState.hxx>
#include <TPIDState.hxx>
#include "TEventFolder.hxx"
#include <TrackTruthInfo.hxx>

#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

#include <TG4VHit.hxx>
#include <TG4HitSegment.hxx>
#include <retrieveHitTruthInfo.hxx>
#include <TG4Trajectory.hxx>

#include <TString.h>
#include <TVector3.h>

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

ClassImp(ND::TSFGReconModule::TSFGAlgoRes);
ND::TSFGReconModule::TSFGAlgoRes::~TSFGAlgoRes() {}

ClassImp(ND::TSFGReconModule::TSFGVertex);
ND::TSFGReconModule::TSFGVertex::~TSFGVertex() {}

ClassImp(ND::TSFGReconModule::TSFGParticle);
ND::TSFGReconModule::TSFGParticle::~TSFGParticle() {}

ClassImp(ND::TSFGReconModule::TSFGTrack);
ND::TSFGReconModule::TSFGTrack::~TSFGTrack() {}

ClassImp(ND::TSFGReconModule::TSFGShower);
ND::TSFGReconModule::TSFGShower::~TSFGShower() {}

ClassImp(ND::TSFGReconModule::TSFGNode);
ND::TSFGReconModule::TSFGNode::~TSFGNode() {}

ClassImp(ND::TSFGReconModule::TSFGCluster);
ND::TSFGReconModule::TSFGCluster::~TSFGCluster() {}

ClassImp(ND::TSFGReconModule::TSFGHit);
ND::TSFGReconModule::TSFGHit::~TSFGHit() {}

ND::TSFGReconModule::TSFGReconModule(const char *name,
                                     const char *title) {

    SetNameTitle(name, title);

    // Enable this module by default:
    fIsEnabled = kTRUE;
    fDescription = "SFG Recon Output";
    fCVSTagName = "FromGit";
    fCVSID = "FromGit";

    fAlgoResults = new TClonesArray("ND::TSFGReconModule::TSFGAlgoRes");
    fVertices = new TClonesArray("ND::TSFGReconModule::TSFGVertex");
    fParticles = new TClonesArray("ND::TSFGReconModule::TSFGParticle");
    fTracks = new TClonesArray("ND::TSFGReconModule::TSFGTrack");
    fShowers = new TClonesArray("ND::TSFGReconModule::TSFGShower");
    fClusters = new TClonesArray("ND::TSFGReconModule::TSFGCluster");
    fNodes = new TClonesArray("ND::TSFGReconModule::TSFGNode");
    fHits = new TClonesArray("ND::TSFGReconModule::TSFGHit");
    fTrueHits = new TClonesArray("ND::TSFGReconModule::TSFGHit");
    fFibers = new TClonesArray("ND::TSFGReconModule::TSFGHit");

    fRejectAlgoResultList.clear();
    fRejectAlgoResultList.push_back(TPRegexp("TSFGBogusExample$"));
}


ND::TSFGReconModule::~TSFGReconModule() {}

void ND::TSFGReconModule::InitializeBranches() {

    TBranch* branch = fOutputTree->Branch("NAlgoResults", &fNAlgoResults,
                                          "NAlgoResults/I", fBufferSize);
    branch ->SetTitle("The number of Algorithm Results ");

    branch = fOutputTree->Branch("AlgoResults", &fAlgoResults,
                                 fBufferSize, fSplitLevel);
    branch ->SetTitle("The TSFGAlgoRes vector of Algorithm Results");

    branch = fOutputTree->Branch("NVertices", &fNVertices,
                                 "NVertices/I", fBufferSize);
    branch->SetTitle(" The number of added vertices");

    branch = fOutputTree->Branch("Vertices", &fVertices,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("  The vector of vertices");

    branch = fOutputTree->Branch("NParticles", &fNParticles,
                                 "NParticles/I", fBufferSize);
    branch->SetTitle(" The number of particles ");

    branch = fOutputTree->Branch("Particles", &fParticles,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of particles");

    branch = fOutputTree->Branch("NTracks", &fNTracks,
                                 "NTracks/I", fBufferSize);
    branch->SetTitle("  The number of tracks");

    branch = fOutputTree->Branch("Tracks", &fTracks,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle(" The vector of tracks");

    branch = fOutputTree->Branch("NShowers", &fNShowers,
                                 "NShowers/I", fBufferSize);
    branch->SetTitle(" The number of showers");

    branch = fOutputTree->Branch("Showers", &fShowers,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of showers");

    branch = fOutputTree->Branch("NClusters", &fNClusters,
                                 "NClusters/I", fBufferSize);
    branch ->SetTitle("The number of clusters ");

    branch = fOutputTree->Branch("Clusters", &fClusters,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of clusters");

    branch = fOutputTree->Branch("NNodes", &fNNodes,
                                 "NNodes/I", fBufferSize);
    branch->SetTitle(" The number of track nodes");

    branch = fOutputTree->Branch("Nodes", &fNodes,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of nodes");

    branch = fOutputTree->Branch("NHits", &fNHits,
                                 "NHits/I", fBufferSize);
    branch->SetTitle("The number of hits that are saved ");

    branch = fOutputTree->Branch("Hits", &fHits,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of reconstructed hits");

    branch = fOutputTree->Branch("NFibers", &fNFibers,
                                 "NFiber/I", fBufferSize);
    branch->SetTitle("The number of fiber hits");

    branch = fOutputTree->Branch("Fibers", &fFibers,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of fiber hits");

    branch = fOutputTree->Branch("NTrueHits", &fNTrueHits,
                                 "NTrueHits/I", fBufferSize);
    branch->SetTitle("The energy deposition from the simulation");

    branch = fOutputTree->Branch("TrueHits", &fTrueHits,
                                 fBufferSize, fSplitLevel);
    branch->SetTitle("The vector of true hits");

}


bool ND::TSFGReconModule::FillTree(ND::TND280Event& event) {

    fNAlgoResults = 0;
    fNVertices = 0;
    fNParticles = 0;
    fNTracks = 0;
    fNShowers = 0;
    fNClusters = 0;
    fNNodes = 0;
    fNHits = 0;
    fNFibers = 0;
    fNTrueHits = 0;

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

    ND::THandle<ND::THitSelection> sfgFibers(event.GetHitSelection("sfg"));
    if (sfgFibers) {
        for (ND::THitSelection::iterator h = sfgFibers->begin();
             h != sfgFibers->end(); ++h) {
            AddHit(*h,fNFibers,fFibers);
        }
    }
    else {
        ND280Info("Missing fiber hits");
    }

    ND::THandle<ND::THitSelection> sfgTruth(event.GetHitSelection("sfg_truth"));
    if (sfgTruth) {
        for (ND::THitSelection::iterator h = sfgTruth->begin();
             h != sfgTruth->end(); ++h) {
            AddHit(*h,fNTrueHits,fTrueHits);
        }
    }
    else {
        ND280Info("Missing truth hits");
    }

    // Get the true trajectories
    fTrueTraj = event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
    if(!fTrueTraj){
      ND280Log("Missing truth trajectories");
    }

    // Get the main result to save.
    ND::THandle<ND::TAlgorithmResult> sfgRecon = event.GetFit("TSFGRecon");
    if (!sfgRecon) {
        sfgRecon = event.GetFit("sfgRecon");
        if (!sfgRecon) {
            ND280Info("No sfgRecon in Event");
            return true;
        }
    }

    // Fill the trees.
    fTempHitMap.clear();
    FillAlgorithmResult(sfgRecon, -1);
    fTempHitMap.clear();

    return true;
}

template <class T>
void ND::TSFGReconModule::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->NHits = (int)baseHits->size();

        if (saveHits) {
            for (ND::THitSelection::const_iterator hit = baseHits->begin();
                 hit != baseHits->end(); ++hit) {
                int hitId = FillHit(*hit);
                if (hitId < 0) continue;
                basePtr->Hits.push_back(hitId);
            }
        }
    }
    else {
        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);
    }

};

int
ND::TSFGReconModule::FillAlgorithmResult(
    const ND::THandle<ND::TAlgorithmResult> algoRes,
    int 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 level
    // 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: " << baseName
                         << " (" << cut->GetPattern() << ")");
            return -1;
        }
    }

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

    int algoResId = fNAlgoResults;
    TSFGAlgoRes *sfgAlgoResult
        = new ((*fAlgoResults)[fNAlgoResults++]) TSFGAlgoRes;

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

    // Get the final recon object container and save all of the objects.  This
    // needs to be first so that the objects for the main algorithm are the
    // first thing in all of the vectors.
    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) {
                sfgAlgoResult->Vertices.push_back(objId.id);
            }
            else if (objId.type == ObjectId::kParticle) {
                sfgAlgoResult->Particles.push_back(objId.id);
            }
            else if (objId.type == ObjectId::kTrack) {
                sfgAlgoResult->Tracks.push_back(objId.id);
            }
            else if (objId.type == ObjectId::kShower) {
                sfgAlgoResult->Showers.push_back(objId.id);
            }
            else if (objId.type == ObjectId::kCluster) {
                sfgAlgoResult->Clusters.push_back(objId.id);
            }
        }
    }

    // Get the used hits from the algo result, and save them as a cluster.
    // This ends up right after the reconstructed clusters in the cluster
    // vector.
    ND::THandle<ND::THitSelection> usedHits = algoRes->GetHitSelection("used");
    if (usedHits && (!usedHits->empty())) {
        ND::THitSelection::const_iterator hit;

        for (hit = usedHits->begin();
             hit != usedHits->end(); ++hit) {
            int hitId = FillHit(*hit);
            if (hitId < 0) continue;
            sfgAlgoResult->Hits.push_back(hitId);
        }
        ND::THandle<ND::TReconCluster> usedCluster(new ND::TReconCluster);
        usedCluster->FillFromHits((baseName+"Used").Data(), *usedHits);

        ObjectId clusterId = FillReconObject(usedCluster, true);
        sfgAlgoResult->UsedHitCluster = clusterId.id;
    }
    else {
        sfgAlgoResult->UsedHitCluster     = -1;
    }

    // Get the unused hits from the algo result, and save them as a cluster.
    // This ends up right after the used hit cluster in the cluster
    // vector.
    ND::THandle<ND::THitSelection> unusedHits
        = algoRes->GetHitSelection("unused");
    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);
        sfgAlgoResult->UnusedHitCluster = clusterId.id;
    }
    else {
        sfgAlgoResult->UnusedHitCluster     = -1;
    }

    // Handle any daughter algorithms.  They are processed in reverse order so
    // that the "most important one" is first.
    for (TDataVector::reverse_iterator d = algoRes->rbegin();
         d != algoRes->rend(); ++d) {
        if (!(*d)) continue;

        std::string class_name = (*d)->ClassName();
        if (class_name.find("ND::TAlgorithmResult")
            == std::string::npos) continue;

        ND::THandle<ND::TAlgorithmResult>
            subResult = (*d)->Get<ND::TAlgorithmResult>(".");
        if (!subResult) continue;

        int daughterId = FillAlgorithmResult(subResult, algoResId);

        if (daughterId < 0) continue;

        sfgAlgoResult->AlgoResults.push_back(daughterId);
    }

    return algoResId;
}

ND::TSFGReconModule::ObjectId
ND::TSFGReconModule::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) {
        return FillVertexObject(reconObject, saveHits);
    }
    else if (reco_class_name.find("ND::TReconPID") != std::string::npos) {
        return FillParticleObject(reconObject, saveHits);
    }
    else if (reco_class_name.find("ND::TReconTrack") != std::string::npos) {
        return FillTrackObject(reconObject, saveHits);
    }
    else if (reco_class_name.find("ND::TReconShower") != std::string::npos) {
        return FillShowerObject(reconObject, saveHits);
    }
    else if (reco_class_name.find("ND::TReconCluster") != std::string::npos) {
        return FillClusterObject(reconObject, saveHits);
    }

    ND280Error("Unknown object passed to TSFGReconModule");
    ObjectId blank;
    return blank;
}

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

    ObjectId vertexId;
    vertexId.id   = fNVertices;
    vertexId.type = ObjectId::kVertex;
    TSFGVertex* sfgVertex = new ((*fVertices)[fNVertices++]) TSFGVertex;

    FillBaseObject(sfgVertex, vertex, saveHits);

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

    sfgVertex->Status        = vertex->GetStatus();
    sfgVertex->Quality       = vertex->GetQuality();
    sfgVertex->NDOF          = vertex->GetNDOF();
    sfgVertex->Position      = vertexState->GetPosition();
    sfgVertex->PosVariance   = vertexState->GetPositionVariance();

    sfgVertex->Truth_TrajIds.clear();
    sfgVertex->Truth_HitCount.clear();
    sfgVertex->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) {
            sfgVertex->Truth_TrajIds.push_back((*parentId).first);
            sfgVertex->Truth_HitCount.push_back((*parentId).second.first);
            sfgVertex->Truth_ChargeShare.push_back((*parentId).second.second);
        }

        sfgVertex->Truth_PrimaryTrajIds = HitTruthPrimaryInfo(vertexHits);
    }

    return vertexId;
}

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

    ObjectId particleId;
    particleId.id   = fNParticles;
    particleId.type = ObjectId::kParticle;
    TSFGParticle* sfgParticle = new ((*fParticles)[fNParticles++]) TSFGParticle;

    FillBaseObject(sfgParticle, particle, saveHits);

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

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

    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) {
            sfgParticle->Truth_TrajIds.push_back((*parentId).first);
            sfgParticle->Truth_HitCount.push_back((*parentId).second.first);
            sfgParticle->Truth_ChargeShare.push_back((*parentId).second.second);
        }

        sfgParticle->Truth_PrimaryTrajIds = HitTruthPrimaryInfo(particleHits);
    }

    // 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) {
            sfgParticle->PID.push_back(alt->GetParticleId());
            sfgParticle->PID_weight.push_back(alt->GetPIDWeight());
        }
        else {
            ND280Warn("TReconPID alternate is not a TReconPID - Ignored");
        }
    }

    // Fill the true particle when it is found
    double pur, eff;
    ND::THandle<ND::TG4Trajectory> G4track = GetG4Trajectory(*particle, pur, eff);
    if(G4track){
      sfgParticle->TrueParticleID = G4track->GetTrackId();
      sfgParticle->TrueParticlePur = pur;
      sfgParticle->TrueParticleEff = eff;
    }
    else{
      sfgParticle->TrueParticleID = -1;
      sfgParticle->TrueParticlePur = -1;
      sfgParticle->TrueParticleEff = -1;
    }

    return particleId;
}

ND::TSFGReconModule::ObjectId
ND::TSFGReconModule::FillTrackObject(const ND::THandle<ND::TReconTrack> track,
                                     bool saveHits = true) {
    ObjectId trackId;
    trackId.id   = fNTracks;
    trackId.type = ObjectId::kTrack;
    TSFGTrack* sfgTrack = new ((*fTracks)[fNTracks++]) TSFGTrack;

    FillBaseObject(sfgTrack, track, saveHits);

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

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

    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) {
            sfgTrack->Truth_TrajIds.push_back((*parentId).first);
            sfgTrack->Truth_HitCount.push_back((*parentId).second.first);
            sfgTrack->Truth_ChargeShare.push_back((*parentId).second.second);
        }

        sfgTrack->Truth_PrimaryTrajIds = HitTruthPrimaryInfo(trackHits);
    }

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

    // Get the track length from the nodes
    sfgTrack->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();

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

        ++node;
    }

    // Fill the true particle when it is found
    double pur, eff;
    ND::THandle<ND::TG4Trajectory> G4track = GetG4Trajectory(*track, pur, eff);
    if(G4track){
      sfgTrack->TrueParticleID = G4track->GetTrackId();
      sfgTrack->TrueParticlePur = pur;
      sfgTrack->TrueParticleEff = eff;
    }
    else{
      sfgTrack->TrueParticleID = -1;
      sfgTrack->TrueParticlePur = -1;
      sfgTrack->TrueParticleEff = -1;
    }

    return trackId;
}

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

    ObjectId showerId;
    showerId.id   = fNShowers;
    showerId.type = ObjectId::kShower;
    TSFGShower* sfgShower = new ((*fShowers)[fNShowers++]) TSFGShower;

    FillBaseObject(sfgShower, 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();
    for (ND::TReconNodeContainer::const_iterator node = nodeCont.begin();
         node != nodeCont.end(); ++node) {
        ND::THandle<ND::TReconCluster> cluster = (*node)->GetObject();
        ObjectId clusterId = FillClusterObject(cluster, saveHits);
        sfgShower->Clusters.push_back(clusterId.id);
    }

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

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

    sfgShower->Width         = 0.0;
    sfgShower->Length        = 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);

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

        sfgShower->Truth_PrimaryTrajIds = HitTruthPrimaryInfo(showerHits);


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

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

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

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

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

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

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

        }
    }

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

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

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

    return showerId;
}

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

    ObjectId clusterId;
    clusterId.id   = fNClusters;
    clusterId.type = ObjectId::kCluster;
    TSFGCluster* sfgCluster = new ((*fClusters)[fNClusters++]) TSFGCluster;

    FillBaseObject(sfgCluster, cluster, saveHits);

    sfgCluster->EDeposit      = cluster->GetEDeposit();
    sfgCluster->Position      = cluster->GetPosition();
    sfgCluster->PosVariance   = cluster->GetPositionVariance();

    TMatrixF moments = cluster->GetMoments();
    if (sfgCluster->arraySize ==
        (moments.GetNrows() * moments.GetNcols()) ) {
        std::copy(moments.GetMatrixArray(),
                  moments.GetMatrixArray()+sfgCluster->arraySize,
                  sfgCluster->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) {
            sfgCluster->Truth_TrajIds.push_back((*parentId).first);
            sfgCluster->Truth_HitCount.push_back((*parentId).second.first);
            sfgCluster->Truth_ChargeShare.push_back((*parentId).second.second);
        }

        sfgCluster->Truth_PrimaryTrajIds = HitTruthPrimaryInfo(clusterHits);
    }

    // Fill the true particle when it is found
    double pur, eff;
    ND::THandle<ND::TG4Trajectory> G4track = GetG4Trajectory(*cluster, pur, eff);
    if(G4track){
      sfgCluster->TrueParticleID = G4track->GetTrackId();
      sfgCluster->TrueParticlePur = pur;
      sfgCluster->TrueParticleEff = eff;
    }
    else{
      sfgCluster->TrueParticleID = -1;
      sfgCluster->TrueParticlePur = -1;
      sfgCluster->TrueParticleEff = -1;
    }

    return clusterId;
}

int ND::TSFGReconModule::FillNode(
    const ND::THandle<ND::TReconNode> node,
    bool saveHits) {

    int nodeId = fNNodes;
    TSFGNode* sfgNode = new ((*fNodes)[fNNodes++]) TSFGNode;

   
    


    // 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){
      sfgNode->Position      = posDirState->GetPosition();
      sfgNode->PosVariance   = posDirState->GetPositionVariance();
      sfgNode->Direction     = posDirState->GetDirection();
      sfgNode->DirVariance   = posDirState->GetDirectionVariance();
      
      /*const ND::TMEDepositState* eDepState = NULL;
      eDepState = dynamic_cast<const ND::TMEDepositState*>(ND::GetPointer(s1));
      if(eDepState){
        sfgNode->EDeposit_fit = eDepState->GetEDeposit();
      }*/
    }
    #else
      sfgNode->Position      = GetPosition(s1);
      sfgNode->PosVariance   = GetPositionVariance(s1);
      sfgNode->Direction     = ND::TPIDState::GetStateDirection(s1);
      sfgNode->DirVariance   = ND::TPIDState::GetStateDirectionVariance(s1);
      
      //sfgNode->EDeposit_fit = ND::TPIDState::GetStateEDeposit(s1);
    #endif
   
    }

    // Get fitted energy deposit from node state
    ND::THandle<ND::TTrackState> s2 = node->GetState();
    if(s2){
      sfgNode->EDeposit_fit = s2->GetEDeposit();
    }

    // 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) {
                sfgNode->EDeposit = eDepState->GetEDeposit();
            }

        }
     #else 
        sfgNode->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) {

                int hitId = FillHit(*hit);
                if (hitId >= 0) {
                    sfgNode->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) {
                sfgNode->Truth_TrajIds.push_back((*parentId).first);
                sfgNode->Truth_HitCount.push_back((*parentId).second.first);
                sfgNode->Truth_ChargeShare.push_back((*parentId).second.second);
            }

            sfgNode->Truth_PrimaryTrajIds = HitTruthPrimaryInfo(nodeHits);

        }
    }

    return nodeId;
}

int
ND::TSFGReconModule::FillHit(const ND::THandle<ND::THit> hit) {
    if (!hit) {
        ND280Error("Trying to handle illegal hit");
        return -1;
    }

    // Check if this Hit has already been filled, if so pass back the
    // index, instead of creating a new one.
    HitMap::const_iterator mapIter;
    mapIter = fTempHitMap.find(hit);
    if ( mapIter != fTempHitMap.end() ) {
        return mapIter->second;
    }

    // Save the index of the location for the new hit.
    int hitId = AddHit(hit,fNHits,fHits);
    if (hitId < 0) return hitId;
    fTempHitMap[hit] = hitId;

    return hitId;
}

int
ND::TSFGReconModule::AddHit(const ND::THandle<ND::THit> hit,
                            int& nDest,
                            ND::TSFGReconModule::TSFGHitContainer* dest) {
    if (!hit) {
        ND280Error("Trying to handle illegal hit");
        return -1;
    }

    int hitId = nDest;
    TSFGHit *sfgHit = new ((*dest)[nDest++]) TSFGHit;

    sfgHit->GeomId = hit->GetGeomId().AsInt();
    sfgHit->Charge = hit->GetCharge();
    sfgHit->Time   = hit->GetTime();
    sfgHit->Position = hit->GetPosition();

    // Add information that is needed for neural network study
    // For the moment need track hit tagging and hit segment
    FillHitInfoNeededForNN(hit,sfgHit);
   
    return hitId;
}

std::map< int, std::pair<int, float> >
ND::TSFGReconModule::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::TSFGReconModule::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;
}

// ---- For SFGD neural network study ----

void ND::TSFGReconModule::FillHitInfoNeededForNN(const ND::THandle<ND::THit> hit, TSFGHit *sfgHit){

    // Clean the vectors
    sfgHit->HitSegPosition.clear();
    sfgHit->HitSegDirection.clear();
    sfgHit->HitSegTruePDG.clear();
    sfgHit->HitSegTrueP.clear();
    sfgHit->HitSegTrueEdepo.clear();

    // First get number of contributing MPPCs to this hit
    int NumContrMPPCs = hit->GetContributorCount();
    
    // Hit segments contributing to MPPC readout (3 2D planes)
    std::vector<ND::TG4HitSegment*> HitSegXYPlane;
    std::vector<ND::TG4HitSegment*> HitSegXZPlane;
    std::vector<ND::TG4HitSegment*> HitSegYZPlane;
    
    // One vector containing all of hit segments
    std::vector<ND::TG4HitSegment*> HitSegAllPlanes;
    
    // Then loop over all MPPCs
    for (int MPPCHit = 0; MPPCHit < NumContrMPPCs; MPPCHit++) {
 
      // Get the true hits contributing to the hit segments
      ND::THandle<ND::THit> currContrMPPC = hit->GetContributor(MPPCHit);
      std::vector<ND::TG4VHit*> TrueG4Hit = HitTruthInfo::GetHitTruthInfo(currContrMPPC);
      
      // Loop over all G4VHits
      for (auto G4Hit : TrueG4Hit) {
      
        // Cast G4VHit to hit segment
        ND::TG4HitSegment* G4HitSeg = dynamic_cast<ND::TG4HitSegment*>(G4Hit);
        if (G4HitSeg == NULL) {
          ND280Log("Cast of TG4Hit to TG4HitSegments unsuccessful");
          continue;
        }
 
        if (MPPCHit == 1) HitSegXZPlane.push_back(G4HitSeg);  
        if (MPPCHit == 2) HitSegYZPlane.push_back(G4HitSeg);
        if (MPPCHit == 3) HitSegXYPlane.push_back(G4HitSeg);
 
      } // End of G4VHits loop
      
    } // End of MPPC loop
    
    // The number of true tracks contributed to the hit
    std::vector<int> TrueTrkID;
    std::vector<int> TrueTrkParentID;
    std::vector<int> InCubeTag;

    std::vector<ND::TG4HitSegment*> XZandYZIntersection; // intermediate result
    std::vector<ND::TG4HitSegment*> XZandXYIntersection; // intermediate result
    std::vector<ND::TG4HitSegment*> XYandYZIntersection; // intermediate result
    
    // intersection of hit segments from all three planes 
    // i.e. the hit segments contributing to the current cube
    std::vector<ND::TG4HitSegment*> currContrHitSegV;

    // compare the three sets of hit segments contributing to the MPPCs, and 
    // find the intersection - these are the hit segments contributing to current cube
    std::sort(HitSegXZPlane.begin(), HitSegXZPlane.end());
    std::sort(HitSegYZPlane.begin(), HitSegYZPlane.end());
    std::sort(HitSegXYPlane.begin(), HitSegXYPlane.end());

    std::set_intersection(HitSegXZPlane.begin(), HitSegXZPlane.end(), HitSegYZPlane.begin(), HitSegYZPlane.end(), back_inserter(XZandYZIntersection));
    std::sort(XZandYZIntersection.begin(), XZandYZIntersection.end());
    std::set_intersection(XZandYZIntersection.begin(), XZandYZIntersection.end(), HitSegXYPlane.begin(), HitSegXYPlane.end(), back_inserter(currContrHitSegV));
    
    // loop over all hit segments contributing to the current cube (true hit + xtalk)
    // i.e. those hit segments in the intersection of all three planes
    for (auto HitSeg : currContrHitSegV) {

        // for true hits (hitSegDistanceI < 5 mm for I = X, Y, and Z), 
        // sum their energy deposit in the current cube
        if (checkHitSegInCube(hit, HitSeg) == true) {
            // get primary track ID
            ND::THandle<ND::TG4Trajectory> traj_temp = fTrueTraj->GetTrajectory(HitSeg->GetPrimaryId());
            if(!traj_temp) continue;
            TrueTrkID.push_back(traj_temp->GetTrackId());
            TrueTrkParentID.push_back(traj_temp->GetParentId());
            InCubeTag.push_back(1);
            
            // Fill true trajectory information
            FillTrueTrajInfoForHit(HitSeg,sfgHit);
        }
        else{
          ND::THandle<ND::TG4Trajectory> traj_temp = fTrueTraj->GetTrajectory(HitSeg->GetPrimaryId());
          if(!traj_temp) continue;
          TrueTrkID.push_back(traj_temp->GetTrackId());
          TrueTrkParentID.push_back(traj_temp->GetParentId());
          InCubeTag.push_back(0);
        }

    }  // end hit segment loop 
    
    // Check the length of TrueTrkID, if >= 1 then determine the single or multiple tag
    if(TrueTrkID.size()==0){
      sfgHit->TrkHitTag = 3;
    }
    else{
      // Only select tracks with ParentID = 0
      std::vector<int> PassTrkID;

      for(size_t i = 0; i < TrueTrkID.size(); i++){
        if(TrueTrkParentID[i]==0) PassTrkID.push_back(TrueTrkID[i]);
      }

      if(PassTrkID.size()>=1){
        std::sort(PassTrkID.begin(),PassTrkID.end());
        int NElement = PassTrkID.size();

        if(PassTrkID[0]!=PassTrkID[NElement-1]) sfgHit->TrkHitTag = 1;
        else{
          sfgHit->TrkHitTag = 2;
          bool trackhit = false;
          for(size_t i = 0; i < TrueTrkID.size(); i++){
            if(InCubeTag[i]==1) trackhit = true;
          }
          if(trackhit==false) sfgHit->TrkHitTag = 3;
        }
      }
      else{
        sfgHit->TrkHitTag = 2;
        bool trackhit = false;
        for(size_t i = 0; i < TrueTrkID.size(); i++){
          if(InCubeTag[i]==1) trackhit = true;
        }
        if(trackhit==false) sfgHit->TrkHitTag = 3;
      }
      
    }
    
    // Then check the nearby cubes if there are multiple tracks
    if(checkNearbyMultiTrk(hit,HitSegAllPlanes,fTrueTraj)==true) sfgHit->TrkHitTag = 1;        

}

bool ND::TSFGReconModule::checkHitSegInCube(ND::THandle<ND::THit> hit, ND::TG4HitSegment* HitSeg) {
// check if a hit segment is inside the cube corresponding to the hit

    //i.e. check if the hit segment midpoint is +/- 5 mm away from the cube position
    int dist = 5;
    int negdist = -5;
    bool isInside = false;

    // define position of hit segment as the midpoint
    double HitSegPosX = 0.5*(HitSeg->GetStartX() + HitSeg->GetStopX());
    double HitSegPosY = 0.5*(HitSeg->GetStartY() + HitSeg->GetStopY());
    double HitSegPosZ = 0.5*(HitSeg->GetStartZ() + HitSeg->GetStopZ());

    // distance of hit segment position to cube center (= position of reco hit) 
    // (if < 5 mm from cube center, it is a true hit)
    double HitSegDistanceX = hit->GetPosition().X() - HitSegPosX;
    double HitSegDistanceY = hit->GetPosition().Y() - HitSegPosY;
    double HitSegDistanceZ = hit->GetPosition().Z() - HitSegPosZ;

    bool HitSegXIn = HitSegDistanceX > negdist && HitSegDistanceX < dist;
    bool HitSegYIn = HitSegDistanceY > negdist && HitSegDistanceY < dist;
    bool HitSegZIn = HitSegDistanceZ > negdist && HitSegDistanceZ < dist;

    // for true hits (hitSegDistanceI < 5 mm for I = X, Y, and Z), 
    // sum their energy deposit in the current cube
    if (HitSegXIn && HitSegYIn && HitSegZIn) {
        isInside = true;
    }
    return isInside;
}

bool ND::TSFGReconModule::checkNearbyMultiTrk(ND::THandle<ND::THit> hit, std::vector<ND::TG4HitSegment*> HitSegAllPlanes, ND::THandle<ND::TG4TrajectoryContainer> true_traj){

  // Get the center position of the cube
  double cube_posx = hit->GetPosition().X();
  double cube_posy = hit->GetPosition().Y();
  double cube_posz = hit->GetPosition().Z();

  // A vector to save the hit segment primary ID in adjacent cubes
  std::vector<int> nearby_trkid;
  std::vector<int> nearby_trkparentid;

  bool InAdjacent;

  // Loop over all hit segments
  for (auto HitSeg : HitSegAllPlanes) {

    // Estimate hit segment position
    double hitseg_posx = 0.5*(HitSeg->GetStartX() + HitSeg->GetStopX());
    double hitseg_posy = 0.5*(HitSeg->GetStartY() + HitSeg->GetStopY());
    double hitseg_posz = 0.5*(HitSeg->GetStartZ() + HitSeg->GetStopZ());

    // Check if the hit segment is in the adjacent cube 
    // Along X direction
    if(TMath::Abs(cube_posx-hitseg_posx)>5 && TMath::Abs(cube_posx-hitseg_posx)<15 &&
       TMath::Abs(cube_posy-hitseg_posy)<5 && TMath::Abs(cube_posz-hitseg_posz)<5){
      InAdjacent = true;
    }
    // Along Y direction
    else if(TMath::Abs(cube_posy-hitseg_posy)>5 && TMath::Abs(cube_posy-hitseg_posy)<15 &&
            TMath::Abs(cube_posx-hitseg_posx)<5 && TMath::Abs(cube_posz-hitseg_posz)<5){
      InAdjacent = true;
    }
    // Along Z direction
    else if(TMath::Abs(cube_posz-hitseg_posz)>5 && TMath::Abs(cube_posz-hitseg_posz)<15 &&
            TMath::Abs(cube_posx-hitseg_posx)<5 && TMath::Abs(cube_posy-hitseg_posy)<5){
      InAdjacent = true;
    }
    else{
      InAdjacent = false;
    }

    if(InAdjacent==true){
      ND::THandle<ND::TG4Trajectory> traj_temp = true_traj->GetTrajectory(HitSeg->GetPrimaryId());
      if(!traj_temp) continue;
      nearby_trkid.push_back(traj_temp->GetTrackId());
      nearby_trkparentid.push_back(traj_temp->GetParentId());
    }

  }

  // Check the track ID
  int N = nearby_trkid.size();

  if(N<1) return false;
  else{

    // Only select tracks with ParentID = 0
    std::vector<int> pass_trkid;
    for(int i = 0; i < N; i++){
      if(nearby_trkparentid[i]==0) pass_trkid.push_back(nearby_trkid[i]);
    }

    int n = pass_trkid.size();

    if(n<1) return false;
    else{
      std::sort(pass_trkid.begin(),pass_trkid.end());

      if(pass_trkid[0]!=pass_trkid[n-1]) return true;
      else return false;
    }

  }

}

void ND::TSFGReconModule::FillTrueTrajInfoForHit(ND::TG4HitSegment *HitSeg, TSFGHit *sfgHit){

  // Compute hit segment center position
  TLorentzVector hitseg_pos;
  hitseg_pos.SetXYZT(0.5*(HitSeg->GetStartX() + HitSeg->GetStopX()),
                     0.5*(HitSeg->GetStartY() + HitSeg->GetStopY()),
                     0.5*(HitSeg->GetStartZ() + HitSeg->GetStopZ()),
                     0.5*(HitSeg->GetStartT() + HitSeg->GetStopT()));

  // Compute hit segment direction
  TVector3 hitseg_dir;
  hitseg_dir.SetXYZ(HitSeg->GetStopX() - HitSeg->GetStartX(),
                    HitSeg->GetStopY() - HitSeg->GetStartY(),
                    HitSeg->GetStopZ() - HitSeg->GetStartZ());

  if(hitseg_dir.Mag() > 0) hitseg_dir *= (1. / hitseg_dir.Mag());

  // Get the IDs of the true trajectory contributing to this hit segment
  std::vector<int> traj_ids = HitSeg->GetContributors();
  
  // Loop over true trajectory IDs to save their information
  int curr_id = traj_ids[0];
  for(size_t i = 0; i < traj_ids.size(); i++){
  
    if(i != 0 && traj_ids[i] == curr_id) continue;
    curr_id = traj_ids[i];
  
    ND::THandle<ND::TG4Trajectory> traj_temp = fTrueTraj->GetTrajectory(traj_ids[i]);
    if(!traj_temp) continue;
  
    // Fill the position (same for the same hit segment)
    sfgHit->HitSegPosition.push_back(hitseg_pos);
    
    // Fill the direction
    sfgHit->HitSegDirection.push_back(hitseg_dir);
    
    // True particle PDG of the trajectory
    sfgHit->HitSegTruePDG.push_back(traj_temp->GetPDGEncoding());
    
    // True momentum of the trajectory
    int traj_mom = (traj_temp->GetInitialMomentum().Vect()).Mag();
    sfgHit->HitSegTrueP.push_back(traj_mom);
  
    // True energy loss from this hit segment
    sfgHit->HitSegTrueEdepo.push_back(HitSeg->GetSecondaryDeposit());
  
  }

}

template <class T>
ND::THandle<ND::TG4Trajectory> ND::TSFGReconModule::GetG4Trajectory(
    const T& object, double& purity, double& eff){

    purity = -0xABCDEF;
    eff = -0xABCDEF;
    
    ND::TG4Trajectory traj;
    ND::THandle<ND::THitSelection> hits = object.GetHits();
    if (!hits) {
        return ND::THandle<ND::TG4Trajectory>();
    }
    if (hits->size() < 1) {
        return ND::THandle<ND::TG4Trajectory>();
    }

    int G4TrkId = TrackTruthInfo::GetG4TrajIDHits(*hits, eff, purity);
    const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
    ND::THandle<ND::TG4TrajectoryContainer> G4trajectories =
    event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
    if (!G4trajectories) {
        return ND::THandle<ND::TG4Trajectory>();
    }
    ND::TG4TrajectoryContainer::iterator trajIter = fTrueTraj->begin();
    ND::TG4TrajectoryContainer::iterator trajEnd = fTrueTraj->end();
    for (; trajIter != trajEnd; ++trajIter) {
        ND::TG4Trajectory* trajectory = &(trajIter->second);
        if (trajectory->GetTrackId() == G4TrkId) {
            ND::THandle<ND::TG4Trajectory> traj(new ND::TG4Trajectory((*trajectory)));
            return traj;
        }
    }
    
    return ND::THandle<ND::TG4Trajectory>();

}