TGlobalReconModule.cxx

Go to the documentation of this file.

#include <algorithm>
#include <functional>
#include <vector>

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

#include "TDatum.hxx"
#include "TRealDatum.hxx"
#include "TReconNode.hxx"
#include "TReconShower.hxx"
#include "TReconTrack.hxx"

#include "TGeomInfo.hxx"

#include "TComboHit.hxx"
#include "TG4HitSegment.hxx"
#include "TReconHit.hxx"
#include "TTPCmcDigit.hxx"

#include "TEventFolder.hxx"
#include "TRecPackManager.hxx"
#include "TValidationUtils.hxx"

#include "ReconObjectUtils.hxx"
#include "ReconPrintout.hxx"
#include "TIntegerDatum.hxx"
#include "TrackTruthInfo.hxx"
#include "TrackingUtils.hxx"
#include "fgdUtils.hxx"
#include "pidUtils.hxx"

#include "TReconStateUtilities.hxx"
#include "TTrackerReconUtils.hxx"

#include "TTrackerECALReconModule.hxx"

#include "TGlobalReconModule.hxx"

const int NMAXALTERNATES = 6;
const int NMAXECAL       = 5;
const int NMAXP0DECAL    = 2;
const int NMAXTPC        = 6;
const int NMAXFGD        = 2;
const int NMAXSFG        = 2;
const int NMAXP0D        = 1;
const int NMAXSMRD       = 4;
const int NMAXTRACKER    = 10;
const int NMAXTPCOTHER   = 20;

const bool debug_kine = false;
const bool debug_pids = false;
const bool debug_broken = false;
const bool debug_extrap = false;
const bool debug_vertex_extrap = false;

const int DEFAULT_MIN = 99999;
const int DEFAULT_MAX = -99999;
using std::cout;
using std::endl;
ClassImp(ND::TGlobalReconModule::TTpcPID);
ND::TGlobalReconModule::TTpcPID::~TTpcPID() {}

ClassImp(ND::TGlobalReconModule::TGlobalPID);
ND::TGlobalReconModule::TGlobalPID::~TGlobalPID() { /*if (TrueParticle) delete
    TrueParticle;*/
}

ND::TGlobalReconModule::TGlobalPID::TGlobalPID() {
    Alternates = new TClonesArray("ND::TGlobalReconModule::TGlobalPIDAlternate",
        NMAXALTERNATES);
    ECAL = new TClonesArray("ND::TGlobalReconModule::TECALObject", NMAXECAL);
    P0DECAL =
    new TClonesArray("ND::TGlobalReconModule::TECALObject", NMAXP0DECAL);
    TPC = new TClonesArray("ND::TGlobalReconModule::TTPCObject", NMAXTPC);
    FGD = new TClonesArray("ND::TGlobalReconModule::TFGDObject", NMAXFGD);
    SFG = new TClonesArray("ND::TGlobalReconModule::TSFGObject", NMAXSFG);
    P0D = new TClonesArray("ND::TGlobalReconModule::TP0DObject", NMAXP0D);
    SMRD = new TClonesArray("ND::TGlobalReconModule::TSMRDObject", NMAXSMRD);
    TRACKER =
    new TClonesArray("ND::TGlobalReconModule::TTrackerObject", NMAXTRACKER);
    HitsSaved = new TClonesArray("ND::TGlobalReconModule::TGlobalHit", 4);
    
    // TPC digit decoding is a little too verbose.
    // Disable it here.
    ND::TND280Log::SetLogLevel("eventUnpack", ND::TND280Log::QuietLevel);
    
    NAlternates = 0;
    Alternates->Clear();
    
    NDsECALs = 0;
    NTrECALs = 0;
    NECALs = 0;
    ECAL->Clear();
    
    NP0DECALs = 0;
    P0DECAL->Clear();
    
    NTPCs = 0;
    TPC->Clear();
    
    NFGDs = 0;
    FGD->Clear();
    
    NSFGs = 0;
    SFG->Clear();
    
    NSMRDs = 0;
    SMRD->Clear();
    
    NP0Ds = 0;
    P0D->Clear();
    
    NTRACKERs = 0;
    TRACKER->Clear();
    
    NHitsSaved = 0;
    HitsSaved->Clear();
}

ClassImp(ND::TGlobalReconModule::TVertexConstituent);
ND::TGlobalReconModule::TVertexConstituent::~TVertexConstituent() {}

ClassImp(ND::TGlobalReconModule::TGlobalVertex);
ND::TGlobalReconModule::TGlobalVertex::~TGlobalVertex() {}
ND::TGlobalReconModule::TGlobalVertex::TGlobalVertex() {
    Constituents = new TClonesArray("ND::TGlobalReconModule::TVertexConstituent",
        NCONSTITUENTS);
    TrueVertices = new TClonesArray("ND::TTrueVertex", NCONSTITUENTS);
    NConstituents = 0;
    Constituents->Clear();
    NTrueVertices = 0;
    TrueVertices->Clear();
}

ClassImp(ND::TGlobalReconModule::TP0DObject);
ND::TGlobalReconModule::TP0DObject::~TP0DObject() {}

ClassImp(ND::TGlobalReconModule::TECALObject);
ND::TGlobalReconModule::TECALObject::~TECALObject() {}

ClassImp(ND::TGlobalReconModule::TTPCObject);
ND::TGlobalReconModule::TTPCObject::~TTPCObject() {}

ClassImp(ND::TGlobalReconModule::TTPCOtherObject);
ND::TGlobalReconModule::TTPCOtherObject::~TTPCOtherObject() {}

ClassImp(ND::TGlobalReconModule::TFGDObject);
ND::TGlobalReconModule::TFGDObject::~TFGDObject() {}

ClassImp(ND::TGlobalReconModule::TSFGObject);
ND::TGlobalReconModule::TSFGObject::~TSFGObject() {}

ClassImp(ND::TGlobalReconModule::TSMRDObject);
ND::TGlobalReconModule::TSMRDObject::~TSMRDObject() {}

ClassImp(ND::TGlobalReconModule::TTrackerObject);
ND::TGlobalReconModule::TTrackerObject::~TTrackerObject() {}

ClassImp(ND::TGlobalReconModule::TFgdTimeBin);
ND::TGlobalReconModule::TFgdTimeBin::~TFgdTimeBin() {}

ClassImp(ND::TGlobalReconModule::TSfgTimeBin);
ND::TGlobalReconModule::TSfgTimeBin::~TSfgTimeBin() {}

ND::TGlobalReconModule::TFgdTimeBin::TFgdTimeBin() {
    FGD1Unused = new TClonesArray("ND::TGlobalReconModule::TGlobalHit", 20);
    FGD2Unused = new TClonesArray("ND::TGlobalReconModule::TGlobalHit", 20);
    
    NFGD1Unused = 0;
    FGD1Unused->Clear();
    
    NFGD2Unused = 0;
    FGD2Unused->Clear();
}

ND::TGlobalReconModule::TSfgTimeBin::TSfgTimeBin() {
    SFGUnused = new TClonesArray("ND::TGlobalReconModule::TGlobalHit", 20);
    NSFGUnused = 0;
    SFGUnused->Clear();
}

ClassImp(ND::TGlobalReconModule::TGlobalHit);
ND::TGlobalReconModule::TGlobalHit::~TGlobalHit() {}

ClassImp(ND::TGlobalReconModule::TSMRDHit);
ND::TGlobalReconModule::TSMRDHit::~TSMRDHit() {}

ClassImp(ND::TGlobalReconModule::TOutermostHits);
ND::TGlobalReconModule::TOutermostHits::~TOutermostHits() {}

bool SortHitsInX(const ND::THandle<ND::THit>& h1,
    const ND::THandle<ND::THit>& h2) 
{
    return h1->GetPosition().X() < h2->GetPosition().X();
}
bool SortHitsInY(const ND::THandle<ND::THit>& h1,
    const ND::THandle<ND::THit>& h2) 
{
    return h1->GetPosition().Y() < h2->GetPosition().Y();
}
bool SortHitsInZ(const ND::THandle<ND::THit>& h1,
    const ND::THandle<ND::THit>& h2) 
{
    return h1->GetPosition().Z() < h2->GetPosition().Z();
}

bool SortNodesInX(const ND::THandle<ND::TReconNode>& n1,
    const ND::THandle<ND::TReconNode>& n2) 
{
    return TrackingUtils::GetPosition(n1->GetState()).X() <
    TrackingUtils::GetPosition(n2->GetState()).X();
}
bool SortNodesInY(const ND::THandle<ND::TReconNode>& n1,
    const ND::THandle<ND::TReconNode>& n2) 
{
    return TrackingUtils::GetPosition(n1->GetState()).Y() <
    TrackingUtils::GetPosition(n2->GetState()).Y();
}
bool SortNodesInZ(const ND::THandle<ND::TReconNode>& n1,
    const ND::THandle<ND::TReconNode>& n2)
{
    return TrackingUtils::GetPosition(n1->GetState()).Z() <
    TrackingUtils::GetPosition(n2->GetState()).Z();
}
bool SortNodesInXReverse(const ND::THandle<ND::TReconNode>& n1,
    const ND::THandle<ND::TReconNode>& n2) 

{
    return TrackingUtils::GetPosition(n1->GetState()).X() >=
    TrackingUtils::GetPosition(n2->GetState()).X();
}
bool SortNodesInYReverse(const ND::THandle<ND::TReconNode>& n1,
    const ND::THandle<ND::TReconNode>& n2) 
{
    return TrackingUtils::GetPosition(n1->GetState()).Y() >=
    TrackingUtils::GetPosition(n2->GetState()).Y();
}
bool SortNodesInZReverse(const ND::THandle<ND::TReconNode>& n1,
    const ND::THandle<ND::TReconNode>& n2) 
{
    return TrackingUtils::GetPosition(n1->GetState()).Z() >=
    TrackingUtils::GetPosition(n2->GetState()).Z();
}
#include "../cvstags/TGlobalReconModule.cxx"


//*************************************************************
ND::TGlobalReconModule::TGlobalReconModule(const char* name,
    const char* title) 
{
    //*************************************************************
    
    fTestTPCInfo = false;
    fTestTPCInfo = true;
    
    SetNameTitle(name, title);
    
    // Enable this module by default:
    fIsEnabled = kTRUE;
    fDescription = "Global Recon Output";
    fCVSTagName = CVSTAG;
    fCVSID = CVSID;
    
    // Tree variables
    fNPIDs = 0;
    fPIDs = new TClonesArray("ND::TGlobalReconModule::TGlobalPID", 20);
    
    fNTPCOthers = 0;
    fTPCOthers =
    new TClonesArray("ND::TGlobalReconModule::TTPCOtherObject", NMAXTPCOTHER);
    
    fNVertices = 0;
    fVertices = new TClonesArray("ND::TGlobalReconModule::TGlobalVertex", 20);
    fPVInd = 0;
    fNSVertices = 0;
    
    fNFgdTimeBins = 0;
    fFgdTimeBins = new TClonesArray("ND::TGlobalReconModule::TFgdTimeBin", 20);
    
    fNSfgTimeBins = 0;
    fSfgTimeBins = new TClonesArray("ND::TGlobalReconModule::TSfgTimeBin", 20);
    
    fSMRDUnused = new TClonesArray("ND::TGlobalReconModule::TSMRDHit", 100);
    
    fP0DUnused = new TClonesArray("ND::TGlobalReconModule::TGlobalHit", 100);
    
    fNTPCPIDs = 0;
    if (fTestTPCInfo) {
        fTPCPIDs = new TClonesArray("ND::TGlobalReconModule::TTpcPID", 20);
    } else {
        fTPCPIDs = NULL;
    }
    
    fUseECalVertices = false;
    fRecPackInitialized = false;
}

void ND::TGlobalReconModule::FillConfigTree(TTree* configTree)
{
    if (fFilledConfigTree) {
        ND280NamedWarn(
            "TGlobalReconModule",
            "Module: " << this->GetName()
            << " has already written the config tree this run.");
    }
    configTree->Branch("AllModules", &fALLMODULES);
    configTree->Branch("UseECalVertices", &fUseECalVertices);
    
    TAnalysisModuleBase::FillConfigTree(configTree);
}

//*************************************************************
ND::TGlobalReconModule::~TGlobalReconModule() {}
//*************************************************************

Bool_t ND::TGlobalReconModule::Configure(std::string& opt) 
{
    if (opt == "UseECalVertices") {
        fUseECalVertices = true;
        return true;
    }
    return false;
}

//*************************************************************
Bool_t ND::TGlobalReconModule::ProcessFirstEvent(ND::TND280Event& event) 
{
    //*************************************************************
    return true;
}

//*************************************************************
void ND::TGlobalReconModule::InitializeModule() {}
//*************************************************************

//*************************************************************
void ND::TGlobalReconModule::InitializeBranches() 
{
    //*************************************************************
    SetSplitLevel(1);
    
    if (fTestTPCInfo) 
    {
        fOutputTree->Branch("NTPCPIDs", &fNTPCPIDs, "NTPCPIDs/I", fBufferSize);
        
        fOutputTree->Branch("TPCPIDs", "TClonesArray", &fTPCPIDs, fBufferSize,
            fSplitLevel);
    }
    
    // The number of vertex objects that have been saved.
    fOutputTree->Branch("NVertices", &fNVertices, "NVertices/I", fBufferSize)
    ->SetTitle("The number of added primary vertices. ");
    
    // Vertex branch
    fOutputTree->Branch("Vertices", &fVertices, fBufferSize, fSplitLevel)
    ->SetTitle("The TGlobalVertex vector of vertices. ");
    
    // FGD hit time bins
    fOutputTree->Branch("NFgdTimeBins", &fNFgdTimeBins, "NFgdTimeBins/I",
        fBufferSize)
    ->SetTitle("  Number of hit time bins in the FGD as determined by fgdRecon. ");
    fOutputTree->Branch("FgdTimeBins", &fFgdTimeBins, fBufferSize, fSplitLevel)
    ->SetTitle(" Information for each hit time bin (TFgdTimeBins). ");
    // SFG hit time bins
    fOutputTree->Branch("NSfgTimeBins", &fNSfgTimeBins, "NSfgTimeBins/I",
        fBufferSize)
    ->SetTitle("  Number of hit time bins in the SFG as determined by sfgRecon. ");
    fOutputTree->Branch("SfgTimeBins", &fSfgTimeBins, fBufferSize, fSplitLevel)
    ->SetTitle(" Information for each hit time bin (TSfgTimeBins). ");
    // The number of objects that have been saved.
    fOutputTree->Branch("NPIDs", &fNPIDs, "NPIDs/I", fBufferSize)
    ->SetTitle(" The number of global objects.");
    // A branch for the clone array of Global object information.
    fOutputTree->Branch("PIDs", "TClonesArray", &fPIDs, fBufferSize, fSplitLevel)
    ->SetTitle(" The vector of TGlobalPID.");
    // The number of objects that have been saved.
    fOutputTree->Branch("NTPCOthers", &fNTPCOthers, "NTPCOthers/I", fBufferSize)
    ->SetTitle("The number of TPC other objects. ");
    // A branch for the clone array of Global object information.
    fOutputTree->Branch("TPCOthers", "TClonesArray", &fTPCOthers, fBufferSize,
        fSplitLevel)
    ->SetTitle("The vector of TPC other objects (TTPCOtherObject). ");
    // number of unused hits in FGD1, FGD2, P0D and SMRD
    fOutputTree->Branch("NFGD1Unused", &fNFGD1Unused, "NFGD1Unused/I", fBufferSize)
    ->SetTitle(" Number of hits unused in FGD1");
    fOutputTree->Branch("NFGD2Unused", &fNFGD2Unused, "NFGD2Unused/I", fBufferSize)
    ->SetTitle(" Number of hits unused in FGD2");
    fOutputTree->Branch("NSFGUnused", &fNSFGUnused, "NSFGUnused/I", fBufferSize)
    ->SetTitle(" Number of hits unused in SFG");
    fOutputTree->Branch("NP0DUnused", &fNP0DUnused, "NP0DUnused/I", fBufferSize)
    ->SetTitle("  Number of hits unused in the P0D");
    fOutputTree->Branch("P0DUnused", "TClonesArray", &fP0DUnused, fBufferSize,
        fSplitLevel)
    ->SetTitle("The vector unused hits in the P0D (TGlobalHit). ");
    fOutputTree->Branch("SMRDUnused", "TClonesArray", &fSMRDUnused, fBufferSize,
        fSplitLevel)
    ->SetTitle(" The vector unused hits in the SMRD (TSMRDHit). ");
    fOutputTree->Branch("NSMRDUnused", &fNSMRDUnused, "NSMRDU nused/I",
        fBufferSize)
    ->SetTitle("  Number of all hits unused in the SMRD"); 
    fOutputTree->Branch("NSMRDTopUnused", &fNSMRDTopUnused, "NSMRDTopUnused/I",
        fBufferSize)
    ->SetTitle(" Number of hits unused in the SMRDTop ");
    fOutputTree->Branch("NSMRDBottomUnused", &fNSMRDBottomUnused,
        "NSMRDBottomUnused/I", fBufferSize)
    ->SetTitle(" Number of hits unused in the SMRDBottom.");
    fOutputTree->Branch("NSMRDLeftUnused", &fNSMRDLeftUnused, "NSMRDLeftUnused/I",
        fBufferSize)
    ->SetTitle("Number of hits unused in the SMRDLeft. ");
    fOutputTree->Branch("NSMRDRightUnused", &fNSMRDRightUnused,
        "NSMRDRightUnused/I", fBufferSize)
    ->SetTitle("  Number of hits unused in the SMRDRight.");
    
    fOutputTree->Branch("NTPCUnused", &fNTPCUnused, "NTPCUnused/I", fBufferSize)
    ->SetTitle(" Number of hits unused in the TPC. ");
    
    // Save P0D outermost hits
    fOutputTree->Branch("P0DOutermostHits", "TOutermostHits", &fP0DOutermostHits,
        fBufferSize, fSplitLevel)
    ->SetTitle("Outermost hits in the P0D. ");
    // time of the median hit of the earliest track
    fOutputTree->Branch("EarliestTrackTime", &fEarliestTrackMedianHitTime,
        fBufferSize)
    ->SetTitle("  Median hit time of the earliest track. ");
    // Counters for each detector
    fOutputTree->Branch("NDsECAL", &fNDsECAL, "NDsECAL/I", fBufferSize)
    ->SetTitle("Number of objects containing the DsECAL. "); 
    fOutputTree->Branch("NTrECAL", &fNTrECAL, "NTrECAL/I", fBufferSize)
    ->SetTitle("  Number of objects containing the TrECAL."); 
    fOutputTree->Branch("NP0DECAL", &fNP0DECAL, "NP0DECAL/I", fBufferSize)
    ->SetTitle(" Number of objects containing the P0DECAL. "); 
    fOutputTree->Branch("NTPC", &fNTPC, "NTPC/I", fBufferSize)
    ->SetTitle("  Number of objects containing the TPC.");
    fOutputTree->Branch("NFGD", &fNFGD, "NFGD/I", fBufferSize)
    ->SetTitle("  Number of objects containing the FGD.");
    fOutputTree->Branch("NSFG", &fNSFG, "NSFG/I", fBufferSize)
    ->SetTitle("  Number of objects containing the SFG.");
    fOutputTree->Branch("NP0D", &fNP0D, "NP0D/I", fBufferSize)
    ->SetTitle(" Number of objects containing the P0D.");
    fOutputTree->Branch("NSMRD", &fNSMRD, "NSMRD/I", fBufferSize)
    ->SetTitle("   Number of objects containing the SMRD.");
    fOutputTree->Branch("NDsECALIso", &fNDsECALIso, "NDsECALIso/I", fBufferSize)
    ->SetTitle("Number of objects in the DsECAL only. "); 
    fOutputTree->Branch("NTrECALIso", &fNTrECALIso, "NTrECALIso/I", fBufferSize)
    ->SetTitle(" Number of objects in the TrECAL only. ");
    fOutputTree->Branch("NP0DECALIso", &fNP0DECALIso, "NP0DECALIso/I",
        fBufferSize)
    ->SetTitle("  Number of objects in the P0DECAL only.");
    fOutputTree->Branch("NTPCIso", &fNTPCIso, "NTPCIso/I", fBufferSize)
    ->SetTitle(" Number of objects in the TPC only. ");
    fOutputTree->Branch("NFGDIso", &fNFGDIso, "NFGDIso/I", fBufferSize)
    ->SetTitle(" Number of objects in the FGD only. ");
    fOutputTree->Branch("NSFGIso", &fNSFGIso, "NSFGIso/I", fBufferSize)
    ->SetTitle(" Number of objects in the SFG only. ");
    fOutputTree->Branch("NP0DIso", &fNP0DIso, "NP0DIso/I", fBufferSize)
    ->SetTitle(" Number of objects in the P0D only.");
    fOutputTree->Branch("NSMRDIso", &fNSMRDIso, "NSMRDIso/I", fBufferSize)
    ->SetTitle("Number on objects in the SMRD only. ");                 
}

//*************************************************************
bool ND::TGlobalReconModule::FillTree(ND::TND280Event& event)
{
    //*************************************************************
    
    if (!fRecPackInitialized){
        // add a new RecpackManager
        ND::tman().add_manager("globalAnalysis");
        ND::tman().select_manager("globalAnalysis");
        
        // only surfaces with measurements are tried inside a volume
        ND::rpman("globalAnalysis")
        .navigation_svc()
        .navigator(RP::particle_helix)
        .set_unique_surface(true);
        
        // enable energy loss correction
        ND::rpman("globalAnalysis")
        .model_svc()
        .enable_correction(RP::particle_helix, RP::eloss, true);
        
        // enable MS
        ND::rpman("globalAnalysis")
        .model_svc()
        .enable_noiser(RP::particle_helix, RP::ms, true);
        
        fEquation = ND::rpman()
        .model_svc()
        .model(RP::particle_helix)
        .retrieve_tool<HelixEquation>("equation");
        
        if (debug_extrap) {
            ND::tman("globalAnalysis").SetVerbosity(0, 2, 0, 0);
            ND::rpman("globalAnalysis")
            .model_svc()
            .model(RP::particle_helix)
            .tool("correction/eloss")
            .set_verbosity(Messenger::VVERBOSE);
            ND::rpman("globalAnalysis")
            .navigation_svc()
            .navigator(RP::particle_helix)
            .master_inspector()
            .set_verbosity(Messenger::VVERBOSE);
            ND::rpman("globalAnalysis")
            .navigation_svc()
            .inspector("eloss")
            .set_verbosity(Messenger::VVERBOSE);
        }
        InitializeExtrapolationToDetectors();
        
        fRecPackInitialized = true;
    }
    
    if (fTestTPCInfo) 
    {
        fNTPCPIDs = 0;
        fTPCPIDs->Clear();
    }
    
    fNVertices = 0;
    fVertices->Clear();
    fNPIDs = 0;
    fNSVertices = 0;
    
    fNPIDs = 0;
    fPIDs->Clear();
    
    fNTPCOthers = 0;
    fTPCOthers->Clear();
    
    fNFgdTimeBins = 0;
    fFgdTimeBins->Clear();
    
    fNSfgTimeBins = 0;
    fSfgTimeBins->Clear();
    
    fNFGD1Unused = 0;
    fNFGD2Unused = 0;
    
    fNSFGUnused = 0;
    
    fNP0DUnused = 0;
    fP0DUnused->Clear();
    
    fSMRDUnused->Clear();
    fNSMRDUnused = 0;
    fNSMRDTopUnused = 0;
    fNSMRDBottomUnused = 0;
    fNSMRDLeftUnused = 0;
    fNSMRDRightUnused = 0;
    
    fNTPCUnused = 0;
    
    fEarliestTrackMedianHitTime = DEFAULT_MIN;
    
    fNDsECAL = 0;
    fNTrECAL = 0;
    fNP0DECAL = 0;
    fNTPC = 0;
    fNFGD = 0;
    fNSFG = 0;
    fNP0D = 0;
    fNSMRD = 0;
    
    fNDsECALIso = 0;
    fNTrECALIso = 0;
    fNP0DECALIso = 0;
    fNTPCIso = 0;
    fNFGDIso = 0;
    fNSFGIso = 0;
    fNP0DIso = 0;
    fNSMRDIso = 0;
    
    // Get the main result to save.
    ND::THandle<ND::TAlgorithmResult> globalRecon = event.GetFit("oaRecon");
    if (!globalRecon) 
    {
        return true;
    }
    
    // Get the results from reconstruction
    ND::THandle<ND::TReconObjectContainer> allObjects =
    globalRecon->GetResultsContainer("final");
    
    if (!allObjects) return true;
    
    // This is only for tests
    if (fTestTPCInfo) 
    {
        ND::THandle<ND::TAlgorithmResult> tpcRecon = event.GetFit("tpcRecon");
        if (tpcRecon)
        {
            // Get the results from reconstruction
            ND::THandle<ND::TReconObjectContainer> tpcObjects =
            tpcRecon->GetResultsContainer("TPCPids");
            
            if (tpcObjects)
            {
                // Now loop over objects.
                for (ND::TReconObjectContainer::iterator tt = tpcObjects->begin();
                    tt != tpcObjects->end(); tt++)
                {
                    ND::THandle<ND::TReconBase> object = *tt;
                    if (object) FillTPCPID(object);
                }
            }
        }
    }
    
    // fill the map of broken tracks (to be used later)
    FillBrokenTracksMap(*allObjects);
    
    // Now loop over object and fill the global analysis objects
    fGlobalIndexMap.clear();
    for (ND::TReconObjectContainer::iterator tt = allObjects->begin();
         tt != allObjects->end(); tt++)
    {
        ND::THandle<ND::TReconBase> object = *tt;
        if (!object) continue;
        
        // The way to be done in future
        //   ND::THandle<ND::TReconVertex> vertex = *tt;
        //     if (vertex)
        //         FillVertex(object);
        //     else//Fill track/pid informtaion
        
        FillGlobalPIDs(event, object);
    }
    
    FillVertexInfo(event, allObjects);
    
    FillSfgTimeBins(event);
    
    FillFgdTimeBins(event);
    
    FillUnusedHits(event);
    
    FillTPCOther(event);
    
    FillSFGOther(event);
    
    // Fill P0D outermost hits
    double p0d_charge_cut = 7;
    ND::THandle<ND::THitSelection> p0dHits = event.GetHitSelection("p0d");
    if (p0dHits) FillOutermostHits(*p0dHits, p0d_charge_cut, fP0DOutermostHits);
    
    fGlobalIndexMap.clear();
    fBrokenIndexMap.clear();
    
    return true;
}


//*************************************************************
void ND::TGlobalReconModule::FillGlobalPIDs(
    ND::TND280Event& event, ND::THandle<ND::TReconBase> object)
{
//*************************************************************
    
  if (debug_pids) std::cout << "Fill Global PID: " << object << std::endl;
  
  // Get the compact detector number and number of used detectors
  unsigned long dets = GetDetectorNumber(object);
  int      det_count = CountDetector(object);
  
  // update counters for specfic detectors
  if (object->UsesDetector(ND::TReconBase::kDSECal)) 
  {
    fNDsECAL++;
    if (det_count == 1) fNDsECALIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kECal)) 
  {
    fNTrECAL++;
    if (det_count == 1) fNTrECALIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kPECal))
  {
    fNP0DECAL++;
    if (det_count == 1) fNP0DECALIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kFGD)) 
  {
    fNFGD++;
    if (det_count == 1) fNFGDIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kSFG)) 
  {
    fNSFG++;
    if (det_count == 1) fNSFGIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kTPC))
  {
    fNTPC++;
    if (det_count == 1) fNTPCIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kP0D))
  {
    fNP0D++;
    if (det_count == 1) fNP0DIso++;
  }
  if (object->UsesDetector(ND::TReconBase::kSMRD))
  {
    fNSMRD++;
    if (det_count == 1) fNSMRDIso++;
  }
  
  // __________________________________________________________
  // Fill the reconstructed Global PIDs.
  
  // try to get the state of the object. If it dosen't exists skip the object.
  try 
  {
    ND::THandle<ND::TReconState> tstate = object->GetState();
  } catch (ND::EReconObject const& e)
  {
    std::cout << "TReconBase with no state. Skip object !!" << std::endl;
    return;
  }
  
  // Fill the map of indexes used for vertex link
  fGlobalIndexMap[object] = fNPIDs;
  
  // create a new TGlobalPID
  TGlobalPID* globalObject = new ((*fPIDs)[fNPIDs++]) TGlobalPID;
  
  // the unique ID of the global object
  globalObject->UniqueID = object->GetUniqueID();
  
  //--------- Fill the vector of unique IDs of associated tracks
  if (fBrokenIndexMap.find(object) != fBrokenIndexMap.end())
    globalObject->BrokenUniqueIDs = fBrokenIndexMap[object];
  
  //--------- Information about the detectors that used in this object
  globalObject->Detectors = dets;
  FillDetectorUsed(object, globalObject->DetectorUsed);
  
  //------- General information about the track --------
  
  globalObject->isForward = (TrackingUtils::GetDirection(object).Z() > 0);
  globalObject->SenseOK = TrackingUtils::GetSenseOK(*object);
  globalObject->IsCurvBack = TrackingUtils::GetCurvBackTrack(*object);
  globalObject->Charge = TrackingUtils::GetCharge(object);
  globalObject->EDeposit = TrackingUtils::GetEDeposit(object);
  globalObject->Width =
  TrackingUtils::GetWidth(object).X();  //! take only the first component
  globalObject->Cone = TrackingUtils::GetCone(object);
  globalObject->Length = ComputeTrackLength(object);
  
  globalObject->AlgorithmName = object->GetAlgorithmName();
  globalObject->Status = object->CheckStatus(object->kSuccess);
  globalObject->Chi2 = object->GetQuality();
  globalObject->NDOF = object->GetNDOF();
  try 
  {
    globalObject->NNodes = object->GetNodes().size();
  } catch (ND::EReconObject const& e)
  {
    globalObject->NNodes = 0;
  }
  
  globalObject->NHits = DEFAULT_MAX;
  if (object->GetHits()) globalObject->NHits = object->GetHits()->size();
  
  // The number of constituents. For composite objects
  globalObject->NConstituents = 0;
  globalObject->IsFgdECalIncremental = false;
  if (object->GetConstituents()) 
  {
    globalObject->NConstituents = object->GetConstituents()->size();
    
    ND::TReconObjectContainer t1const =
    ReconObjectUtils::GetConstituentsRecursive(object, true);
    
    if (t1const.size() > 0) 
    {
      ND::TReconObjectContainer::const_iterator it1;
      
      for (it1 = t1const.begin(); it1 != t1const.end(); it1++)
      {
        ND::THandle<ND::TReconBase> obj = *it1;
        if (obj) 
        {
          if (obj->GetAlgorithmName() == "FgdECalIncremental") 
          {
            globalObject->IsFgdECalIncremental = true;
            break;
          }
        }
      }
    }
  }
  
  //---------- Fill info about first and last hits
  if (object->GetHits()) 
  {
//#define FILL_ALL_GLOBAL_HITS
#ifdef FILL_ALL_GLOBAL_HITS
    FillAllHits(*(object->GetHits()), *globalObject);
#else
    FillFirstLastHits(*(object->GetHits()), *globalObject);
#endif
  }
  
  // ------------ Fill the node position info
  globalObject->NodePosition.clear();
  const ND::TReconNodeContainer& nodes = object->GetNodes();
  ND::TReconNodeContainer::const_iterator it;
  for(it = nodes.begin(); it != nodes.end(); it++){
    ND::THandle<ND::TReconState> nodeState = (*it)->GetState();
    if(!nodeState) continue;
    
    TLorentzVector node_pos = ND::TReconStateUtilities::GetPosition(nodeState);
    globalObject->NodePosition.push_back(node_pos);
  }
  
  // ------------ Fill the hit position info
  globalObject->HitPosition.clear();
  const ND::THandle<ND::THitSelection> hits = object->GetHits();
  ND::THitSelection::const_iterator it_h;
  for(it_h = hits->begin(); it_h != hits->end(); it_h++){
    ND::THandle<ND::THit> hit = *it_h;
    if(!hit) continue;
    
    TVector3 hit_pos = hit->GetPosition();
    TLorentzVector pos_temp(hit_pos.X(),hit_pos.Y(),hit_pos.Z(),
                            hit->GetTime());
    globalObject->HitPosition.push_back(pos_temp);
  }  
  
  //---------- Fill time info: for each group of nodes in a sub-detector we
  // store the
  // times of two edge nodes
  
  // do not apply for TPC iso
  
  if (object->GetDetectors() != ND::TReconBase::kTPC1 &&
      object->GetDetectors() != ND::TReconBase::kTPC2 &&
      object->GetDetectors() != ND::TReconBase::kTPC3)
  {
    const ND::TReconNodeContainer& tnodes = object->GetNodes();
    unsigned long det_prev = 0;
    double first_time = DEFAULT_MIN;
    double last_time = DEFAULT_MAX;
    bool   first = true;
    
    ND::TReconNodeContainer::const_iterator in;
    for (in = tnodes.begin(); in != tnodes.end(); in++)
    {
      ND::THandle<ND::TReconBase> nodeObject = (*in)->GetObject();
      if (!nodeObject) continue;
      
      ND::THandle<ND::TReconState> nodeState = (*in)->GetState();
      if (!nodeState) continue;
      
      unsigned long det_current = GetDetectorNumber(nodeObject);
      
      double nodeTime = TrackingUtils::GetPosition(nodeState).T();
      
      if (det_current != det_prev){ 
        if (!first) {
          globalObject->NodeTimes.push_back(
            std::make_pair(det_prev, TVector2(first_time, last_time)));
        }
        first = false;
        
        det_prev = det_current;
        first_time = nodeTime;
        last_time = nodeTime;
      } else{
        last_time = nodeTime;
      }
    }
    
    // last group of nodes
    if (det_prev != 0){
      globalObject->NodeTimes.push_back(
        std::make_pair(det_prev, TVector2(first_time, last_time)));
    }
  }
  
  //--------- Fill the true particle when it is found
  double pur, eff;
  ND::THandle<ND::TG4Trajectory> G4track = GetG4Trajectory(*object, pur, eff);
  FillTrueParticle(G4track, pur, eff, globalObject->TrueParticle);
  
  //--------- Fill the kinematic extrapolated to the true vertex
  FillKinematicsAtTrueVertex(
    G4track, object, globalObject->PositionAtTrueVertex,
    globalObject->DirectionAtTrueVertex, globalObject->MomentumAtTrueVertex,
    globalObject->MomentumErrorAtTrueVertex, globalObject->ChargeAtTrueVertex,
    globalObject->StateCovAtTrueVertex);
  
  //------- Front and Back kinematics
  
  // Get the edge states
  ND::THandle<ND::TReconState> firstState, lastState;
  
  ND::tman().GetFirstAndLastStates(object, firstState, lastState);
  
  FillKinematics(firstState, globalObject->FrontPosition,
    globalObject->FrontDirection, globalObject->FrontMomentum,
    globalObject->FrontMomentumError, globalObject->FrontCharge,
    globalObject->FrontStateCov);
  
  FillKinematics(lastState, globalObject->BackPosition,
    globalObject->BackDirection, globalObject->BackMomentum,
    globalObject->BackMomentumError, globalObject->BackCharge,
    globalObject->BackStateCov);
  
  //------- Momentumby range
  
  // globalObject->RangeMomentum  = TrackingUtils::GetMomentumByRange(*object);
