TReconPerformanceEvalModule.cxx

Go to the documentation of this file.

#include "TReconPerformanceEvalModule.hxx"

#include <TParticlePDG.h>
#include <TString.h>
#include <HelixModel.hxx>
#include <ND280GeomId.hxx>
#include <ReconObjectUtils.hxx>
#include <ReconPrintout.hxx>
#include <TG4Trajectory.hxx>
#include <TG4TrajectoryPoint.hxx>
#include <TGeomIdManager.hxx>
#include <TEventReconUtils.hxx>
#include <TrackTruthInfo.hxx>
#include <TrackingUtils.hxx>
#include <cmath>
#include "TDatum.hxx"
#include "TGeomInfo.hxx"
#include "TOARuntimeParams.hxx"
#include "TPrincipal.h"
#include "TRealDatum.hxx"
#include "TReconNode.hxx"

ClassImp(ND::TReconPerformanceEvalModule::TGlobalReconObject);
ClassImp(ND::TReconPerformanceEvalModule::TGlobalReconConstituent);
ClassImp(ND::TReconPerformanceEvalModule::TGlobalReconStateInfo);
ClassImp(ND::TReconPerformanceEvalModule::TGlobalReconNodeInfo);
ClassImp(ND::TReconPerformanceEvalModule::TReconStateInfo);
ClassImp(ND::TReconPerformanceEvalModule::TDownToTrackerInfo);
ClassImp(ND::TReconPerformanceEvalModule::THelixModelInfo);
ClassImp(ND::TReconPerformanceEvalModule::TMatchingChi2Info);
ClassImp(ND::TReconPerformanceEvalModule::TGlobalTruthInfo);
ClassImp(ND::TReconPerformanceEvalModule::TTruthInfo);

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

ND::TReconPerformanceEvalModule::TReconPerformanceEvalModule(
    const char* name, const char* title) {
  SetNameTitle(name, title);
  fIsEnabled = kTRUE;
  fDescription = "Global Recon Inputs and Outputs";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;

  // Tree variables
  fGlobalReconObjects = new TClonesArray(
      "ND::TReconPerformanceEvalModule::TGlobalReconObject", 200);

  // Whether to save info on all global nodes or just the first and last ones
  fSaveAllGlobalNodes = true;
}

ND::TReconPerformanceEvalModule::~TReconPerformanceEvalModule() {}

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

bool ND::TReconPerformanceEvalModule::IsFullEventWorthSaving(
    ND::TND280Event& event) {
  return (fNGlobalReconObjects ? true : false);
}

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

void ND::TReconPerformanceEvalModule::InitializeBranches() {
  // Number of global recon objects
  fOutputTree->Branch("NGlobalReconObject", &fNGlobalReconObjects,
                      "NGlobalReconObject/I", fBufferSize);

  // TClonesArray of TGlobalReconObject objects
  fOutputTree->Branch("GlobalReconObject", &fGlobalReconObjects, fBufferSize,
                      fSplitLevel);

  // Number of hits in each subdetector
  fOutputTree->Branch("HitInfo", &fHitInfo, fBufferSize, fSplitLevel);
}

void ND::TReconPerformanceEvalModule::FillConfigTree(TTree* configTree) {
  if (fFilledConfigTree) {
    ND280NamedWarn(
        "TReconPerformanceEvalModule",
        "Module: " << this->GetName()
                   << " has already written the config tree this run.");
  }

  configTree->Branch("SaveAllGlobalNodes", &fSaveAllGlobalNodes,
                     "SaveAllGlobalNodes/O");

  TAnalysisModuleBase::FillConfigTree(configTree);
}

bool ND::TReconPerformanceEvalModule::FillTree(ND::TND280Event& event) {
  // Initialise
  fNGlobalReconObjects = 0;
  fGlobalReconObjects->Clear();
  fHitInfo.clear();

  ND::THandle<ND::TAlgorithmResult> globalResult = event.GetFit("oaRecon");

  if (globalResult) {
    ND280Verbose("Got global result for event ID " << event.GetEventId());
    for (ND::TDataVector::iterator dataIter = globalResult->begin();
         dataIter != globalResult->end(); ++dataIter) {
      // Looking for ObjectContainers that are labelled correctly
      std::string class_name = (*dataIter)->ClassName();
      TString object_name = (*dataIter)->GetName();

      if (class_name.find("ND::TReconObjectContainer") != std::string::npos &&
          object_name.EqualTo("final")) {
        ND::TReconObjectContainer* reconContainer =
            dynamic_cast<ND::TReconObjectContainer*>(*dataIter);

        if (!reconContainer) {
          ND280Warn("Dynamic cast of ReconObjectContainer didn't work");
        }

        // Loop looking for objects in the top recon container
        for (ND::TReconObjectContainer::iterator reconIter =
                 reconContainer->begin();
             reconIter != reconContainer->end(); ++reconIter) {
          // Get the global object and fill information related to it
          ND::THandle<ND::TReconBase> globalObject = (*reconIter);

          TGlobalReconObject* globalReconObject =
              new ((*fGlobalReconObjects)[fNGlobalReconObjects])
                  TGlobalReconObject;

          FillGlobalInfo(event, globalObject, globalReconObject);
          fNGlobalReconObjects++;

          ND::TReconObjectContainer conContainer =
              GetConstituentsOfInterest(globalObject);
          ND280Verbose("  Next PID - got " << conContainer.size()
                                           << " constituents, "
                                           << globalObject->ConvertDetector());

          // Loop over the constituents and fill information related to each
          for (ND::TReconObjectContainer::iterator conIter =
                   conContainer.begin();
               conIter != conContainer.end(); ++conIter) {
            ND::THandle<ND::TReconBase> base = (*conIter);

            ND280Verbose("    " << base->ClassName() << " from "
                                << base->ConvertDetector());
            TGlobalReconConstituent* constituent =
                new ((*(globalReconObject
                            ->Constituents))[globalReconObject->NConstituents])
                    TGlobalReconConstituent;
            FillConstituentInfo(event, base, globalObject, constituent);

            // Fill counters of how many objects we have in total, and in each
            // module
            globalReconObject->NConstituents++;

            std::string det = base->ConvertDetector();
            (*(globalReconObject->NModuleConstituents))[det]++;
          }

          // For PIDs that are made of tracker + "something", see what effect
          // adding the "something" has
          if (ContainsTracker(globalObject) && !IsOnlyTracker(globalObject)) {
            ND::TReconObjectContainer trackerConContainer =
                GetDownToTracker(globalObject);
            ND280Verbose("  Next PID - got "
                         << trackerConContainer.size()
                         << " constituents on way to tracker, "
                         << globalObject->ConvertDetector());

            for (ND::TReconObjectContainer::iterator conIter =
                     trackerConContainer.begin();
                 conIter != trackerConContainer.end(); ++conIter) {
              ND::THandle<ND::TReconBase> base = (*conIter);

              ND280Verbose("    " << base->ClassName() << " from "
                                  << base->ConvertDetector());
              TDownToTrackerInfo* constituent =
                  new ((*(globalReconObject->DownToTrackerConstituents))
                           [globalReconObject->NDownToTrackerConstituents])
                      TDownToTrackerInfo;
              FillDownToTrackerInfo(base, constituent);
              globalReconObject->NDownToTrackerConstituents++;
            }
          }

          // Fill information on the matching chi^2s that went in to this object
          ND::TReconObjectContainer mConContainer =
              ReconObjectUtils::GetConstituentsRecursive(globalObject, true);
          ND280Verbose("  Got " << mConContainer.size()
                                << " constituents for matching chi^2");

          for (ND::TReconObjectContainer::iterator conIter =
                   mConContainer.begin();
               conIter != mConContainer.end(); ++conIter) {
            ND::THandle<ND::TReconBase> base = (*conIter);

            if (base->Has<ND::TRealDatum>("matching_chi2ndf")) {
              TMatchingChi2Info* matchInfo = new (
                  (*(globalReconObject->MatchingChi2Info))
                      [globalReconObject->NMatchingChi2Info]) TMatchingChi2Info;
              double matchVal = TrackingUtils::GetMatchingChi2(*base);
              ND::THandle<ND::TReconObjectContainer> cons =
                  base->GetConstituents();

              if (cons && cons->size() == 2) {
                ND::THandle<ND::TReconBase> con1 = cons->at(0);
                ND::THandle<ND::TReconBase> con2 = cons->at(1);
                matchInfo->Det1 = con1->ConvertDetector();
                matchInfo->Det2 = con2->ConvertDetector();
                matchInfo->IsShower1 = TrackingUtils::IsShower(con1);
                matchInfo->IsShower2 = TrackingUtils::IsShower(con2);
                matchInfo->MatchingChi2 = matchVal;
                matchInfo->SetOK = true;
              }

              globalReconObject->NMatchingChi2Info++;
            }
          }
        }
      }
    }
  }

  // Fill hit information
  ND::THandle<ND::TDataVector> hitsData = event.Get<ND::TDataVector>("hits");
  for (ND::TDataVector::iterator dataIter = hitsData->begin();
       dataIter != hitsData->end(); ++dataIter) {
    TString object_name = (*dataIter)->GetName();
    ND::THitSelection* hits = dynamic_cast<ND::THitSelection*>(*dataIter);
    fHitInfo[object_name.Data()] = hits->size();
  }

  return true;
}