// globalObject->RangeMomentum  = TrackingUtils::GetMomentumByRange(*object);
    globalObject->RangeMomentumMuon =
    TrackingUtils::GetTRealDatumValue(*object, "prange_Muon", DEFAULT_MAX); 
    globalObject->RangeMomentumMuonFlip = TrackingUtils::GetTRealDatumValue(
        *object, "prange_Muon_flipped", DEFAULT_MAX);
    globalObject->RangeMomentumEndToTPCMuon =
    TrackingUtils::GetTRealDatumValue(*object, "prange_end_to_tpc_Muon", DEFAULT_MAX);
   
    globalObject->RangeMomentumPion =
    TrackingUtils::GetTRealDatumValue(*object, "prange_Pion", DEFAULT_MAX); 
    globalObject->RangeMomentumPionFlip = TrackingUtils::GetTRealDatumValue(
        *object, "prange_Pion_flipped", DEFAULT_MAX);
    globalObject->RangeMomentumEndToTPCPion =
    TrackingUtils::GetTRealDatumValue(*object, "prange_end_to_tpc_Pion", DEFAULT_MAX);
   
    globalObject->RangeMomentumElectron = TrackingUtils::GetTRealDatumValue(
        *object, "prange_Electron", DEFAULT_MAX);
    globalObject->RangeMomentumElectronFlip = TrackingUtils::GetTRealDatumValue(
        *object, "prange_Electron_flipped", DEFAULT_MAX);
    globalObject->RangeMomentumEndToTPCElectron =
    TrackingUtils::GetTRealDatumValue(*object, "prange_end_to_tpc_Electron", DEFAULT_MAX);
    
    globalObject->RangeMomentumProton =
    TrackingUtils::GetTRealDatumValue(*object, "prange_Proton", DEFAULT_MAX);
    globalObject->RangeMomentumProtonFlip = TrackingUtils::GetTRealDatumValue(
        *object, "prange_Proton_flipped", DEFAULT_MAX);
    globalObject->RangeMomentumEndToTPCProton =
    TrackingUtils::GetTRealDatumValue(*object, "prange_end_to_tpc_Proton", DEFAULT_MAX);
    
  
  //----------  Set the PID information
  ND::THandle<ND::TReconPID> PID = object;
  if (PID) {
    globalObject->PIDWeights.push_back(PID->GetPIDWeight());
    globalObject->ParticleIds.push_back(ComputeParticleId(PID));
    
    ND::TReconObjectContainer::const_iterator it;
    for (it = PID->GetAlternates().begin(); it != PID->GetAlternates().end();
        it++) {
      ND::THandle<ND::TReconPID> alter = *it;
      
      // check that this particle Id is not yet present (alternate for the main
      // hypothesis, sense alternates)
      int pid = ComputeParticleId(alter);
      if (std::find(globalObject->ParticleIds.begin(),
          globalObject->ParticleIds.end(),
          pid) == globalObject->ParticleIds.end()) {
        globalObject->PIDWeights.push_back(alter->GetPIDWeight());
        globalObject->ParticleIds.push_back(ComputeParticleId(alter));
      }
    }
  }
  
  //--------- Fill the information needed to do the association between the TREx
  // output and global recon ---------
  globalObject->TPCGasPatternIDs.clear();
  globalObject->TPCGasPathIDs.clear();
  globalObject->TPCGasJunctionIDs.clear();
  
  std::set<int> IDs;
  
  // Fill the patterns
  IDs.clear();
  ReconObjectUtils::GetTRExPatternIDs(object, IDs, kGas);
  
  globalObject->TPCGasPatternIDs.reserve(IDs.size());
  std::copy(IDs.begin(), IDs.end(),
    std::back_inserter(globalObject->TPCGasPatternIDs));
  
  // Fill the paths
  IDs.clear();
  ReconObjectUtils::GetTRExPathIDs(object, IDs, kGas);
  
  globalObject->TPCGasPathIDs.reserve(IDs.size());
  std::copy(IDs.begin(), IDs.end(),
    std::back_inserter(globalObject->TPCGasPathIDs));
  
  // Fill the junctions
  IDs.clear();
  ReconObjectUtils::GetTRExJunctionIDs(object, IDs, kGas);
  
  globalObject->TPCGasJunctionIDs.reserve(IDs.size());
  std::copy(IDs.begin(), IDs.end(),
    std::back_inserter(globalObject->TPCGasJunctionIDs));
  
  //---------- Extrapolate to the entrance and exit of each subdetector and save
  // position, direction and momentum
  FillExtrapolationToDetectors(object, *globalObject);
  
  //---------- get FGD vertex activity before calling FillFGDInfo but after
  // front/back position defined
  // Claculate Vertex Activity for Front position of reconstructed track
  if (ND::TGeomInfo::Get().FGD().IsInFGD1(globalObject->FrontPosition.Vect()) ||
      ND::TGeomInfo::Get().FGD().IsInFGD2(globalObject->FrontPosition.Vect()))
  GetFGDSimpleVA(event, object, globalObject->FrontPosition, 1, 1);
  
  // Claculate End Activity for Back position of reconstructed track
  if (ND::TGeomInfo::Get().FGD().IsInFGD1(globalObject->BackPosition.Vect()) ||
      ND::TGeomInfo::Get().FGD().IsInFGD2(globalObject->BackPosition.Vect()))
  GetFGDSimpleVA(event, object, globalObject->BackPosition, 1, 2);
  
  // Calculate Vertex and End Activities for Entering and Exiting positions to
  // FGD
  // Active volume of FGD1 and FGD2 as in oaGeomInfo TFGDGeom.cxx
  TVector3 fFGD1ActiveMin;
  TVector3 fFGD1ActiveMax;
  TVector3 fFGD2ActiveMin;
  TVector3 fFGD2ActiveMax;
  fFGD1ActiveMin.SetXYZ(-932.17, -877.17, 116.045);
  fFGD1ActiveMax.SetXYZ(932.17, 987.17, 446.955);
  fFGD2ActiveMin.SetXYZ(-932.17, -877.17, 1474.04);
  fFGD2ActiveMax.SetXYZ(932.17, 987.17, 1806.96);
  
  // FrontPosition in FGD1
  //-----------------------------
  if (ND::TGeomInfo::Get().FGD().IsInFGD1(
    globalObject->FrontPosition.Vect())) {  // Front in FGD1
    
    if (ND::TGeomInfo::Get().FGD().IsInFGD1(
      globalObject->BackPosition.Vect())) {  // Back in FGD1
      // std::cout<<"Front FGD1 Back FGD1"<<std::endl;
    } else if (ND::TGeomInfo::Get().FGD().IsInFGD2(
               globalObject->BackPosition.Vect())) {  // Back in FGD2
      
      if (globalObject->ExitPosition[3].Z() != -99999) {
          TLorentzVector ExitPosFGD1;
          ExitPosFGD1.SetX(globalObject->ExitPosition[3].X());
          ExitPosFGD1.SetY(globalObject->ExitPosition[3].Y());
          ExitPosFGD1.SetZ(globalObject->ExitPosition[3].Z());
          ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD1, fFGD1ActiveMin,
            fFGD1ActiveMax);
          GetFGDSimpleVA(event, object, ExitPosFGD1, 2, 2);
      }
      
      if (globalObject->EntrancePosition[4].Z() != -99999) {
          TLorentzVector EntrancePosFGD2;
          EntrancePosFGD2.SetX(globalObject->EntrancePosition[4].X());
          EntrancePosFGD2.SetY(globalObject->EntrancePosition[4].Y());
          EntrancePosFGD2.SetZ(globalObject->EntrancePosition[4].Z());
          ND::TGlobalReconModule::CheckFGDdist(EntrancePosFGD2, fFGD2ActiveMin,
            fFGD2ActiveMax);
          GetFGDSimpleVA(event, object, EntrancePosFGD2, 2, 1);
      }
      
    } else {  // Back !FGD1 and !FGD2
        
      if (globalObject->ExitPosition[3].Z() != -99999) {
          TLorentzVector ExitPosFGD1;
          ExitPosFGD1.SetX(globalObject->ExitPosition[3].X());
          ExitPosFGD1.SetY(globalObject->ExitPosition[3].Y());
          ExitPosFGD1.SetZ(globalObject->ExitPosition[3].Z());
          
          ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD1, fFGD1ActiveMin,
            fFGD1ActiveMax);
          GetFGDSimpleVA(event, object, ExitPosFGD1, 2, 2);
      }
    
      // check if track passed through FGD2
      if ((globalObject->EntrancePosition[4].Z() != -99999) &&
          (globalObject->ExitPosition[4].Z() != -99999)) {  
        // track Passed through FGD2
        TLorentzVector EntrancePosFGD2;
        EntrancePosFGD2.SetX(globalObject->EntrancePosition[4].X());
        EntrancePosFGD2.SetY(globalObject->EntrancePosition[4].Y());
        EntrancePosFGD2.SetZ(globalObject->EntrancePosition[4].Z());
        TLorentzVector ExitPosFGD2;
        ExitPosFGD2.SetX(globalObject->ExitPosition[4].X());
        ExitPosFGD2.SetY(globalObject->ExitPosition[4].Y());
        ExitPosFGD2.SetZ(globalObject->ExitPosition[4].Z());
        
        // check if distance between Entrance / ExitPosition is outside of FGD2
        // volume
        // if out-of-active distance is less than 2mm then assume edge of active
        // volume being exit position
        ND::TGlobalReconModule::CheckFGDdist(EntrancePosFGD2, fFGD2ActiveMin,
            fFGD2ActiveMax);
        ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD2, fFGD2ActiveMin,
            fFGD2ActiveMax);
        GetFGDSimpleVA(event, object, EntrancePosFGD2, 2, 1);
        GetFGDSimpleVA(event, object, ExitPosFGD2, 2, 2);
      }
    }
  }
  // FrontPosition in FGD2
  //-----------------------------
  else if (ND::TGeomInfo::Get().FGD().IsInFGD2(
      globalObject->FrontPosition.Vect())) {  // Front in FGD2
    if (ND::TGeomInfo::Get().FGD().IsInFGD2(
      globalObject->BackPosition.Vect())) {  // Back in FGD2
      // std::cout<<"Front FGD2 Back FGD2"<<std::endl;
    } else if (ND::TGeomInfo::Get().FGD().IsInFGD1(
      globalObject->BackPosition.Vect())) {  // Back in FGD1
      
      if (globalObject->ExitPosition[4].Z() != -99999) {
        // track exiting FGD2
        TLorentzVector ExitPosFGD2;
        ExitPosFGD2.SetX(globalObject->ExitPosition[4].X());
        ExitPosFGD2.SetY(globalObject->ExitPosition[4].Y());
        ExitPosFGD2.SetZ(globalObject->ExitPosition[4].Z());
        ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD2, fFGD2ActiveMin,
            fFGD2ActiveMax);
        GetFGDSimpleVA(event, object, ExitPosFGD2, 2, 2);
      }
      
      if (globalObject->EntrancePosition[3].Z() != -99999) {
        // track entering FGD1
        TLorentzVector EntrancePosFGD1;
        EntrancePosFGD1.SetX(globalObject->EntrancePosition[3].X());
        EntrancePosFGD1.SetY(globalObject->EntrancePosition[3].Y());
        EntrancePosFGD1.SetZ(globalObject->EntrancePosition[3].Z());
        ND::TGlobalReconModule::CheckFGDdist(EntrancePosFGD1, fFGD1ActiveMin,
            fFGD1ActiveMax);
        GetFGDSimpleVA(event, object, EntrancePosFGD1, 2, 1);
      }
    } else {  // Back !FGD1 and !FGD2
         
      if (globalObject->ExitPosition[4].Z() != -99999) {
        // track leaving FGD2
        TLorentzVector ExitPosFGD2;
        ExitPosFGD2.SetX(globalObject->ExitPosition[4].X());
        ExitPosFGD2.SetY(globalObject->ExitPosition[4].Y());
        ExitPosFGD2.SetZ(globalObject->ExitPosition[4].Z());
        ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD2, fFGD2ActiveMin,
            fFGD2ActiveMax);
        GetFGDSimpleVA(event, object, ExitPosFGD2, 2, 2);
      }
     
      // check if track passed through FGD1
      if ((globalObject->EntrancePosition[3].Z() != -99999) &&
          (globalObject->ExitPosition[3].Z() != -99999)) {  
        // Passed through FGD1
        TLorentzVector EntrancePosFGD1;
        EntrancePosFGD1.SetX(globalObject->EntrancePosition[3].X());
        EntrancePosFGD1.SetY(globalObject->EntrancePosition[3].Y());
        EntrancePosFGD1.SetZ(globalObject->EntrancePosition[3].Z());
        ND::TGlobalReconModule::CheckFGDdist(EntrancePosFGD1, fFGD1ActiveMin,
          fFGD1ActiveMax);
        
        TLorentzVector ExitPosFGD1;
        ExitPosFGD1.SetX(globalObject->ExitPosition[3].X());
        ExitPosFGD1.SetY(globalObject->ExitPosition[3].Y());
        ExitPosFGD1.SetZ(globalObject->ExitPosition[3].Z());
        ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD1, fFGD1ActiveMin,
          fFGD1ActiveMax);
        
        GetFGDSimpleVA(event, object, EntrancePosFGD1, 2, 1);
        GetFGDSimpleVA(event, object, ExitPosFGD1, 2, 2);
      }
    }
  }
      
  // FrontPosition outside of FGD1 and FGD2
  //-----------------------------
  else {  // Front !FGD1 and !FGD2
    // check if track enters FDG1
    if (globalObject->EntrancePosition[3].Z() != -99999) {
        TLorentzVector EntrancePosFGD1;
        EntrancePosFGD1.SetX(globalObject->EntrancePosition[3].X());
        EntrancePosFGD1.SetY(globalObject->EntrancePosition[3].Y());
        EntrancePosFGD1.SetZ(globalObject->EntrancePosition[3].Z());
        ND::TGlobalReconModule::CheckFGDdist(EntrancePosFGD1, fFGD1ActiveMin,
            fFGD1ActiveMax);
        GetFGDSimpleVA(event, object, EntrancePosFGD1, 2, 1);
    }
    
    if (globalObject->ExitPosition[3].Z() !=
        -99999) {  // check if track leaves FGD1
      TLorentzVector ExitPosFGD1;
      ExitPosFGD1.SetX(globalObject->ExitPosition[3].X());
      ExitPosFGD1.SetY(globalObject->ExitPosition[3].Y());
      ExitPosFGD1.SetZ(globalObject->ExitPosition[3].Z());
      ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD1, fFGD1ActiveMin,
          fFGD1ActiveMax);
      GetFGDSimpleVA(event, object, ExitPosFGD1, 2, 2);
    }
        
    if (globalObject->EntrancePosition[4].Z() !=
      -99999) {  // check if track enters FGD2
      TLorentzVector EntrancePosFGD2;
      EntrancePosFGD2.SetX(globalObject->EntrancePosition[4].X());
      EntrancePosFGD2.SetY(globalObject->EntrancePosition[4].Y());
      EntrancePosFGD2.SetZ(globalObject->EntrancePosition[4].Z());
      ND::TGlobalReconModule::CheckFGDdist(EntrancePosFGD2, fFGD2ActiveMin,
        fFGD2ActiveMax);
      GetFGDSimpleVA(event, object, EntrancePosFGD2, 2, 1);
    }
        
    if (globalObject->ExitPosition[4].Z() !=
          -99999) {  // check if track leaves FGD2
      TLorentzVector ExitPosFGD2;
      ExitPosFGD2.SetX(globalObject->ExitPosition[4].X());
      ExitPosFGD2.SetY(globalObject->ExitPosition[4].Y());
      ExitPosFGD2.SetZ(globalObject->ExitPosition[4].Z());
      ND::TGlobalReconModule::CheckFGDdist(ExitPosFGD2, fFGD2ActiveMin,
          fFGD2ActiveMax);
      GetFGDSimpleVA(event, object, ExitPosFGD2, 2, 2);
    }
  }
  // end calulation of Entering and Exiting FGD activities
  
  //---------- Fill subdetector information
  FillECALInfo   (object, *globalObject);
  FillP0DInfo    (object, *globalObject);
  FillTPCInfo    (object, *globalObject);
  FillFGDInfo    (object, *globalObject);
  FillSFGInfo    (object, *globalObject);
  FillSMRDInfo   (object, *globalObject);
  FillTrackerInfo(object, *globalObject);
  
  //---------- Fill the alternates information
  FillGlobalPIDAlternates(G4track, object, *globalObject);
  
  // Add subdetector recomputation of some variables for isolated objects
  // For ECAL
  if (det_count == 1 && globalObject->NECALs == 1) {
    TECALObject* sub = (TECALObject*)(*(globalObject->ECAL))[0];
    globalObject->Cone = sub->Cone;
    globalObject->EDeposit = sub->EDeposit;
    globalObject->Width = sub->Width.X();
    if (globalObject->Length == 0) globalObject->Length = sub->Length;
  }
  // For P0D
  else if (det_count == 1 && globalObject->NP0Ds == 1) {
    TP0DObject* sub = (TP0DObject*)(*(globalObject->P0D))[0];
    globalObject->Cone = sub->Cone;
    globalObject->EDeposit = sub->EDeposit;
    globalObject->Width = sub->Width;
    if (globalObject->Length == 0) globalObject->Length = sub->Length;
  }
}

//*************************************************************
void ND::TGlobalReconModule::CheckFGDdist(TLorentzVector& FGDPos,
    TVector3 fFGDActiveMin,
    TVector3 fFGDActiveMax) 
{
    //*************************************************************
    TVector3 dist;
    
    dist.SetXYZ(0.0, 0.0, 0.0);
    if (FGDPos.X() < fFGDActiveMin.X()) {
        dist.SetX(fFGDActiveMin.X() - FGDPos.X());
        if (dist.X() != 0.0 && abs(dist.X()) < 2.0)
            FGDPos.SetX(fFGDActiveMin.X() + 0.1);
    }
    if (FGDPos.X() > fFGDActiveMax.X()) {
        dist.SetX(fFGDActiveMax.X() - FGDPos.X());
        if (dist.X() != 0.0 && abs(dist.X()) < 2.0)
            FGDPos.SetX(fFGDActiveMax.X() - 0.1);
    }
    if (FGDPos.Y() < fFGDActiveMin.Y()) {
        dist.SetY(fFGDActiveMin.Y() - FGDPos.Y());
        if (dist.Y() != 0.0 && abs(dist.Y()) < 2.0)
            FGDPos.SetY(fFGDActiveMin.Y() + 0.1);
    }
    if (FGDPos.Y() > fFGDActiveMax.Y()) {
        dist.SetY(fFGDActiveMax.Y() - FGDPos.Y());
        if (dist.Y() != 0.0 && abs(dist.Y()) < 2.0)
            FGDPos.SetY(fFGDActiveMax.Y() - 0.1);
    }
    if (FGDPos.Z() < fFGDActiveMin.Z()) {
        dist.SetZ(fFGDActiveMin.Z() - FGDPos.Z());
        if (dist.Z() != 0.0 && abs(dist.Z()) < 2.0)
            FGDPos.SetZ(fFGDActiveMin.Z() + 0.1);
    }
    if (FGDPos.Z() > fFGDActiveMax.Z()) {
        dist.SetZ(fFGDActiveMax.Z() - FGDPos.Z());
        if (dist.Z() != 0.0 && abs(dist.Z()) < 2.0)
            FGDPos.SetZ(fFGDActiveMax.Z() - 0.1);
    }
}

//*************************************************************
int ND::TGlobalReconModule::ComputeParticleId(ND::THandle<ND::TReconPID> PID) {
    //*************************************************************
    
    // only three hypothesis considered at recon level
    
    int pdg = 0;
    
    if (PID->GetParticleId() == ND::TReconPID::kElectron)
        pdg = 11;
    else if (PID->GetParticleId() == ND::TReconPID::kMuon)
        pdg = 13;
    else if (PID->GetParticleId() == ND::TReconPID::kProton)
        pdg = 2212;
    
    return pdg;
}