/// Get the constituents of a global PID that we wish to investigate further.
/// Gets the list of unique highest-level constituents in each subdetector.
ND::TReconObjectContainer
ND::TReconPerformanceEvalModule::GetConstituentsOfInterest(
    ND::THandle<ND::TReconBase> input) {
  ND::TReconObjectContainer allCons;

  static const int NDETS = 19;
  ND::TReconBase::Status dets[NDETS] = {
      ND::TReconBase::kTPC1,        ND::TReconBase::kFGD1,
      ND::TReconBase::kTPC2,        ND::TReconBase::kFGD2,
      ND::TReconBase::kTPC3,        ND::TReconBase::kDSECal,
      ND::TReconBase::kRightTECal,  ND::TReconBase::kLeftTECal,
      ND::TReconBase::kBottomTECal, ND::TReconBase::kTopTECal,
      ND::TReconBase::kRightPECal,  ND::TReconBase::kLeftPECal,
      ND::TReconBase::kBottomPECal, ND::TReconBase::kTopPECal,
      ND::TReconBase::kP0D,         ND::TReconBase::kRightSMRD,
      ND::TReconBase::kLeftSMRD,    ND::TReconBase::kBottomSMRD,
      ND::TReconBase::kTopSMRD};

  for (int ii = 0; ii < NDETS; ii++) {
    ND::TReconObjectContainer cons =
        ReconObjectUtils::GetAllConstituentsInDetector(input, dets[ii]);
    for (ND::TReconObjectContainer::const_iterator it1 = cons.begin();
         it1 != cons.end(); it1++) {
      // Ensure the constituent is unique
      bool unique = true;
      for (ND::TReconObjectContainer::const_iterator it2 = allCons.begin();
           it2 != allCons.end(); it2++) {
        if ((*it1) == (*it2)) {
          unique = false;
          break;
        }
      }
      if (unique) {
        allCons.push_back(*it1);
      }
    }
  }

  return allCons;
}

/// Get the list of constituents that are found on the way down to getting the
/// tracker-only constituent
ND::TReconObjectContainer ND::TReconPerformanceEvalModule::GetDownToTracker(
    ND::THandle<ND::TReconBase> input) {
  ND::TReconObjectContainer allCons;

  if (IsOnlyTracker(input)) {
    allCons.push_back(input);
  } else if (ContainsTracker(input)) {
    allCons.push_back(input);
    ND::THandle<ND::TReconObjectContainer> cons = input->GetConstituents();

    if (cons) {
      for (ND::TReconObjectContainer::const_iterator it = cons->begin();
           it != cons->end(); it++) {
        ND::TReconObjectContainer moreCons = GetDownToTracker(*it);

        for (ND::TReconObjectContainer::const_iterator it2 = moreCons.begin();
             it2 != moreCons.end(); it2++) {
          allCons.push_back(*it2);
        }
      }
    }
  }

  return allCons;
}

bool ND::TReconPerformanceEvalModule::IsOnlyTracker(
    ND::THandle<ND::TReconBase> input) {
  ND::TReconBase::Status notTrackerDets[14] = {
      ND::TReconBase::kDSECal,     ND::TReconBase::kRightTECal,
      ND::TReconBase::kLeftTECal,  ND::TReconBase::kBottomTECal,
      ND::TReconBase::kTopTECal,   ND::TReconBase::kRightPECal,
      ND::TReconBase::kLeftPECal,  ND::TReconBase::kBottomPECal,
      ND::TReconBase::kTopPECal,   ND::TReconBase::kP0D,
      ND::TReconBase::kRightSMRD,  ND::TReconBase::kLeftSMRD,
      ND::TReconBase::kBottomSMRD, ND::TReconBase::kTopSMRD};

  // Ensure that the object doesn't use non-tracker modules
  for (int ii = 0; ii < 14; ii++) {
    if (input->UsesDetector(notTrackerDets[ii])) {
      return false;
    }
  }

  // Ensure that it does use tracker modules
  return ContainsTracker(input);
}

bool ND::TReconPerformanceEvalModule::ContainsTracker(
    ND::THandle<ND::TReconBase> input) {
  return (input->UsesDetector(ND::TReconBase::kTPC) ||
          input->UsesDetector(ND::TReconBase::kFGD));
}

/// Fill information in the output tree relating to the global PID
void ND::TReconPerformanceEvalModule::FillGlobalInfo(
    ND::TND280Event& event, ND::THandle<ND::TReconBase> reconObject,
    TGlobalReconObject* output) {
  ND::THandle<ND::TReconState> globalState;

  try {
    globalState = reconObject->GetState();
  } catch (ND::EReconObject e) {
    ND280Warn("TReconBase with no state. Skip object!");
    return;
  }

  output->UniqueID = reconObject->GetUniqueID();
  output->Momentum = GetStateMomentum(globalState);
  //  output->MomentumByRange = TrackingUtils::GetMomentumByRange(*reconObject);
  const double MomentumByRangeDefault = -10000;
  output->MomentumByRangeMuon = TrackingUtils::GetTRealDatumValue(
      *reconObject, "prange_muon", MomentumByRangeDefault);
  output->MomentumByRangeMuonFlip = TrackingUtils::GetTRealDatumValue(
      *reconObject, "prange_muon_flipped", MomentumByRangeDefault);
  output->MomentumByRangeElectron = TrackingUtils::GetTRealDatumValue(
      *reconObject, "prange_electron", MomentumByRangeDefault);
  output->MomentumByRangeElectronFlip = TrackingUtils::GetTRealDatumValue(
      *reconObject, "prange_electron_flipped", MomentumByRangeDefault);
  output->MomentumByRangeProton = TrackingUtils::GetTRealDatumValue(
      *reconObject, "prange_proton", MomentumByRangeDefault);
  output->MomentumByRangeProtonFlip = TrackingUtils::GetTRealDatumValue(
      *reconObject, "prange_proton_flipped", MomentumByRangeDefault);

  output->Position = GetStatePosition(globalState);
  output->Direction = GetStateDirection(globalState);
  output->Quality = reconObject->GetQuality();
  output->NDOF = reconObject->GetNDOF();
  output->SubdetectorString = reconObject->ConvertDetector();
  output->StatusString = reconObject->ConvertStatus();

  ND::THandle<ND::TReconPID> pid = reconObject;
  if (pid) {
    output->ParticleID = pid->ConvertParticleId();
    output->PIDWeight = pid->GetPIDWeight();
    output->Charge = pid->GetCharge();
  }

  ND::TReconNodeContainer nodes = reconObject->GetNodes();
  for (ND::TReconNodeContainer::iterator it = nodes.begin(); it != nodes.end();
       it++) {
    // Only save information about the first and last nodes, unless the flag to
    // save them all has
    // been switched on.
    if (fSaveAllGlobalNodes || output->NGlobalNodes == 0 ||
        output->NGlobalNodes == ((int)nodes.size() - 1)) {
      TGlobalReconNodeInfo* nodeInfo =
          new ((*(output->GlobalNodes))[output->NGlobalNodesSaved])
              TGlobalReconNodeInfo;
      ND::THandle<ND::TReconState> state = (*it)->GetState();
      FillGlobalStateInfo(state, nodeInfo->NodeState);

      ND::THandle<ND::TReconBase> object = (*it)->GetObject();
      if (!object) {
        ND280Warn("Unable to get object of global node");
      } else {
        ND::THandle<ND::TReconState> nodeState;
        try {
          nodeState = object->GetState();
        } catch (ND::EReconObject e) {
          ND280Warn("TReconBase with no state. Skip node!");
          continue;
        }

        FillGlobalStateInfo(nodeState, nodeInfo->ObjectState);
      }

      output->NGlobalNodesSaved++;
    }
    output->NGlobalNodes++;
  }

  FillGlobalTruthInfo(event, reconObject, output->Truth);

  output->SetOK = true;

  return;
}

/// Accessor for getting position from TReconState. Supported by all TxxxState
/// that inherit from TReconState.
TLorentzVector ND::TReconPerformanceEvalModule::GetStatePosition(
    ND::THandle<ND::TReconState> state) {
  TLorentzVector pos(-10000, -10000, -10000, -10000);

  ND::THandle<ND::TPIDState> pid = state;
  if (pid) {
    pos = pid->GetPosition();
  }
  ND::THandle<ND::TTrackState> tr = state;
  if (tr) {
    pos = tr->GetPosition();
  }
  ND::THandle<ND::TShowerState> sh = state;
  if (sh) {
    pos = sh->GetPosition();
  }
  ND::THandle<ND::TClusterState> cl = state;
  if (cl) {
    pos = cl->GetPosition();
  }
  ND::THandle<ND::TVertexState> ve = state;
  if (ve) {
    pos = ve->GetPosition();
  }

  return pos;
}

/// Accessor for getting position variance from TReconState. Supported by all
/// TxxxState that inherit from TReconState.
TLorentzVector ND::TReconPerformanceEvalModule::GetStatePositionVariance(
    ND::THandle<ND::TReconState> state) {
  TLorentzVector pos(-10000, -10000, -10000, -10000);

  ND::THandle<ND::TPIDState> pid = state;
  if (pid) {
    pos = pid->GetPositionVariance();
  }
  ND::THandle<ND::TTrackState> tr = state;
  if (tr) {
    pos = tr->GetPositionVariance();
  }
  ND::THandle<ND::TShowerState> sh = state;
  if (sh) {
    pos = sh->GetPositionVariance();
  }
  ND::THandle<ND::TClusterState> cl = state;
  if (cl) {
    pos = cl->GetPositionVariance();
  }
  ND::THandle<ND::TVertexState> ve = state;
  if (ve) {
    pos = ve->GetPositionVariance();
  }

  return pos;
}

/// Accessor for getting position from TReconState. Only supported by TPIDState.
double ND::TReconPerformanceEvalModule::GetStateMomentum(
    ND::THandle<ND::TReconState> state) {
  double mom(-10000);

  ND::THandle<ND::TPIDState> pid = state;
  if (pid) {
    mom = pid->GetMomentum();
  }

  return mom;
}

/// Accessor for getting direction from TReconState. Only supported by
/// TPIDState, TTrackState and TShowerState.
TVector3 ND::TReconPerformanceEvalModule::GetStateDirection(
    ND::THandle<ND::TReconState> state) {
  TVector3 dir(-10000, -10000, -10000);

  ND::THandle<ND::TPIDState> pid = state;
  if (pid) {
    dir = pid->GetDirection();
  }
  ND::THandle<ND::TTrackState> tr = state;
  if (tr) {
    dir = tr->GetDirection();
  }
  ND::THandle<ND::TShowerState> sh = state;
  if (sh) {
    dir = sh->GetDirection();
  }

  return dir;
}