//*************************************************************
void ND::TGlobalReconModule::FillKinematics(ND::THandle<ND::TReconState> state,
    TLorentzVector& pos, TVector3& dir,
    double& mom, double& momErr,
    double& charge,
    double (&cov)[7][7]) 
{
//*************************************************************
    
    if (state) {
        pos = TrackingUtils::GetPosition(state);
        dir = TrackingUtils::GetDirection(state);
        mom = TrackingUtils::GetMomentum(state);
        momErr = TrackingUtils::GetMomentumError(state);
        charge = TrackingUtils::GetCharge(state);
        
        // get the cov matrix
        RP::State rp_state;
        bool ok = ND::converter().TReconState_to_State(*state, rp_state);
        size_t dim = 0;
        
        if (ok) {
            RP::rep().convert(rp_state, RP::pos_dir_curv);
            dim = rp_state.vector().num_row();
        }
        
        // should have vector of 7: pos,  dir, curv
        if (dim == 7) {
            const EMatrix& C = rp_state.matrix();
            
            for (size_t i = 0; i < 7; i++)
                for (size_t j = 0; j < 7; j++) cov[i][j] = C[i][j];
            
        } else {
            std::fill(&cov[0][0], &cov[0][0] + 7 * 7, DEFAULT_MAX);
        }
        
    }  // state is ok
    else {
        pos = TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        dir = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        mom = DEFAULT_MAX;
        momErr = DEFAULT_MAX;
        charge = DEFAULT_MAX;
        
        std::fill(&cov[0][0], &cov[0][0] + 7 * 7, DEFAULT_MAX);
    }
    
    if (debug_kine) {
        std::cout << "position:  " << pos.X() << " " << pos.Y() << " " << pos.Z()
        << std::endl;
        std::cout << "direction: " << dir.X() << " " << dir.Y() << " " << dir.Z()
        << std::endl;
        std::cout << "momentum:  " << mom << " " << momErr << std::endl;
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillKinematics(ND::THandle<ND::TReconState> state,
    TLorentzVector& pos,
    TLorentzVector& posVar,
    TVector3& dir, TVector3& dirVar,
    double& mom, double& momErr)
{
    //*************************************************************
    
    if (state) {
        pos = TrackingUtils::GetPosition(state);
        posVar = TrackingUtils::GetPositionVariance(state);
        dir = TrackingUtils::GetDirection(state);
        dirVar = TrackingUtils::GetDirectionVariance(state);
        mom = TrackingUtils::GetMomentum(state);
        momErr = TrackingUtils::GetMomentumError(state);
    } else {
        pos = TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        posVar = TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        dir = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        dirVar = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        mom = DEFAULT_MAX;
        momErr = DEFAULT_MAX;
    }
    
    if (debug_kine) {
        std::cout << "position:  " << pos.X() << " " << pos.Y() << " " << pos.Z()
        << std::endl;
        std::cout << "direction: " << dir.X() << " " << dir.Y() << " " << dir.Z()
        << std::endl;
        std::cout << "momentum:  " << mom << " " << momErr << std::endl;
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillKinematicsAtTrueVertex(
    ND::THandle<ND::TG4Trajectory> G4track, ND::THandle<ND::TReconBase> object,
    TLorentzVector& pos, TVector3& dir, double& mom, double& momErr,
    double& charge, double (&cov)[7][7]) 
{
    //*************************************************************
    
    // Propagate the first state in the track to a surface normal to z at the
    // vertex position.
    
    ND::THandle<ND::TReconState> tstate = object->GetState();
    
    ND::THandle<ND::TReconState> newstate = ND::THandle<ND::TReconState>();
    
    // The G4 track should exist and the track should have TPC, FGD, SFG or P0D
    // components
    if (tstate && G4track && (object->UsesDetector(ND::TReconBase::kTPC) ||
          object->UsesDetector(ND::TReconBase::kFGD) ||
          object->UsesDetector(ND::TReconBase::kSFG) ||
        object->UsesDetector(ND::TReconBase::kP0D))) {
      TVector3 snorm;
      double dir_max = -100;
      int idir_max = 2;
      TVector3 dirS = TrackingUtils::GetDirection(object);
      for (int i = 0; i < 3; i++) {
        if (fabs(dirS(i)) > dir_max) {
            dir_max = fabs(dirS(i));
            idir_max = i;
        }
      }
      if (idir_max == 0)
        snorm = TVector3(1, 0, 0);
      else if (idir_max == 1)
        snorm = TVector3(0, 1, 0);
      else if (idir_max == 2)
        snorm = TVector3(0, 0, 1);
      
      if (debug_vertex_extrap){
        std::cout << "Extrapolation to true vertex: dir = (" << dirS.X() << ", "
                  << dirS.Y() << ", " << dirS.Z() << ")" << std::endl;
      }
      if (debug_vertex_extrap)
        ND::tman("globalAnalysis").SetVerbosity(0, 6, 0, 0);
      
      newstate = ND::THandle<ND::TReconState>(new ND::TTrackState);
      ND::tman("globalAnalysis").PropagateToSurface(
          *tstate, G4track->GetInitialPosition().Vect(),
          snorm, *newstate);
      
      if (debug_vertex_extrap)
          ND::tman("globalAnalysis").SetVerbosity(0, 0, 0, 0);
    }
    
    FillKinematics(newstate, pos, dir, mom, momErr, charge, cov);
}
//*************************************************************
void ND::TGlobalReconModule::FillKinematicsAtTrueVertex(
    ND::THandle<ND::TG4Trajectory> G4track, ND::THandle<ND::TReconBase> object,
    TLorentzVector& pos, TLorentzVector& posVar, TVector3& dir,
    TVector3& dirVar, double& mom, double& momErr) 
{
    //*************************************************************
    
    // Propagate the first state in the track to a surface normal to z at the
    // vertex position.
    
    ND::THandle<ND::TReconState> tstate = object->GetState();
    
    ND::THandle<ND::TReconState> newstate = ND::THandle<ND::TReconState>();
    
    // The G4 track should exist and the track should have TPC, FGD, SFG or P0D
    // components
    if (tstate && G4track && (object->UsesDetector(ND::TReconBase::kTPC) ||
          object->UsesDetector(ND::TReconBase::kFGD) ||
          object->UsesDetector(ND::TReconBase::kSFG) ||
        object->UsesDetector(ND::TReconBase::kP0D))) {

    TVector3 snorm;
    double dir_max = -100;
    int idir_max = 2;
    TVector3 dirS = TrackingUtils::GetDirection(object);
    for (int i = 0; i < 3; i++) {
        if (fabs(dirS(i)) > dir_max) {
            dir_max = fabs(dirS(i));
            idir_max = i;
        }
    }
    if (idir_max == 0)
        snorm = TVector3(1, 0, 0);
    else if (idir_max == 1)
        snorm = TVector3(0, 1, 0);
    else if (idir_max == 2)
        snorm = TVector3(0, 0, 1);
    
    if (debug_vertex_extrap)
        std::cout << "Extrapolation to true vertex: dir = (" << dirS.X() << ", "
    << dirS.Y() << ", " << dirS.Z() << ")" << std::endl;
    
    if (debug_vertex_extrap)
        ND::tman("globalAnalysis").SetVerbosity(0, 6, 0, 0);
    
    newstate = ND::THandle<ND::TReconState>(new ND::TTrackState);
    ND::tman("globalAnalysis")
    .PropagateToSurface(*tstate, G4track->GetInitialPosition().Vect(),
        snorm, *newstate);
    
    if (debug_vertex_extrap)
        ND::tman("globalAnalysis").SetVerbosity(0, 0, 0, 0);
    }
    
    FillKinematics(newstate, pos, posVar, dir, dirVar, mom, momErr);
}

//*************************************************************
void ND::TGlobalReconModule::FillGlobalPIDAlternates(
    ND::THandle<ND::TG4Trajectory> G4track, ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject)
{
    //*************************************************************
    
    // __________________________________________________________
    // Fill the Alternates information for Global PIDs.
    
    // make sure it is a PID
    ND::THandle<ND::TReconPID> pid = object;
    if (!pid) return;
    
    if (!(object->UsesDetector(ND::TReconBase::kTPC)) ||
        object->CheckStatus(ND::TReconBase::kLikelihoodFit))
    return;
    
    // do nothing for TPC-only objects
    if (object->GetDetectors() == ND::TReconBase::kTPC1 ||
        object->GetDetectors() == ND::TReconBase::kTPC2 ||
    object->GetDetectors() == ND::TReconBase::kTPC3)
    return;
    
    // do nothing for curving back tracks
    if (TrackingUtils::GetCurvBackTrack(*object)) return;
    
    const ND::TReconObjectContainer& alter = pid->GetAlternates();
    ND::TReconObjectContainer::const_iterator it;
    for (it = alter.begin(); it != alter.end(); it++) {
        // create a new TGlobalPIDAlternate
        TGlobalPIDAlternate* PIDAlternate =
        new ((*globalObject.Alternates)[globalObject.NAlternates++])
        TGlobalPIDAlternate;
        FillGlobalPIDAlternate(G4track, *it, *PIDAlternate);
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillGlobalPIDAlternate(
    ND::THandle<ND::TG4Trajectory> G4track, ND::THandle<ND::TReconBase> object,
    TGlobalPIDAlternate& PIDAlternate)
{
    //*************************************************************
    
    // __________________________________________________________
    // Fill the Alternates information for inner (P0D+Tracker) Global PIDs.
    
    // try to get the state of the object. If it dosen't exists skip the object.
    try {
        ND::THandle<ND::TReconState> tstate = object->GetState();
    } catch (ND::EReconObject e) {
        std::cout << "TReconBase with no state. Skip object !!" << std::endl;
        return;
    }
    
    //--------- Information about the detectors that used in this object
    PIDAlternate.Detectors = GetDetectorNumber(object);
    FillDetectorUsed(object, PIDAlternate.DetectorUsed);
    
    //-------  General information
    
    PIDAlternate.Status = object->CheckStatus(object->kSuccess);
    PIDAlternate.Chi2 = object->GetQuality();
    PIDAlternate.NDOF = object->GetNDOF();
    PIDAlternate.Charge = TrackingUtils::GetCharge(object);
    PIDAlternate.Length = ComputeTrackLength(object);
    
    //------- Front and Back kinematics
    ND::THandle<ND::TReconState> firstState =
    TrackingUtils::GetFirstState(*object);
    FillKinematics(firstState, PIDAlternate.FrontPosition,
        PIDAlternate.FrontPositionVar, PIDAlternate.FrontDirection,
        PIDAlternate.FrontDirectionVar, PIDAlternate.FrontMomentum,
        PIDAlternate.FrontMomentumError);
    
    ND::THandle<ND::TReconState> lastState = TrackingUtils::GetLastState(*object);
    FillKinematics(lastState, PIDAlternate.BackPosition,
        PIDAlternate.BackPositionVar, PIDAlternate.BackDirection,
        PIDAlternate.BackDirectionVar, PIDAlternate.BackMomentum,
        PIDAlternate.BackMomentumError);
    
    //-------- Fill the kinematic extrapolated to the true vertex
    FillKinematicsAtTrueVertex(
        G4track, object, PIDAlternate.PositionAtTrueVertex,
        PIDAlternate.PositionVarAtTrueVertex, PIDAlternate.DirectionAtTrueVertex,
        PIDAlternate.DirectionVarAtTrueVertex, PIDAlternate.MomentumAtTrueVertex,
        PIDAlternate.MomentumErrorAtTrueVertex);
    
    //---------- Set the PID information
    ND::THandle<ND::TReconPID> PID = object;
    if (PID) {
        PIDAlternate.PIDWeight = PID->GetPIDWeight();
        PIDAlternate.ParticleId = ComputeParticleId(PID);
    }
}

//*************************************************************
double ND::TGlobalReconModule::ComputeTrackLength(
    ND::THandle<ND::TReconBase> object) 
{
    //*************************************************************
    
    // Default for showers is 0
    double length = 0;
    
    double length_tmp;
    
    length = DEFAULT_MAX;
    
    if (!object) return length;
    
    if (TrackingUtils::IsShower(object)) return 0.;
    
    ND::THandle<ND::TReconPID> pid = object;
    ND::THandle<ND::TReconTrack> track = object;
    
    if (!pid && !track) return length;
    
    bool ok = ND::tman().ComputeObjectLength(object, length_tmp);
    
    if (ok) length = length_tmp;
    
    return length;
}

//*************************************************************
void ND::TGlobalReconModule::FillVertexInfo(
    ND::TND280Event& event,
    ND::THandle<ND::TReconObjectContainer> globalObjects)
{
//*************************************************************
    // Container names
    // containers are "trakerPrimaryVertex" and "trackerVertices"
    TString tname  = "trackerRecon";
    TString tvname = "tracker";  
    TString gname  = "globalVertex";
    TString gvname = "globalVertex";
    
    // We assume it's global vertices first
    bool global = true;
    TString reconame = gname;
    TString vname = gvname;
    
    if (fUseECalVertices) {
        ND::THandle<ND::TAlgorithmResult> ResultBin = event.GetFit("globalVertex");
        
        ND::THandle<ND::TReconObjectContainer> pvertices =
        ResultBin->GetResultsContainer("vertices");
        if (!pvertices) {
            ND280NamedInfo("ECalVertexInTGRM",
            "Leaving FillVertexInfo as we couldn't find the container.");
            return;
        }
        
        ND::THandle<ND::TG4PrimaryVertexContainer> g4PrimVerts;
        if (fIsMC) {
            g4PrimVerts =
            event.Get<ND::TG4PrimaryVertexContainer>("truth/G4PrimVertex00");
        }
        for (size_t prim_vert_it = 0; prim_vert_it < pvertices->size();
            ++prim_vert_it) {
            ND::THandle<ND::TReconVertex> vert = pvertices->at(prim_vert_it);
            if (!vert) {
                continue;
            }
            
            // Fill the primary
            FillVertex(event,vert, true);
            
            if (fIsMC && bool(g4PrimVerts)) {
                MatchTrueVertex(g4PrimVerts);
            }
            
            if (vert->GetConstituents()) {
                for (size_t sec_vert_it = 0;
                    sec_vert_it < vert->GetConstituents()->size(); ++sec_vert_it) {
                    ND::THandle<ND::TReconVertex> sVert =
                    vert->GetConstituents()->at(sec_vert_it);
                    
                    if (!sVert) {
                        continue;
                    }
                    bool cool_beans = FillVertex(event,sVert, false);
                    if (cool_beans) {
                        ND280NamedInfo("ECalVertexInTGRM", "Added secondary vertex.");
                        if (fIsMC && bool(g4PrimVerts)) {
                            MatchTrueVertex(g4PrimVerts);
                        }
                    }
                }
            }
            
            ND280NamedInfo("ECalVertexInTGRM", "Added an ecalVertex.");
        }
        ND280NamedInfo("ECalVertexInTGRM", "Leaving FillVertexInfo having examined "
                       << pvertices->size() << " vertices.");
    }  // End us ECal Vertices
    else {
        // Fill primary and secondary vertex info
        //**********************************************************
        for (int ibin = 0; ibin < NTIMEBINS; ibin++) {
            TString bin_result_name = reconame + "_timebin";
            bin_result_name += ibin;
            
            ND::THandle<ND::TAlgorithmResult> ResultBin =
            event.GetFit(bin_result_name);
            // Check presence of global vertexing, if not try tracker vertices
            if (ibin == 0 && !ResultBin) {
                global = false;
                reconame = tname;
                vname = tvname;
                ResultBin = event.GetFit("trackerRecon_timebin0");
            }
            
            // If we didn't find a bin, stop looking
            if (!ResultBin) break;
            
            // Fill tracker to global tracks map
            if (!global) {
                TString container_name = "finalTrackerObjects";
                ND::THandle<ND::TReconObjectContainer> finalTracks =
                ResultBin->GetResultsContainer(container_name);
                if (!finalTracks) continue;
                if (finalTracks->size() == 0) continue;
                fTrackerGlobalIndexMap.clear();
                DoAssociationBetweenTrackerAndGlobalObjects(*globalObjects,
                    *finalTracks);
            }
            
            TString primary_name = vname + "PrimaryVertex";
            ND::THandle<ND::TReconObjectContainer> pvertices =
            ResultBin->GetResultsContainer(primary_name);
            if (!pvertices) continue;
            if (pvertices->size() == 0) continue;
            
            ND::THandle<ND::TReconVertex> pv = pvertices->front();
            if (!pv) continue;
            
            // Fill the primary vertex, if failed,  do not proceed
            if (!FillVertex(event,pv, true)) continue;
            
            ND::THandle<ND::TG4PrimaryVertexContainer> g4PrimVert =
            event.Get<ND::TG4PrimaryVertexContainer>("truth/G4PrimVertex00");
            if (g4PrimVert) MatchTrueVertex(g4PrimVert);
            
            TString vertices_name = vname + "Vertices";
            ND::THandle<ND::TReconObjectContainer> svertices =
            ResultBin->GetResultsContainer(vertices_name);
            if (!svertices) continue;
            
            for (ND::TReconObjectContainer::iterator siter = svertices->begin();
                 siter != svertices->end(); siter++) {
                ND::THandle<ND::TReconVertex> sv = *siter;
                if (!sv) continue;
                
                // Skip the primary vertex
                if (sv == pv) continue;
                
                bool vertexCreated = false;
                
                // Fill secondary vertex
                vertexCreated = FillVertex(event,sv, false);
                
                if (g4PrimVert && vertexCreated) MatchTrueVertex(g4PrimVert);
            }
        }
    }
    
    fTrackerGlobalIndexMap.clear();
}

//*************************************************************
void ND::TGlobalReconModule::FillGlobalHit(ND::THandle<ND::THit> hit,
    TGlobalHit& gHit) 
{
    //*************************************************************
    
    if (hit) {
        gHit.Position = hit->GetPosition();
        gHit.Time = hit->GetTime();
        gHit.PositionError = hit->GetUncertainty();
        gHit.Charge = hit->GetCharge();
        gHit.Type = hit->IsXHit() + (hit->IsYHit()) * 10 + (hit->IsZHit()) * 100;
    } else {
        gHit.Position = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        gHit.Time = DEFAULT_MAX;
        gHit.PositionError = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        gHit.Charge = DEFAULT_MAX;
        gHit.Type = DEFAULT_MAX;
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillGlobalHit(ND::THandle<ND::THit> hit,
    TSMRDHit& smrdHit) 
{
    //*************************************************************
    
    if (hit) {
        smrdHit.Position = hit->GetPosition();
        smrdHit.Time = hit->GetTime();
        smrdHit.PositionError = hit->GetUncertainty();
        smrdHit.Charge = hit->GetCharge();
        smrdHit.Type = hit->IsXHit() + (hit->IsYHit()) * 10 + (hit->IsZHit()) * 100;
        
    } else {
        smrdHit.Position = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        smrdHit.Time = DEFAULT_MAX;
        smrdHit.PositionError = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        smrdHit.Charge = DEFAULT_MAX;
        smrdHit.Type = DEFAULT_MAX;
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillSmrdHit(ND::THandle<ND::THit> hit,
    TSMRDHit& smrdHit) 
{
    //*************************************************************
    enum ESMRDWall { kTop = 0, kBottom, kRight, kLeft };
    // fill general info
    FillGlobalHit(hit, smrdHit);
    if (hit) {
        bool isTopWall = ND::TGeomInfo::SMRD().IsSMRDTopWall(hit);
        bool isBottomWall = ND::TGeomInfo::SMRD().IsSMRDBottomWall(hit);
        bool isRightWall = ND::TGeomInfo::SMRD().IsSMRDRightWall(hit);
        bool isLeftWall = ND::TGeomInfo::SMRD().IsSMRDLeftWall(hit);
        
        if (isTopWall) smrdHit.Wall = kTop;
        if (isBottomWall) smrdHit.Wall = kBottom;
        if (isRightWall) smrdHit.Wall = kRight;
        if (isLeftWall) smrdHit.Wall = kLeft;
        smrdHit.Yoke = ND::TGeomInfo::SMRD().GetYokeRingNumber(hit);
        smrdHit.Layer = ND::TGeomInfo::SMRD().GetSMRDRingNumber(hit);
        smrdHit.Tower = ND::TGeomInfo::SMRD().GetSMRDTowerNumber(hit);
        smrdHit.Counter = ND::TGeomInfo::SMRD().GetSMRDScintNumber(hit);
        
    } else {
        // fill with unmeaningful values
        smrdHit.Wall = -1;
        smrdHit.Yoke = -1;
        smrdHit.Layer = -1;
        smrdHit.Tower = -1;
        smrdHit.Counter = -1;
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillSubBaseObject(
    ND::THandle<ND::TReconBase> object, TSubBaseObject& sub, int det)
{
    //*************************************************************
    
    if (ReconObjectUtils::GetOriginalObject(object))
        sub.UniqueID = ReconObjectUtils::GetOriginalObject(object)->GetUniqueID();
    else
        sub.UniqueID = object->GetUniqueID();
    
    sub.Detector = det;
    sub.Status = object->CheckStatus(object->kSuccess);
    sub.NDOF = object->GetNDOF();
    try {
        if (object->Get<ND::TIntegerDatum>("NNodes"))
            sub.NNodes = object->Get<ND::TIntegerDatum>("NNodes")
        ->GetValue();  // If this is a TPC object,should use the
        // new (P6) method of getting the number
        // of nodes
        else
            sub.NNodes = object->GetNodes().size();
    } catch (ND::EReconObject const& e) {
        sub.NNodes = 0;
    }
    
    sub.Chi2 = object->GetQuality();
    sub.EDeposit = TrackingUtils::GetEDeposit(object);
    
    // The number of constituents. For composite objects
    sub.NConstituents = 0;
    if (object->GetConstituents())
        sub.NConstituents = object->GetConstituents()->size();
    
    sub.NHits = DEFAULT_MAX;
    if (object->GetHits()) sub.NHits = object->GetHits()->size();
    
    // the length of the object for track-like objects
    sub.Length = ComputeTrackLength(object);
    
    //--------- Front and back Kinematics
    // Get the edge states
    ND::THandle<ND::TReconState> firstState, lastState;
    
    ND::tman().GetFirstAndLastStates(object, firstState, lastState);
    
    FillKinematics(firstState, sub.FrontPosition, sub.FrontPositionVar,
        sub.FrontDirection, sub.FrontDirectionVar, sub.FrontMomentum,
        sub.FrontMomentumError);
    
    FillKinematics(lastState, sub.BackPosition, sub.BackPositionVar,
        sub.BackDirection, sub.BackDirectionVar, sub.BackMomentum,
        sub.BackMomentumError);
    
    // Fill the true particle when it is found
    double pur, eff;
    ND::THandle<ND::TG4Trajectory> G4track = GetG4Trajectory(*object, pur, eff);
    FillTrueParticle(G4track, pur, eff, sub.TrueParticle);
}

namespace {
    template <typename T, size_t N>
    size_t ArrSize(T const (&arr)[N]) {
        return N;
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillECALInfo(ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject)
{
    //*************************************************************
    
    static const ND::TReconBase::Status ECalDetBits[] = {
        ND::TReconBase::kLeftPECal, ND::TReconBase::kRightPECal,
        ND::TReconBase::kTopPECal,  ND::TReconBase::kBottomPECal,
        ND::TReconBase::kLeftTECal, ND::TReconBase::kRightTECal,
        ND::TReconBase::kTopTECal,  ND::TReconBase::kBottomTECal,
    ND::TReconBase::kDSECal};
    
    // Get constituents in any of the calorimeters
    for (size_t det_bit = 0; det_bit < ArrSize(ECalDetBits); det_bit++) {
        if (!object->UsesDetector(ECalDetBits[det_bit])) {
            continue;
        }
        
        // one needs to take into account that ecalRecon doesn`t combine tracks
        // that
        // cross
        // Left/Right sub-volumes for TECal Top and Bottom
        ND::THandle<ND::TReconObjectContainer> ecalObjects(
            new ND::TReconObjectContainer());
        
        ReconObjectUtils::GetConstituentsInDetector(object, ECalDetBits[det_bit],
            *ecalObjects);
        
        for (size_t eo_it = 0; eo_it < ecalObjects->size(); eo_it++) {
            ND::THandle<ND::TReconBase> ecalObject = ecalObjects->at(eo_it);
            
            if (!ecalObject) {
                continue;
            }
            
            ND280NamedVerbose(
                "TGlobalReconModule",
                "Transposing ECal object from: " << ecalObject->ConvertDetector());
            
            TECALObject* ECalObjToFill;
            // P0DECAL constituents
            if (det_bit < 4) {
                if (globalObject.NP0DECALs >= NMAXP0DECAL) {
                    ND280NamedWarn("TGlobalReconModule",
                        "Ignoring as already have too many P0DECal objects.");
                    continue;
                }
                ECalObjToFill =
                new ((*(globalObject.P0DECAL))[globalObject.NP0DECALs++])
                TECALObject;
            }
            // TrECAL constituents
            else if (det_bit < 8) {
                if (globalObject.NECALs >= NMAXECAL) {
                    ND280NamedWarn("TGlobalReconModule",
                        "Ignoring as already have too many ECal objects.");
                    continue;
                }
                ECalObjToFill =
                new ((*(globalObject.ECAL))[globalObject.NECALs++]) TECALObject;
                globalObject.NTrECALs++;
            }
            // DsECAL constituents
            else {
                if (globalObject.NECALs >= NMAXECAL) {
                    ND280NamedWarn("TGlobalReconModule",
                        "Ignoring as already have too many ECal objects.");
                    continue;
                }
                ECalObjToFill =
                new ((*(globalObject.ECAL))[globalObject.NECALs++]) TECALObject;
                globalObject.NDsECALs++;
            }
            
            // fill general information
            FillSubBaseObject(ecalObject, *ECalObjToFill, det_bit + 1);
            // Still need to fill TSubBaseShowerObject::Cone
            
            // Set defaults:
            ECalObjToFill->LLR_Quality = -9999.;      
            ECalObjToFill->LLR_MIP_EM = -9999.;
            ECalObjToFill->LLR_MIP_Pion = -9999.;
            ECalObjToFill->LLR_EM_HIP = -9999.;
            ECalObjToFill->LLR_MIP_EM_LowMomentum = -9999.;
            ECalObjToFill->LLR_MIP_HIP_VACut = -9999.;
            ECalObjToFill->AverageHitTime = DEFAULT_MAX;
            ECalObjToFill->AverageZPosition = DEFAULT_MAX;
            ECalObjToFill->EMEnergyFit_Result = DEFAULT_MAX;
            ECalObjToFill->EMEnergyFit_Uncertainty = DEFAULT_MAX;
            ECalObjToFill->FilledAsType = -1;
            ECalObjToFill->PID_Circularity = DEFAULT_MAX;
            ECalObjToFill->PID_Angle = DEFAULT_MAX;
            ECalObjToFill->PID_ShowerAngle = DEFAULT_MAX;
            ECalObjToFill->PID_Asymmetry = DEFAULT_MAX;
            ECalObjToFill->PID_TruncatedMaxRatio = DEFAULT_MAX;
            ECalObjToFill->PID_TransverseChargeRatio = DEFAULT_MAX;
            ECalObjToFill->PID_FrontBackRatio = DEFAULT_MAX;
            ECalObjToFill->Containment = DEFAULT_MAX;
            ECalObjToFill->Cone = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            ECalObjToFill->EDeposit = DEFAULT_MAX;
            ECalObjToFill->Thrust = DEFAULT_MAX;
            ECalObjToFill->ThrustOrigin =
            TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            ECalObjToFill->ThrustAxis =
            TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            ECalObjToFill->InnerMostLayerHit = -1;
            ECalObjToFill->OuterMostLayerHit = -1;
            ECalObjToFill->Width = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            
            if ((det_bit < 4) || !fUseECalVertices) {
                if (!ecalObject->GetConstituents()) {
                    ND280NamedWarn("TGlobalReconModule",
                        "ECal PID contained no constituents.");
                    continue;
                }
                // Unlike TTrackerECALReconModule, we're just going to assume the first
                // constituent as thats what global did before.
                ND::THandle<ND::TReconTrack> track_constit =
                (*ecalObject->GetConstituents()->begin());
                ND::THandle<ND::TReconShower> shower_constit =
                (*ecalObject->GetConstituents()->begin());
                ND::TTrackerECALReconModule::TECALReconObject ECalIsoObjToFill;
                // Only the implicit dynamic cast that is valid will result in a new
                // object.
                if (track_constit) {
                    if (!track_constit->GetEDeposit()) {
                        ND280NamedWarn("TGlobalReconModule",
                            "ECal PID-constituent track deposited no energy.");
                        continue;
                    }
                    // Making a TECALReconObject ClonesArray entry...
                    ND::TTrackerECALReconModule::TECALReconObjectFactory(
                        track_constit, ecalObject, false, &ECalIsoObjToFill);
                    
                    ECalObjToFill->LLR_Quality = ECalIsoObjToFill.PID_LLR_Quality;
                    if (ECalObjToFill->LLR_Quality == 0) {
                        // only fill the PID variables if the PID has good quality.
                        ECalObjToFill->LLR_MIP_EM = ECalIsoObjToFill.PID_LLR_MIP_EM;
                        ECalObjToFill->LLR_MIP_Pion = ECalIsoObjToFill.PID_LLR_MIP_Pion;
                        ECalObjToFill->LLR_EM_HIP = ECalIsoObjToFill.PID_LLR_EM_HIP;
                        ECalObjToFill->LLR_MIP_EM_LowMomentum =
                        ECalIsoObjToFill.PID_LLR_MIP_EM_LowMomentum;
                    }
                    ECalObjToFill->AverageHitTime = ECalIsoObjToFill.AverageHitTime;
                    ECalObjToFill->AverageZPosition = ECalIsoObjToFill.AverageZPosition;
                    ECalObjToFill->EMEnergyFit_Result =
                    ECalIsoObjToFill.EMEnergyFit_Result;
                    ECalObjToFill->EMEnergyFit_Uncertainty =
                    ECalIsoObjToFill.EMEnergyFit_Uncertainty;
                    ECalObjToFill->FilledAsType = 0;
                    ECalObjToFill->PID_Circularity = ECalIsoObjToFill.PID_Circularity;
                    ECalObjToFill->PID_Angle = ECalIsoObjToFill.PID_Angle;
                    ECalObjToFill->PID_ShowerAngle = ECalIsoObjToFill.PID_ShowerAngle;
                    ECalObjToFill->PID_Asymmetry = ECalIsoObjToFill.PID_Asymmetry;
                    ECalObjToFill->PID_TruncatedMaxRatio =
                    ECalIsoObjToFill.PID_TruncatedMaxRatio;
                    ECalObjToFill->PID_TransverseChargeRatio =
                    ECalIsoObjToFill.PID_TransverseChargeRatio;
                    ECalObjToFill->PID_FrontBackRatio =
                    ECalIsoObjToFill.PID_FrontBackRatio;
                    ECalObjToFill->Containment = ECalIsoObjToFill.Containment;
                    
                    ECalObjToFill->Cone = ECalIsoObjToFill.Shower.ConeAngle;
                    ECalObjToFill->EDeposit = ECalIsoObjToFill.Track.EDeposit;
                    
                    ECalObjToFill->Thrust = ECalIsoObjToFill.Thrust;
                    ECalObjToFill->ThrustOrigin = ECalIsoObjToFill.ThrustOrigin;
                    ECalObjToFill->ThrustAxis = ECalIsoObjToFill.ThrustAxis;
                    ECalObjToFill->InnerMostLayerHit =
                    ECalIsoObjToFill.MostUpStreamLayerHit;
                    ECalObjToFill->OuterMostLayerHit =
                    ECalIsoObjToFill.MostDownStreamLayerHit;
                    
                    ECalObjToFill->Width = ECalIsoObjToFill.Track.Width;
                    // The length has already been set, but it is the base object variable 
                    ECalObjToFill->Length = ECalObjToFill->TSubBaseObject::Length;

                    ND280NamedVerbose(
                        "TGlobalReconModule",
                        "Added information about an ECal Track PID constituent.");
                }  // End fill as track
                else if (shower_constit) {
                    if (!shower_constit->GetEDeposit()) {
                        ND280NamedWarn("TGlobalReconModule",
                            "ECal PID-constituent shower deposited no energy.");
                        continue;
                    }
                    
                    ND::TTrackerECALReconModule::TECALReconObjectFactory(
                        shower_constit, ecalObject, false, &ECalIsoObjToFill);
                    ECalObjToFill->LLR_Quality = ECalIsoObjToFill.PID_LLR_Quality;
                    ECalObjToFill->LLR_MIP_EM = ECalIsoObjToFill.PID_LLR_MIP_EM;
                    ECalObjToFill->LLR_MIP_Pion = ECalIsoObjToFill.PID_LLR_MIP_Pion;
                    ECalObjToFill->LLR_EM_HIP = ECalIsoObjToFill.PID_LLR_EM_HIP;
                    ECalObjToFill->LLR_MIP_EM_LowMomentum =
                    ECalIsoObjToFill.PID_LLR_MIP_EM_LowMomentum;
                    ECalObjToFill->AverageHitTime = ECalIsoObjToFill.AverageHitTime;
                    ECalObjToFill->EMEnergyFit_Result =
                    ECalIsoObjToFill.EMEnergyFit_Result;
                    ECalObjToFill->EMEnergyFit_Uncertainty =
                    ECalIsoObjToFill.EMEnergyFit_Uncertainty;
                    ECalObjToFill->FilledAsType = 1;
                    ECalObjToFill->PID_Circularity = ECalIsoObjToFill.PID_Circularity;
                    ECalObjToFill->PID_Angle = ECalIsoObjToFill.PID_Angle;
                    ECalObjToFill->PID_ShowerAngle = ECalIsoObjToFill.PID_ShowerAngle;
                    ECalObjToFill->PID_Asymmetry = ECalIsoObjToFill.PID_Asymmetry;
                    ECalObjToFill->PID_TruncatedMaxRatio =
                    ECalIsoObjToFill.PID_TruncatedMaxRatio;
                    ECalObjToFill->PID_TransverseChargeRatio =
                    ECalIsoObjToFill.PID_TransverseChargeRatio;
                    ECalObjToFill->PID_FrontBackRatio =
                    ECalIsoObjToFill.PID_FrontBackRatio;
                    ECalObjToFill->Containment = ECalIsoObjToFill.Containment;
                    
                    ECalObjToFill->Cone = ECalIsoObjToFill.Shower.ConeAngle;
                    ECalObjToFill->Length = ECalIsoObjToFill.ObjectLength;
                    ECalObjToFill->EDeposit = ECalIsoObjToFill.Shower.EDeposit;
                    
                    ECalObjToFill->Thrust = ECalIsoObjToFill.Thrust;
                    ECalObjToFill->ThrustOrigin = ECalIsoObjToFill.ThrustOrigin;
                    ECalObjToFill->ThrustAxis = ECalIsoObjToFill.ThrustAxis;
                    ECalObjToFill->InnerMostLayerHit =
                    ECalIsoObjToFill.MostUpStreamLayerHit;
                    ECalObjToFill->OuterMostLayerHit =
                    ECalIsoObjToFill.MostDownStreamLayerHit;
                    
                    Double_t widthHolder =
                    ECalObjToFill->Length * tan(ECalObjToFill->PID_ShowerAngle);
                    ECalObjToFill->Width =
                    TVector3(widthHolder, widthHolder, widthHolder);
                    
                    ND280NamedVerbose(
                        "TGlobalReconModule",
                        "Added information about an ECal Shower PID constituent.");
                }  // End fill as shower
                else {
                    ND280NamedWarn(
                        "TGlobalReconModule",
                        "When transposing an ECal object from: "
                        << ecalObject->ConvertDetector()
                        << " found a constituent of type: "
                        << (*ecalObject->GetConstituents()->begin())->ClassName()
                        << ", expected TReconTrack or TReconShower.");
                }
            }  // End Transpose P0DECal and non ECal-Vtx tracks
            
            else {  // If we are a Tr/Ds object and we are using ECal vertexing.
                
                ND::TTrackerECALReconModule::TECALReconVertexTrack
                ECalIsoVtxTrackToFill;
                ND::TTrackerECALReconModule::TECALReconVertexTrackFactory(
                    ecalObject, false, &ECalIsoVtxTrackToFill,
                    ND::THandle<ND::TG4TrajectoryContainer>(NULL), false);
                
                // Combine these two to build MIP_HIP
                ECalObjToFill->LLR_MIP_EM = ECalIsoVtxTrackToFill.PID_LLR_MIP_EM;
                ECalObjToFill->LLR_EM_HIP = ECalIsoVtxTrackToFill.PID_LLR_EM_HIP;
                
                ECalObjToFill->AverageHitTime = ECalIsoVtxTrackToFill.AverageHitTime;
                
                ECalObjToFill->InnerMostLayerHit =
                ECalIsoVtxTrackToFill.FrontLayerNumber;
                ECalObjToFill->OuterMostLayerHit =
                ECalIsoVtxTrackToFill.BackLayerNumber;
                
                ECalObjToFill->EDeposit = ECalIsoVtxTrackToFill.TotalCharge;
                
                ECalObjToFill->LLR_MIP_HIP_VACut =
                ECalIsoVtxTrackToFill.PID_LLR_MIP_HIP_VACut;
                
                ECalObjToFill->Containment = ECalIsoVtxTrackToFill.Containment;
                
                ECalObjToFill->FilledAsType = 2;
                
                ND280NamedVerbose(
                    "TGlobalReconModule",
                    "Added information about an ECal Vertex Track constituent.");
                // Most non-standard members are unfilled.
            }
        }  // End loop over ECal objects in module
    }    // End loop over ECal modules
}

//*************************************************************
void ND::TGlobalReconModule::FillTrackerInfo(ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject) 
{
//*************************************************************
    
    // We created dummy objects for the 18 possible Tracker combinations:
    
    // No detector missing
    // TPC1-FGD1
    // TPC1-FGD1-TPC2
    // TPC1-FGD1-TPC2-FGD2
    // TPC1-FGD1-TPC2-FGD2-TPC3
    // FGD1-TPC2
    // FGD1-TPC2-FGD2
    // FGD1-TPC2-FGD2-TPC3
    // TPC2-FGD2
    // TPC2-FGD2-TPC3
    // FGD2-TPC3
    
    ND::TReconTrack dummy[18];
    dummy[0].AddDetector(ND::TReconBase::kTPC1);
    dummy[0].AddDetector(ND::TReconBase::kFGD1);
    
    dummy[1].AddDetector(ND::TReconBase::kTPC1);
    dummy[1].AddDetector(ND::TReconBase::kFGD1);
    dummy[1].AddDetector(ND::TReconBase::kTPC2);
    
    dummy[2].AddDetector(ND::TReconBase::kTPC1);
    dummy[2].AddDetector(ND::TReconBase::kFGD1);
    dummy[2].AddDetector(ND::TReconBase::kTPC2);
    dummy[2].AddDetector(ND::TReconBase::kFGD2);
    
    dummy[3].AddDetector(ND::TReconBase::kTPC1);
    dummy[3].AddDetector(ND::TReconBase::kFGD1);
    dummy[3].AddDetector(ND::TReconBase::kTPC2);
    dummy[3].AddDetector(ND::TReconBase::kFGD2);
    dummy[3].AddDetector(ND::TReconBase::kTPC3);
    
    dummy[4].AddDetector(ND::TReconBase::kFGD1);
    dummy[4].AddDetector(ND::TReconBase::kTPC2);
    
    dummy[5].AddDetector(ND::TReconBase::kFGD1);
    dummy[5].AddDetector(ND::TReconBase::kTPC2);
    dummy[5].AddDetector(ND::TReconBase::kFGD2);
    
    dummy[6].AddDetector(ND::TReconBase::kFGD1);
    dummy[6].AddDetector(ND::TReconBase::kTPC2);
    dummy[6].AddDetector(ND::TReconBase::kFGD2);
    dummy[6].AddDetector(ND::TReconBase::kTPC3);
    
    dummy[7].AddDetector(ND::TReconBase::kTPC2);
    dummy[7].AddDetector(ND::TReconBase::kFGD2);
    
    dummy[8].AddDetector(ND::TReconBase::kTPC2);
    dummy[8].AddDetector(ND::TReconBase::kFGD2);
    dummy[8].AddDetector(ND::TReconBase::kTPC3);
    
    dummy[9].AddDetector(ND::TReconBase::kFGD2);
    dummy[9].AddDetector(ND::TReconBase::kTPC3);
    
    // one FGD missing
    // TPC1-FGD1-TPC2-TPC3
    // TPC1-FGD2-TPC2-TPC3
    // TPC1-TPC2
    // TPC2-TPC3
    
    dummy[10].AddDetector(ND::TReconBase::kTPC1);
    dummy[10].AddDetector(ND::TReconBase::kFGD1);
    dummy[10].AddDetector(ND::TReconBase::kTPC2);
    dummy[10].AddDetector(ND::TReconBase::kTPC3);
    
    dummy[11].AddDetector(ND::TReconBase::kTPC1);
    dummy[11].AddDetector(ND::TReconBase::kTPC2);
    dummy[11].AddDetector(ND::TReconBase::kFGD2);
    dummy[11].AddDetector(ND::TReconBase::kTPC3);
    
    dummy[12].AddDetector(ND::TReconBase::kTPC1);
    dummy[12].AddDetector(ND::TReconBase::kTPC2);
    
    dummy[13].AddDetector(ND::TReconBase::kTPC2);
    dummy[13].AddDetector(ND::TReconBase::kTPC3);
    
    // two FGDs missing
    // TPC1-TPC2-TPC3
    
    dummy[14].AddDetector(ND::TReconBase::kTPC1);
    dummy[14].AddDetector(ND::TReconBase::kTPC2);
    dummy[14].AddDetector(ND::TReconBase::kTPC3);
    
    // one TPC missing
    // TPC1-FGD1-TPC3
    // TPC1-FGD2-TPC3
    
    dummy[15].AddDetector(ND::TReconBase::kTPC1);
    dummy[15].AddDetector(ND::TReconBase::kFGD1);
    dummy[15].AddDetector(ND::TReconBase::kTPC3);
    
    dummy[16].AddDetector(ND::TReconBase::kTPC1);
    dummy[16].AddDetector(ND::TReconBase::kFGD2);
    dummy[16].AddDetector(ND::TReconBase::kTPC3);
    
    // two FGDs and one TPC missing
    // TPC1-TPC3
    
    dummy[17].AddDetector(ND::TReconBase::kTPC1);
    dummy[17].AddDetector(ND::TReconBase::kTPC3);
    
    for (int i = 0; i < 18; i++) {
        if (!object->UsesDetector(dummy[i].GetDetectors())) continue;
        
        ND::THandle<ND::TReconBase> trackerObject =
        ReconObjectUtils::GetConstituentInDetector(object,
            dummy[i].GetDetectors());
        
        if (trackerObject && globalObject.NTRACKERs < NMAXTRACKER) {
            TTrackerObject* tracker =
            new ((*(globalObject.TRACKER))[globalObject.NTRACKERs++])
            TTrackerObject;
            
            // fill general information 
            int  det = GetTrackerDetectorNumber(trackerObject);
            FillSubBaseObject(trackerObject, *tracker, det);
            
            // Fill charge for object.
            tracker->Charge = TrackingUtils::GetCharge(trackerObject);
            CheckMatchingFailure(tracker, trackerObject);
        }
    }
}

//*****************************************************
int ND::TGlobalReconModule::GetTrackerDetectorNumber(
    ND::THandle<ND::TReconBase> object)
{
//*****************************************************
    
    int det = 0;
    
    if (object->UsesDetector(ND::TReconBase::kTPC1)) det = 1;
    if (object->UsesDetector(ND::TReconBase::kTPC2)) det = det * 10 + 2;
    if (object->UsesDetector(ND::TReconBase::kTPC3)) det = det * 10 + 3;
    if (object->UsesDetector(ND::TReconBase::kFGD1)) det = det * 10 + 4;
    if (object->UsesDetector(ND::TReconBase::kFGD2)) det = det * 10 + 5;
    if (object->UsesDetector(ND::TReconBase::kSFG))  det = det * 10 + 6;
    // What about the HATs?????  (AJF 19/3/24) .. added these 3 
    if (object->UsesDetector(ND::TReconBase::kTopHAT))  det = det * 10 + 7;
    if (object->UsesDetector(ND::TReconBase::kBottomHAT))  det = det * 10 + 8;
    if (object->UsesDetector(ND::TReconBase::kTOF))  det = det * 10 + 9;

    return det;
}

//*************************************************************
void ND::TGlobalReconModule::FillTPCInfo(ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject) 
{
    //*************************************************************
    
    ND::TReconBase::Status dets[3] = {
    ND::TReconBase::kTPC1, ND::TReconBase::kTPC2, ND::TReconBase::kTPC3};
    
    for (int i = 0; i < 3; i++) {
        if (!object->UsesDetector(dets[i])) continue;
        
        // account for possible >1 segments in one TPC
        ND::TReconObjectContainer tpcObjects =
        ReconObjectUtils::GetAllConstituentsInDetector(object, dets[i]);
        
        if (tpcObjects.size() == 0) continue;
        
        for (ND::TReconObjectContainer::iterator iter = tpcObjects.begin();
            iter != tpcObjects.end(); iter++) {
        ND::THandle<ND::TReconBase> tpcObject = *iter;
        
        if (tpcObject && globalObject.NTPCs < NMAXTPC) {
            TTPCObject* tpc =
            new ((*(globalObject.TPC))[globalObject.NTPCs++]) TTPCObject;
            FillTPCObject(tpcObject, *tpc, i + 1);
        }
            }
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillTPCObject(
    ND::THandle<ND::TReconBase> tpcObject, TTPCObject& tpc, int det) 

{
    //*************************************************************
    
    // fill general information
    FillSubBaseObject(tpcObject, tpc, det);
    
    tpc.Charge = TrackingUtils::GetCharge(tpcObject);
    
    tpc.NTrun = 0;
    tpc.Ccorr = DEFAULT_MAX;
    tpc.PullEle = DEFAULT_MAX;
    tpc.PullMuon = DEFAULT_MAX;
    tpc.PullPion = DEFAULT_MAX;
    tpc.PullKaon = DEFAULT_MAX;
    tpc.PullProton = DEFAULT_MAX;
    
    tpc.SigmaEle = DEFAULT_MAX;
    tpc.SigmaMuon = DEFAULT_MAX;
    tpc.SigmaPion = DEFAULT_MAX;
    tpc.SigmaKaon = DEFAULT_MAX;
    tpc.SigmaProton = DEFAULT_MAX;
    
    tpc.dEdxexpEle = DEFAULT_MAX;
    tpc.dEdxexpMuon = DEFAULT_MAX;
    tpc.dEdxexpPion = DEFAULT_MAX;
    tpc.dEdxexpKaon = DEFAULT_MAX;
    tpc.dEdxexpProton = DEFAULT_MAX;
    
    tpc.PDist = DEFAULT_MAX;
    tpc.PEField = DEFAULT_MAX;
    
    tpc.TPCLength = DEFAULT_MAX;
    tpc.TrackType = DEFAULT_MAX;
    
    // True hit segment information
    int n_hitseg = 0;
    double edepo_sum = 0;
    double edepo_sec_sum = 0;
    double dedx_ave = 0;
    double dedx_sec_ave = 0;
    
    ND::THandle<ND::THitSelection> tpc_hits = tpcObject->GetHits();
    if(tpc_hits){
      // Loop over hits
      for(ND::THitSelection::iterator hit = tpc_hits->begin();
          hit != tpc_hits->end(); hit++){
      
        // Get true contributed G4VHits
        std::vector<ND::TG4VHit*> g4vhit = HitTruthInfo::GetHitTruthInfo(*hit);
        
        // Loop over G4VHits
        for(std::vector<ND::TG4VHit*>::iterator g4Hit = g4vhit.begin();
            g4Hit != g4vhit.end(); g4Hit++){
            
          // Cast the G4VHit into hit segment
          ND::TG4HitSegment* G4HitSeg = dynamic_cast<ND::TG4HitSegment*>(*g4Hit);
          if(G4HitSeg == NULL) continue;
          
          // Get the energy deposition of this hit segment
          double edepo_temp = G4HitSeg->GetEnergyDeposit();
          double edepo_sec_temp = G4HitSeg->GetSecondaryDeposit();
          
          // Get the "length" of this hit segment
          //TVector3 pos_temp;
          //pos_temp.SetXYZ(G4HitSeg->GetStopX() - G4HitSeg->GetStartX(),
          //                G4HitSeg->GetStopY() - G4HitSeg->GetStartY(),
          //                G4HitSeg->GetStopZ() - G4HitSeg->GetStartZ());
          //double len_temp = pos_temp.Mag();
          double len_temp = G4HitSeg->GetTrackLength();
          
          // Record information
          n_hitseg += 1;
          edepo_sum += edepo_temp;
          edepo_sec_sum += edepo_sec_temp;
          if(len_temp > 0){
            dedx_ave += (edepo_temp / len_temp);
            dedx_sec_ave += (edepo_sec_temp / len_temp);
          }
          
        } // End of G4VHit loop
      
      } // End of hit loop
    }

    tpc.TrueEDeposit = edepo_sum;
    tpc.TrueEDeposit_sec = edepo_sec_sum;
    if(n_hitseg > 0){
      tpc.TruedEdx = dedx_ave / n_hitseg;
      tpc.TruedEdx_sec = dedx_sec_ave / n_hitseg;
    }
    else{
      tpc.TruedEdx = 0;
      tpc.TruedEdx_sec = 0;
    }
       
    // Fill PID information when exists
    ND::THandle<ND::TReconPID> pid = tpcObject;
    
    if (pid) {
                        
        try {
            GetTPCDatum<ND::TRealDatum, int>(pid, "NVerRows", tpc.NVerRows);
        } catch (ND::EBadConversion& ex) {
            GetTPCDatum<ND::TIntegerDatum, int>(pid, "NVerRows", tpc.NVerRows);
        }
        
        try {
            GetTPCDatum<ND::TRealDatum, int>(pid, "NHorRows", tpc.NVerRows);
        } catch (ND::EBadConversion& ex) {
            GetTPCDatum<ND::TIntegerDatum, int>(pid, "NHorRows", tpc.NVerRows);
        }
        
        try {
            GetTPCDatum<ND::TRealDatum, int>(pid, "T0Source", tpc.T0Source);
        } catch (ND::EBadConversion& ex) {
            GetTPCDatum<ND::TIntegerDatum, int>(pid, "T0Source", tpc.T0Source);
        }
        
        // Changed from tpcPid_PullEle to NTrun to fix bug 569 - no reason why
        // this
        // should work,
        // but it did, so the bug still needs investigating.
        ND::THandle<ND::TRealDatum> datum_ele =
        (pid)->Get<ND::TRealDatum>("tpcPid_NTrun");
        if (datum_ele) {
            tpc.NTrun = pid->Get<ND::TRealDatum>("tpcPid_NTrun")->GetValue();
            tpc.Craw = pid->Get<ND::TRealDatum>("tpcPid_Craw")->GetValue();
            tpc.Ccorr = pid->Get<ND::TRealDatum>("tpcPid_Ccorr")->GetValue();
            tpc.PullEle = pid->Get<ND::TRealDatum>("tpcPid_PullEle")->GetValue();
            tpc.PullMuon = pid->Get<ND::TRealDatum>("tpcPid_PullMuon")->GetValue();
            tpc.PullPion = pid->Get<ND::TRealDatum>("tpcPid_PullPion")->GetValue();
            tpc.PullKaon = pid->Get<ND::TRealDatum>("tpcPid_PullKaon")->GetValue();
            tpc.PullProton =
            pid->Get<ND::TRealDatum>("tpcPid_PullProton")->GetValue();
            tpc.dEdxexpEle =
            pid->Get<ND::TRealDatum>("tpcPid_dEdxexpEle")->GetValue();
            tpc.dEdxexpMuon =
            pid->Get<ND::TRealDatum>("tpcPid_dEdxexpMuon")->GetValue();
            tpc.dEdxexpPion =
            pid->Get<ND::TRealDatum>("tpcPid_dEdxexpPion")->GetValue();
            tpc.dEdxexpKaon =
            pid->Get<ND::TRealDatum>("tpcPid_dEdxexpKaon")->GetValue();
            tpc.dEdxexpProton =
            pid->Get<ND::TRealDatum>("tpcPid_dEdxexpProton")->GetValue();
            tpc.SigmaEle = pid->Get<ND::TRealDatum>("tpcPid_SigmaEle")->GetValue();
            tpc.SigmaMuon = pid->Get<ND::TRealDatum>("tpcPid_SigmaMuon")->GetValue();
            tpc.SigmaPion = pid->Get<ND::TRealDatum>("tpcPid_SigmaPion")->GetValue();
            tpc.SigmaKaon = pid->Get<ND::TRealDatum>("tpcPid_SigmaKaon")->GetValue();
            tpc.SigmaProton =
            pid->Get<ND::TRealDatum>("tpcPid_SigmaProton")->GetValue();
        }
        
        if (pid->Get<ND::TIntegerDatum>("TrackType")) {
            tpc.TrackType = pid->Get<ND::TIntegerDatum>("TrackType")->GetValue();
        }
        ND::tman("globalAnalysis").GetTPCLength(tpcObject, tpc.TPCLength);
    }  // pid
    
    // Get the number of fitted vertical and horizontal cluster, relevant from
    // TREx era
    tpc.NbFittedVerticalClusters = 0;
    tpc.NbFittedHorizontalClusters = 0;
    if (tpcObject->Get<ND::TIntegerDatum>("NbFittedVerticalClusters"))
        tpc.NbFittedVerticalClusters =
    tpcObject->Get<ND::TIntegerDatum>("NbFittedVerticalClusters")
    ->GetValue();
    
    if (tpcObject->Get<ND::TIntegerDatum>("NbFittedHorizontalClusters"))
        tpc.NbFittedHorizontalClusters =
    tpcObject->Get<ND::TIntegerDatum>("NbFittedHorizontalClusters")
    ->GetValue();
    
    std::vector<double> vect_tmp;
    
    if (GetTPCDatum<ND::TRealDatum, double>(tpcObject, "RefitDistortionCorrMom",
        vect_tmp) &&
        vect_tmp.size() == 1)
    tpc.PDist = vect_tmp[0];
    
    if (GetTPCDatum<ND::TRealDatum, double>(tpcObject, "RefitEFieldDistortionMom",
        vect_tmp) &&
        vect_tmp.size() == 1)
    tpc.PEField = vect_tmp[0];
}

//*************************************************************
template <class T, typename I>
bool ND::TGlobalReconModule::GetTPCDatum(ND::THandle<ND::TReconBase> object,
    const std::string& datum,
    std::vector<I>& vvector) const 
{
    //*************************************************************
    
    // Clear the vector
    vvector.clear();
    
    bool ok = false;
    if (!object) return false;
    
    // Check whether the object has the datum
    if (object->Get<T>(datum.c_str())) {
        vvector.push_back(object->Get<T>(datum.c_str())->GetValue());
        ok = true;
    } else {
        // Try to get first layer constituents
        ND::THandle<ND::TReconObjectContainer> constituents1 =
        object->GetConstituents();
        
        if (!constituents1) return false;
        
        // The datum is expected either in each constituent itself or in its first
        // contributor
        for (ND::TReconObjectContainer::iterator iter = constituents1->begin();
            iter != constituents1->end(); iter++) {
            ND::THandle<ND::TReconBase> object_tmp1 = *iter;
            if (!object_tmp1) continue;
            
            if (object_tmp1->Get<T>(datum.c_str())) {
                vvector.push_back(object_tmp1->Get<T>(datum.c_str())->GetValue());
                ok = true;
                continue;
            }
            
            ND::THandle<ND::TReconObjectContainer> constituents2 =
            object_tmp1->GetConstituents();
            if (!constituents2 || constituents2->size() < 1) continue;
            
            ND::THandle<ND::TReconBase> object_tmp2 = *(constituents2->begin());
            if (!object_tmp2) continue;
            
            if (object_tmp2->Get<T>(datum.c_str())) {
                vvector.push_back(object_tmp2->Get<T>(datum.c_str())->GetValue());
                ok = true;
                continue;
            }
        }
    }
    
    return ok;
}

//*************************************************************
void ND::TGlobalReconModule::FillSFGInfo(
    ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject) 
{
//*************************************************************
  
  ND::TReconBase::Status det  = ND::TReconBase::kSFG;
  
  do {
    if(!object->UsesDetector(det)) break;
    
    ND::TReconObjectContainer sfgObjects;
    
    // for curving-back tracks allow using two constituents in one SFG
    if (TrackingUtils::GetCurvBackTrack(*object)) {
      sfgObjects =
        ReconObjectUtils::GetAllConstituentsInDetector(object, det);
    } else {
      //!!! WHICH SFG CONSTITUENT ONE WANTS TO USE WHEN THERE ARE TO OR MORE
      ND::THandle<ND::TReconBase> sfgConst =
        ReconObjectUtils::GetConstituentInDetector(object, det);
      if (sfgConst) sfgObjects.push_back(sfgConst);
    }
    
    if (sfgObjects.size() == 0) break;
    
    for (ND::TReconObjectContainer::iterator iter = sfgObjects.begin();
        iter != sfgObjects.end(); iter++) {
      ND::THandle<ND::TReconBase> sfgObject = *iter;
      std::map<ND::THandle<THit>,int> hitIndices;
      
      if (sfgObject && globalObject.NSFGs < NMAXSFG) {
        TSFGObject* sfg =
          new ((*(globalObject.SFG))[globalObject.NSFGs++]) TSFGObject;
        // Fill general information 
        // (for SFG, detector index should always be = 1)
        FillSubBaseObject(sfgObject, *sfg, 1 );
        
        // Fill information of individual SFG hits
        sfg->NHits  = 0;
        ND::THandle<ND::THitSelection> tbHits = sfgObject->GetHits();
        if (tbHits) {
          int ihit = 0;
          for (ND::THitSelection::iterator hit = tbHits->begin();
               hit != tbHits->end(); hit++) {
            // Fill index map (will be used for filling the node hits)
            hitIndices.insert(std::make_pair(*hit, ihit));
            ++ihit;
            // Fill the information about individual SFG hit
            TVector3 pos = (*hit)->GetPosition();
            TLorentzVector hitpos = TLorentzVector(
              pos[0],pos[1],pos[2],(*hit)->GetTime());
            sfg->HitPosition.push_back(hitpos);
            sfg->HitCharge  .push_back((*hit)->GetCharge());
            ND::THitSelection thisHit;
            thisHit.push_back(*hit);
            sfg->HitTrajID.push_back(
              TrackTruthInfo::GetG4TrajIDHits(thisHit));
          }  // End loop over hits
          sfg->NHits = ihit;
        } // Finished filling information on individual SFG hits

        double edep_sum = 0;
        // Fill information on individual SFG nodes
        // (probably only used by tracks and particles, but
        // technically valid for all TReconBase objects)
        sfg->NNodes = 0;
        const ND::TReconNodeContainer& nodeCont = sfgObject->GetNodes();
        ND::TReconNodeContainer::const_iterator it_node;
        for (it_node = nodeCont.begin();
             it_node != nodeCont.end(); ++it_node) {
          if(!(*it_node)) continue;
          // Variables to fill node information
          TLorentzVector   pos    = TLorentzVector(-999,-999,-999,-999);
          TLorentzVector   posvar = TLorentzVector(-999,-999,-999,-999);
          TVector3         dir    = TVector3(-999,-999,-999);
          TVector3         dirvar = TVector3(-999,-999,-999);
          float            edep   = 0;
          double edep_fit = 0;
          int              nhits  = 0;
          std::vector<int> nodeHits;
          // Start
          const ND::THandle<ND::TReconNode> node = *it_node;
          ND::THandle<ND::TReconState> s1 = node->GetState();
          if (s1) {
            pos    = ND::TReconStateUtilities::GetPosition(s1);
            posvar = ND::TReconStateUtilities::GetPositionVariance(s1);
            dir    = ND::TPIDState::GetStateDirection(s1);
            dirvar = ND::TPIDState::GetStateDirectionVariance(s1);
          }
          ND::THandle<ND::TTrackState> s2 = node->GetState();
          if(s2){
            edep_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) {
              edep =   ND::TPIDState::GetStateEDeposit(objState);
              edep_sum += ND::TPIDState::GetStateEDeposit(objState);
            }
            // Node hits
            const ND::THandle<ND::THitSelection> 
              nodeObjectHits = nodeObject->GetHits();
            if (nodeObjectHits) {
              ND::THitSelection::const_iterator hit;
              for (hit = nodeObjectHits->begin();
                   hit != nodeObjectHits->end(); ++hit) {
                auto find = hitIndices.find(*hit);
                if(find != hitIndices.end()){
                    nodeHits.push_back(hitIndices[*hit]);
                    ++nhits;
                }
              }
            }
          }
          // Fill variables
          sfg->NodePosition   .push_back(pos);
          sfg->NodePosVariance.push_back(posvar);
          sfg->NodeDirection  .push_back(dir);
          sfg->NodeDirVariance.push_back(dirvar);
          sfg->NodeEDeposit   .push_back(edep);
          sfg->NodeEDeposit_fit.push_back(edep_fit);
          sfg->NodeNHits      .push_back(nhits);
          sfg->NodeHits       .push_back(nodeHits);
          ++sfg->NNodes;
        } // Finished filling information on individual SFG nodes
        // Fill length
        double length_temp = 0;
        it_node = nodeCont.begin();
        ND::THandle<ND::TReconState> nodeState = (*it_node)->GetState();
        TVector3 frontPos = ND::TReconStateUtilities::GetPosition(nodeState).Vect();
        ++it_node;
        while(it_node != nodeCont.end()){
          if(!(*it_node)){
            ++it_node;
            continue;
          }
          nodeState = (*it_node)->GetState();
          TVector3 currPos = ND::TReconStateUtilities::GetPosition(nodeState).Vect();
          length_temp += (currPos - frontPos).Mag();
          frontPos = currPos;
          ++it_node;
        }  
        sfg->Length = length_temp;
        // Fill energy deposit
        ND::THandle<ND::TReconPID>     sfgPID     = sfgObject;
        ND::THandle<ND::TReconTrack>   sfgTrack   = sfgObject;
        ND::THandle<ND::TReconCluster> sfgCluster = sfgObject;
        if(sfgPID){
          sfg->Type = 0; // particle
          // If the object is TReconPID, fill position/direction/EDeposit.
          ND::THandle<ND::TPIDState> pidState = sfgPID->GetState();
          sfg->Position      = pidState->GetPosition();
          sfg->PosVariance   = pidState->GetPositionVariance();
          sfg->Direction     = pidState->GetDirection();
          sfg->DirVariance   = pidState->GetDirectionVariance();
          sfg->EDeposit      = edep_sum; // Sum of all the node Edep
          // Get track properties
          // (Set value to -1 if not available)
          int contain = TrackingUtils::GetTRealDatumValue(*sfgObject, "Containment", -1);
          sfg->ProtonBraggCut = TrackingUtils::GetTRealDatumValue(*sfgObject, "ProtonBraggCut", -1);
          sfg->PionBraggCut = TrackingUtils::GetTRealDatumValue(*sfgObject, "PionBraggCut", -1);
          sfg->NegPionBraggCut = TrackingUtils::GetTRealDatumValue(*sfgObject, "NegPionBraggCut", -1);
          
          // Fill PID weights for each hypothesis
          std::vector<int> ParticleIds;
          std::vector<double> PIDWeights;
          // Loop over all PID alternates
          // For SFG all PID hypothesis are saved in alternates
          ND::TReconObjectContainer::const_iterator it;
          for(it = sfgPID->GetAlternates().begin(); 
              it != sfgPID->GetAlternates().end(); it++){
            ND::THandle<ND::TReconPID> alter = *it;
            if(contain == 1){ // SFG contained
              if(alter->GetParticleId() == ND::TReconPID::kProton){
                ParticleIds.push_back(2212);
                PIDWeights.push_back(alter->GetPIDWeight());
              }
              else if(alter->GetParticleId() == ND::TReconPID::kPion){
                ParticleIds.push_back(211);
                PIDWeights.push_back(alter->GetPIDWeight());
              }              
              else if(alter->GetParticleId() == ND::TReconPID::kMuon){
                ParticleIds.push_back(13);
                PIDWeights.push_back(alter->GetPIDWeight());
              }
              else if(alter->GetParticleId() == ND::TReconPID::kEM){
                ParticleIds.push_back(11);
                PIDWeights.push_back(alter->GetPIDWeight());
              }               
            }
            else{ // SFG escaped
              if(alter->GetParticleId() == ND::TReconPID::kHeavyTrack){
                ParticleIds.push_back(2212);
                PIDWeights.push_back(alter->GetPIDWeight());                
              }
              else if(alter->GetParticleId() == ND::TReconPID::kLightTrack){
                ParticleIds.push_back(211);
                PIDWeights.push_back(alter->GetPIDWeight());  
                ParticleIds.push_back(13);
                PIDWeights.push_back(alter->GetPIDWeight());               
              }
              else if(alter->GetParticleId() == ND::TReconPID::kEM){
                ParticleIds.push_back(11);
                PIDWeights.push_back(alter->GetPIDWeight());                
              }              
            }
          } // End of PID alternates loop                        
          // Save PID weights for different hypothesis
          for(unsigned int hypo = 0; hypo < ParticleIds.size(); hypo++){
            if(ParticleIds[hypo] == 2212)
              sfg->PIDWeightProt = PIDWeights[hypo];
            else if(ParticleIds[hypo] == 211)
              sfg->PIDWeightPion = PIDWeights[hypo];
            else if(ParticleIds[hypo] == 13)
              sfg->PIDWeightMuon = PIDWeights[hypo];
            else if(ParticleIds[hypo] == 11)
              sfg->PIDWeightElec = PIDWeights[hypo];
          }
        }else if(sfgTrack){
          sfg->Type = 1; // track
          // If the object is TReconTrack, fill position/direction/EDeposit.
          ND::THandle<ND::TTrackState> trackState = sfgTrack->GetState();
          sfg->Position      = trackState->GetPosition();
          sfg->PosVariance   = trackState->GetPositionVariance();
          sfg->Direction     = trackState->GetDirection();
          sfg->DirVariance   = trackState->GetDirectionVariance();
          sfg->EDeposit      = trackState->GetEDeposit();
        }else if(sfgCluster){
          sfg->Type = 2; // cluster
          // If the object is TReconCluster, fill position/direction/EDeposit.
          sfg->Position      = sfgCluster->GetPosition();
          sfg->PosVariance   = sfgCluster->GetPositionVariance();
          sfg->EDeposit      = sfgCluster->GetEDeposit();
        }
        // fill average hit time
        ND::THandle<ND::TRealDatum> sfgTimeRD =
        sfgObject->Get<ND::TRealDatum>("averageHitTime");
        if (sfgTimeRD) {
            sfg->avgtime = sfgTimeRD->GetValue();
        } else {
            sfg->avgtime = DEFAULT_MAX;
        }
      }
    }  // loop through SFG costituents
  } while (false);  // always stops...
}


//*************************************************************
void ND::TGlobalReconModule::FillSFGOther(ND::TND280Event& event) 
{
//*************************************************************
 
  // Get the main result to save.
  ND::THandle<ND::TAlgorithmResult> sfgRecon = event.GetFit("sfgRecon");
  if (!sfgRecon) return;
  // Get the results from reconstruction
  ND::THandle<ND::TReconObjectContainer> 
    sfgObjects = sfgRecon->GetResultsContainer("final");
  if (!sfgObjects) return;
  
  // Loop over all the objects and save clusters and vertices.
  ND::TReconObjectContainer::iterator it;
  for (it = sfgObjects->begin();
       it != sfgObjects->end(); it++){
    ND::THandle<ND::TReconBase> object = *it;
    if (!object) continue;
    ND::THandle<ND::TReconCluster> cluster = object;
    ND::THandle<ND::TReconVertex>  vertex  = object;
    if(!cluster && ! vertex) continue; 
    if(cluster){
      // Fill the cluster as GlobalPID      
      // Get the compact detector number and number of used detectors
      object->AddDetector(ND::TReconBase::kSFG);
      unsigned long dets = GetDetectorNumber(object);
      int      det_count = CountDetector(object);
      if (object->UsesDetector(ND::TReconBase::kSFG)){
        fNSFG++;
        if (det_count == 1) fNSFGIso++;
      }
      // Try to get the state of the object. 
      // If it dosen't exists skip the object.
      try 
      {
        ND::THandle<ND::TReconState> tstate = object->GetState();
      } catch (ND::EReconObject const& e)
      {
        std::cout << "TReconBase with no state. Skip object !!" << std::endl;
        return;
      }
      // Fill the map of indexes used for vertex link
      fGlobalIndexMap[object] = fNPIDs;
      // create a new TGlobalPID
      TGlobalPID* globalObject = new ((*fPIDs)[fNPIDs++]) TGlobalPID;
      // the unique ID of the global object
      globalObject->UniqueID = object->GetUniqueID();
      // Fill the vector of unique IDs of associated tracks
      if (fBrokenIndexMap.find(object) != fBrokenIndexMap.end())
        globalObject->BrokenUniqueIDs = fBrokenIndexMap[object];
      // Information about the detectors that used in this object
      globalObject->Detectors = dets;
      FillDetectorUsed(object, globalObject->DetectorUsed);
      globalObject->Charge   = TrackingUtils::GetCharge(object);
      globalObject->EDeposit = TrackingUtils::GetEDeposit(object);
      globalObject->Length   = ComputeTrackLength(object);
      globalObject->Status   = object->CheckStatus(object->kSuccess);
      globalObject->Chi2     = object->GetQuality();
      globalObject->NDOF     = object->GetNDOF();
      globalObject->NNodes   = 0;
      globalObject->AlgorithmName = object->GetAlgorithmName();
      globalObject->NHits = DEFAULT_MAX;
      if (object->GetHits()) 
        globalObject->NHits = object->GetHits()->size();
      globalObject->NConstituents = 0;
      //---------- Fill info about first and last hits
      if (object->GetHits()){
#ifdef FILL_ALL_GLOBAL_HITS
        FillAllHits(*(object->GetHits()), *globalObject);
#else
        FillFirstLastHits(*(object->GetHits()), *globalObject);
#endif
      }
      // Fill SFG Info
      TSFGObject* sfg =
        new ((*(globalObject->SFG))[globalObject->NSFGs++]) TSFGObject;
      sfg->Type = 2; // cluster
      // Fill information of individual SFG hits
      sfg->NHits = 0;
      ND::THandle<ND::THitSelection> tbHits = object->GetHits();
      if (tbHits) {
        for (ND::THitSelection::iterator hit = tbHits->begin();
             hit != tbHits->end(); hit++) {
          // Fill the information about individual SFG hit
          TVector3 pos = (*hit)->GetPosition();
          TLorentzVector hitpos = TLorentzVector(
            pos[0],pos[1],pos[2],(*hit)->GetTime());
          sfg->HitPosition.push_back(hitpos);
          sfg->HitCharge  .push_back((*hit)->GetCharge());
          ND::THitSelection thisHit;
          thisHit.push_back(*hit);
          sfg->HitTrajID.push_back(
            TrackTruthInfo::GetG4TrajIDHits(thisHit));
          ++sfg->NHits;
        }  // End loop over hits
      } // Finished filling information on individual SFG hits
      ND::THandle<ND::TReconCluster> sfgCluster = object;
      if(sfgCluster){
        sfg->Position      = sfgCluster->GetPosition();
        sfg->PosVariance   = sfgCluster->GetPositionVariance();
        sfg->EDeposit      = sfgCluster->GetEDeposit();
      }
      // fill general information 
      // (for SFG, detector index should always be = 1)
      FillSubBaseObject(object, *sfg, 1 );
      // fill average hit time
      ND::THandle<ND::TRealDatum> sfgTimeRD =
        object->Get<ND::TRealDatum>("averageHitTime");
      if (sfgTimeRD) {
        sfg->avgtime = sfgTimeRD->GetValue();
      } else {
        sfg->avgtime = DEFAULT_MAX;
      }
      // This is not the ideal way to save the cluster position, but 
      // is necessary to identify the bunch number.
      globalObject->FrontPosition = sfg->Position;
      globalObject->BackPosition  = sfg->Position;
    }else if(vertex){
      ////////////////////////////////////////////////////////////////
      // Check that the vertex has a state
      ND::THandle<ND::TVertexState> vstate = vertex->GetState();
      if (!vstate) continue;
      // Create a global vertex
      TGlobalVertex* gVertex = new ((*fVertices)[fNVertices++]) TGlobalVertex;
      gVertex->Type          = 1; // label as SFG vertex
      gVertex->NTrueVertices = 0;
      // Index of last added primary
      fPVInd                 = fNVertices - 1; 
      gVertex->PrimaryIndex  = -1;
      // Basic info
      gVertex->Position      = vstate->GetPosition();
      gVertex->Variance      = vstate->GetPositionVariance();
      gVertex->AlgorithmName = vertex->GetAlgorithmName();
      gVertex->Status        = vertex->GetStatus();
      gVertex->Quality       = vertex->GetQuality();
      gVertex->NDOF          = vertex->GetNDOF();
      // Associated tracks
      ND::THandle<ND::TReconObjectContainer> tracksAssociated 
        = vertex->GetConstituents();
      if (tracksAssociated->size() == 0) continue;
      std::map<ND::THandle<ND::TReconBase>,
      ND::THandle<ND::TReconObjectContainer>>::iterator cit;
      // Loop over associated tracks
      ND::TReconObjectContainer::iterator assocIter;
      for (assocIter = tracksAssociated->begin();
           assocIter != tracksAssociated->end(); assocIter++){
        // Skip vertex constituents
        if (ND::THandle<ND::TReconVertex> vtxDummy = (*assocIter)) continue;
        // Fill track info
        if (gVertex->NConstituents == NCONSTITUENTS) break;
        ND::THandle<ND::TReconBase> track = *assocIter;
        if (!track) continue;
        ND::THandle<ND::TReconState> frontState, backState;
        ND::tman().GetFirstAndLastStates(track, frontState, backState);
        if (!frontState || !backState) continue;
        TLorentzVector fpos = TrackingUtils::GetPosition(frontState);
        TLorentzVector bpos = TrackingUtils::GetPosition(backState);
        // Get index of the global track in the overall fPIDs list
        bool inMap = false;
        std::map<ND::THandle<ND::TReconBase>, int>::iterator mapIt;
        for (mapIt=fGlobalIndexMap.begin();mapIt!=fGlobalIndexMap.end();mapIt++){
          if (!mapIt->first) continue;
          // It should have SFG segment.
          ND::TReconBase::Status det = ND::TReconBase::kSFG;
          if(!((mapIt->first)->UsesDetector(det))) continue;
          ND::THandle<ND::TReconBase> sfgSeg =
            ReconObjectUtils::GetConstituentInDetector(mapIt->first, det);
          ND::THandle<ND::TReconState> frontTemp, backTemp;
          ND::tman().GetFirstAndLastStates(sfgSeg, frontTemp, backTemp);
          if (!frontTemp || !backTemp) continue;
          TLorentzVector fpos_tmp = TrackingUtils::GetPosition(frontTemp);
          TLorentzVector bpos_tmp = TrackingUtils::GetPosition(backTemp);
          // Compare the front and back position
          if((fpos.Vect()==fpos_tmp.Vect() && bpos.Vect()==bpos_tmp.Vect()) ||
             (fpos.Vect()==bpos_tmp.Vect() && bpos.Vect()==fpos_tmp.Vect())){
            inMap = true; break;
          }
        }
        if(!inMap) continue;
        // Set PID
        int GlobalPIDIndex = mapIt->second;
        // Fill constituent track info.
        TVertexConstituent* vConst =
          new ((*(gVertex->Constituents))[gVertex->NConstituents++])
          TVertexConstituent;
        vConst->Charge = (int)TrackingUtils::GetCharge(track);
        vConst->Quality = track->GetQuality();
        double   mom  = TrackingUtils::GetMomentum(frontState);
        TVector3 mdir = TrackingUtils::GetDirection(frontState);
        mdir *= mom;
        vConst->Momentum = mdir;
        vConst->PID = GlobalPIDIndex;
      } // loop over associated tracks
    }
  }
}

//*************************************************************
void ND::TGlobalReconModule::FillFGDInfo(ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject) 
{
    //*************************************************************
    
    // Get the tracker recon version of the this object, if it is a FGD
    // iso-recon
    // object.
    // The global recon version of the object lacks the PID information.
    ND::THandle<ND::TReconBase> newFgdIso =
    GetTrackerReconVersionOfFGDIsoTrack(object);
    if (newFgdIso) {
        object = newFgdIso;
    }
    
    ND::TReconBase::Status dets[2] = {ND::TReconBase::kFGD1,
    ND::TReconBase::kFGD2};
    
    for (int i = 0; i < 2; i++) {
        if (!object->UsesDetector(dets[i])) continue;
        
        ND::TReconObjectContainer fgdObjects;
        
        // for curving-back tracks allow using two constituents in one FGD
        if (TrackingUtils::GetCurvBackTrack(*object)) {
            fgdObjects =
            ReconObjectUtils::GetAllConstituentsInDetector(object, dets[i]);
        } else {
            //!!! WHICH FGD CONSTITUENT ONE WANTS TO USE FOR THE THROUGH(FGD)-GOING
            //! TRACKS?
            ND::THandle<ND::TReconBase> fgdConst =
            ReconObjectUtils::GetConstituentInDetector(object, dets[i]);
            if (fgdConst) fgdObjects.push_back(fgdConst);
        }
        
        if (fgdObjects.size() == 0) continue;
        
        for (ND::TReconObjectContainer::iterator iter = fgdObjects.begin();
            iter != fgdObjects.end(); iter++) {
        ND::THandle<ND::TReconBase> fgdObject = *iter;
        
        if (fgdObject && globalObject.NFGDs < NMAXFGD) {
            TFGDObject* fgd =
            new ((*(globalObject.FGD))[globalObject.NFGDs++]) TFGDObject;          
            
            // Fill information on individual FGD hits
            ND::THandle<ND::THitSelection> tbHits = fgdObject->GetHits();
            if (tbHits) {
                for (ND::THitSelection::iterator hit = tbHits->begin();
                    hit != tbHits->end(); hit++) {
                // Fill the information about individual FGD hit
                TVector3 pos = (*hit)->GetPosition();
                int layer = ND::TGeomInfo::FGD().ActivePlane(pos.Z());
                fgd->fgdhits_Layer.push_back(layer);
                fgd->fgdhits_Position.push_back(pos);
                fgd->fgdhits_Time.push_back((*hit)->GetTime());
                fgd->fgdhits_Charge.push_back((*hit)->GetCharge());
                ND::THitSelection thisHit;
                thisHit.push_back(*hit);
                fgd->fgdhits_TrajID.push_back(
                    TrackTruthInfo::GetG4TrajIDHits(thisHit));
                
                    }  // End loop over hits
            }
                // Finished filling information on individual FGD hits
            
            // Fill node information
            fgd->NodePosition.clear();
            fgd->NodeEDeposit.clear();
            
            const ND::TReconNodeContainer& nodeCont = fgdObject->GetNodes();
            ND::TReconNodeContainer::const_iterator it_node;
            for(it_node = nodeCont.begin(); it_node != nodeCont.end(); ++it_node){

              if(!(*it_node)) continue;
              
              // Variables to fill node information
              TLorentzVector pos_temp(-999,-999,-999,-999);
              double edepo_temp = -999;
              
              const ND::THandle<ND::TReconNode> node = *it_node;
              // Get the node state
              ND::THandle<ND::TReconState> node_state = node->GetState();
              if(node_state){
                // Retrieve the position
                pos_temp = ND::TReconStateUtilities::GetPosition(node_state);
              }
              // Get object from node
              ND::THandle<ND::TReconBase> node_object = node->GetObject();
              if(node_object){
                ND::THandle<ND::TReconState> object_state = node_object->GetState();
                if(object_state){
                  edepo_temp = ND::TPIDState::GetStateEDeposit(object_state);
                }
              }
              // (If energy deposition can't get from object->state, then from state)
              if(edepo_temp < 0){
                if(node_state){
                  edepo_temp = node_state->GetValue(node_state->GetEDepositIndex());
                }
              }
              
              fgd->NodePosition.push_back(pos_temp);
              fgd->NodeEDeposit.push_back(edepo_temp);

            } // End of node loop
        
            // fill general information
            FillSubBaseObject(fgdObject, *fgd, i + 1);
            
            // fill average hit time
            ND::THandle<ND::TRealDatum> fgdTimeRD =
            fgdObject->Get<ND::TRealDatum>("averageHitTime");
            if (fgdTimeRD) {
                fgd->avgtime = fgdTimeRD->GetValue();
            } else {
                fgd->avgtime = DEFAULT_MAX;
            }

            // Fill true energy deposit
            fgd->TrueE = -999.;
            ND::THandle<ND::TReconPID>   pid   = fgdObject;
            ND::THandle<ND::TReconTrack> track = fgdObject; 
            if (pid){  
                fgd->TrueE = pidUtils::GetTruthEDep(pid);
            }
            else if (track){
                fgd->TrueE = pidUtils::GetTruthEDep(track);  
            }
          
            
            // Fill charge per layer information
            for (int ilayer = 0; ilayer < 30; ilayer++) {
                fgd->chargePerLayer[ilayer] = 0.;
                fgd->chargePerLayerAttenCorr[ilayer] = 0.;
            }
            std::map<int, double> chargeMap;
            if (fgdUtils::GetChargeByLayers(fgdObject->GetHits(), chargeMap)) {
                for (std::map<int, double>::iterator it = chargeMap.begin();
                    it != chargeMap.end(); it++) {
                int layer = (*it).first;
                if (i == 1) {
                    layer -= 30;
                }
                fgd->chargePerLayer[layer] = (*it).second;
                    }
            }
            std::map<int, double> chargeMapAC;
            if (fgdUtils::GetFibreCorrectedChargeByLayersReconFit(
                fgdObject->GetHits(), fgdObject, chargeMapAC)) {
            for (std::map<int, double>::iterator it = chargeMapAC.begin();
                it != chargeMapAC.end(); it++) {
            int layer = (*it).first;
            if (i == 1) {
                layer -= 30;
            }
            fgd->chargePerLayerAttenCorr[layer] = (*it).second;
                }
                }
                // fill PID Information
                
                // fully contained info
                ND::THandle<ND::TIntegerDatum> RD_geomFC =
                fgdObject->Get<ND::TIntegerDatum>("fgdPidPar_geomFC");
                fgd->fgdContainment = 0;
                if (RD_geomFC) {
                    fgd->fgdContainment = RD_geomFC->GetValue();
                }
                
                // measured E and x
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_E =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_E");
                fgd->E = -1.;
                if (RD_fgdPidPar_E) {
                    fgd->E = RD_fgdPidPar_E->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_x =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_x");
                fgd->x = -1.;
                if (RD_fgdPidPar_x) {
                    fgd->x = RD_fgdPidPar_x->GetValue();
                }
                
                // expected E at the measured x
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Muon =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_E_exp_EvsXPull_Muon");
                fgd->E_exp_muon = -1.;
                if (RD_fgdPidPar_E_exp_EvsXPull_Muon) {
                    fgd->E_exp_muon = RD_fgdPidPar_E_exp_EvsXPull_Muon->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Pion =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_E_exp_EvsXPull_Pion");
                fgd->E_exp_pion = -1.;
                if (RD_fgdPidPar_E_exp_EvsXPull_Pion) {
                    fgd->E_exp_pion = RD_fgdPidPar_E_exp_EvsXPull_Pion->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Proton =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_E_exp_EvsXPull_Proton");
                fgd->E_exp_proton = -1.;
                if (RD_fgdPidPar_E_exp_EvsXPull_Proton) {
                    fgd->E_exp_proton = RD_fgdPidPar_E_exp_EvsXPull_Proton->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_E_exp_EvsXPull_Electron =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_E_exp_EvsXPull_Electron");
                fgd->E_exp_electron = -1.;
                if (RD_fgdPidPar_E_exp_EvsXPull_Electron) {
                    fgd->E_exp_electron =
                    RD_fgdPidPar_E_exp_EvsXPull_Electron->GetValue();
                }
                
                // sigma E at the measured x
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_error_EvsXPull_Muon =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_error_EvsXPull_Muon");
                fgd->sigmaE_muon = -1.;
                if (RD_fgdPidPar_error_EvsXPull_Muon) {
                    fgd->sigmaE_muon = RD_fgdPidPar_error_EvsXPull_Muon->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_error_EvsXPull_Pion =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_error_EvsXPull_Pion");
                fgd->sigmaE_pion = -1.;
                if (RD_fgdPidPar_error_EvsXPull_Pion) {
                    fgd->sigmaE_pion = RD_fgdPidPar_error_EvsXPull_Pion->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_error_EvsXPull_Proton =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_error_EvsXPull_Proton");
                fgd->sigmaE_proton = -1.;
                if (RD_fgdPidPar_error_EvsXPull_Proton) {
                    fgd->sigmaE_proton = RD_fgdPidPar_error_EvsXPull_Proton->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidPar_error_EvsXPull_Electron =
                fgdObject->Get<ND::TRealDatum>("fgdPidPar_error_EvsXPull_Electron");
                fgd->sigmaE_electron = -1.;
                if (RD_fgdPidPar_error_EvsXPull_Electron) {
                    fgd->sigmaE_electron =
                    RD_fgdPidPar_error_EvsXPull_Electron->GetValue();
                }
                
                // Pull values
                ND::THandle<ND::TRealDatum> RD_fgdPidWgt_EvsXPull_Muon =
                fgdObject->Get<ND::TRealDatum>("fgdPidWgt_EvsXPull_Muon");
                fgd->PullMuon = -1.e10;
                if (RD_fgdPidWgt_EvsXPull_Muon) {
                    fgd->PullMuon = RD_fgdPidWgt_EvsXPull_Muon->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidWgt_EvsXPull_Pion =
                fgdObject->Get<ND::TRealDatum>("fgdPidWgt_EvsXPull_Pion");
                fgd->PullPion = -1.e10;
                if (RD_fgdPidWgt_EvsXPull_Pion) {
                    fgd->PullPion = RD_fgdPidWgt_EvsXPull_Pion->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidWgt_EvsXPull_Proton =
                fgdObject->Get<ND::TRealDatum>("fgdPidWgt_EvsXPull_Proton");
                fgd->PullProton = -1.e10;
                if (RD_fgdPidWgt_EvsXPull_Proton) {
                    fgd->PullProton = RD_fgdPidWgt_EvsXPull_Proton->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidWgt_EvsXPull_Electron =
                fgdObject->Get<ND::TRealDatum>("fgdPidWgt_EvsXPull_Electron");
                fgd->PullElectron = -1.e10;
                if (RD_fgdPidWgt_EvsXPull_Electron) {
                    fgd->PullElectron = RD_fgdPidWgt_EvsXPull_Electron->GetValue();
                }
                ND::THandle<ND::TRealDatum> RD_fgdPidWgt_EvsXPull_NotSet =
                fgdObject->Get<ND::TRealDatum>("fgdPidWgt_EvsXPull_NotSet");
                fgd->PullNotSet = -1.e10;
                if (RD_fgdPidWgt_EvsXPull_NotSet) {
                    fgd->PullNotSet = RD_fgdPidWgt_EvsXPull_NotSet->GetValue();
                }
                
                int fgdNum;
                // FGD1 VA/EA variables
                // -------------------------
                if (i == 0) {
                    // Vertex activity values
                    int fgdNum = 2;
                    ND::THandle<ND::TRealDatum> RD_fgdVA_FgdN;
                    RD_fgdVA_FgdN = object->Get<ND::TRealDatum>("fgdVA_FgdN");
                    if (RD_fgdVA_FgdN) {
                        fgdNum = int(RD_fgdVA_FgdN->GetValue());
                    }
                    
                    fgd->isFgdVA = 0;
                    if (fgdNum == i) fgd->isFgdVA = 1;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_flag;
                    RD_fgdVA_flag = object->Get<ND::TRealDatum>("fgdVA_flag");
                    if (RD_fgdVA_flag && i == fgdNum) {
                        fgd->isFgdVA_flag = int(RD_fgdVA_flag->GetValue());
                    } else
                        fgd->isFgdVA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_globpassVA_flag;
                    RD_globpassVA_flag = object->Get<ND::TRealDatum>("globpassVA_flag");
                    if (RD_globpassVA_flag && i == fgdNum) {
                        fgd->globpassVA_flag = int(RD_globpassVA_flag->GetValue());
                    } else
                        fgd->globpassVA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_upMinZ;
                    RD_fgdVA_upMinZ = object->Get<ND::TRealDatum>("fgdVA_upMinZ");
                    if (RD_fgdVA_upMinZ && i == fgdNum) {
                        fgd->fgdVA_upMinZ = RD_fgdVA_upMinZ->GetValue();
                    } else
                        fgd->fgdVA_upMinZ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_downMaxZ;
                    RD_fgdVA_downMaxZ = object->Get<ND::TRealDatum>("fgdVA_downMaxZ");
                    if (RD_fgdVA_downMaxZ && i == fgdNum) {
                        fgd->fgdVA_downMaxZ = RD_fgdVA_downMaxZ->GetValue();
                    } else
                        fgd->fgdVA_downMaxZ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_upCosTheta;
                    RD_fgdVA_upCosTheta = object->Get<ND::TRealDatum>("fgdVA_upCosTheta");
                    if (RD_fgdVA_upCosTheta && i == fgdNum) {
                        fgd->fgdVA_upCosTheta = RD_fgdVA_upCosTheta->GetValue();
                    } else
                        fgd->fgdVA_upCosTheta = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_downCosTheta;
                    RD_fgdVA_downCosTheta =
                    object->Get<ND::TRealDatum>("fgdVA_downCosTheta");
                    if (RD_fgdVA_downCosTheta && i == fgdNum) {
                        fgd->fgdVA_downCosTheta = RD_fgdVA_downCosTheta->GetValue();
                    } else
                        fgd->fgdVA_downCosTheta = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_otherUpQ;
                    RD_fgdVA_otherUpQ = object->Get<ND::TRealDatum>("fgdVA_otherUpQ");
                    if (RD_fgdVA_otherUpQ && i == fgdNum) {
                        fgd->fgdVA_otherUpQ = RD_fgdVA_otherUpQ->GetValue();
                    } else
                        fgd->fgdVA_otherUpQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_otherDownQ;
                    RD_fgdVA_otherDownQ = object->Get<ND::TRealDatum>("fgdVA_otherDownQ");
                    if (RD_fgdVA_otherDownQ && i == fgdNum) {
                        fgd->fgdVA_otherDownQ = RD_fgdVA_otherDownQ->GetValue();
                    } else
                        fgd->fgdVA_otherDownQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verQ;
                    RD_fgdVA_verQ = object->Get<ND::TRealDatum>("fgdVA_verQ");
                    if (RD_fgdVA_verQ && i == fgdNum) {
                        fgd->fgdVA_verQ = RD_fgdVA_verQ->GetValue();
                    } else
                        fgd->fgdVA_verQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verLayQ;
                    RD_fgdVA_verLayQ = object->Get<ND::TRealDatum>("fgdVA_verLayQ");
                    if (RD_fgdVA_verLayQ && i == fgdNum) {
                        fgd->fgdVA_verLayQ = RD_fgdVA_verLayQ->GetValue();
                    } else
                        fgd->fgdVA_verLayQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNearQ;
                    RD_fgdVA_verNearQ = object->Get<ND::TRealDatum>("fgdVA_verNearQ");
                    if (RD_fgdVA_verNearQ && i == fgdNum) {
                        fgd->fgdVA_verNearQ = RD_fgdVA_verNearQ->GetValue();
                    } else
                        fgd->fgdVA_verNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNearQ;
                    RD_fgdVA_verNextNearQ =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNearQ");
                    if (RD_fgdVA_verNextNearQ && i == fgdNum) {
                        fgd->fgdVA_verNextNearQ = RD_fgdVA_verNextNearQ->GetValue();
                    } else
                        fgd->fgdVA_verNextNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNextNearQ;
                    RD_fgdVA_verNextNextNearQ =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNextNearQ");
                    if (RD_fgdVA_verNextNextNearQ && i == fgdNum) {
                        fgd->fgdVA_verNextNextNearQ = RD_fgdVA_verNextNextNearQ->GetValue();
                    } else
                        fgd->fgdVA_verNextNextNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_totalQ;
                    RD_fgdVA_totalQ = object->Get<ND::TRealDatum>("fgdVA_totalQ");
                    if (RD_fgdVA_totalQ && i == fgdNum) {
                        fgd->fgdVA_totalQ = RD_fgdVA_totalQ->GetValue();
                    } else
                        fgd->fgdVA_totalQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_totalCorrE;
                    RD_fgdVA_totalCorrE = object->Get<ND::TRealDatum>("fgdVA_totalCorrE");
                    if (RD_fgdVA_totalCorrE && i == fgdNum) {
                        fgd->fgdVA_totalCorrE = RD_fgdVA_totalCorrE->GetValue();
                    } else
                        fgd->fgdVA_totalCorrE = 0;
                    
                    // New Vertex Activity for FGD1
                    
                    // Vertex Activity in FGD1 2x3
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNearQ_rect;
                    RD_fgdVA_verNearQ_rect =
                    object->Get<ND::TRealDatum>("fgdVA_verNearQ_rect");
                    if (RD_fgdVA_verNearQ_rect && i == fgdNum) {
                        fgd->fgdVA_verNearQ_rect = RD_fgdVA_verNearQ_rect->GetValue();
                    } else
                        fgd->fgdVA_verNearQ_rect = 0;
                    
                    // Vertex activity in FGD1 2x5
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNearQ_rect;
                    RD_fgdVA_verNextNearQ_rect =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNearQ_rect");
                    if (RD_fgdVA_verNextNearQ_rect && i == fgdNum) {
                        fgd->fgdVA_verNextNearQ_rect =
                        RD_fgdVA_verNextNearQ_rect->GetValue();
                    } else
                        fgd->fgdVA_verNextNearQ_rect = 0;
                    
                    // Vertex activity in FGD1 2x7
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNextNearQ_rect;
                    RD_fgdVA_verNextNextNearQ_rect =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNextNearQ_rect");
                    if (RD_fgdVA_verNextNextNearQ_rect && i == fgdNum) {
                        fgd->fgdVA_verNextNextNearQ_rect =
                        RD_fgdVA_verNextNextNearQ_rect->GetValue();
                    } else
                        fgd->fgdVA_verNextNextNearQ_rect = 0;
                    
                    // End Activity values
                    // ----------------------------
                    fgdNum = 2;
                    ND::THandle<ND::TRealDatum> RD_fgdEA_FgdN;
                    RD_fgdEA_FgdN = object->Get<ND::TRealDatum>("fgdEA_FgdN");
                    if (RD_fgdEA_FgdN) {
                        fgdNum = int(RD_fgdEA_FgdN->GetValue());
                    }
                    
                    fgd->isFgdEA = 0;
                    if (fgdNum == i) fgd->isFgdEA = 1;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_flag;
                    RD_fgdEA_flag = object->Get<ND::TRealDatum>("fgdEA_flag");
                    if (RD_fgdEA_flag && i == fgdNum) {
                        fgd->isFgdEA_flag = int(RD_fgdEA_flag->GetValue());
                    } else
                        fgd->isFgdEA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_globpassEA_flag;
                    RD_globpassEA_flag = object->Get<ND::TRealDatum>("globpassEA_flag");
                    if (RD_globpassEA_flag && i == fgdNum) {
                        fgd->globpassEA_flag = int(RD_globpassEA_flag->GetValue());
                    } else
                        fgd->globpassEA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_otherUpQ;
                    RD_fgdEA_otherUpQ = object->Get<ND::TRealDatum>("fgdEA_otherUpQ");
                    if (RD_fgdEA_otherUpQ && i == fgdNum) {
                        fgd->fgdEA_otherUpQ = RD_fgdEA_otherUpQ->GetValue();
                    } else
                        fgd->fgdEA_otherUpQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_otherDownQ;
                    RD_fgdEA_otherDownQ = object->Get<ND::TRealDatum>("fgdEA_otherDownQ");
                    if (RD_fgdEA_otherDownQ && i == fgdNum) {
                        fgd->fgdEA_otherDownQ = RD_fgdEA_otherDownQ->GetValue();
                    } else
                        fgd->fgdEA_otherDownQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verQ;
                    RD_fgdEA_verQ = object->Get<ND::TRealDatum>("fgdEA_verQ");
                    if (RD_fgdEA_verQ && i == fgdNum) {
                        fgd->fgdEA_verQ = RD_fgdEA_verQ->GetValue();
                    } else
                        fgd->fgdEA_verQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verLayQ;
                    RD_fgdEA_verLayQ = object->Get<ND::TRealDatum>("fgdEA_verLayQ");
                    if (RD_fgdEA_verLayQ && i == fgdNum) {
                        fgd->fgdEA_verLayQ = RD_fgdEA_verLayQ->GetValue();
                    } else
                        fgd->fgdEA_verLayQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNearQ;
                    RD_fgdEA_verNearQ = object->Get<ND::TRealDatum>("fgdEA_verNearQ");
                    if (RD_fgdEA_verNearQ && i == fgdNum) {
                        fgd->fgdEA_verNearQ = RD_fgdEA_verNearQ->GetValue();
                    } else
                        fgd->fgdEA_verNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNearQ;
                    RD_fgdEA_verNextNearQ =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNearQ");
                    if (RD_fgdEA_verNextNearQ && i == fgdNum) {
                        fgd->fgdEA_verNextNearQ = RD_fgdEA_verNextNearQ->GetValue();
                    } else
                        fgd->fgdEA_verNextNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNextNearQ;
                    RD_fgdEA_verNextNextNearQ =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNextNearQ");
                    if (RD_fgdEA_verNextNextNearQ && i == fgdNum) {
                        fgd->fgdEA_verNextNextNearQ = RD_fgdEA_verNextNextNearQ->GetValue();
                    } else
                        fgd->fgdEA_verNextNextNearQ = 0;
                    
                    // End Activity in FGD1 2x3
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNearQ_rect;
                    RD_fgdEA_verNearQ_rect =
                    object->Get<ND::TRealDatum>("fgdEA_verNearQ_rect");
                    if (RD_fgdEA_verNearQ_rect && i == fgdNum) {
                        fgd->fgdEA_verNearQ_rect = RD_fgdEA_verNearQ_rect->GetValue();
                    } else
                        fgd->fgdEA_verNearQ_rect = 0;
                    
                    // End activity in FGD1 2x5
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNearQ_rect;
                    RD_fgdEA_verNextNearQ_rect =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNearQ_rect");
                    if (RD_fgdEA_verNextNearQ_rect && i == fgdNum) {
                        fgd->fgdEA_verNextNearQ_rect =
                        RD_fgdEA_verNextNearQ_rect->GetValue();
                    } else
                        fgd->fgdEA_verNextNearQ_rect = 0;
                    
                    // End activity in FGD1 2x7
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNextNearQ_rect;
                    RD_fgdEA_verNextNextNearQ_rect =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNextNearQ_rect");
                    if (RD_fgdEA_verNextNextNearQ_rect && i == fgdNum) {
                        fgd->fgdEA_verNextNextNearQ_rect =
                        RD_fgdEA_verNextNextNearQ_rect->GetValue();
                    } else
                        fgd->fgdEA_verNextNextNearQ_rect = 0;
                    
                }
                //--------------------------------------------------
                // FGD2 VA/EA variables
                else if (i == 1) {
                    // Vertex activity values
                    fgdNum = 2;
                    ND::THandle<ND::TRealDatum> RD_fgdVA_FgdN_FGD2;
                    RD_fgdVA_FgdN_FGD2 = object->Get<ND::TRealDatum>("fgdVA_FgdN_FGD2");
                    if (RD_fgdVA_FgdN_FGD2) {
                        fgdNum = int(RD_fgdVA_FgdN_FGD2->GetValue());
                    }
                    
                    fgd->isFgdVA = 0;
                    if (fgdNum == i) fgd->isFgdVA = 1;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_flag_FGD2;
                    RD_fgdVA_flag_FGD2 = object->Get<ND::TRealDatum>("fgdVA_flag_FGD2");
                    if (RD_fgdVA_flag_FGD2 && i == fgdNum) {
                        fgd->isFgdVA_flag = int(RD_fgdVA_flag_FGD2->GetValue());
                    } else
                        fgd->isFgdVA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_globpassVA_flag_FGD2;
                    RD_globpassVA_flag_FGD2 =
                    object->Get<ND::TRealDatum>("globpassVA_flag_FGD2");
                    if (RD_globpassVA_flag_FGD2 && i == fgdNum) {
                        fgd->globpassVA_flag = int(RD_globpassVA_flag_FGD2->GetValue());
                    } else
                        fgd->globpassVA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_upMinZ_FGD2;
                    RD_fgdVA_upMinZ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_upMinZ_FGD2");
                    if (RD_fgdVA_upMinZ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_upMinZ = RD_fgdVA_upMinZ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_upMinZ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_downMaxZ_FGD2;
                    RD_fgdVA_downMaxZ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_downMaxZ_FGD2");
                    if (RD_fgdVA_downMaxZ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_downMaxZ = RD_fgdVA_downMaxZ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_downMaxZ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_upCosTheta_FGD2;
                    RD_fgdVA_upCosTheta_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_upCosTheta_FGD2");
                    if (RD_fgdVA_upCosTheta_FGD2 && i == fgdNum) {
                        fgd->fgdVA_upCosTheta = RD_fgdVA_upCosTheta_FGD2->GetValue();
                    } else
                        fgd->fgdVA_upCosTheta = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_downCosTheta_FGD2;
                    RD_fgdVA_downCosTheta_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_downCosTheta_FGD2");
                    if (RD_fgdVA_downCosTheta_FGD2 && i == fgdNum) {
                        fgd->fgdVA_downCosTheta = RD_fgdVA_downCosTheta_FGD2->GetValue();
                    } else
                        fgd->fgdVA_downCosTheta = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_otherUpQ_FGD2;
                    RD_fgdVA_otherUpQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_otherUpQ_FGD2");
                    if (RD_fgdVA_otherUpQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_otherUpQ = RD_fgdVA_otherUpQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_otherUpQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_otherDownQ_FGD2;
                    RD_fgdVA_otherDownQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_otherDownQ_FGD2");
                    if (RD_fgdVA_otherDownQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_otherDownQ = RD_fgdVA_otherDownQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_otherDownQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verQ_FGD2;
                    RD_fgdVA_verQ_FGD2 = object->Get<ND::TRealDatum>("fgdVA_verQ_FGD2");
                    if (RD_fgdVA_verQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_verQ = RD_fgdVA_verQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_verQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verLayQ_FGD2;
                    RD_fgdVA_verLayQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_verLayQ_FGD2");
                    if (RD_fgdVA_verLayQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_verLayQ = RD_fgdVA_verLayQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_verLayQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNearQ_FGD2;
                    RD_fgdVA_verNearQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_verNearQ_FGD2");
                    if (RD_fgdVA_verNearQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_verNearQ = RD_fgdVA_verNearQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_verNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNearQ_FGD2;
                    RD_fgdVA_verNextNearQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNearQ_FGD2");
                    if (RD_fgdVA_verNextNearQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_verNextNearQ = RD_fgdVA_verNextNearQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_verNextNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNextNearQ_FGD2;
                    RD_fgdVA_verNextNextNearQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNextNearQ_FGD2");
                    if (RD_fgdVA_verNextNextNearQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_verNextNextNearQ =
                        RD_fgdVA_verNextNextNearQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_verNextNextNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_totalQ_FGD2;
                    RD_fgdVA_totalQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_totalQ_FGD2");
                    if (RD_fgdVA_totalQ_FGD2 && i == fgdNum) {
                        fgd->fgdVA_totalQ = RD_fgdVA_totalQ_FGD2->GetValue();
                    } else
                        fgd->fgdVA_totalQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdVA_totalCorrE_FGD2;
                    RD_fgdVA_totalCorrE_FGD2 =
                    object->Get<ND::TRealDatum>("fgdVA_totalCorrE_FGD2");
                    if (RD_fgdVA_totalCorrE_FGD2 && i == fgdNum) {
                        fgd->fgdVA_totalCorrE = RD_fgdVA_totalCorrE_FGD2->GetValue();
                    } else
                        fgd->fgdVA_totalCorrE = 0;
                    
                    // New Vertex Activity for FGD2 (avoid taking into account the scint
                    // layers after the water ones)
                    
                    // Vertex Activity in FGD2 2x3
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNearQ_FGD2_rect;
                    RD_fgdVA_verNearQ_FGD2_rect =
                    object->Get<ND::TRealDatum>("fgdVA_verNearQ_FGD2_rect");
                    if (RD_fgdVA_verNearQ_FGD2_rect && i == fgdNum) 
                    {
                        fgd->fgdVA_verNearQ_rect = RD_fgdVA_verNearQ_FGD2_rect->GetValue();
                    } else
                        fgd->fgdVA_verNearQ_rect = 0;
                    
                    // Vertex activity in FGD2 2x5
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNearQ_FGD2_rect;
                    RD_fgdVA_verNextNearQ_FGD2_rect =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNearQ_FGD2_rect");
                    if (RD_fgdVA_verNextNearQ_FGD2_rect && i == fgdNum) {
                        fgd->fgdVA_verNextNearQ_rect =
                        RD_fgdVA_verNextNearQ_FGD2_rect->GetValue();
                    } else
                        fgd->fgdVA_verNextNearQ_rect = 0;
                    
                    // Vertex activity in FGD2 2x7
                    ND::THandle<ND::TRealDatum> RD_fgdVA_verNextNextNearQ_FGD2_rect;
                    RD_fgdVA_verNextNextNearQ_FGD2_rect =
                    object->Get<ND::TRealDatum>("fgdVA_verNextNextNearQ_FGD2_rect");
                    if (RD_fgdVA_verNextNextNearQ_FGD2_rect && i == fgdNum) {
                        fgd->fgdVA_verNextNextNearQ_rect =
                        RD_fgdVA_verNextNextNearQ_FGD2_rect->GetValue();
                    } else
                        fgd->fgdVA_verNextNextNearQ_rect = 0;
                    
                    // End Activity values
                    // ----------------------------
                    fgdNum = 2;
                    ND::THandle<ND::TRealDatum> RD_fgdEA_FgdN_FGD2;
                    RD_fgdEA_FgdN_FGD2 = object->Get<ND::TRealDatum>("fgdEA_FgdN_FGD2");
                    if (RD_fgdEA_FgdN_FGD2) {
                        fgdNum = int(RD_fgdEA_FgdN_FGD2->GetValue());
                    }
                    
                    fgd->isFgdEA = 0;
                    if (fgdNum == i) fgd->isFgdEA = 1;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_flag_FGD2;
                    RD_fgdEA_flag_FGD2 = object->Get<ND::TRealDatum>("fgdEA_flag_FGD2");
                    if (RD_fgdEA_flag_FGD2 && i == fgdNum) {
                        fgd->isFgdEA_flag = int(RD_fgdEA_flag_FGD2->GetValue());
                    } else
                        fgd->isFgdEA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_globpassEA_flag_FGD2;
                    RD_globpassEA_flag_FGD2 =
                    object->Get<ND::TRealDatum>("globpassEA_flag_FGD2");
                    if (RD_globpassEA_flag_FGD2 && i == fgdNum) {
                        fgd->globpassEA_flag = int(RD_globpassEA_flag_FGD2->GetValue());
                    } else
                        fgd->globpassEA_flag = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_otherUpQ_FGD2;
                    RD_fgdEA_otherUpQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdEA_otherUpQ_FGD2");
                    if (RD_fgdEA_otherUpQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_otherUpQ = RD_fgdEA_otherUpQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_otherUpQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_otherDownQ_FGD2;
                    RD_fgdEA_otherDownQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdEA_otherDownQ_FGD2");
                    if (RD_fgdEA_otherDownQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_otherDownQ = RD_fgdEA_otherDownQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_otherDownQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verQ_FGD2;
                    RD_fgdEA_verQ_FGD2 = object->Get<ND::TRealDatum>("fgdEA_verQ_FGD2");
                    if (RD_fgdEA_verQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_verQ = RD_fgdEA_verQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_verQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verLayQ_FGD2;
                    RD_fgdEA_verLayQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdEA_verLayQ_FGD2");
                    if (RD_fgdEA_verLayQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_verLayQ = RD_fgdEA_verLayQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_verLayQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNearQ_FGD2;
                    RD_fgdEA_verNearQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdEA_verNearQ_FGD2");
                    if (RD_fgdEA_verNearQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_verNearQ = RD_fgdEA_verNearQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_verNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNearQ_FGD2;
                    RD_fgdEA_verNextNearQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNearQ_FGD2");
                    if (RD_fgdEA_verNextNearQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_verNextNearQ = RD_fgdEA_verNextNearQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_verNextNearQ = 0;
                    
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNextNearQ_FGD2;
                    RD_fgdEA_verNextNextNearQ_FGD2 =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNextNearQ_FGD2");
                    if (RD_fgdEA_verNextNextNearQ_FGD2 && i == fgdNum) {
                        fgd->fgdEA_verNextNextNearQ =
                        RD_fgdEA_verNextNextNearQ_FGD2->GetValue();
                    } else
                        fgd->fgdEA_verNextNextNearQ = 0;
                    
                    // End Activity in FGD2 2x3
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNearQ_FGD2_rect;
                    RD_fgdEA_verNearQ_FGD2_rect =
                    object->Get<ND::TRealDatum>("fgdEA_verNearQ_FGD2_rect");
                    if (RD_fgdEA_verNearQ_FGD2_rect && i == fgdNum) {
                        fgd->fgdEA_verNearQ_rect = RD_fgdEA_verNearQ_FGD2_rect->GetValue();
                    } else
                        fgd->fgdEA_verNearQ_rect = 0;
                    
                    // End activity in FGD2 2x5
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNearQ_FGD2_rect;
                    RD_fgdEA_verNextNearQ_FGD2_rect =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNearQ_FGD2_rect");
                    if (RD_fgdEA_verNextNearQ_FGD2_rect && i == fgdNum) {
                        fgd->fgdEA_verNextNearQ_rect =
                        RD_fgdEA_verNextNearQ_FGD2_rect->GetValue();
                    } else
                        fgd->fgdEA_verNextNearQ_rect = 0;
                    
                    // End activity in FGD2 2x7
                    ND::THandle<ND::TRealDatum> RD_fgdEA_verNextNextNearQ_FGD2_rect;
                    RD_fgdEA_verNextNextNearQ_FGD2_rect =
                    object->Get<ND::TRealDatum>("fgdEA_verNextNextNearQ_FGD2_rect");
                    if (RD_fgdEA_verNextNextNearQ_FGD2_rect && i == fgdNum) {
                        fgd->fgdEA_verNextNextNearQ_rect =
                        RD_fgdEA_verNextNextNearQ_FGD2_rect->GetValue();
                    } else
                        fgd->fgdEA_verNextNextNearQ_rect = 0;
                }
                // -----------------------------------------
        }
            }  // loop through FGD costituents
    }    // uses particular FGD
}

//*************************************************************
void ND::TGlobalReconModule::FillP0DInfo(ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject) 
{
    //*************************************************************
    
    if (!object->UsesDetector(ND::TReconBase::kP0D)) return;
    
    ND::THandle<ND::TReconBase> p0dObject =
    ReconObjectUtils::GetConstituentInDetector(object, ND::TReconBase::kP0D);
    
    if (p0dObject && globalObject.NP0Ds < NMAXP0D) {
        TP0DObject* p0d =
        new ((*(globalObject.P0D))[globalObject.NP0Ds++]) TP0DObject;
        
        // fill general information
        FillSubBaseObject(p0dObject, *p0d, 1);
        
        p0d->Cone = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        p0d->Width = DEFAULT_MAX;
        p0d->EDeposit = DEFAULT_MAX;
        
        // fill average node time
        ND::THandle<ND::TRealDatum> p0dTimeRD =
        p0dObject->Get<ND::TRealDatum>("averageHitTime");
        if (p0dTimeRD) {
            p0d->avgtime = p0dTimeRD->GetValue();
        } else {
            p0d->avgtime = DEFAULT_MAX;
        }
        
        // 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;
        
        ND::THandle<ND::TReconPID> pid = p0dObject;
        
        for (alternate = pid->GetAlternates().begin();
            alternate != pid->GetAlternates().end(); alternate++) {
        ND::THandle<ND::TReconPID> alt = (*alternate);
        if (alt) {
            p0d->ParticleId.push_back(alt->GetParticleId());
            p0d->PIDWeight.push_back(alt->GetPIDWeight());
        } else
            ND280Warn("TReconPID alternate is not a TReconPID - Ignored");
            }
            
            ND::THandle<ND::TReconObjectContainer> conContainer =
            p0dObject->GetConstituents();
            if (conContainer) {
                // The  PID variables are not saved in the main object but in its
                // constituent
                ND::THandle<ND::TReconBase> p0dObject2 = *(conContainer->begin());
                if (p0dObject2) {
                    // The cone angle is only accesible form the shower object
                    p0d->Cone = TrackingUtils::GetCone(p0dObject2);
                    p0d->EDeposit = TrackingUtils::GetEDeposit(p0dObject2);
                    
                    if (GetObjectType(p0dObject).find("ND::TReconShower") ==
                        std::string::npos)
                    return;
                    
                    // Compute the length and the width for showers
                    ND::THandle<ND::THitSelection> hits = p0dObject2->GetHits();
                    if (hits) {
                        double totalCharge = 0.0;
                        double avgLen = 0.0;
                        double avgLen2 = 0.0;
                        
                        for (ND::THitSelection::iterator h = hits->begin(); h != hits->end();
                            ++h) {
                        // Find the width.
                        TVector3 diff = (*h)->GetPosition() - p0d->FrontPosition.Vect();
                        TVector3 perp =
                        diff - (diff * p0d->FrontDirection) * p0d->FrontDirection;
                        if ((*h)->IsXHit()) {
                            diff.SetY(0.0);
                            perp.SetY(0.0);
                        } else {
                            diff.SetX(0.0);
                            perp.SetX(0.0);
                        }
                        double q = std::max(3.0, (*h)->GetCharge());
                        p0d->Width += q * perp.Mag();
                        double len = diff.Mag();
                        avgLen += q * len;
                        avgLen2 += q * len * len;
                        totalCharge += q;
                            }
                            p0d->Width /= std::max(totalCharge, 1.0);
                            avgLen /= std::max(totalCharge, 1.0);
                            avgLen2 /= std::max(totalCharge, 1.0);
                            p0d->Length = avgLen2 - avgLen * avgLen;
                            p0d->Length = std::sqrt(std::max(0.0, p0d->Length));
                    }
                }
            }
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillSMRDInfo(ND::THandle<ND::TReconBase> object,
    TGlobalPID& globalObject) 
{
    //*************************************************************
    
    ND::TReconBase::Status dets[4] = {
        ND::TReconBase::kTopSMRD, ND::TReconBase::kBottomSMRD,
    ND::TReconBase::kLeftSMRD, ND::TReconBase::kRightSMRD};
    
    for (int i = 0; i < 4; i++) {
        if (!object->UsesDetector(dets[i])) continue;
        
        // account for possible >1 segments in one SMRD quadrant
        ND::TReconObjectContainer smrdObjects =
        ReconObjectUtils::GetAllConstituentsInDetector(object, dets[i]);
        
        if (smrdObjects.size() == 0) continue;
        
        for (ND::TReconObjectContainer::iterator iter = smrdObjects.begin();
            iter != smrdObjects.end(); iter++) {
        ND::THandle<ND::TReconBase> smrdObject = *iter;
        
        if (smrdObject && globalObject.NSMRDs < NMAXSMRD) {
            TSMRDObject* smrd =
            new ((*(globalObject.SMRD))[globalObject.NSMRDs++]) TSMRDObject;
            
            // fill general information
            FillSubBaseObject(smrdObject, *smrd, i + 1);
            
            // explicitly set SMRD EDeposit as a some of the charge of all hits
            smrd->EDeposit = DEFAULT_MAX;
            ND::THandle<ND::TRealDatum> totalCharge =
            smrdObject->Get<ND::TRealDatum>("totalCharge");
            if (totalCharge)
                smrd->EDeposit =
            smrdObject->Get<ND::TRealDatum>("totalCharge")->GetValue();
            
            // fill average  hit time
            smrd->avgtime = DEFAULT_MAX;
            ND::THandle<ND::TRealDatum> smrdTimeRD =
            smrdObject->Get<ND::TRealDatum>("averageHitTime");
            if (smrdTimeRD) smrd->avgtime = smrdTimeRD->GetValue();
        }
            }
    }
}

//*************************************************************
void ND::TGlobalReconModule::InitializeExtrapolationToDetectors() 
{
    //*************************************************************
    
    std::vector<ND::TGeomModuleBase*>::const_iterator it;
    
    fALLMODULES.resize(NDETSEXTRAP);
    
    // initialize all names
    for (int i = 0; i < NDETSEXTRAP; i++) fALLMODULES[i] = "";
    
    // Get TPC volumes
    std::vector<ND::TGeomModuleBase*> TPCMODULES =
    ND::TGeomInfo::Get().TPC().GetModules();
    int index = 0;
    for (it = TPCMODULES.begin(); it != TPCMODULES.end(); it++) {
        std::string name = (*it)->GetPath() + "/GasGap_0/Drift_0";
        fALLMODULES[index] = name;
        index++;
    }
    
    // Get FGD volumes
    std::vector<ND::TGeomModuleBase*> FGDMODULES =
    ND::TGeomInfo::Get().FGD().GetModules();
    index = 3;
    for (it = FGDMODULES.begin(); it != FGDMODULES.end(); it++) {
        fALLMODULES[index] = (*it)->GetPath() + "/Active_0";
        index++;
    }
    
    // Get P0D volumes
    std::vector<ND::TGeomModuleBase*> P0DMODULES =
    ND::TGeomInfo::Get().P0D().GetModules();
    index = 5;
    for (it = P0DMODULES.begin(); it != P0DMODULES.end(); it++) {
        fALLMODULES[index] = (*it)->GetPath();
        index++;
    }
    
    // Get DsECAL volume
    const std::vector<ND::TGeomModuleBase*> ECALMODULES =
    ND::TGeomInfo::Get().ECAL().GetModules();
    index = 6;
    for (it = ECALMODULES.begin(); it != ECALMODULES.end(); it++) {
        if ((*it)->GetPath().find("DsECal") != std::string::npos) {
            fALLMODULES[index] = (*it)->GetPath() + "/Active_0";
            index++;
        }
    }
    
    // Get SMRD volumes
    std::string SMRD_name[6];
    std::string magnet_name = ND::tman().GetSetup().GetVolumePath("Magnet_0");
    if (magnet_name != "") {
        SMRD_name[0] = magnet_name + "/RightClam_0/SMRD_0/Side_0";
        SMRD_name[1] = magnet_name + "/RightClam_0/SMRD_0/Top_0";
        SMRD_name[2] = magnet_name + "/RightClam_0/SMRD_0/Bottom_0";
        SMRD_name[3] = magnet_name + "/LeftClam_0/SMRD_0/Side_0";
        SMRD_name[4] = magnet_name + "/LeftClam_0/SMRD_0/Top_0";
        SMRD_name[5] = magnet_name + "/LeftClam_0/SMRD_0/Bottom_0";
        
        index = 7;
        for (int i = 0; i < 6; i++) {
            fALLMODULES[index] = SMRD_name[i];
            index++;
        }
    }
    
    // Get Barrel ECAL volumes
    index = 13;
    for (it = ECALMODULES.begin(); it != ECALMODULES.end(); it++) {
        if ((*it)->GetPath().find("DsECal") == std::string::npos) {
            fALLMODULES[index] = (*it)->GetPath() + "/Active_0";
            index++;
        }
    }
    
    //  const int NDETS = fALLMODULES.size();
    
    for (int i = 0; i < NDETSEXTRAP; i++) {
        std::string volname = fALLMODULES[i];
        if (volname == "") continue;
        const Volume& vol = ND::gman().GetSetup().volume(volname);
        const surface_vector& surfs = vol.surfaces("outer");
        stc_tools::extend_vector<Surface*>(surfs, fDetSurfaces[i]);
        
        fDetIndex[vol.name("setup_name")] = i;
    }
    
    if (debug_extrap) {
        std::cout << "list of volumes: " << std::endl;
        for (int i = 0; i < NDETSEXTRAP; i++) {
            std::cout << i << " --> " << fALLMODULES[i] << std::endl;
        }
        
        std::cout << "fDetIndex" << std::endl;
        //  std::map<std::string, int>::const_iterator it2;
        //  dict::dictionary<int>::const_iterator it2;
        for (unsigned int it2 = 0; it2 < fDetIndex.keys().size(); it2++) {
            std::cout << fDetIndex.keys()[it2] << " "
            << fDetIndex[fDetIndex.keys()[it2]] << std::endl;
            int i = fDetIndex[fDetIndex.keys()[it2]];
            for (unsigned int j = 0; j < fDetSurfaces[i].size(); j++)
                std::cout << "    " << j << ":  "
            << (fDetSurfaces[i])[j]->name("setup_name") << std::endl;
        }
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillExtrapolationToDetectors(
    ND::THandle<ND::TReconBase> object, TGlobalPID& globalObject) 
{
    //*************************************************************
    
    if (debug_extrap) {
        std::cout << "--------------------------------------------------"
        << std::endl;
        std::cout << "TGlobalReconModule:: Extrapolate to subdetectors: "
        << std::endl;
        std::cout << "--------------------------------------------------"
        << std::endl;
    }
    
    // Initialize variables
    for (int i = 0; i < NDETSEXTRAP; i++) {
        fPassedDetector[i] = false;
        globalObject.EntranceOK[i] = 0;
        globalObject.ExitOK[i] = 0;
        
        globalObject.EntrancePosition[i] =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        globalObject.ExitPosition[i] =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        
        globalObject.EntranceDirection[i] =
        TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        globalObject.ExitDirection[i] =
        TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        
        globalObject.EntranceMomentum[i] = DEFAULT_MAX;
        globalObject.ExitMomentum[i] = DEFAULT_MAX;
        
        globalObject.EntrancePositionErr[i] =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        globalObject.ExitPositionErr[i] =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        
        globalObject.EntranceDirectionErr[i] =
        TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        globalObject.ExitDirectionErr[i] =
        TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        
        globalObject.EntranceMomentumErr[i] = DEFAULT_MAX;
        globalObject.ExitMomentumErr[i] = DEFAULT_MAX;
    }
    
    // convert the track to a RecPack trajectory
    Trajectory traj;
    bool ok = ND::converter().TReconBase_to_Trajectory(object, traj);
    if (!ok) return;
    
    if (!traj.status("fitted")) {
        return;
    }
    
    if (traj.state(traj.first_fitted_node()).vector()[5] < 0) {
        traj.sort_nodes(RP::z, -1);
        if (debug_extrap)
            std::cout << "Negative sense: sort nodes in reverse z order" << std::endl;
    }
    
    if (debug_extrap) std::cout << "Trajectory: " << traj << std::endl;
    
    // extrapolate forwards
    FillExtrapolationToDetectors(object, globalObject, traj, 1);
    
    // extrapolate backwards
    FillExtrapolationToDetectors(object, globalObject, traj, -1);
    
    if (debug_extrap) {
        std::cout << "---------------------------------------------------------"
        << std::endl;
        std::cout << "TGlobalReconModule:: Extrapolate to subdetectors results:"
        << std::endl;
        std::cout << "---------------------------------------------------------"
        << std::endl;
        
        std::cout << object << std::endl;
        
        for (int i = 0; i < NDETSEXTRAP; i++) {
            if (globalObject.EntranceOK[i])
                std::cout << i << "   x = " << globalObject.EntrancePosition[i].X()
            << "\t y = " << globalObject.EntrancePosition[i].Y()
            << "\t z = " << globalObject.EntrancePosition[i].Z()
            << "\t ex = " << globalObject.EntrancePositionErr[i].X()
            << "\t ey = " << globalObject.EntrancePositionErr[i].Y()
            << "\t ez = " << globalObject.EntrancePositionErr[i].Z()
            << std::endl;
            if (globalObject.ExitOK[i])
                std::cout << i << "   x = " << globalObject.ExitPosition[i].X()
            << "\t y = " << globalObject.ExitPosition[i].Y()
            << "\t z = " << globalObject.ExitPosition[i].Z()
            << "\t ex = " << globalObject.ExitPositionErr[i].X()
            << "\t ey = " << globalObject.ExitPositionErr[i].Y()
            << "\t ez = " << globalObject.ExitPositionErr[i].Z()
            << std::endl;
            
            if (globalObject.EntranceOK[i])
                std::cout << i << "   ux = " << globalObject.EntranceDirection[i].X()
            << "\t uy = " << globalObject.EntranceDirection[i].Y()
            << "\t uz = " << globalObject.EntranceDirection[i].Z()
            << "\t ex = " << globalObject.EntranceDirectionErr[i].X()
            << "\t ey = " << globalObject.EntranceDirectionErr[i].Y()
            << "\t ez = " << globalObject.EntranceDirectionErr[i].Z()
            << std::endl;
            if (globalObject.ExitOK[i])
                std::cout << i << "   ux = " << globalObject.ExitDirection[i].X()
            << "\t uy = " << globalObject.ExitDirection[i].Y()
            << "\t uz = " << globalObject.ExitDirection[i].Z()
            << "\t ex = " << globalObject.ExitDirectionErr[i].X()
            << "\t ey = " << globalObject.ExitDirectionErr[i].Y()
            << "\t ez = " << globalObject.ExitDirectionErr[i].Z()
            << std::endl;
            
            if (globalObject.EntranceOK[i])
                std::cout << i << "   P0 = " << (globalObject.EntranceMomentum[i])
            << "   eP0 = " << (globalObject.EntranceMomentumErr[i])
            << std::endl;
            if (globalObject.ExitOK[i])
                std::cout << i << "   P1 = " << (globalObject.ExitMomentum[i])
            << "   eP1 = " << (globalObject.ExitMomentumErr[i])
            << std::endl;
        }
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillExtrapolationToDetectors(
    ND::THandle<ND::TReconBase> object, TGlobalPID& globalObject,
    Trajectory& traj, int sense) 
{
    //*************************************************************
    
    // 1. Take the first state in the trajectory (taking into account the sense)
    int ifirst = traj.first_fitted_node();
    State state = traj.state(ifirst);
    
    if (debug_extrap) {
        std::cout << "Extrapolate to subdetectors: " << std::endl;
        std::cout << "1. First state: " << state << std::endl;
    }
    
    // 2. Find the volume associated to the first state
    
    const Volume* volume = NULL;
    EVector position = fEquation.position(state, 0.);
    bool ok = ND::gman().GetSetup().volume(position, volume);
    if (!ok) {
        if (debug_extrap)
            std::cout << "Volume not found. Position:  " << position << std::endl;
        return;
    }
    
    int current_det = -1;
    for (unsigned int i = 0; i < fALLMODULES.size(); i++) {
        if (volume->name("setup_name").find(fALLMODULES[i]) != std::string::npos)
            current_det = i;
    }
    
    if (current_det < 0) {
        if (debug_extrap) {
            std::cout << "Extrapolate to subdetectors: " << std::endl;
            std::cout << "2. Current volume not in list, name = "
            << volume->name("setup_name") << std::endl;
        }
        return;
    }
    
    if (debug_extrap) {
        std::cout << "Extrapolate to subdetectors: " << std::endl;
        std::cout << "2. Current volume: index = " << current_det
        << ", name = " << volume->name("setup_name") << std::endl;
    }
    
    fPassedDetector[current_det] = false;
    int iexit = -1;
    bool in_found = true;
    int nLoops{0};
    while (volume) {
        nLoops++;

        if (debug_extrap) {
           std::cout<<" Loop: "<<nLoops<<std::endl;
           std::cout << "Extrapolate to subdetectors: " << std::endl;
           std::cout << "3. Exit from volume: " << volume->name("setup_name")
           << std::endl;
           std::cout << "   ifirst: " << ifirst << ", iexit = " << iexit
           << ", in found = " << in_found << std::endl;
       }
        
        // 3. Loop over states until we exit from the volume. Only when there is
        // an
        // state inside the current volume (in_found=true)
        if (in_found) {
            // by default the exit state is the first one
            iexit = ifirst;
            if (sense == 1) {
                int iprev = ifirst;
                for (unsigned int i = ifirst + 1; i < traj.nodes().size(); i++) {
                    if (traj.node(i).status("fitted")) {
                        position = fEquation.position(traj.state(i), 0.);
                        bool ok = volume->is_inside(position);
                        if (!ok) {
                            // this is the first state outside the volume. Take the previous
                            // as exit state
                            iexit = iprev;
                            break;
                        }
                        iprev = i;
                    }
                }
            }
            if (debug_extrap)
                std::cout << "   Last state in volume: " << iexit << std::endl;
            
            // take the state
            state = traj.state(iexit);
            
            if (debug_extrap) std::cout << state << std::endl;
        } else {
            if (debug_extrap)
                std::cout << "   No state inside volume. Propagate to exit surface "
            "from the entrance state "
            << std::endl;
        }
        
        // 4. Find the closest surface to the last state in the volume (or the
        // previous extrapolation if there is no state in the volume)
        surface_vector exit_surfaces;
        stc_tools::extend_vector<Surface*>(fDetSurfaces[current_det],
            exit_surfaces);
        
        double end_length = 1000000 * sense;
        double delta_length;
        const Surface* exit_surf;
        
        ok = ND::rpman("globalAnalysis")
        .navigation_svc()
        .next_surface(state, end_length, exit_surfaces, delta_length,
            exit_surf);
        
        if (ok) {
            if (debug_extrap) {
                std::cout << "Extrapolate to subdetectors: " << std::endl;
                std::cout << "4. Exit surface. Length: " << delta_length << std::endl;
                std::cout << "                 name:   "
                << exit_surf->name("setup_name") << std::endl;
                std::cout << "                 volume: " << volume->name("setup_name")
                << std::endl;
            }
        } else {
            if (debug_extrap) {
                std::cout << "Extrapolate to subdetectors: " << std::endl;
                std::cout << "4. Exit surface. Fail in finding exit surface. stop !!! "
                << std::endl;
            }
            return;
        }
        
        // 5. propagate to the surface
        ok = ND::rpman("globalAnalysis")
        .navigation_svc()
        .propagate(0.1 * sense, state);
        ok = ND::rpman("globalAnalysis")
        .navigation_svc()
        .propagate(*exit_surf, state, delta_length);
        
        if (ok) {
            // take the position as the exit position
            EVector exitPosition = fEquation.position(state, 0.);
            EVector exitDirection = fEquation.direction(state, 0.);
            double exitMomentum = DEFAULT_MAX;
            if (state.vector()[6] != 0) exitMomentum = fabs(1 / state.vector()[6]);
            
            TLorentzVector posErr(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            if (state.matrix()[0][0] > 0 && state.matrix()[1][1] > 0 &&
                state.matrix()[2][2] > 0)
            posErr = TLorentzVector(sqrt(state.matrix()[0][0]),
                sqrt(state.matrix()[1][1]),
                sqrt(state.matrix()[2][2]), 0);
            
            TVector3 dirErr(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            if (state.matrix()[3][3] > 0 && state.matrix()[4][4] > 0 &&
                state.matrix()[5][5] > 0)
            dirErr =
            TVector3(sqrt(state.matrix()[3][3]), sqrt(state.matrix()[4][4]),
                sqrt(state.matrix()[5][5]));
            
            double exitMomentumErr = DEFAULT_MAX;
            if (state.matrix()[6][6] > 0 && state.vector()[6] != 0)
                exitMomentumErr =
            fabs(1 / pow(state.vector()[6], 2)) * sqrt(state.matrix()[6][6]);
            
            if (debug_extrap) {
                std::cout << "Extrapolate to subdetectors: " << std::endl;
                std::cout << "5. Exit position:  " << print(exitPosition) << std::endl;
                std::cout << "   Exit direction: " << print(exitDirection) << std::endl;
                std::cout << "   Exit momentum:  " << exitMomentum << std::endl;
            }
            
            TLorentzVector pos;
            ok = ND::converter().EVector3_to_TVector(exitPosition, pos);
            
            TLorentzVector dir4;
            ok = ND::converter().EVector3_to_TVector(exitDirection, dir4);
            TVector3 dir = dir4.Vect();
            
            if (sense == 1) {
                globalObject.ExitOK[current_det] = 1;
                globalObject.ExitPosition[current_det] = pos;
                globalObject.ExitDirection[current_det] = dir;
                globalObject.ExitMomentum[current_det] = exitMomentum;
                globalObject.ExitPositionErr[current_det] = posErr;
                globalObject.ExitDirectionErr[current_det] = dirErr;
                globalObject.ExitMomentumErr[current_det] = exitMomentumErr;
            } else {
                globalObject.EntranceOK[current_det] = 1;
                globalObject.EntrancePosition[current_det] = pos;
                globalObject.EntranceDirection[current_det] = dir;
                globalObject.EntranceMomentum[current_det] = exitMomentum;
                globalObject.EntrancePositionErr[current_det] = posErr;
                globalObject.EntranceDirectionErr[current_det] = dirErr;
                globalObject.EntranceMomentumErr[current_det] = exitMomentumErr;
            }
            
        } else {
            if (debug_extrap) {
                std::cout << "Extrapolate to subdetectors: " << std::endl;
                std::cout << "5. Exit position: "
                << "surface not intersected !!!" << std::endl;
            }
        }
        
        fPassedDetector[current_det] = true;
        
        // 6. Find the closest surface to the exit state among all other volumes
        
        // list of possible entrance surfaces
        surface_vector entrance_surfaces;
        for (int det = 0; det < NDETSEXTRAP; det++) {
            if (!fPassedDetector[det]) {
                stc_tools::extend_vector<Surface*>(fDetSurfaces[det],
                    entrance_surfaces);
            }
        }
        
        const Surface* entrance_surf;
        ok = ND::rpman("globalAnalysis")
        .navigation_svc()
        .next_surface(state, end_length, entrance_surfaces, delta_length,
            entrance_surf);
        
        if (!ok) {
            if (debug_extrap) std::cout << "No suface intersected " << std::endl;
            break;
        }
        
        // find the parent volume
        const Volume* volume_try;
        ok = ND::gman().GetSetup().volume(*entrance_surf, volume_try);
        
        if (debug_extrap) {
            std::cout << "Extrapolate to subdetectors: " << std::endl;
            std::cout << "6. Entrance surface. Length: " << delta_length << std::endl;
            std::cout << "                     name:   "
            << entrance_surf->name("setup_name") << std::endl;
            std::cout << "                     volume: "
            << volume_try->name("setup_name") << std::endl;
        }
        
        // 7.  Check if the next state in the trajectory is in that volume
        in_found = false;
        if (sense == 1) {
            if ((unsigned int)(iexit + 1) < traj.nodes().size()) {
                for (unsigned int i = iexit + 1; i < traj.nodes().size(); i++) {
                    if (traj.node(i).status("fitted")) {
                        position = fEquation.position(traj.state(i), 0.);
                        bool ok = volume_try->is_inside(position);
                        if (ok) {
                            ifirst = i;
                            in_found = true;
                        }
                        break;
                    }
                }
            }
        }
        
        if (debug_extrap) {
            std::cout << "Extrapolate to subdetectors: " << std::endl;
            std::cout << "7. Enter in volume: " << volume_try->name("setup_name")
            << std::endl;
            std::cout << "   ifirst: " << ifirst << ", iexit = " << iexit
            << ", in found = " << in_found << std::endl;
        }
        
        //    iexit=ifirst
        
        if (debug_extrap) {
            if (in_found)
                std::cout << "   Next state inside volume: " << ifirst << std::endl;
            else
                std::cout << "   No state inside volume. Propagate to entrance surface "
            "from the exit of the previous volume "
            << std::endl;
        }
        
        // 8.  propagate to the surface
        
        if (in_found) {
            // from the first state in the current volume
            state = traj.state(ifirst);
        }
        
        ok = ND::rpman("globalAnalysis")
        .navigation_svc()
        .propagate(0.1 * sense, state);
        ok = ND::rpman("globalAnalysis")
        .navigation_svc()
        .propagate(*entrance_surf, state, delta_length);
        
        if (ok) {
            volume = volume_try;
            current_det = -1;
            for (unsigned int i = 0; i < fALLMODULES.size(); i++) {
                if (volume->name("setup_name").find(fALLMODULES[i]) !=
                    std::string::npos)
                current_det = i;
            }
            
            // take the position as the entrance position
            EVector entrancePosition = fEquation.position(state, 0.);
            EVector entranceDirection = fEquation.direction(state, 0.);
            double entranceMomentum = DEFAULT_MAX;
            
            if (state.vector()[6] != 0)
                entranceMomentum = fabs(1 / state.vector()[6]);
            
            TLorentzVector posErr(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            if (state.matrix()[0][0] > 0 && state.matrix()[1][1] > 0 &&
                state.matrix()[2][2] > 0)
            posErr = TLorentzVector(sqrt(state.matrix()[0][0]),
                sqrt(state.matrix()[1][1]),
                sqrt(state.matrix()[2][2]), 0);
            
            TVector3 dirErr(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
            if (state.matrix()[3][3] > 0 && state.matrix()[4][4] > 0 &&
                state.matrix()[5][5] > 0)
            dirErr =
            TVector3(sqrt(state.matrix()[3][3]), sqrt(state.matrix()[4][4]),
                sqrt(state.matrix()[5][5]));
            
            double entranceMomentumErr = DEFAULT_MAX;
            if (state.matrix()[6][6] > 0 && state.vector()[6] != 0)
                entranceMomentumErr =
            fabs(1 / pow(state.vector()[6], 2)) * sqrt(state.matrix()[6][6]);
            
            if (debug_extrap) {
                std::cout << "8. Entrance position:  " << print(entrancePosition)
                << std::endl;
                std::cout << "   Entrance direction: " << print(entranceDirection)
                << std::endl;
                std::cout << "   Entrance momentum:  " << entranceMomentum << std::endl;
            }
            
            TLorentzVector pos;
            ok = ND::converter().EVector3_to_TVector(entrancePosition, pos);
            
            TLorentzVector dir4;
            ok = ND::converter().EVector3_to_TVector(entranceDirection, dir4);
            TVector3 dir = dir4.Vect();
            
            if (sense == 1) {
                globalObject.EntranceOK[current_det] = 1;
                globalObject.EntrancePosition[current_det] = pos;
                globalObject.EntranceDirection[current_det] = dir;
                globalObject.EntranceMomentum[current_det] = entranceMomentum;
                globalObject.EntrancePositionErr[current_det] = posErr;
                globalObject.EntranceDirectionErr[current_det] = dirErr;
                globalObject.EntranceMomentumErr[current_det] = entranceMomentumErr;
            } else {
                globalObject.ExitOK[current_det] = 1;
                globalObject.ExitPosition[current_det] = pos;
                globalObject.ExitDirection[current_det] = dir;
                globalObject.ExitMomentum[current_det] = entranceMomentum;
                globalObject.ExitPositionErr[current_det] = posErr;
                globalObject.ExitDirectionErr[current_det] = dirErr;
                globalObject.ExitMomentumErr[current_det] = entranceMomentumErr;
            }
        } else {
            if (debug_extrap) {
                std::cout << "Extrapolate to subdetectors: " << std::endl;
                std::cout << "8. Entrance position: "
                << "surface not intersected !!!" << std::endl;
            }
             if(nLoops>12)
            {
             ND280Log("ND::TGlobalReconModule::FillExtrapolationToDetectors : Extrapolation to subdetectors failed, and "<<nLoops<<" volumes have been attempted so abandoning extrapolation to avoid infinite loop.");
             return;
            }
        }
        
        // 9. Repit 3 to 8 for all other volumes.
    }
}

//*****************************************************************************
bool ND::TGlobalReconModule::FillVertex(
    ND::TND280Event& event, ND::THandle<ND::TReconVertex> vertex, bool primary)
{
//*****************************************************************************
  // Fill a reconstructed vertex and associated tracks
  
  // Maximum number of secondary vertices in spill
  if (!primary && fNSVertices == NSVERTICES) return false;
  // Check that the vertex has a state
  ND::THandle<ND::TVertexState> vs = vertex->GetState();
  if (!vs) return false;
  // Create a global vertex
  TGlobalVertex* gVertex = new ((*fVertices)[fNVertices++]) TGlobalVertex;
  gVertex->Type          = 0;
  gVertex->NTrueVertices = 0;
  
  // Primary
  if (primary) {
    fPVInd = fNVertices - 1;     // index of last added primary
    gVertex->PrimaryIndex = -1;  //
  } else {
    // For secondaries set the index of the corresponding primary vertex
    gVertex->PrimaryIndex = fPVInd;
    fNSVertices++;
  }
  // Basic info
  gVertex->Position      = vs->GetPosition();
  gVertex->Variance      = vs->GetPositionVariance();
  gVertex->AlgorithmName = vertex->GetAlgorithmName();
  gVertex->Status        = vertex->GetStatus();
  gVertex->Quality       = vertex->GetQuality();
  gVertex->NDOF          = vertex->GetNDOF();
  
  // Associated tracks
  ND::THandle<ND::TReconObjectContainer> tracksAssociated 
    = vertex->GetConstituents();
  if (tracksAssociated->size() == 0) {
    ND280NamedWarn("ECalVertexInTGRMFV",
        "Leaving as vertex has no associated constituents.");
    return true;
  }
  std::map<ND::THandle<ND::TReconBase>,
  ND::THandle<ND::TReconObjectContainer> >::iterator cit;
  // Loop over associated tracks
  for (ND::TReconObjectContainer::iterator assocIter =
       tracksAssociated->begin();
       assocIter != tracksAssociated->end(); assocIter++){
    if (ND::THandle<ND::TReconVertex> vtxDummy = (*assocIter)) {
      ND280NamedDebug("ECalVertexInTGRMFV", "Skipping vertex constituent: "
        << vtxDummy->GetUniqueID());
      continue;
    }
    
    // Fill track info
    if (gVertex->NConstituents == NCONSTITUENTS) break;
    ND::THandle<ND::TReconBase> track = *assocIter;
    if (!track) continue;
    ND::THandle<ND::TReconState> state = track->GetState();
    if (!state) continue;
    
    // The original track stored as constituent for Kalman vertices
    ND::THandle<ND::TReconBase> otrack;
    if (gVertex->AlgorithmName == "Kalman") {
      ND::THandle<ND::TReconObjectContainer> trackC = track->GetConstituents();
      if (trackC->size() == 0) continue;
      otrack = trackC->front();
      if (!otrack) continue;
    }
    // Get index of the global track in the overall fPIDs list
    int GlobalPIDIndex = -1;
    bool inMap = false;
    std::map<ND::THandle<ND::TReconBase>, int>::iterator it, miter;
    
    TLorentzVector pos_track = TrackingUtils::GetPosition(state);
    double mom_track = TrackingUtils::GetMomentum(state);
    
    // Use original track for Kalman vertices
    // First look for ND::TIntegerDatum>("PIDID") in the track itself
    int PIDID=0;
    if(track->Has<ND::TIntegerDatum>("PIDID")){
      PIDID=track->Get<ND::TIntegerDatum>("PIDID")->at(0);
    }
    // failing that look in contained track
    if (otrack) {
      ND::THandle<ND::TReconState> ostate = otrack->GetState();
      mom_track = TrackingUtils::GetMomentum(ostate);
      pos_track = TrackingUtils::GetPosition(ostate);
      if(PIDID<1&&otrack->Has<ND::TIntegerDatum>("PIDID")){
        PIDID=otrack->Get<ND::TIntegerDatum>("PIDID")->at(0);
      }
    }
    if(PIDID!=0){
      for (it = fGlobalIndexMap.begin(); it != fGlobalIndexMap.end(); it++) {
        if (!it->first) continue;
        if((it->first)->Has<ND::TIntegerDatum>("PIDID")){
          if( (it->first)->Get<ND::TIntegerDatum>("PIDID")->at(0)==PIDID){
            inMap = true;
            miter = it;
            break;
          }
        }
      }
    }
    /*}
    }
    else
    { // no contained copy , look directly at the track itself
    if(track->Has<ND::TIntegerDatum>("PIDID"))
    {
    int PIDID=track->Get<ND::TIntegerDatum>("PIDID")->at(0);
    if(PIDID!=0)
    {
    for (it = fGlobalIndexMap.begin(); it != fGlobalIndexMap.end(); it++) {
    if (!it->first) continue;
    if((it->first)->Has<ND::TIntegerDatum>("PIDID"))
    {
    if( (it->first)->Get<ND::TIntegerDatum>("PIDID")->at(0)==PIDID)
    {
    inMap = true;
    miter = it;
    break;
    }
    }
    }
    }
    }
    }
    */
    // Global vertex#
    if(!inMap){
      for (it = fGlobalIndexMap.begin(); it != fGlobalIndexMap.end(); it++) {
        if (!it->first) continue;
        ND::THandle<ND::TReconState> state_tmp = (it->first)->GetState();
        if (!state_tmp) continue;
        double mom_tmp = TrackingUtils::GetMomentum(state_tmp);
        TLorentzVector pos_tmp = TrackingUtils::GetPosition(state_tmp);
        if (pos_track == pos_tmp && mom_track == mom_tmp){
          miter = it;
          inMap = true;
          break;
        }
      }
    }
    
    // Tracker vertex
    if (!inMap) {
      if (otrack)
        miter = fTrackerGlobalIndexMap.find(otrack);
      else
        miter = fTrackerGlobalIndexMap.find(track);
      if (miter != fTrackerGlobalIndexMap.end()) inMap = true;
    }
    
    if (fUseECalVertices) {
      if (!inMap) {
        ND280NamedWarn(
            "ECalVertexInTGRMFV",
            "Didn't find this constituent("
            << track->GetUniqueID() << ", " << track->ConvertDetector()
            << " " << track->ClassName() << ") in the global map of PIDs.");
        gVertex->Quality = 0;
      } else {
        ND280NamedVerbose("ECalVertexInTGRMFV",
            "Adding constituent("
            << track->GetUniqueID() << ", "
            << track->ConvertDetector()
            << ") which was in the global map of PIDs.");
        gVertex->Quality = 1;
      }
    }
    
    // Set PID
    if (inMap) {
      GlobalPIDIndex = miter->second;
    } else {  
      // Dont add the constituent and issue a warning
      // for the moment log it anyway,  since this indicates a special case
      // and better to pay attention
      std::cout << " WARNING: TGlobalReconModule::FillVertex(): GV constituent "
                << "not found among recon tracks! Event ID is "
                << event.GetEventId() << std::endl;
      continue;
    }
    
    TVertexConstituent* vConst =
    new ((*(gVertex->Constituents))[gVertex->NConstituents++])
    TVertexConstituent;
    
    vConst->Charge = (int)TrackingUtils::GetCharge(track);
    vConst->Quality = track->GetQuality();
    
    double mom = TrackingUtils::GetMomentum(state);
    TVector3 mdir = TrackingUtils::GetDirection(state);
    mdir *= mom;
    
    vConst->Momentum = mdir;
    vConst->PID = GlobalPIDIndex;
  }
  return true;
}

//*************************************************************
void ND::TGlobalReconModule::MatchTrueVertex(
    ND::THandle<ND::TG4PrimaryVertexContainer> g4PrimVert) 
{
    //*************************************************************
    
    TGlobalVertex* gVertex = (TGlobalVertex*)(*fVertices)[fNVertices - 1];
    std::vector<TLorentzVector> vertexPos;
    std::vector<std::vector<int> > vertexPIDs;
    
    int nConst = gVertex->NConstituents;
    
    // Loop over constituents and create a list of trueVertices with their
    // constituents for this global vertex
    for (int i = 0; i < nConst; i++) {
        if (gVertex->NTrueVertices == NCONSTITUENTS)
            break;  // why is needed this limit? (max 6 globalV->trueV)
        TVertexConstituent* vConst =
        (TVertexConstituent*)(*(gVertex->Constituents))[i];
        int this_PID = vConst->PID;
        TGlobalPID* globalObject = (TGlobalPID*)(*fPIDs)[this_PID];
        TTrueParticle globalTrue = globalObject->TrueParticle;
        TLorentzVector vPos = globalTrue.Vertex.Position;
        bool inList = false;
        for (int j = 0; j < (int)vertexPos.size(); j++) {
            if (vertexPos[j] == vPos) inList = true;
        }
        if (inList) continue;
        vertexPos.push_back(vPos);
        
        // Fill globalV->trueV
        TTrueVertex* tVertex =
        new ((*(gVertex->TrueVertices))[gVertex->NTrueVertices++]) TTrueVertex;
        *tVertex = globalTrue.Vertex;
        
        // ************** the following part of this function, that evaluates
        // Cleanliness and Completeness
        // ************** has been replaced with in highland, where objects are
        // "bunched"
        
        // Loop over all PIDs, and store those that come from this true vertex
        std::vector<int> IDs;
        for (int k = 0; k < fNPIDs; k++) {
            TGlobalPID* thisGlobalObject = (TGlobalPID*)(*(fPIDs))[k];
            
            if (thisGlobalObject->NNodes > 0)
                ;
            else
                continue;  // track not used in vertexing algorithm
            if (thisGlobalObject->Status != 1)
                continue;  // track not used in vertexing algorithm
            if (thisGlobalObject->DetectorUsed[0] != 1 &&
                thisGlobalObject->DetectorUsed[1] != 1 &&
                thisGlobalObject->DetectorUsed[2] != 1 &&
                thisGlobalObject->DetectorUsed[3] != 1 &&
                thisGlobalObject->DetectorUsed[4] != 1)
            continue;  // track not used in vertexing algorithm (only FGDs or TPCs
            // tracks)
            
            TTrueParticle thisGlobalTrue = thisGlobalObject->TrueParticle;
            if (thisGlobalTrue.ID == -1) continue;  // track without true
            
            // here we should also select only primary true tracks and only in the
            // same bunch (that's why we do this in highland)
            
            TLorentzVector thisVPos = thisGlobalTrue.Vertex.Position;
            // broken tracks are double counted
            if (thisVPos == vPos) IDs.push_back(k);
        }
        vertexPIDs.push_back(IDs);
    }
    
    // Set efficiencies and purities FOR EACH TRUEV associated to the GlbV under
    // study
    // Loop over the found true vertices
    for (int i = 0; i < gVertex->NTrueVertices; i++) {
        std::vector<int> IDs = vertexPIDs[i];
        int nIDs = (int)IDs.size();
        int nMatch = 0;
        
        // Loop over vertex constituents
        for (int j = 0; j < nConst; j++) {
            TVertexConstituent* vConst =
            (TVertexConstituent*)(*(gVertex->Constituents))[j];
            int this_PID = vConst->PID;
            // Check if they match true vertex PID list
            for (int k = 0; k < nIDs; k++) {
                if (this_PID == IDs[k]) nMatch++;
            }
        }
        
        TTrueVertex* tVertex = (TTrueVertex*)(*gVertex->TrueVertices)[i];
        
        tVertex->Pur =
        (double)nMatch / nConst;  // nConstituents of the GlbV under study
        tVertex->Eff =
        (double)nMatch /
        nIDs;  // nIDS = n of true particle coming from this true vertex
    }
}

//*************************************************************
void ND::TGlobalReconModule::FillTrueParticle(
    ND::THandle<ND::TG4Trajectory> G4track, double pur, double eff,
    ND::TTrueParticle& trueParticle) 
{
    //*************************************************************
    
    if (G4track) {
        trueParticle.ID = G4track->GetTrackId();
        
        trueParticle.Pur = pur;
        trueParticle.Eff = eff;
        
        ND::TG4PrimaryVertex G4vertex;
        bool ok = GetG4Vertex(G4track, G4vertex);
        FillTrueVertex(ok, G4vertex, DEFAULT_MAX, DEFAULT_MAX, trueParticle.Vertex);
        
    } else {
        trueParticle.ID = -1;
        
        trueParticle.Pur = DEFAULT_MAX;
        trueParticle.Eff = DEFAULT_MAX;
        
        ND::TG4PrimaryVertex G4vertex;
        FillTrueVertex(false, G4vertex, DEFAULT_MAX, DEFAULT_MAX,
            trueParticle.Vertex);
    }
}

//*************************************************************
ND::THandle<ND::TG4Trajectory> ND::TGlobalReconModule::GetParent(
    ND::THandle<ND::TG4Trajectory> G4track) 
{
    //*************************************************************
    
    const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
    ND::THandle<ND::TG4TrajectoryContainer> trajectories =
    event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
    
    if (G4track->GetParentId() != 0) {
        ND::THandle<ND::TG4Trajectory> traj(
            new ND::TG4Trajectory(((*trajectories)[G4track->GetParentId()])));
        return traj;
    } else
        return ND::THandle<ND::TG4Trajectory>();
    
    /*
    // create a map with trajectories ordered in descending kinetic energy
    if(trajectories){
    ND::TG4TrajectoryContainer::iterator trajIter;
    for(trajIter = trajectories->begin();trajIter != trajectories->end();
    ++trajIter) {
    if ((*trajIter)->GetTrackId() == G4track->GetParentId()) return *trajIter;
    }
    }
    
    return ND::THandle<ND::TG4Trajectory>();
    */
}

//*************************************************************
void ND::TGlobalReconModule::FillTrueVertex(
    bool found, const ND::TG4PrimaryVertex& G4vertex, double pur, double eff,
    ND::TTrueVertex& vertex) 
{
    //*************************************************************
    
    if (found) {
        vertex.Position = G4vertex.GetPosition();
        vertex.ID = G4vertex.GetInteractionNumber();
        vertex.Pur = pur;
        vertex.Eff = eff;
    } else {
        vertex.Position =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        vertex.ID = -1;
        vertex.Pur = DEFAULT_MAX;
        vertex.Eff = DEFAULT_MAX;
    }
}

//*************************************************************
ND::THandle<ND::TG4PrimaryVertex> ND::TGlobalReconModule::GetG4Vertex(
    const ND::TReconBase&, double& pur, double& eff)
{
    //*************************************************************
    
    return ND::THandle<ND::TG4PrimaryVertex>();
}

//*************************************************************
bool ND::TGlobalReconModule::GetG4Vertex(ND::THandle<ND::TG4Trajectory> G4track,
    ND::TG4PrimaryVertex& G4vertex) 
{
    //*************************************************************
    
    // 1. Find the primary particle associated to this G4 trajectory
    
    const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
    ND::THandle<ND::TG4TrajectoryContainer> trajectories =
    event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
    if (!trajectories) return false;
    
    int nextparent;
    int count = 0;
    int primaryID = G4track->GetTrackId();
    
    while ((primaryID > 0) &&
        ((nextparent = (*trajectories)[primaryID].GetParentId()) > 0)) {
    primaryID = nextparent;
    if (count++ > 100000000) {
        // just to make sure it doesn't go into an infinite loop somehow
        ND280Error(
            "TGlobalReconModule: Particle has more than 100000000 levels of "
            "parents!");
        throw ND::EeventAnalysisInfiniteLoop2();
    }
        }
        
        // 2. Find the primary vertex associated to the primary particle
        ND::THandle<ND::TG4PrimaryVertexContainer> vertices =
        event.Get<ND::TG4PrimaryVertexContainer>("truth/G4PrimVertex00");
        if (!vertices) return false;
        
        std::vector<ND::TG4PrimaryVertex>::iterator it1;
        for (it1 = vertices->begin(); it1 != vertices->end(); it1++) {
            const ND::TG4PrimaryParticleContainer& particles =
            it1->GetPrimaryParticles();
            
            std::vector<ND::TG4PrimaryParticle>::const_iterator it2;
            for (it2 = particles.begin(); it2 != particles.end(); it2++) {
                if (it2->GetTrackId() == primaryID) {
                    G4vertex = *it1;
                    return true;
                }
            }
        }
        
        return false;
}

//*************************************************************
bool ND::TGlobalReconModule::GetIncomingParticle(
    const ND::TG4PrimaryVertex& G4vertex, ND::TG4PrimaryParticle& incoming)
{
    //*************************************************************
    
    for (std::vector<ND::TG4PrimaryVertex>::const_iterator infoVtxIter =
        G4vertex.GetInfoVertex().begin();
        infoVtxIter != G4vertex.GetInfoVertex().end(); ++infoVtxIter) {
    std::vector<ND::TG4PrimaryParticle>::const_iterator particleIter;
    for (particleIter = infoVtxIter->GetPrimaryParticles().begin();
        particleIter != infoVtxIter->GetPrimaryParticles().end();
        ++particleIter) {
    Int_t pdg = particleIter->GetPDGCode();
    if (particleIter->GetTrackId() == -1) {
        // is an incoming particle (including target)
        if (TMath::Abs(pdg) == 12 || TMath::Abs(pdg) == 14 ||
            TMath::Abs(pdg) == 16) {
        incoming = *particleIter;
        return true;
            }
    }
        }
        }
        
        return false;
}

//*************************************************************
int ND::TGlobalReconModule::GetNumberOfHits(
    ND::THandle<ND::TReconBase> object)
{
    //*************************************************************
    
    int nHits = 0;
    if (object->GetHits())
        nHits += object->GetHits()->size();
    else if (object->GetConstituents()) {
        ND::TReconObjectContainer::const_iterator it;
        for (it = object->GetConstituents()->begin();
            it != object->GetConstituents()->end(); it++) {
        nHits += GetNumberOfHits(*it);
            }
    }
    
    return nHits;
}

//*****************************************************
unsigned long ND::TGlobalReconModule::GetDetectorNumber(
    ND::THandle<ND::TReconBase> object)
{
//*****************************************************
    
    unsigned long det = 0;
    
    if (object->UsesDetector(ND::TReconBase::kFGD1)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kFGD1);
    if (object->UsesDetector(ND::TReconBase::kFGD2)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kFGD2);
    if (object->UsesDetector(ND::TReconBase::kTPC1)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kTPC1);
    if (object->UsesDetector(ND::TReconBase::kTPC2)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kTPC2);
    if (object->UsesDetector(ND::TReconBase::kTPC3))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kTPC3);
    if (object->UsesDetector(ND::TReconBase::kSFG)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kSFG);
    if (object->UsesDetector(ND::TReconBase::kP0D))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kP0D);
    if (object->UsesDetector(ND::TReconBase::kDSECal)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kDsECal);
    if (object->UsesDetector(ND::TReconBase::kTopPECal))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kP0DECalTop);
    if (object->UsesDetector(ND::TReconBase::kBottomPECal))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kP0DECalBottom);
    if (object->UsesDetector(ND::TReconBase::kLeftPECal))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kP0DECalLeft);
    if (object->UsesDetector(ND::TReconBase::kRightPECal)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kP0DECalRight);
    if (object->UsesDetector(ND::TReconBase::kTopTECal))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kBrlECalTop);
    if (object->UsesDetector(ND::TReconBase::kBottomTECal))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kBrlECalBottom);
    if (object->UsesDetector(ND::TReconBase::kLeftTECal))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kBrlECalLeft);
    if (object->UsesDetector(ND::TReconBase::kRightTECal)) 
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kBrlECalRight);
    if (object->UsesDetector(ND::TReconBase::kSMRD))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kSMRD);
    if (object->UsesDetector(ND::TReconBase::kTopHAT))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kTopHAT);
    if (object->UsesDetector(ND::TReconBase::kBottomHAT))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kBottomHAT);
    if (object->UsesDetector(ND::TReconBase::kTOF))
      eventAnalysisEnums::SetDetectorUsed(det,eventAnalysisEnums::kTOF);
    return det;
}

int ND::TGlobalReconModule::CountDetector(
    ND::THandle<ND::TReconBase> object){ 
  int count = 0;
  if (object->UsesDetector(ND::TReconBase::kTPC1))        ++count;
  if (object->UsesDetector(ND::TReconBase::kTPC2))        ++count;
  if (object->UsesDetector(ND::TReconBase::kTPC3))        ++count;
  if (object->UsesDetector(ND::TReconBase::kFGD1))        ++count;
  if (object->UsesDetector(ND::TReconBase::kFGD2))        ++count;
  if (object->UsesDetector(ND::TReconBase::kSFG))         ++count;
  if (object->UsesDetector(ND::TReconBase::kTopHAT))      ++count;
  if (object->UsesDetector(ND::TReconBase::kBottomHAT))   ++count;
  if (object->UsesDetector(ND::TReconBase::kTOF))         ++count;
  if (object->UsesDetector(ND::TReconBase::kP0D))         ++count;
  if (object->UsesDetector(ND::TReconBase::kDSECal))      ++count;
  if (object->UsesDetector(ND::TReconBase::kTopPECal))    ++count;
  if (object->UsesDetector(ND::TReconBase::kBottomPECal)) ++count;
  if (object->UsesDetector(ND::TReconBase::kLeftPECal))   ++count;
  if (object->UsesDetector(ND::TReconBase::kRightPECal))  ++count;
  if (object->UsesDetector(ND::TReconBase::kTopTECal))    ++count;
  if (object->UsesDetector(ND::TReconBase::kBottomTECal)) ++count;
  if (object->UsesDetector(ND::TReconBase::kLeftTECal))   ++count;
  if (object->UsesDetector(ND::TReconBase::kRightTECal))  ++count;
  if (object->UsesDetector(ND::TReconBase::kSMRD))        ++count;
  return count;
}


//*****************************************************
void ND::TGlobalReconModule::FillDetectorUsed(
    ND::THandle<ND::TReconBase> object, bool dets[]) 
{
//*****************************************************
    
    // Fill the used detectors. This should be consistent with 
    // psycheND280Utils/TreeConverterUtils::SetDetectorUsedField.
    
    for (int i = 0; i < NDETSUSED; i++) dets[i] = false;
    
    if (object->UsesDetector(ND::TReconBase::kFGD1))        dets[0]  = true;
    if (object->UsesDetector(ND::TReconBase::kFGD2))        dets[1]  = true;
    if (object->UsesDetector(ND::TReconBase::kTPC1))        dets[2]  = true;
    if (object->UsesDetector(ND::TReconBase::kTPC2))        dets[3]  = true;
    if (object->UsesDetector(ND::TReconBase::kTPC3))        dets[4]  = true;
    if (object->UsesDetector(ND::TReconBase::kSFG))         dets[5]  = true;
    if (object->UsesDetector(ND::TReconBase::kTopHAT))      dets[6]  = true;
    if (object->UsesDetector(ND::TReconBase::kBottomHAT))   dets[7]  = true;
    if (object->UsesDetector(ND::TReconBase::kTOF))         dets[8]  = true;
    if (object->UsesDetector(ND::TReconBase::kP0D))         dets[9]  = true;
    if (object->UsesDetector(ND::TReconBase::kDSECal))      dets[10] = true;
    if (object->UsesDetector(ND::TReconBase::kTopTECal))    dets[11] = true;
    if (object->UsesDetector(ND::TReconBase::kBottomTECal)) dets[12] = true;
    if (object->UsesDetector(ND::TReconBase::kLeftTECal))   dets[13] = true;
    if (object->UsesDetector(ND::TReconBase::kRightTECal))  dets[14] = true;
    if (object->UsesDetector(ND::TReconBase::kTopPECal))    dets[15] = true;
    if (object->UsesDetector(ND::TReconBase::kBottomPECal)) dets[16] = true;
    if (object->UsesDetector(ND::TReconBase::kLeftPECal))   dets[17] = true;
    if (object->UsesDetector(ND::TReconBase::kRightPECal))  dets[18] = true;
    if (object->UsesDetector(ND::TReconBase::kTopSMRD))     dets[19] = true;
    if (object->UsesDetector(ND::TReconBase::kBottomSMRD))  dets[20] = true;
    if (object->UsesDetector(ND::TReconBase::kLeftSMRD))    dets[21] = true;
    if (object->UsesDetector(ND::TReconBase::kRightSMRD))   dets[22] = true;
    
}

//*****************************************************
void ND::TGlobalReconModule::FillTPCOther(ND::TND280Event& event) 
{
    //*****************************************************
    
    ND::THandle<ND::TAlgorithmResult> tpcRecon = event.GetFit("tpcRecon");
    if (tpcRecon) {
        ND::THandle<ND::TReconObjectContainer> otherObjects =
        tpcRecon->GetResultsContainer("TPCother");
        if (otherObjects) {
            for (ND::TReconObjectContainer::iterator tt = otherObjects->begin();
                tt != otherObjects->end(); tt++) {
            ND::THandle<ND::TReconBase> other = *tt;
            
            if (other && fNTPCOthers < NMAXTPCOTHER) {
                TTPCOtherObject* tpcOther =
                new ((*fTPCOthers)[fNTPCOthers++]) TTPCOtherObject;
                
                tpcOther->Detector = GetDetectorNumber(other);
                tpcOther->Chi2 = other->GetQuality();
                tpcOther->EDeposit = TrackingUtils::GetEDeposit(other);
                tpcOther->Charge = TrackingUtils::GetCharge(other);
                tpcOther->NHits = DEFAULT_MAX;
                if (other->GetHits()) tpcOther->NHits = other->GetHits()->size();
                
                // Set the back and front position, variance and direction
                ND::THandle<ND::TReconState> firstState =
                TrackingUtils::GetFirstState(*other);
                if (firstState) {
                    tpcOther->FrontPosition = TrackingUtils::GetPosition(firstState);
                    tpcOther->FrontDirection = TrackingUtils::GetDirection(firstState);
                    tpcOther->Momentum = TrackingUtils::GetMomentum(firstState);
                } else {
                    tpcOther->FrontPosition = TLorentzVector(DEFAULT_MAX, DEFAULT_MAX,
                        DEFAULT_MAX, DEFAULT_MAX);
                    tpcOther->FrontDirection =
                    TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
                    tpcOther->Momentum = DEFAULT_MAX;
                }
                
                ND::THandle<ND::TReconState> lastState =
                TrackingUtils::GetLastState(*other);
                if (lastState) {
                    tpcOther->BackPosition = TrackingUtils::GetPosition(lastState);
                    tpcOther->BackDirection = TrackingUtils::GetDirection(lastState);
                } else {
                    tpcOther->BackPosition = TLorentzVector(DEFAULT_MAX, DEFAULT_MAX,
                        DEFAULT_MAX, DEFAULT_MAX);
                    tpcOther->BackDirection =
                    TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
                }
                
                // Fill the true particle when it is found
                double pur, eff;
                ND::THandle<ND::TG4Trajectory> G4track =
                GetG4Trajectory(*other, pur, eff);
                FillTrueParticle(G4track, pur, eff, tpcOther->TrueParticle);
            }
                }
        }
    }
}

//*****************************************************
void ND::TGlobalReconModule::FillUnusedHits(ND::TND280Event& event) 
{
    //*****************************************************
    
    // Now add TPC hits
    ND::THandle<ND::TAlgorithmResult> tpcRecon = event.GetFit("tpcRecon");
    if (tpcRecon) {
        // Get TPC unused hits
        ND::THandle<ND::THitSelection> tpcHits =
        tpcRecon->GetHitSelection("TPCunused");
        if (tpcHits) {
            fNTPCUnused = tpcHits->size();
        }
    }
    
    // Now add P0D hits
    ND::THandle<ND::TAlgorithmResult> p0dRecon = event.GetFit("p0dRecon");
    if (p0dRecon) {
        // Get P0D unused hits
        ND::THandle<ND::THitSelection> p0dHits =
        p0dRecon->GetHitSelection("unused");
        if (p0dHits) {
            for (ND::THitSelection::iterator hit = p0dHits->begin();
                hit != p0dHits->end(); hit++) {
            TGlobalHit* gHit = new ((*fP0DUnused)[fNP0DUnused++]) TGlobalHit;
            FillGlobalHit(*hit, *gHit);
                }
        }
    }
    
    ND::THandle<ND::TAlgorithmResult> smrdRecon = event.GetFit("smrdRecon");
    if (smrdRecon) {
        // Get SMRD unused hits
        ND::THandle<ND::THitSelection> smrdHits =
        smrdRecon->GetHitSelection("unused");
        if (smrdHits) {
            for (ND::THitSelection::iterator h = smrdHits->begin();
                h != smrdHits->end(); h++) {
            // fill info for smrd unused hit
            TSMRDHit* smrdHit = NULL;
            smrdHit = new ((*fSMRDUnused)[fNSMRDUnused++]) TSMRDHit;
            FillSmrdHit(*h, *smrdHit);
            
            if (ND::TGeomInfo::SMRD().IsSMRDTopWall((*h))) {
                fNSMRDTopUnused++;
            } else if (ND::TGeomInfo::SMRD().IsSMRDBottomWall((*h))) {
                fNSMRDBottomUnused++;
            } else if (ND::TGeomInfo::SMRD().IsSMRDLeftWall((*h))) {
                fNSMRDLeftUnused++;
            } else if (ND::TGeomInfo::SMRD().IsSMRDRightWall((*h))) {
                fNSMRDRightUnused++;
            }
                }
        }
    }
}

//*****************************************************
void ND::TGlobalReconModule::FillFirstLastHits(ND::THitSelection& hits,
    TGlobalPID& globalObject) 
{
    //*****************************************************
    
    // sort hits in increasing Z
    std::sort(hits.begin(), hits.end(), SortHitsInZ);
    
    TGlobalHit* hitSaved[10];
    bool firstX = false;
    bool lastX = false;
    bool firstY = false;
    bool lastY = false;
    
    for (ND::THitSelection::iterator hit = hits.begin(); hit != hits.end();
        hit++) {
    if (globalObject.NHitsSaved >= 10 || (firstX && firstY)) break;
    
    if ((*hit)->IsXHit() && !firstX) {
        hitSaved[globalObject.NHitsSaved] =
        new ((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) TGlobalHit;
        FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]);
        globalObject.NHitsSaved++;
        firstX = true;
        if ((*hit)->IsYHit()) firstY = true;
    } else if ((*hit)->IsYHit() && !firstY) {
        hitSaved[globalObject.NHitsSaved] =
        new ((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) TGlobalHit;
        FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]);
        globalObject.NHitsSaved++;
        firstY = true;
        if ((*hit)->IsXHit()) firstX = true;
    }
        }
        
        for (ND::THitSelection::reverse_iterator hit = hits.rbegin();
            hit != hits.rend(); hit++) {
        if (globalObject.NHitsSaved >= 10 || (lastX && lastY)) break;
        
        if ((*hit)->IsXHit() && !lastX) {
            hitSaved[globalObject.NHitsSaved] =
            new ((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) TGlobalHit;
            FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]);
            globalObject.NHitsSaved++;
            lastX = true;
            if ((*hit)->IsYHit()) lastY = true;
        } else if ((*hit)->IsYHit() && !lastY) {
            hitSaved[globalObject.NHitsSaved] =
            new ((*(globalObject.HitsSaved))[globalObject.NHitsSaved]) TGlobalHit;
            FillGlobalHit(*hit, *hitSaved[globalObject.NHitsSaved]);
            globalObject.NHitsSaved++;
            lastY = true;
            if ((*hit)->IsXHit()) lastX = true;
        }
            }
}

//*****************************************************
void ND::TGlobalReconModule::FillAllHits(
    ND::THitSelection& hits,
    TGlobalPID&        globalObject) {
//*****************************************************
  for(ND::THitSelection::iterator hit = hits.begin(); 
      hit != hits.end(); hit++) {
    TGlobalHit* gHit =
      new ((*(globalObject.HitsSaved))[globalObject.NHitsSaved++])TGlobalHit;
    FillGlobalHit(*hit, *gHit);
  }
}

//*****************************************************
void ND::TGlobalReconModule::FillOutermostHits(ND::THitSelection& hits,
    double charge_cut,
    TOutermostHits& outer) 
{
    //*****************************************************
    
    double minX = 1000000;
    double maxX = -1000000;
    double minY = 1000000;
    double maxY = -1000000;
    double minZ = 1000000;
    double maxZ = -1000000;
    
    ND::THandle<ND::THit> hitSave[6];
    
    for (ND::THitSelection::iterator hit = hits.begin(); hit != hits.end();
        hit++) {
    if ((*hit)->GetCharge() < charge_cut) continue;
    
    if ((*hit)->IsXHit()) {
        if ((*hit)->GetPosition().X() > maxX) {
            maxX = (*hit)->GetPosition().X();
            hitSave[0] = *hit;
        }
        if ((*hit)->GetPosition().X() < minX) {
            minX = (*hit)->GetPosition().X();
            hitSave[1] = *hit;
        }
    }
    if ((*hit)->IsYHit()) {
        if ((*hit)->GetPosition().Y() > maxY) {
            maxY = (*hit)->GetPosition().Y();
            hitSave[2] = *hit;
        }
        if ((*hit)->GetPosition().Y() < minY) {
            minY = (*hit)->GetPosition().Y();
            hitSave[3] = *hit;
        }
    }
    if ((*hit)->IsZHit()) {
        if ((*hit)->GetPosition().Z() > maxZ) {
            maxZ = (*hit)->GetPosition().Z();
            hitSave[4] = *hit;
        }
        if ((*hit)->GetPosition().Z() < minZ) {
            minZ = (*hit)->GetPosition().Z();
            hitSave[5] = *hit;
        }
    }
        }
        
        FillGlobalHit(hitSave[0], outer.hitMaxX);
        FillGlobalHit(hitSave[1], outer.hitMinX);
        FillGlobalHit(hitSave[2], outer.hitMaxY);
        FillGlobalHit(hitSave[3], outer.hitMinY);
        FillGlobalHit(hitSave[4], outer.hitMaxZ);
        FillGlobalHit(hitSave[5], outer.hitMinZ);
}


//*****************************************************
void ND::TGlobalReconModule::FillSfgTimeBins(ND::TND280Event& event) 
{
//*****************************************************
    
  
  TTrackerReconUtils trackerUtils;
  fNSfgTimeBins = trackerUtils.GetSFGTimeBinNumber(event);
  for (int ibin = 0; ibin < fNSfgTimeBins; ibin++) {
    TSfgTimeBin* bin = new ((*fSfgTimeBins)[ibin]) TSfgTimeBin;
    ND::THandle<ND::TAlgorithmResult> tbResult
      = trackerUtils.GetSFGTimeBinResult(event, ibin);
    
    // Initialize the bin info
    bin->minTime           = 2.e10;
    bin->maxTime           = -2.e10;
    bin->nHits             = 0;
    bin->rawChargeSum      = 0.;
    bin->chargeWeightedPos = TVector3(0., 0., 0.);
    bin->NSFGUnused        = 0;
    bin->g4ID              = 0;
    bin->SFGUnused->Clear();
    bin->SFGHitPositionX.clear();
    bin->SFGHitPositionY.clear();
    bin->SFGHitPositionZ.clear();
    bin->SFGHitT        .clear();
    bin->SFGHitQ        .clear();
    bin->SFGHitTrajID   .clear();
    
    if (tbResult) {
      std::tie(bin->minTime, bin->maxTime) 
        = trackerUtils.GetTimeBinBoundary(tbResult);
      
      // Fill the SFG unused hits
      ND::THandle<ND::THitSelection> sfgUnused =
      tbResult->GetHitSelection("unused");
      if (sfgUnused) {
        for (ND::THitSelection::iterator hit = sfgUnused->begin();
            hit != sfgUnused->end(); hit++) {
          TGlobalHit* gHit = NULL;
          if (ND::GeomId::SFG::IsCube((*hit)->GetGeomId()))
              gHit = new ((*bin->SFGUnused)[bin->NSFGUnused++]) TGlobalHit;
          else
              continue;
          
          FillGlobalHit(*hit, *gHit);
        }
      }
      // Cumulative result for all time bins
      fNSFGUnused += bin->NSFGUnused;
      
      
      // Fill information on individual SFG hits
      ND::THandle<ND::THitSelection> tbHits = tbResult->GetHitSelection("used");
      if (tbHits) {
        bin->nHits = tbHits->size();
        for (ND::THitSelection::iterator hit = tbHits->begin();
             hit != tbHits->end(); hit++) {
          // Fill the information about individual SFG hit
          TVector3 pos = (*hit)->GetPosition();
          bin->SFGHitPositionX.push_back(pos.X());
          bin->SFGHitPositionY.push_back(pos.Y());
          bin->SFGHitPositionZ.push_back(pos.Z());
          bin->SFGHitT        .push_back((*hit)->GetTime());
          bin->SFGHitQ        .push_back((*hit)->GetCharge());
          bin->rawChargeSum += (*hit)->GetCharge();
            ND::THitSelection thisHit;
            thisHit.push_back(*hit);
            bin->SFGHitTrajID.push_back(TrackTruthInfo::GetG4TrajIDHits(thisHit));
        }  // End loop over hits
      } // Finished filling information on individual SFG hits
    }
  }
}


//*****************************************************
void ND::TGlobalReconModule::FillFgdTimeBins(ND::TND280Event& event) 
{
    //*****************************************************
    
    bool isMC = event.GetContext().IsMC();
    
    fNFgdTimeBins = fgdUtils::getNumberTimeBins(event);
    for (int ibin = 0; ibin < fNFgdTimeBins; ibin++) {
        TFgdTimeBin* bin = new ((*fFgdTimeBins)[ibin]) TFgdTimeBin;
        ND::THandle<ND::TAlgorithmResult> tbResult =
        fgdUtils::getTimeBinResult(event, ibin);
        
        bin->minTime = 2.e10;
        bin->maxTime = -2.e10;
        
        for (int ifgd = 0; ifgd < 2; ifgd++) {
            bin->nHits[ifgd] = 0;
            bin->rawChargeSum[ifgd] = 0.;
            bin->chargeWeightedPos[ifgd] = TVector3(0., 0., 0.);
            for (int ilayer = 0; ilayer < 30; ilayer++) {
                bin->chargePerLayer[ifgd][ilayer] = 0.;
            }
        }
        
        bin->NFGD1Unused = 0;
        bin->NFGD2Unused = 0;
        
        bin->FGD1Unused->Clear();
        bin->FGD2Unused->Clear();
        
        if (tbResult) {
            ND::THandle<ND::TRealDatum> binStartTime =
            tbResult->Get<ND::TRealDatum>("binStartTime");
            if (binStartTime) {
                bin->minTime = binStartTime->GetValue();
            }
            
            ND::THandle<ND::TRealDatum> binEndTime =
            tbResult->Get<ND::TRealDatum>("binEndTime");
            if (binEndTime) {
                bin->maxTime = binEndTime->GetValue();
            }
            
            // Fill the FGD unused hits
            ND::THandle<ND::THitSelection> fgdUnused =
            tbResult->GetHitSelection("unused");
            if (fgdUnused) {
                for (ND::THitSelection::iterator hit = fgdUnused->begin();
                    hit != fgdUnused->end(); hit++) {
                TGlobalHit* gHit = NULL;
                if (ND::TGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()))
                    gHit = new ((*bin->FGD1Unused)[bin->NFGD1Unused++]) TGlobalHit;
                else if (ND::TGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))
                    gHit = new ((*bin->FGD2Unused)[bin->NFGD2Unused++]) TGlobalHit;
                else
                    continue;
                
                FillGlobalHit(*hit, *gHit);
                    }
            }
            
            // Also include hits from 2D fgdIsoRecon tracks for now
            ND::THandle<ND::TReconObjectContainer> xz =
            tbResult->GetResultsContainer("xzUnmatchedFgdReconTracks");
            if (xz) {
                for (ND::TReconObjectContainer::iterator tt = xz->begin();
                    tt != xz->end(); tt++) {
                ND::THandle<ND::TReconBase> object = *tt;
                if (object) {
                    ND::THandle<ND::THitSelection> hits = object->GetHits();
                    if (hits) {
                        for (ND::THitSelection::iterator hit = hits->begin();
                            hit != hits->end(); hit++) {
                        TGlobalHit* gHit = NULL;
                        if (ND::TGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()))
                            gHit =
                        new ((*bin->FGD1Unused)[bin->NFGD1Unused++]) TGlobalHit;
                        else if (ND::TGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))
                            gHit =
                        new ((*bin->FGD2Unused)[bin->NFGD2Unused++]) TGlobalHit;
                        else
                            continue;
                        
                        FillGlobalHit(*hit, *gHit);
                            }
                    }
                }
                    }
            }
            ND::THandle<ND::TReconObjectContainer> yz =
            tbResult->GetResultsContainer("yzUnmatchedFgdReconTracks");
            if (yz) {
                for (ND::TReconObjectContainer::iterator tt = yz->begin();
                    tt != yz->end(); tt++) {
                ND::THandle<ND::TReconBase> object = *tt;
                if (object) {
                    ND::THandle<ND::THitSelection> hits = object->GetHits();
                    if (hits) {
                        for (ND::THitSelection::iterator hit = hits->begin();
                            hit != hits->end(); hit++) {
                        TGlobalHit* gHit = NULL;
                        if (ND::TGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()))
                            gHit =
                        new ((*bin->FGD1Unused)[bin->NFGD1Unused++]) TGlobalHit;
                        else if (ND::TGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))
                            gHit =
                        new ((*bin->FGD2Unused)[bin->NFGD2Unused++]) TGlobalHit;
                        else
                            continue;
                        
                        FillGlobalHit(*hit, *gHit);
                            }
                    }
                }
                    }
            }
            
            // Cumulative result for all time bins
            fNFGD1Unused += bin->NFGD1Unused;
            fNFGD2Unused += bin->NFGD2Unused;
            
            ND::THandle<ND::THitSelection> tbHits = tbResult->GetHitSelection("fgd");
            if (tbHits) {
                // Extract FGD1 and FGD2 hits
                ND::THitSelection hits[2];
                
                TVector3 chargeSum[2];
                for (ND::THitSelection::iterator hit = tbHits->begin();
                    hit != tbHits->end(); hit++) {
                // Fill the information about individual FGD hit
                double charge = (*hit)->GetCharge();
                TVector3 pos = (*hit)->GetPosition();
                int layer = ND::TGeomInfo::FGD().ActivePlane(pos.Z());
                float hit_pos = 0;
                
                // Only need to store on coordinate, since these are 2D hits
                if ((*hit)->IsXHit()) {
                    hit_pos = pos.X();
                } else {
                    hit_pos = pos.Y();
                }
                
                float time = (*hit)->GetTime();
                
                // Save separate information for FGD1 and FGD2
                if (ND::TGeomInfo::FGD().IsInFGD1(pos)) {
                    bin->FGD1HitLayer.push_back(layer);
                    bin->FGD1HitPosition.push_back(hit_pos);
                    bin->FGD1HitQ.push_back(charge);
                    bin->FGD1HitT.push_back(time);
                    ND::THitSelection thisHit;
                    thisHit.push_back(*hit);
                    bin->FGD1HitTrajID.push_back(
                        TrackTruthInfo::GetG4TrajIDHits(thisHit));
                } else {
                    bin->FGD2HitLayer.push_back(layer);
                    bin->FGD2HitPosition.push_back(hit_pos);
                    bin->FGD2HitQ.push_back(charge);
                    bin->FGD2HitT.push_back(time);
                    ND::THitSelection thisHit;
                    thisHit.push_back(*hit);
                    bin->FGD2HitTrajID.push_back(
                        TrackTruthInfo::GetG4TrajIDHits(thisHit));
                }
                // Finished filling information on individual FGD hits
                
                TVector3 isHit((*hit)->IsXHit(), (*hit)->IsYHit(), (*hit)->IsZHit());
                //    std::cout << "isHit = (" << isHit.X() << "," << isHit.Y() <<
                //"," << isHit.Z() << ")" << std::endl;
                bool isInFgd[2];
                isInFgd[0] = ND::TGeomInfo::FGD().IsInFGD1((*hit)->GetPosition());
                isInFgd[1] = ND::TGeomInfo::FGD().IsInFGD2((*hit)->GetPosition());
                for (int ifgd = 0; ifgd < 2; ifgd++) {
                    //      std::cout << "isInFgd[" << ifgd << "] = " << isInFgd[ifgd]
                    //<< std::endl;
                    if (isInFgd[ifgd]) {
                        hits[ifgd].push_back(*hit);
                        bin->rawChargeSum[ifgd] += (*hit)->GetCharge();
                        
                        // Fill chargePerLayer
                        ND::TGeometryId geomId = (*hit)->GetGeomId();
                        int xory = ND::TGeomInfo::FGD().GetXorY(geomId);
                        int module = ND::TGeomInfo::FGD().GetModule(geomId);
                        int layer = module * 2 + xory;
                        bin->chargePerLayer[ifgd][layer] += (*hit)->GetCharge();
                        // std::cout << "Z, XorY, Module, Layer = " << pos.Z() << ", "
                        // <<
                        // xory << ", " << module << ", " << layer << std::endl;
                        
                        for (int ix = 0; ix < 3; ix++) {
                            if (isHit(ix)) {
                                //      std::cout << "filling hit for component " << ix <<
                                //" in FGD" << ifgd+1 << std::endl;
                                bin->chargeWeightedPos[ifgd](ix) += charge * pos(ix);
                                chargeSum[ifgd](ix) += charge;
                            }
                        }
                    }
                }
                    }
                    
                    for (int ifgd = 0; ifgd < 2; ifgd++) {
                        for (int ix = 0; ix < 3; ix++) {
                            if (chargeSum[ifgd](ix) > 1.e-3) {
                                bin->chargeWeightedPos[ifgd](ix) *= 1. / chargeSum[ifgd](ix);
                            } else {
                                bin->chargeWeightedPos[ifgd](ix) = -1.e-10;
                            }
                        }
                        bin->nHits[ifgd] = hits[ifgd].size();
                    }
                    
                    // Fill outermost hits for this time bin
                    double charge_cut = 0;
                    FillOutermostHits(hits[0], charge_cut, bin->FGD1OutermostHits);
                    FillOutermostHits(hits[1], charge_cut, bin->FGD2OutermostHits);
                    
                    if (isMC)
                        bin->g4ID = TrackTruthInfo::GetG4TrajIDHits(*tbHits);
                    else
                        bin->g4ID = 0;
            }
        }
    }
}

// Helper method
// Get the median FGD hit time for this track.
//*************************************************************
double ND::TGlobalReconModule::MedianObjectTime(
    ND::THandle<ND::TReconBase> object)
{
    //*************************************************************
    std::vector<double> times;
    ND::THandle<ND::THitSelection> hits = object->GetHits();
    if (hits)
    {
        // make a vector of all the times for the 'good' hits
        for (ND::THitSelection::iterator hit = (hits)->begin();
            hit != (hits)->end(); hit++)
        {
            if ((ND::TGeomInfo::FGD().IsInFGD1((*hit)->GetPosition()) ||
                ND::TGeomInfo::FGD().IsInFGD2((*hit)->GetPosition()))) {
            times.push_back((*hit)->GetTime());
                }
        }
    }
    std::sort(times.begin(), times.end());
    
    unsigned len = times.size();
    
    if (len == 0)
        // return 0.0;
    return DEFAULT_MIN;
    else if (len % 2 != 0) {
        return times[(len - 1) / 2];
    }
    
    ND280Debug(" end median track ");
    
    return (times[len / 2] + times[len / 2 - 1]) / 2.0;
}

//*************************************************************
void ND::TGlobalReconModule::GetFGDSimpleVA(ND::TND280Event& event,
    ND::THandle<ND::TReconBase>& object,
    TLorentzVector& vertexPos,
    int globpass, int vaea) 
{
    //*************************************************************
    
    // time bin loop - match reco. time bin to track to get associated fgd hit
    // selection
    int fgdRecoBinIndex = -1;
    if (!fgdUtils::GetTrackRecoTimeBinIndex(event, object, fgdRecoBinIndex))
        return;
    
    ND::THandle<ND::TAlgorithmResult> binRecon =
    fgdUtils::getTimeBinResult(event, fgdRecoBinIndex);
    if (!binRecon) return;
    
    ND::THandle<ND::THitSelection> fgdHits = binRecon->GetHitSelection("fgd");
    if (!fgdHits) return;
    
    // run vertex activity measurement on track
    fgdUtils::SimpleVertexActivity(0, fgdHits, object, vertexPos, globpass, vaea);
    
    return;
}

//*************************************************************
void ND::TGlobalReconModule::FillTPCPID(ND::THandle<ND::TReconBase> object) 
{
    //*************************************************************
    
    unsigned int dets = GetDetectorNumber(object);
    
    // create a new TGlobalPID
    TTpcPID* tpcObject = new ((*fTPCPIDs)[fNTPCPIDs++]) TTpcPID;
    
    // Fill the true particle when it is found
    double pur, eff;
    ND::THandle<ND::TG4Trajectory> G4track = GetG4Trajectory(*object, pur, eff);
    FillTrueParticle(G4track, pur, eff, tpcObject->TrueParticle);
    
    // extrapolate to the vertex (True vertex for the moment !!!!)
    tpcObject->PositionAtTrueVertex =
    TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
    tpcObject->PositionVarAtTrueVertex =
    TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
    tpcObject->DirectionAtTrueVertex =
    TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
    //------- copy information from TReconPID to TGlobalPID --------
    
    tpcObject->AlgorithmName = object->GetAlgorithmName();
    tpcObject->Status = object->CheckStatus(object->kSuccess);
    tpcObject->Chi2 = object->GetQuality();
    tpcObject->NDOF = object->GetNDOF();
    tpcObject->NNodes = object->GetNodes().size();
    tpcObject->Detectors = dets;
    
    tpcObject->NHits = DEFAULT_MAX;
    if (object->GetHits()) tpcObject->NHits = object->GetHits()->size();
    
    // The number of constituents. For composite objects
    tpcObject->NConstituents = 0;
    if (object->GetConstituents())
        tpcObject->NConstituents = object->GetConstituents()->size();
    
    tpcObject->isForward = (TrackingUtils::GetDirection(object).Z() > 0);
    tpcObject->MomentumAtTrueVertex = TrackingUtils::GetMomentum(object);
    tpcObject->MomentumErrorAtTrueVertex =
    TrackingUtils::GetMomentumError(object);
    tpcObject->Charge = TrackingUtils::GetCharge(object);
    tpcObject->EDeposit = TrackingUtils::GetEDeposit(object);
    tpcObject->Width =
    TrackingUtils::GetWidth(object).X();  //! take only the first component
    tpcObject->Cone = TrackingUtils::GetCone(object);
    tpcObject->Length = DEFAULT_MAX;
    
    // Set the PID back and front position, variance, direction and momentum
    ND::THandle<ND::TReconState> firstState =
    TrackingUtils::GetFirstState(*object);
    if (firstState) {
        tpcObject->FrontPosition = TrackingUtils::GetPosition(firstState);
        tpcObject->FrontPositionVar =
        TrackingUtils::GetPositionVariance(firstState);
        tpcObject->FrontDirection = TrackingUtils::GetDirection(firstState);
        tpcObject->FrontDirectionVar =
        TrackingUtils::GetDirectionVariance(firstState);
        tpcObject->FrontMomentum = TrackingUtils::GetMomentum(firstState);
        tpcObject->FrontMomentumError = TrackingUtils::GetMomentumError(firstState);
    } else {
        tpcObject->FrontPosition =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->FrontPositionVar =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->FrontDirection = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->FrontDirectionVar =
        TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->FrontMomentum = DEFAULT_MAX;
        tpcObject->FrontMomentumError = DEFAULT_MAX;
    }
    
    ND::THandle<ND::TReconState> lastState = TrackingUtils::GetLastState(*object);
    if (lastState) {
        tpcObject->BackPosition = TrackingUtils::GetPosition(lastState);
        tpcObject->BackPositionVar = TrackingUtils::GetPositionVariance(lastState);
        tpcObject->BackDirection = TrackingUtils::GetDirection(lastState);
        tpcObject->BackDirectionVar =
        TrackingUtils::GetDirectionVariance(lastState);
        tpcObject->BackMomentum = TrackingUtils::GetMomentum(lastState);
        tpcObject->BackMomentumError = TrackingUtils::GetMomentumError(lastState);
    } else {
        tpcObject->BackPosition =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->BackPositionVar =
        TLorentzVector(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->BackDirection = TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->BackDirectionVar =
        TVector3(DEFAULT_MAX, DEFAULT_MAX, DEFAULT_MAX);
        tpcObject->BackMomentum = DEFAULT_MAX;
        tpcObject->BackMomentumError = DEFAULT_MAX;
    }
    
    // Set the PID information
    ND::THandle<ND::TReconPID> PID = object;
    if (PID) {
        tpcObject->Likelihoods.push_back(PID->GetPIDWeight());
        ND::TReconObjectContainer::const_iterator it;
        for (it = PID->GetAlternates().begin(); it != PID->GetAlternates().end();
            it++) {
        ND::THandle<ND::TReconPID> alter = *it;
        tpcObject->Likelihoods.push_back(alter->GetPIDWeight());
            }
    }
}

//*****************************************************
ND::THandle<ND::TG4Trajectory> ND::TGlobalReconModule::GetG4Trajectory(
    const ND::TReconBase& object, double& purity, double& eff)
{
    //*****************************************************
    
    purity = -0xABCDEF;
    eff = -0xABCDEF;
    
    //  std::cout<<"****** GET ASSOCIATE TRAJECTORY TO THE TRACK ******
    //  "<<std::endl;
    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 = G4trajectories->begin();
    ND::TG4TrajectoryContainer::iterator trajEnd = G4trajectories->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>();
}

//*****************************************************
void ND::TGlobalReconModule::UpdateCoincidences(ND::THandle<ND::TMCHit> mch,
    ND::TMCDigit* mcd,
    std::vector<int>& coinc,
    int& nhits) 
{
    //*****************************************************
    
    std::vector<TG4VHit*> contrib;
    if (mch)
        contrib = mch->GetContributors();
    else if (mcd)
        contrib = mcd->GetContributors();
    else
        return;
    
    //  std::cout << "loop over hit contribuitors" << std::endl;
    for (std::vector<ND::TG4VHit*>::const_iterator h = contrib.begin();
        h != contrib.end(); ++h) {
    if (!*h) continue;
    ND::TG4HitSegment* g4hitseg = dynamic_cast<ND::TG4HitSegment*>(*h);
    if (g4hitseg != NULL) {
        for (unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
            unsigned int vv = g4hitseg->GetContributors()[i];
            if (vv >= coinc.size()) coinc.resize(vv + 1);
            coinc[vv]++;
            nhits++;
        }
    }
        }
}

//*****************************************************
ND::THandle<ND::TG4Trajectory> ND::TGlobalReconModule::GetMainContributor(
    ND::THandle<ND::TMCHit> mch,
    ND::THandle<ND::TG4TrajectoryContainer> trajectories) 
{
    //*****************************************************
    
    if (!mch) return ND::THandle<ND::TG4Trajectory>();
    if (!trajectories) return ND::THandle<ND::TG4Trajectory>();
    
    std::vector<int> coinc;
    coinc.resize(trajectories->size());
    
    //  std::cout << "loop over hit contribuitors" << std::endl;
    for (std::vector<ND::TG4VHit*>::const_iterator h =
        (mch->GetContributors()).begin();
        h != (mch->GetContributors()).end(); ++h) {
    if (!*h) continue;
    ND::TG4HitSegment* g4hitseg = dynamic_cast<ND::TG4HitSegment*>(*h);
    if (g4hitseg != NULL) {
        for (unsigned int i = 0; i != g4hitseg->GetContributors().size(); ++i) {
            unsigned int vv = g4hitseg->GetContributors()[i];
            if (vv >= coinc.size()) coinc.resize(vv + 1);
            coinc[vv]++;
        }
    }
        }
        
        int maxim = 0;
        int imax = -1;
        for (unsigned int i = 0; i < coinc.size(); i++) {
            if (coinc[i] > maxim) {
                maxim = coinc[i];
                imax = i;
            }
        }
        
        if (imax >= 0) {
            ND::THandle<ND::TG4Trajectory> traj(
                new ND::TG4Trajectory((*trajectories)[imax]));
            return traj;
        }
        
        return ND::THandle<ND::TG4Trajectory>();
}

//*****************************************************
void ND::TGlobalReconModule::GetHitsAssociatedToG4Trajectory(
    const ND::THitSelection& hits, ND::THandle<ND::TG4Trajectory> traj,
    ND::THitSelection& traj_hits) 
{
    //*****************************************************
    
    if (!traj) return;
    
    const ND::TND280Event& event = *(ND::TEventFolder::GetCurrentEvent());
    ND::THandle<ND::TG4TrajectoryContainer> trajectories =
    event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
    
    if (!trajectories) return;
    
    // Loop over reconstructed hits
    for (ND::THitSelection::const_iterator ht = hits.begin(); ht != hits.end();
        ht++) {
    ND::THandle<ND::TComboHit> ch(*ht);
    if (ch) {
        // Loop over hits in Combo
        for (ND::THitSelection::const_iterator ch_member = ch->GetHits().begin();
            ch_member != ch->GetHits().end(); ++ch_member) {
        // Get monte Carlo hit associated.
        ND::THandle<ND::TMCHit> mch;
        ND::THandle<ND::TMultiHit> mhit = *ch_member;
        if (mhit)
            mch = *(mhit->begin());
        else
            mch = *ch_member;
        
        ND::THandle<ND::TG4Trajectory> mainCont =
        GetMainContributor(mch, trajectories);
        if (mainCont)
            if (GetPointer(mainCont) == GetPointer(traj))
            traj_hits.push_back(mch);
            }
    } else {
        ND::THandle<ND::TMultiHit> mhit = *ht;
        ND::THandle<ND::TMCHit> mch;
        
        if (mhit)
            mch = *(mhit->begin());
        else
            mch = *ht;
        
        ND::THandle<ND::TG4Trajectory> mainCont =
        GetMainContributor(mch, trajectories);
        if (mainCont)
            if (GetPointer(mainCont) == GetPointer(traj)) traj_hits.push_back(mch);
    }
        }
}

//*****************************************************
std::string ND::TGlobalReconModule::GetObjectType(
    ND::THandle<ND::TReconBase> object) 
{
    //*****************************************************
    
    std::string reco_class_name = "";
    
    ND::THandle<ND::TReconObjectContainer> conContainer =
    object->GetConstituents();
    if (conContainer) {
        ND::THandle<ND::TReconBase> object2 = *(conContainer->begin());
        if (object2) reco_class_name = object2->ClassName();
    }
    
    return reco_class_name;
}

//*****************************************************
bool ND::TGlobalReconModule::IsTrackLike(ND::THandle<ND::TReconBase> object) 
{
    //*****************************************************
    
    ND::THandle<ND::TReconObjectContainer> conContainer =
    object->GetConstituents();
    
    if (!conContainer) return true;
    
    if (conContainer->size() > 1) return true;
    
    ND::THandle<ND::TReconBase> object2 = *(conContainer->begin());
    if (!object2) return true;
    
    std::string reco_class_name = object2->ClassName();
    
    if (reco_class_name.find("ND::TReconPID") != std::string::npos ||
        reco_class_name.find("ND::TReconTrack") != std::string::npos)
    return true;
    else
        return false;
}

//*****************************************************************************
void ND::TGlobalReconModule::GetConstituentsInTracker(
    ND::THandle<ND::TReconBase> t1, ND::TReconObjectContainer& trackerObjects) 
{
    //*****************************************************************************
    
    // check if this object is tracker only. If so add it
    if (IsTrackerOnly(t1)) {
        trackerObjects.push_back(t1);
        return;
    }
    
    ND::THandle<ND::TReconObjectContainer> t1const = t1->GetConstituents();
    if (!t1const) return;
    
    // loop over constituents
    ND::TReconObjectContainer::const_iterator it1;
    for (it1 = t1const->begin(); it1 != t1const->end(); it1++) {
        ND::THandle<ND::TReconBase> object = *it1;
        if (!object) continue;
        GetConstituentsInTracker(object, trackerObjects);
    }
}

//*****************************************************************************
bool ND::TGlobalReconModule::IsTrackerOnly(ND::THandle<ND::TReconBase> t1)
{
    //*****************************************************************************
    
    if (t1->UsesDetector(ND::TReconBase::kP0D) ||
        t1->UsesDetector(ND::TReconBase::kDSECal) ||
    t1->UsesDetector(ND::TReconBase::kECal) ||
    t1->UsesDetector(ND::TReconBase::kSMRD)) {
    return false;
    }
    
    return true;
}

//*****************************************************************************
void ND::TGlobalReconModule::DoAssociationBetweenTrackerAndGlobalObjects(
    const ND::TReconObjectContainer& globalObjects,
    const ND::TReconObjectContainer& trackerObjects) 
{
    //*****************************************************************************
    
    // 1. loop over all global objects
    for (ND::TReconObjectContainer::const_iterator tt = globalObjects.begin();
        tt != globalObjects.end(); tt++) {
    ND::THandle<ND::TReconBase> gobject = *tt;
    if (!gobject) continue;
    
    // 2. Get the tracker objects in the global object
    ND::TReconObjectContainer trackerObjects2;
    GetConstituentsInTracker(gobject, trackerObjects2);
    
    // 3. loop over tracker objects in the global object
    for (ND::TReconObjectContainer::iterator tt2 = trackerObjects2.begin();
        tt2 != trackerObjects2.end(); tt2++) {
    ND::THandle<ND::TReconBase> tobject2 = *tt2;
    if (!tobject2) continue;
    
    // 4. loop over objects from trackerRecon
    for (ND::TReconObjectContainer::const_iterator tt1 =
        trackerObjects.begin();
        tt1 != trackerObjects.end(); tt1++) {
    ND::THandle<ND::TReconBase> tobject1 = *tt1;
    if (!tobject1) continue;
    
    // 5. Associate the tracker object in the global object with the
    // trackerRecon object
    if (fabs(tobject1->GetQuality() - tobject2->GetQuality()) < 1e-6 &&
        tobject1->GetDetectors() == tobject2->GetDetectors()) {
    // this is the link between a trackerRecon object and a global
    // object
    // index
    fTrackerGlobalIndexMap[tobject1] = fGlobalIndexMap[gobject];
    break;
        }
        }
        }
        }
}



//*****************************************************************************
ND::THandle<ND::TReconBase>
ND::TGlobalReconModule::GetTrackerReconVersionOfFGDIsoTrack(
    ND::THandle<ND::TReconBase> object) 
{
    //*****************************************************************************
    
    ND::THandle<ND::TReconBase> trackerIsoObject;
    if (!object->UsesDetector(ND::TReconBase::kFGD) ||
        object->UsesDetector(ND::TReconBase::kTPC) ||
    object->UsesDetector(ND::TReconBase::kECal) ||
    object->UsesDetector(ND::TReconBase::kDSECal) ||
    object->UsesDetector(ND::TReconBase::kP0D) ||
    object->UsesDetector(ND::TReconBase::kSMRD)) {
    return trackerIsoObject;
    }
    
    bool printStuff = true;
    double objecttime = -1.e10;
    ND::THandle<ND::TReconObjectContainer> cons = object->GetConstituents();
    int ncons = cons->size();
    if (printStuff) {
        std::cout << "Found an FGD-only track" << std::endl;
        if (cons) {
            std::cout << "  and has " << ncons << " constituents" << std::endl;
        } else {
            std::cout << "  and has no constituents" << std::endl;
        }
    }
    ND::THandle<ND::TReconTrack> track = object;
    ND::THandle<ND::TReconPID> pid = object;
    if (track) {
        if (printStuff) std::cout << "   object is a track" << std::endl;
        objecttime = track->GetPosition().T();
    } else if (pid) {
        if (printStuff) std::cout << "   object is a PID" << std::endl;
        objecttime = pid->GetPosition().T();
    } else {
        std::cout << "Not a track or a PID?" << std::endl;
    }
    if (printStuff) {
        std::cout << "object time is " << objecttime << std::endl;
        for (int i = 0; i < ncons; i++) {
            std::cout << "cons " << i << " detectors: P"
            << (*cons)[i]->UsesDetector(ND::TReconBase::kP0D) << " T"
            << (*cons)[i]->UsesDetector(ND::TReconBase::kTPC1)
            << (*cons)[i]->UsesDetector(ND::TReconBase::kFGD1)
            << (*cons)[i]->UsesDetector(ND::TReconBase::kTPC2)
            << (*cons)[i]->UsesDetector(ND::TReconBase::kFGD2)
            << (*cons)[i]->UsesDetector(ND::TReconBase::kTPC3) << " E"
            << (*cons)[i]->UsesDetector(ND::TReconBase::kECal)
            << (*cons)[i]->UsesDetector(ND::TReconBase::kDSECal) << " S"
            << (*cons)[i]->UsesDetector(ND::TReconBase::kSMRD) << std::endl;
        }
    }
    
    // Find corresponding PID in trackerRecon "final" container
    ND::THandle<ND::TAlgorithmResult> trackerRecon =
    ND::TEventFolder::GetCurrentEvent()->GetFit("trackerRecon");
    if (trackerRecon) {
        // Get the results from reconstruction
        ND::THandle<ND::TReconObjectContainer> trackerObjects =
        trackerRecon->GetResultsContainer("final");
        
        if (trackerObjects) {
            for (ND::TReconObjectContainer::iterator tt = trackerObjects->begin();
                tt != trackerObjects->end(); tt++) {
            if ((*tt)->UsesDetector(ND::TReconBase::kFGD) &&
                !(*tt)->UsesDetector(ND::TReconBase::kTPC)) {
            ND::THandle<ND::TReconTrack> track = *tt;
            ND::THandle<ND::TReconPID> pid = *tt;
            double fgdtime = -1.e10;
            if (track) {
                std::cout << "   fgd object is a track?!?!?!?" << std::endl;
                fgdtime = track->GetPosition().T();
            } else if (pid) {
                if (printStuff) std::cout << "   fgd object is a PID" << std::endl;
                fgdtime = pid->GetPosition().T();
            } else {
                std::cout << "fgd not a track or a PID?" << std::endl;
            }
            if (printStuff)
                std::cout << "fgd object time is " << fgdtime << std::endl;
            if (fabs(fgdtime - objecttime) < 1.e-3) {
                if (printStuff) std::cout << "replacing..." << std::endl;
                trackerIsoObject = *tt;
                break;
            }
                } else {
                    if (printStuff) std::cout << "---no TPC component" << std::endl;
                }
                }
        } else {
            std::cout
            << "Why doesn't trackerRecon final results container exist?!?!?!?!?"
            << std::endl;
        }
    } else {
        std::cout << "Why doesn't trackerRecon result exist?!?!?!?!? " << std::endl;
    }
    if (!trackerIsoObject && printStuff) {
        std::cout
        << "!!!!!!!!!!!!!!!! ERROR:  FGD-only object NOT found!!!!!!!!!!!!!!!"
        << std::endl;
    }
    
    return trackerIsoObject;
}

//*****************************************************************************
void ND::TGlobalReconModule::FillBrokenTracksMap(
    const ND::TReconObjectContainer& globalObjects) 
{
    //*****************************************************************************
    
    // This function find the right links between broken tracks
    
    // clear the map from the previous event
    fBrokenIndexMap.clear();
    
    // map of uniqueIDs for the actual broken tracks
    std::map<ND::THandle<ND::TReconBase>, std::vector<UInt_t> > brokenMap1;
    
    // map of uniqueIDs for the tracks associated to a broken track
    std::map<ND::THandle<ND::TReconBase>, std::vector<UInt_t> > brokenMap2;
    
    // loop over global objects
    for (ND::TReconObjectContainer::const_iterator tt = globalObjects.begin();
        tt != globalObjects.end(); tt++) {
    ND::THandle<ND::TReconBase> object = *tt;
    if (!object) continue;
    
    // get the broken IDs recursively
    std::vector<UInt_t> brokenIDs1;
    std::vector<UInt_t> brokenIDs2;
    bool found = GetBrokenIDs(object, brokenIDs1, brokenIDs2);
    
    // if this track is broken fill the maps
    if (found) {
        brokenMap1[object] = brokenIDs1;
        brokenMap2[object] = brokenIDs2;
        
        if (debug_broken) {
            std::cout << " - " << object->GetUniqueID() << " " << object
            << std::endl;
            std::cout << "     --> broken IDs1: ";
            for (unsigned int i = 0; i < brokenIDs1.size(); i++) {
                std::cout << " " << brokenIDs1[i];
            }
            std::cout << std::endl;
            std::cout << "     --> broken IDs2: ";
            for (unsigned int i = 0; i < brokenIDs2.size(); i++) {
                std::cout << " " << brokenIDs2[i];
            }
            std::cout << std::endl;
        }
    }
        }
        
        // now do the actual association
        std::map<ND::THandle<ND::TReconBase>, std::vector<UInt_t> >::const_iterator
        tt1;
        std::map<ND::THandle<ND::TReconBase>, std::vector<UInt_t> >::const_iterator
        tt2;
        for (tt1 = brokenMap1.begin(); tt1 != brokenMap1.end(); tt1++) {
            ND::THandle<ND::TReconBase> object1 = tt1->first;
            if (!object1) continue;
            
            const std::vector<UInt_t>& brokenIDs1 = tt1->second;
            
            std::vector<UInt_t> brokenIDsFinal;
            bool found = false;
            // does any of the brokenIDs of this object exist as a independent global
            // object ?
            for (tt2 = brokenMap2.begin(); tt2 != brokenMap2.end(); tt2++) {
                ND::THandle<ND::TReconBase> object2 = tt2->first;
                if (!object2) continue;
                
                const std::vector<UInt_t>& brokenIDs2 = tt2->second;
                
                for (unsigned int i = 0; i < brokenIDs1.size(); i++) {
                    for (unsigned int j = 0; j < brokenIDs2.size(); j++) {
                        if (brokenIDs1[i] == brokenIDs2[j]) {
                            brokenIDsFinal.push_back(object2->GetUniqueID());
                            found = true;
                            break;
                        }
                    }
                    if (found) break;
                }
                if (found) break;
            }
            
            // fill the map with top level association
            if (found) {
                fBrokenIndexMap[object1] = brokenIDsFinal;
                
                if (debug_broken) {
                    std::cout << " - " << object1->GetUniqueID() << " " << object1
                    << std::endl;
                    std::cout << "     --> broken IDs FINAL: ";
                    for (unsigned int i = 0; i < brokenIDsFinal.size(); i++) {
                        std::cout << " " << brokenIDsFinal[i];
                    }
                    std::cout << std::endl;
                }
            }
        }
}

//*****************************************************************************
bool ND::TGlobalReconModule::GetBrokenIDs(ND::THandle<ND::TReconBase> object,
    std::vector<UInt_t>& brokenIDs1,
    std::vector<UInt_t>& brokenIDs2) 
{
    //*****************************************************************************
    
    bool found = false;
    
    // get the broken IDs for this object
    ND::THandle<ND::TDataVector> broken =
    object->Get<ND::TDataVector>("brokenID");
    if (broken) {
        brokenIDs1.push_back(object->GetUniqueID());
        found = true;
        for (unsigned int i = 0; i < broken->size(); i++) {
            UInt_t ID = broken->At<ND::TIntegerDatum>(i)->GetValue();
            brokenIDs2.push_back(ID);
        }
    }
    
    // Get The BrokenIDs for the original object
    if (ReconObjectUtils::GetOriginalObject(object)) {
        found = GetBrokenIDs(ReconObjectUtils::GetOriginalObject(object),
            brokenIDs1, brokenIDs2);
    } else {
        // get the broken IDs for the constituents (only if the track has no
        // original object)
        ND::THandle<ND::TReconObjectContainer> cons = object->GetConstituents();
        if (!cons) return found;
        
        for (ND::TReconObjectContainer::const_iterator tt = cons->begin();
            tt != cons->end(); tt++) {
        // only for PID constituents
        ND::THandle<ND::TReconPID> object2 = *tt;
        if (!object2) continue;
        
        bool found2 = GetBrokenIDs(object2, brokenIDs1, brokenIDs2);
        if (found2) found = true;
            }
    }
    
    if (found) brokenIDs1.push_back(object->GetUniqueID());
    
    return found;
}

//*************************************************************
void ND::TGlobalReconModule::CheckMatchingFailure(
    TTrackerObject* result, ND::THandle<ND::TReconBase> reco) 

{
    //*************************************************************
    
    bool useFGD1 = reco->UsesDetector(ND::TReconBase::kFGD1);
    bool useFGD2 = reco->UsesDetector(ND::TReconBase::kFGD2);
    
    ND::THandle<ND::TReconBase> fgdCons;
    ND::THandle<ND::TReconPID> fgdPid;
    ND::THandle<ND::TReconTrack> fgdTrack;
    
    if (useFGD1) {
        fgdCons =
        ReconObjectUtils::GetConstituentInDetector(reco, ND::TReconBase::kFGD1);
        
        if (fgdCons) {
            fgdPid = fgdCons;
            if (fgdPid) {
                ND::TReconObjectContainer::const_iterator it =
                fgdPid->GetConstituents()->begin();
                fgdTrack = *it;
            }
        }
        
        if (fgdTrack && fgdTrack->Has<ND::TIntegerDatum>("matchingFailure_flag")) {
            result->matchingFailure_flag =
            fgdTrack->Get<ND::TIntegerDatum>("matchingFailure_flag")->GetValue();
            
        } else
            result->matchingFailure_flag = 0;
    }
    
    else if (useFGD2) {
        fgdCons =
        ReconObjectUtils::GetConstituentInDetector(reco, ND::TReconBase::kFGD2);
        
        if (fgdCons) {
            fgdPid = fgdCons;
            if (fgdPid) {
                ND::TReconObjectContainer::const_iterator it =
                fgdPid->GetConstituents()->begin();
                fgdTrack = *it;
            }
        }
        
        if (fgdTrack && fgdTrack->Has<ND::TIntegerDatum>("matchingFailure_flag"))
            result->matchingFailure_flag =
        fgdTrack->Get<ND::TIntegerDatum>("matchingFailure_flag")->GetValue();
        else
            result->matchingFailure_flag = 0;
    }
    
    else
        result->matchingFailure_flag = 0;
    
    return;
}