/// Get the EDeposit of a TRecon(Track|Shower|Cluster) or a TReconPID that has a
/// single constituent
/// that is a TRecon(Track|Shower|Cluster)
double ND::TReconPerformanceEvalModule::GetEDeposit(
    ND::THandle<ND::TReconBase> object) {
  double dep = -10000;

  ND::THandle<ND::TReconTrack> tr = object;
  ND::THandle<ND::TReconShower> sh = object;
  ND::THandle<ND::TReconCluster> cl = object;
  ND::THandle<ND::TReconPID> pid = object;

  if (pid) {
    ND::THandle<ND::TReconObjectContainer> cons = object->GetConstituents();
    if (cons && cons->size() == 1) {
      ND::THandle<ND::TReconBase> con = cons->front();
      tr = con;
      sh = con;
      cl = con;
    }
  }

  if (tr) {
    dep = tr->GetEDeposit();
  }
  if (sh) {
    dep = sh->GetEDeposit();
  }
  if (cl) {
    dep = cl->GetEDeposit();
  }

  return dep;
}

/// Get the state of the last node in an object
ND::THandle<ND::TReconState> ND::TReconPerformanceEvalModule::GetLastState(
    ND::THandle<ND::TReconBase> object) {
  ND::THandle<ND::TReconState> state;

  for (ND::TReconNodeContainer::const_reverse_iterator it =
           object->GetNodes().rbegin();
       it != object->GetNodes().rend(); it++) {
    state = (*it)->GetState();

    if (state) {
      return state;
    }
  }

  return state;
}

/// Get the state of the first node in an object
ND::THandle<ND::TReconState> ND::TReconPerformanceEvalModule::GetFirstState(
    ND::THandle<ND::TReconBase> object) {
  ND::THandle<ND::TReconState> state;

  for (ND::TReconNodeContainer::const_iterator it = object->GetNodes().begin();
       it != object->GetNodes().end(); it++) {
    state = (*it)->GetState();

    if (state) {
      return state;
    }
  }

  return state;
}

/// Fill information in the output tree relating to a constituent of a global
/// PID.
void ND::TReconPerformanceEvalModule::FillConstituentInfo(
    ND::TND280Event& event, ND::THandle<ND::TReconBase> reconObject,
    ND::THandle<ND::TReconBase> globalObject, TGlobalReconConstituent* output) {
  output->ReconAlgorithm = reconObject->GetAlgorithmName();
  output->DetectorName = reconObject->ConvertDetector();
  output->ClassName = reconObject->ClassName();
  output->StatusString = reconObject->ConvertStatus();
  output->Quality = reconObject->GetQuality();
  output->NDOF = reconObject->GetNDOF();
  output->IsShower = TrackingUtils::IsShower(reconObject);
  output->EDeposit = GetEDeposit(reconObject);

  if (ReconObjectUtils::GetOriginalObject(reconObject))
    output->UniqueID =
        ReconObjectUtils::GetOriginalObject(reconObject)->GetUniqueID();
  else
    output->UniqueID = reconObject->GetUniqueID();

  ND::THandle<ND::TReconPID> pid = reconObject;
  if (pid) {
    output->ParticleID = pid->ConvertParticleId();
    output->PIDWeight = pid->GetPIDWeight();
  }

  if (reconObject->GetHits()) {
    output->NumHits = reconObject->GetHits()->size();
  } else {
    output->NumHits = -1;
  }

  ND::THandle<ND::TRealDatum> averageTime =
      reconObject->Get<ND::TRealDatum>("averageHitTime");
  if (averageTime) {
    output->AverageHitTime = averageTime->GetValue();
  }

  TVector3 refittedStatePos(-10000, -10000, -10000);

  ND::THandle<ND::TReconState> refittedState;
  try {
    refittedState = reconObject->GetState();
  } catch (ND::EReconObject e) {
    ND280Warn("TReconBase with no state.");
  }

  if (!refittedState) {
    ND280Warn("      Unable to get state of refitted object");
  } else {
    ND280Verbose("      Got refitted object state");

    // Overall state of constituent
    refittedStatePos = GetStatePosition(refittedState).Vect();
    FillStateInfo(refittedState, output->GlobalReconConsState);

    ND::THandle<ND::TReconNode> refittedFirstNode =
        TrackingUtils::GetFirstFittedNode(*reconObject);
    if (!refittedFirstNode) {
      ND280Warn("      Unable to get first node of refitted object");
    } else {
      // State of first node of the constituent
      ND::THandle<ND::TReconState> refittedFirstNodeState =
          refittedFirstNode->GetState();
      if (!refittedFirstNodeState) {
        ND280Warn("      Unable to get first node state of refitted object");
      } else {
        FillStateInfo(refittedFirstNodeState,
                      output->GlobalReconConsFirstNodeState);
      }

      // Object of first node of the constituent
      ND::THandle<ND::TReconBase> refittedFirstNodeBase =
          refittedFirstNode->GetObject();
      if (!refittedFirstNodeBase) {
        ND280Warn("      Unable to get first node object of refitted object");
      } else {
        try {
          FillStateInfo(refittedFirstNodeBase->GetState(),
                        output->GlobalReconConsFirstNodeObject);
        } catch (ND::EReconObject e) {
          ND280Warn("TReconBase first node with no state.");
        }
      }
    }

    ND::THandle<ND::TReconNode> refittedLastNode =
        TrackingUtils::GetLastFittedNode(*reconObject);
    if (!refittedLastNode) {
      ND280Warn("      Unable to get last node of refitted object");
    } else {
      // State of last node of the constituent
      ND::THandle<ND::TReconState> refittedLastNodeState =
          refittedLastNode->GetState();
      if (!refittedLastNodeState) {
        ND280Warn("      Unable to get last node state of refitted object");
      } else {
        FillStateInfo(refittedLastNodeState,
                      output->GlobalReconConsLastNodeState);
      }

      // Object of last node of the constituent
      ND::THandle<ND::TReconBase> refittedLastNodeBase =
          refittedLastNode->GetObject();
      if (!refittedLastNodeBase) {
        ND280Warn("      Unable to get last node object of refitted object");
      } else {
        try {
          FillStateInfo(refittedLastNodeBase->GetState(),
                        output->GlobalReconConsLastNodeObject);
        } catch (ND::EReconObject e) {
          ND280Warn("TReconBase last node with no state.");
        }
      }
    }

    // State of closest global node to constituent, extrapolated to the
    // constituent's state position
    ND::THandle<ND::TReconState> closestGlobalNodeState =
        TrackingUtils::GetFirstState(*globalObject);
    ND::TReconNodeContainer globalNodes = globalObject->GetNodes();
    double closestOffset = 1e10;

    ND280Verbose("      Got " << globalNodes.size() << " global nodes");

    for (ND::TReconNodeContainer::iterator it = globalNodes.begin();
         it != globalNodes.end(); it++) {
      ND::THandle<ND::TReconNode> node(*it);
      if (!node) {
        ND280Warn("        Unable to get node from global container");
        continue;
      }

      ND::THandle<ND::TReconState> thisState(node->GetState());
      if (!thisState) {
        ND280Warn("        Unable to get state of global node");
        continue;
      }

      double thisOffset =
          (refittedStatePos - GetStatePosition(thisState).Vect()).Mag();
      if (thisOffset < closestOffset) {
        closestGlobalNodeState = thisState;
        closestOffset = thisOffset;
      }
    }

    if (!closestGlobalNodeState) {
      ND280Warn("      Didn't find closest global node state");
    } else {
      ND280Verbose("      Got closest global node state");
      // Extrapolate the closest global state towards the constituent's state.
      // Extrapolation is done to the surface with normal (0,0,1) in the LOCAL
      // geometry
      // of the module we're in. So, for the Barrel ECals, that means we
      // extrapolate to
      // the surface with normal (1,0,0) or (0,1,0).
      ND::THandle<ND::TReconState> projectedClosestGlobalNodeState(
          new ND::TPIDState);
      TVector3 snorm(0, 0, 1);
      std::string det = reconObject->ConvertDetector();

      if (det.find("LTECal") != std::string::npos ||
          det.find("RTECal") != std::string::npos) {
        snorm = TVector3(1, 0, 0);
      } else if (det.find("TTECal") != std::string::npos ||
                 det.find("BTECal") != std::string::npos) {
        snorm = TVector3(0, 1, 0);
      } else if (det.find("LPECal") != std::string::npos ||
                 det.find("RPECal") != std::string::npos) {
        snorm = TVector3(1, 0, 0);
      } else if (det.find("TPECal") != std::string::npos ||
                 det.find("BPECal") != std::string::npos) {
        snorm = TVector3(0, 1, 0);
      } else if (det.find("LSMRD") != std::string::npos ||
                 det.find("RSMRD") != std::string::npos) {
        snorm = TVector3(1, 0, 0);
      } else if (det.find("TSMRD") != std::string::npos ||
                 det.find("BSMRD") != std::string::npos) {
        snorm = TVector3(0, 1, 0);
      }

      if (fabs(refittedStatePos.X()) < 5000 &&
          fabs(refittedStatePos.Y()) < 5000 &&
          fabs(refittedStatePos.Z()) < 5000) {
        if (ND::tman().PropagateToSurface(*closestGlobalNodeState,
                                          refittedStatePos, snorm,
                                          *projectedClosestGlobalNodeState)) {
          FillStateInfo(projectedClosestGlobalNodeState,
                        output->ClosestGlobalNodeState);
        } else {
          ND280Warn("      Unable to propagate closest global node state");
        }
      } else {
        ND280Warn("      Attempted to propagate to bad position "
                  << refittedStatePos.X() << ", " << refittedStatePos.Y()
                  << ", " << refittedStatePos.Z());
      }
    }

    //
    // Information of first state of original subdetector recon object (before
    // refitting by global recon)
    //

    ND::THandle<ND::TReconBase> original =
        ReconObjectUtils::GetOriginalObject(reconObject);

    if (original) {
      try {
        FillStateInfo(original->GetState(), output->OriginalObjectState);
      } catch (ND::EReconObject e) {
        ND280Warn("      Unable to get state of unfitted object");
      }

      ND::THandle<ND::TReconState> unfittedFirstState = GetFirstState(original);
      if (!unfittedFirstState) {
        ND280Warn("      Unable to get first state of unfitted object");
      } else {
        FillStateInfo(unfittedFirstState, output->OriginalObjectFirstNodeState);
      }

      ND::THandle<ND::TReconState> unfittedLastState = GetLastState(original);
      if (!unfittedLastState) {
        ND280Warn("      Unable to get last state of unfitted object");
      } else {
        FillStateInfo(unfittedLastState, output->OriginalObjectLastNodeState);
      }
    }

    // If the constituent is in a TPC, extrapolate its state towards the
    // ECals/P0D, for matching studies
    if (reconObject->GetDetectors() == ND::TReconBase::kTPC1 ||
        reconObject->GetDetectors() == ND::TReconBase::kTPC2 ||
        reconObject->GetDetectors() == ND::TReconBase::kTPC3) {
      ExtrapolateToModuleAndFill(reconObject, output->StateAtDSECal, "DS");
      ExtrapolateToModuleAndFill(reconObject, output->StateAtTTECal, "TT");
      ExtrapolateToModuleAndFill(reconObject, output->StateAtBTECal, "BT");
      ExtrapolateToModuleAndFill(reconObject, output->StateAtLTECal, "LT");
      ExtrapolateToModuleAndFill(reconObject, output->StateAtRTECal, "RT");
      ExtrapolateToModuleAndFill(reconObject, output->StateAtP0D, "P0D");
      if (reconObject->GetDetectors() == ND::TReconBase::kTPC1 ||
          reconObject->GetDetectors() == ND::TReconBase::kTPC3) {
        ExtrapolateToModuleAndFill(reconObject, output->StateAtTPC2, "TPC2");
      }
    }

    output->SetOK = true;
  }

  //
  // Subdetector truth information
  //
  FillConsTruthInfo(event, reconObject, output->ConsTruth);

  return;
}

/// Extract TReconPIDs from P0DRecon's TReconVertex container.
void ND::TReconPerformanceEvalModule::GetObjectsFromVertices(
    ND::THandle<ND::TReconObjectContainer> input,
    ND::THandle<ND::TReconObjectContainer> output) {
  if (!input || !output) return;

  ND280Verbose("GetObjectsFromVertices: # objects = " << input->size());

  ND::TReconObjectContainer::const_iterator pp, pp1;
  for (pp = input->begin(); pp != input->end(); pp++) {
    ND280Verbose("Next input object: " << *pp);
    ND::THandle<ND::TReconVertex> vertex = *pp;
    if (vertex) {
      ND::THandle<ND::TReconObjectContainer> cons = vertex->GetConstituents();
      if (cons) {
        ND280Verbose("   ---> # cons = " << cons->size());
        for (pp1 = cons->begin(); pp1 != cons->end(); pp1++) {
          output->push_back(*pp1);
        }
      }
    } else {
      output->push_back(*pp);
    }
  }
}

/// Check whether the THandles of the immediate constituents of two objects are
/// equal.
bool ND::TReconPerformanceEvalModule::ConstituentsAreEqual(
    ND::THandle<ND::TReconBase> obj1, ND::THandle<ND::TReconBase> obj2) {
  ND::THandle<ND::TReconObjectContainer> hcon1 = obj1->GetConstituents();
  ND::THandle<ND::TReconObjectContainer> hcon2 = obj2->GetConstituents();

  if (ND::TND280Log::GetLogLevel("eventAnalysis") >= ND::TND280Log::VerboseLevel) {
    std::cout << "obj1 constituents:" << std::endl;
    ReconPrintout::DumpContainerWithHistory(*hcon1);
    std::cout << "obj2 constituents:" << std::endl;
    ReconPrintout::DumpContainerWithHistory(*hcon2);
  }

  if (!hcon1 || !hcon2) {
    ND280Verbose("Unable to get constituents");
    return false;
  }

  ND::TReconObjectContainer con1 = (*hcon1);
  ND::TReconObjectContainer con2 = (*hcon2);

  if (con1.size() != con2.size()) {
    ND280Verbose("Constituent ROCs are different sizes: "
                 << con1.size() << " vs " << con2.size());
    return false;
  }

  for (ND::TReconObjectContainer::iterator it1 = con1.begin();
       it1 != con1.end(); it1++) {
    bool foundMatch = false;

    for (ND::TReconObjectContainer::iterator it2 = con2.begin();
         it2 != con2.end(); it2++) {
      if ((*it1) == (*it2)) {
        foundMatch = true;
        break;
      }
    }

    if (!foundMatch) {
      return false;
    }
  }

  return true;
}

/// Fill information relating to the true trajectory associated to a constituent
/// of a global PID. Fills information
/// from the closest trajectory point to the constituent.
void ND::TReconPerformanceEvalModule::FillConsTruthInfo(
    ND::TND280Event& event, ND::THandle<ND::TReconBase> reconObject,
    TTruthInfo& output) {
  ND::THandle<ND::TReconState> state;
  try {
    state = reconObject->GetState();
  } catch (ND::EReconObject e) {
    return;
  }

  if (!state) {
    return;
  }

  TVector3 statePos = GetStatePosition(state).Vect();
  ND::THandle<ND::TG4TrajectoryContainer> trajectories =
      event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");

  if (trajectories && reconObject->GetHits()) {
    // Use oaUtility to find true particle corresponding to this track, and get
    // the associated trajectory points.
    int G4TrkId = TrackTruthInfo::GetG4TrajIDHits(
        *(reconObject->GetHits()), output.Efficiency, output.Purity);
    ND::THandle<ND::TG4Trajectory> traj = trajectories->GetTrajectory(G4TrkId);
    if (traj) {
      ND280Verbose("      Got trajectory number " << G4TrkId);

      std::vector<ND::TG4TrajectoryPoint> trajPoints =
          traj->GetTrajectoryPoints();

      if (trajPoints.size() > 0) {
        output.Charge = traj->GetParticle()->Charge() / 3.0;

        // Find the momentum / direction of the true particle at the closest
        // point to the position of the recon object's state.
        TVector3 closestOffset(-10000, -10000, -10000);

        for (std::vector<ND::TG4TrajectoryPoint>::iterator trajPointIt =
                 trajPoints.begin();
             trajPointIt != trajPoints.end(); trajPointIt++) {
          ND::TG4TrajectoryPoint trajPoint = (*trajPointIt);

          TVector3 thisOffset = (trajPoint.GetPosition().Vect() - statePos);
          if (fabs(thisOffset.Z()) < fabs(closestOffset.Z())) {
            closestOffset = thisOffset;
            output.Momentum = trajPoint.GetMomentum().Mag();
            output.Direction = trajPoint.GetMomentum().Unit();
            output.Position = trajPoint.GetPosition();
          }
        }

        if (closestOffset.Mag() > 100) {
          ND280Verbose("      Closest MC trajectory point is offset by "
                       << closestOffset.Mag());
        }

        output.SetOK = true;
      }
    } else {
      ND280Warn("      No trajectory number " << G4TrkId);
    }
  } else if (event.GetContext().IsMC()) {
    ND280Warn("      No trajectories, yet is MC!!!");
  }

  return;
}

/// Fill information relating to the true particle associated with a global PID.
/// Fills the initial status of the true trajectory.
void ND::TReconPerformanceEvalModule::FillGlobalTruthInfo(
    ND::TND280Event& event, ND::THandle<ND::TReconBase> reconObject,
    TGlobalTruthInfo& output) {
  ND::THandle<ND::TReconState> state;

  try {
    state = reconObject->GetState();
  } catch (ND::EReconObject) {
    return;
  }

  if (!state) {
    return;
  }

  TVector3 statePos = GetStatePosition(state).Vect();
  ND::THandle<ND::TG4TrajectoryContainer> trajectories =
      event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");

  if (trajectories && reconObject->GetHits()) {
    // Use oaUtility to find true particle corresponding to this track.
    int G4TrkId = TrackTruthInfo::GetG4TrajIDHits(
        *(reconObject->GetHits()), output.Efficiency, output.Purity);
    ND::THandle<ND::TG4Trajectory> traj = trajectories->GetTrajectory(G4TrkId);
    if (traj) {
      ND280Verbose("      Got trajectory number " << G4TrkId);
      output.Charge = traj->GetParticle()->Charge() / 3.0;
      output.Momentum = traj->GetInitialMomentum().Vect().Mag();
      output.Direction = traj->GetInitialMomentum().Vect().Unit();
      output.Position = traj->GetInitialPosition();
      output.SetOK = true;
    } else {
      ND280Warn("      No trajectory number " << G4TrkId);
    }
  } else if (event.GetContext().IsMC()) {
    ND280Warn("      No trajectories, yet is MC!!!");
  }

  return;
}

/// Fill information relating to the state of a global PID constituent.
void ND::TReconPerformanceEvalModule::FillStateInfo(
    ND::THandle<ND::TReconState> state, TReconStateInfo& output) {
  if (!state) {
    return;
  }

  output.Position = GetStatePosition(state);
  output.PositionVariance = GetStatePositionVariance(state);
  output.Momentum = GetStateMomentum(state);
  output.Direction = GetStateDirection(state);

  output.SetOK = true;

  return;
}

/// Fill information relating to the state of a global PID constituent.
void ND::TReconPerformanceEvalModule::FillDownToTrackerInfo(
    ND::THandle<ND::TReconBase> base, TDownToTrackerInfo* output) {
  ND::THandle<ND::TReconState> state;

  try {
    state = base->GetState();
  } catch (ND::EReconObject e) {
    return;
  }

  if (!state) {
    return;
  }

  output->DetectorName = base->ConvertDetector();
  output->MatchingChi2 = TrackingUtils::GetMatchingChi2(*base);
  output->Position = GetStatePosition(state);
  output->Momentum = GetStateMomentum(state);
  output->Direction = GetStateDirection(state);

  output->SetOK = true;

  return;
}

/// Extrapolate an object to the surface of another module. Currently
/// implemented surfaces are:
/// "DS"  - upstream face of DsECal
/// "LT"  - inner face of LTECal
/// "RT"  - inner face of RTECal
/// "TT"  - inner face of TTECal
/// "BT"  - inner face of BTECal
/// "P0D" - downstream face of P0D
/// Extrapolation is done using the tpcRecon HelixModel.
/// The input object must be either a TReconTrack or a TReconPID with a single
/// constituent
/// that is itself a TReconTrack.
void ND::TReconPerformanceEvalModule::ExtrapolateToModuleAndFill(
    ND::THandle<ND::TReconBase> base, THelixModelInfo& output,
    std::string module) {
  if (!base) {
    return;
  }

  // Initialize the tpcRecon HelixModel object.
  HelixModel* helix = new HelixModel();
  ND::THandle<ND::TReconTrack> track = base;

  if (!track) {
    ND::THandle<ND::TReconPID> pid = base;
    ND::THandle<ND::TReconObjectContainer> conContainer =
        base->GetConstituents();
    if (conContainer->size() != 1) {
      ND280Error("Constituent container has > 1 constituent!");
      return;
    }

    for (ND::TReconObjectContainer::iterator conIter = conContainer->begin();
         conIter != conContainer->end(); ++conIter) {
      ND::THandle<ND::TReconBase> base = (*conIter);
      track = base;

      if (!track) {
        ND280Error("Tried to convert non-track to helix");
        return;
      }
    }
  }

  ND::THandle<ND::TTrackState> state;
  try {
    state = track->GetState();
  } catch (ND::EReconObject e) {
    ND280Warn("Couldn't get track state.");
    return;
  }

  double x0 = state->GetPosition().X();
  double y0 = state->GetPosition().Y();
  double z0 = state->GetPosition().Z();

  double tx0 = state->GetDirection().X();
  double ty0 = state->GetDirection().Y() / state->GetDirection().Z();

  double rho = state->GetCurvature();

  helix->SetParameters(x0, y0, z0, tx0, ty0, rho);

  double x, y, z, theta;
  TVector3 minpos, maxpos, pos;

  // Propagate helix to correct co-ordinate. Inner face of TECal modules;
  // upstream face of DsECal; downstream face of P0D.
  if (module == "DS" || module == "P0D") {
    if (module == "DS") {
      GetModuleBox(ND::GeomId::DSECal::Detector(), minpos, maxpos);
      // Slight correction so we extrapolate to ~ first DSECal node
      z = minpos.Z() + 30.5;
    } else {
      GetModuleBox(ND::GeomId::P0D::Detector(), minpos, maxpos);
      // Slight correction so we extrapolate to ~ last P0D node
      z = maxpos.Z() - 42.0;
    }
    x = helix->XCoord(TVector3(0.0, 0.0, z), 2);
    y = helix->YCoord(TVector3(0.0, 0.0, z), 2);
    theta = helix->phiYZ(TVector3(0.0, 0.0, z), 2);
  } else if (module == "TT" || module == "BT") {
    if (module == "TT") {
      GetModuleBox(ND::GeomId::TECal::Module(0, 2), minpos, maxpos);
      // Slight correction so we extrapolate to ~ first TTECal node
      y = minpos.Y() + 45.5;
    } else {
      GetModuleBox(ND::GeomId::TECal::Module(0, 0), minpos, maxpos);
      // Slight correction so we extrapolate to ~ first BTECal node
      y = maxpos.Y() - 45.5;
    }
    x = helix->XCoord(TVector3(0.0, y, 0.0), 1);
    z = helix->ZCoord(TVector3(0.0, y, 0.0), 1);
    theta = helix->phiYZ(TVector3(0.0, 0.0, z), 2);
  } else if (module == "RT" || module == "LT") {
    if (module == "RT") {
      GetModuleBox(ND::GeomId::TECal::Module(0, 1), minpos, maxpos);
      // Slight correction so we extrapolate to ~ first RTECal node
      x = maxpos.X() - 55.5;
    } else {
      GetModuleBox(ND::GeomId::TECal::Module(1, 1), minpos, maxpos);
      // Slight correction so we extrapolate to ~ first LTECal node
      x = minpos.X() + 55.5;
    }
    y = helix->YCoord(TVector3(x, 0.0, 0.0), 0);
    z = helix->ZCoord(TVector3(x, 0.0, 0.0), 0);
    theta = helix->phiYZ(TVector3(0.0, 0.0, z), 2);
  } else if (module == "TPC2") {
    GetModuleBox(ND::GeomId::TPC::TPC2(), minpos, maxpos);
    if (base->GetDetectors() == ND::TReconBase::kTPC1) {
      z = maxpos.Z();
    } else {
      z = minpos.Z();
    }
    x = helix->XCoord(TVector3(0.0, 0.0, z), 2);
    y = helix->YCoord(TVector3(0.0, 0.0, z), 2);
    theta = helix->phiYZ(TVector3(0.0, 0.0, z), 2);
  } else {
    return;
  }

  // Fill the output object with the results of the extrapolation
  output.Position = TVector3(x, y, z);
  output.CosTheta = std::cos(theta);
  output.SetOK = true;

  return;
}

/// Get the minimum and maximum positions of the box defining this TGeometryId.
/// Lovingly stolen from ND::TGeometrySummaryModule::FillBBox().
void ND::TReconPerformanceEvalModule::GetModuleBox(ND::TGeometryId id,
                                                   TVector3& minpos,
                                                   TVector3& maxpos) {
  try {
    TVector3 centre = id.GetPosition();

    ND::TOADatabase::Get().GeomId().CdId(id);
    TGeoShape* vol = gGeoManager->GetCurrentVolume()->GetShape();
    TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(vol);
    if (bbox) {
      TVector3 corner(bbox->GetDX(), bbox->GetDY(), bbox->GetDZ());
      TVector3 globalcorner = ND::TGeomInfo::Get().LocalToMasterVect(
          id, corner.X(), corner.Y(), corner.Z());
      TVector3 min = centre - globalcorner;
      TVector3 max = centre + globalcorner;
      minpos.SetX(std::min(min.X(), max.X()));
      minpos.SetY(std::min(min.Y(), max.Y()));
      minpos.SetZ(std::min(min.Z(), max.Z()));
      maxpos.SetX(std::max(min.X(), max.X()));
      maxpos.SetY(std::max(min.Y(), max.Y()));
      maxpos.SetZ(std::max(min.Z(), max.Z()));
    }
  } catch (ND::EGeomIdInvalid e) {
    // It's possible that we could try passing an invalid geometry id.
    // This happens, for example, if we're looking at run 1 data (or 2010-02
    // MC) and we try to extrapolate to the non-existant Tracker ECal.
    ND280Verbose(" Invalid geometry id passed to GetModuleBox()");
  }

  return;
}

/// Fill information relating to the state of a global PID.
void ND::TReconPerformanceEvalModule::FillGlobalStateInfo(
    ND::THandle<ND::TReconState> state, TGlobalReconStateInfo& output) {
  if (!state) {
    return;
  }

  output.Position = GetStatePosition(state);
  output.PositionVariance = GetStatePositionVariance(state);
  output.Momentum = GetStateMomentum(state);
  output.Direction = GetStateDirection(state);

  output.SetOK = true;

  return;
}

ND::TReconPerformanceEvalModule::TGlobalReconObject::TGlobalReconObject() {
  SetOK = false;

  Constituents = new TClonesArray(
      "ND::TReconPerformanceEvalModule::TGlobalReconConstituent", 200);
  Constituents->Clear();

  NConstituents = 0;
  NModuleConstituents = new std::map<std::string, int>();

  GlobalNodes = new TClonesArray(
      "ND::TReconPerformanceEvalModule::TGlobalReconNodeInfo", 200);
  GlobalNodes->Clear();
  NGlobalNodes = 0;
  NGlobalNodesSaved = 0;

  DownToTrackerConstituents = new TClonesArray(
      "ND::TReconPerformanceEvalModule::TDownToTrackerInfo", 200);
  DownToTrackerConstituents->Clear();
  NDownToTrackerConstituents = 0;

  SubdetectorString = "FAILED";

  ParticleID = "NOT_A_TRECONPID";
  PIDWeight = -10000;

  Momentum = -10000;
  MomentumByRange = -10000;
  Charge = -10000;
  Direction = TVector3(-10000, -10000, -10000);
  UniqueID = 0;

  MatchingChi2Info = new TClonesArray(
      "ND::TReconPerformanceEvalModule::TMatchingChi2Info", 200);
  MatchingChi2Info->Clear();
  NMatchingChi2Info = 0;
}

ND::TReconPerformanceEvalModule::TGlobalReconObject::~TGlobalReconObject() {}

ND::TReconPerformanceEvalModule::TGlobalReconConstituent::
    TGlobalReconConstituent() {
  SetOK = false;
  ReconAlgorithm = "FAILED";
  DetectorName = "FAILED";
  ClassName = "FAILED";
  ParticleID = "NOT_A_TRECONPID";
  PIDWeight = -10000;
  AverageHitTime = -10000;
  UniqueID = 0;
}

ND::TReconPerformanceEvalModule::TGlobalReconConstituent::
    ~TGlobalReconConstituent() {}

ND::TReconPerformanceEvalModule::TGlobalTruthInfo::TGlobalTruthInfo() {
  SetOK = false;
  Position = TLorentzVector(-10000, -10000, -10000, -10000);
  Momentum = -10000;
  Charge = -10000;
  Direction = TVector3(-10000, -10000, -10000);
  Purity = -10000;
  Efficiency = -10000;
}

ND::TReconPerformanceEvalModule::TGlobalTruthInfo::~TGlobalTruthInfo() {}

ND::TReconPerformanceEvalModule::TTruthInfo::TTruthInfo() {
  SetOK = false;
  Position = TLorentzVector(-10000, -10000, -10000, -10000);
  Momentum = -10000;
  Charge = -10000;
  Direction = TVector3(-10000, -10000, -10000);
  Purity = -10000;
  Efficiency = -10000;
}

ND::TReconPerformanceEvalModule::TTruthInfo::~TTruthInfo() {}

ND::TReconPerformanceEvalModule::TGlobalReconNodeInfo::TGlobalReconNodeInfo() {
  SetOK = false;
}

ND::TReconPerformanceEvalModule::TGlobalReconNodeInfo::~TGlobalReconNodeInfo() {
}

ND::TReconPerformanceEvalModule::TGlobalReconStateInfo::
    TGlobalReconStateInfo() {
  SetOK = false;
  Position = TLorentzVector(-10000, -10000, -10000, -10000);
  PositionVariance = TLorentzVector(-10000, -10000, -10000, -10000);
  Momentum = -10000;
  Direction = TVector3(-10000, -10000, -10000);
}

ND::TReconPerformanceEvalModule::TGlobalReconStateInfo::
    ~TGlobalReconStateInfo() {}

ND::TReconPerformanceEvalModule::TReconStateInfo::TReconStateInfo() {
  SetOK = false;
  Position = TLorentzVector(-10000, -10000, -10000, -10000);
  PositionVariance = TLorentzVector(-10000, -10000, -10000, -10000);
  Momentum = -10000;
  Direction = TVector3(-10000, -10000, -10000);
}

ND::TReconPerformanceEvalModule::TReconStateInfo::~TReconStateInfo() {}

ND::TReconPerformanceEvalModule::TDownToTrackerInfo::TDownToTrackerInfo() {
  SetOK = false;
  DetectorName = "FAILED";
  MatchingChi2 = -10000;
  Position = TLorentzVector(-10000, -10000, -10000, -10000);
  Momentum = -10000;
  Direction = TVector3(-10000, -10000, -10000);
}

ND::TReconPerformanceEvalModule::TDownToTrackerInfo::~TDownToTrackerInfo() {}

ND::TReconPerformanceEvalModule::THelixModelInfo::THelixModelInfo() {
  SetOK = false;
  Position = TVector3(-10000, -10000, -10000);
  CosTheta = -10000;
}

ND::TReconPerformanceEvalModule::THelixModelInfo::~THelixModelInfo() {}

ND::TReconPerformanceEvalModule::TMatchingChi2Info::TMatchingChi2Info() {
  SetOK = false;
  Det1 = "FAILED";
  Det2 = "FAILED";
  MatchingChi2 = -10000;
  IsShower1 = false;
  IsShower2 = false;
}

ND::TReconPerformanceEvalModule::TMatchingChi2Info::~TMatchingChi2Info() {}