TSmrdReconModule.cxx

Go to the documentation of this file.

#include <cmath>
#include <string>

#include <TGeomInfo.hxx>
#include <TRecPackManager.hxx>
#include <TReconTrack.hxx>
#include <TTrackState.hxx>

#include <TrackTruthInfo.hxx>

#include "TSmrdReconModule.hxx"

#include "TTFBChannelId.hxx"
#include "TTFBDigit.hxx"

#include <TRealDatum.hxx>

#include <TrackingUtils.hxx>

ClassImp(ND::TSmrdReconModule::TSmrdReconHit);
ND::TSmrdReconModule::TSmrdReconHit::~TSmrdReconHit() {}
ClassImp(ND::TSmrdReconModule::TSmrdIsoTrack);
ND::TSmrdReconModule::TSmrdIsoTrack::~TSmrdIsoTrack() {}

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

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

  fIsEnabled = kTRUE;
  fDescription = "SMRD Recon Output";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;

  fNSmrdReconHits = 0;
  fSmrdReconHits = new TClonesArray("ND::TSmrdReconModule::TSmrdReconHit", 200);
  fNSmrdIsoTracks = 0;
  fSmrdIsoTracks = new TClonesArray("ND::TSmrdReconModule::TSmrdIsoTrack", 200);
}

ND::TSmrdReconModule::~TSmrdReconModule() {}

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

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

void ND::TSmrdReconModule::InitializeBranches() {
  fOutputTree->Branch("NSmrdReconHits", &fNSmrdReconHits, "NSmrdReconHits/I",
      fBufferSize)
  ->SetTitle("The number of SMRD TReconHits.");
  fOutputTree->Branch("NSmrdIsoTracks", &fNSmrdIsoTracks, "NSmrdIsoTracks/I",
                      fBufferSize)
  ->SetTitle("The SMRD TReconHits.");
  fOutputTree->Branch("SmrdReconHits", &fSmrdReconHits, fBufferSize,
                      fSplitLevel)
  ->SetTitle("The number of SMRD isolated tracks.");
  fOutputTree->Branch("SmrdIsoTracks", &fSmrdIsoTracks, fBufferSize,
                      fSplitLevel)
  ->SetTitle("The SMRD isolated tracks.");
}

bool ND::TSmrdReconModule::FillTree(ND::TND280Event &event) {
  fNSmrdReconHits = 0;
  fSmrdReconHits->Clear();

  fNSmrdIsoTracks = 0;
  fSmrdIsoTracks->Clear();

  ND::THandle<ND::TAlgorithmResult> SmrdResult = event.GetFit("smrdRecon");
  ND::THandle<ND::THitSelection> smrdReconHitsMatched;
  ND::THandle<ND::THitSelection> smrdReconHitsUnused;
  ND::THandle<ND::THitSelection> smrdReconHitsUsed;
  ND::THandle<ND::TReconObjectContainer> xyzReconTracks;

  if (SmrdResult) {
    smrdReconHitsUnused = SmrdResult->GetHitSelection("unused");
    smrdReconHitsMatched = SmrdResult->GetHitSelection("matchedSmrdReconHits");
    smrdReconHitsUsed = SmrdResult->GetHitSelection("used");
    xyzReconTracks = SmrdResult->GetResultsContainer("final");
  }

  //----------------------------------------------------------------------------
  // Fill hits information
  //----------------------------------------------------------------------------
  // smrdReconHitsUnused
  if (smrdReconHitsUnused) {
    for (ND::THitSelection::iterator h = smrdReconHitsUnused->begin();
         h != smrdReconHitsUnused->end(); ++h) {
      ND::THandle<ND::THit> hit = *h;
      if (!hit) {
        continue;
      }

      // get hit contributors
      // only deal with "fully-calibrated" two-sided hits
      if ((*h)->GetContributorCount() != 2) {
        ND280Warn("warning: more or less than two smrd single hits contributed"
                  << " to this TReconHit; skipping this TReconHit");
        continue;
      }

      //******************get two single hits-contributors
      ND::THandle<ND::THit> contrhit1 = (*h)->GetContributor(0);
      ND::THandle<ND::THit> contrhit2 = (*h)->GetContributor(1);

      if ((!contrhit1) || (!contrhit2)) {
        ND280Warn("warning: failed get constituents for SMRD TRecon hit; skip "
                  << "this TReconHit");
        continue;
      }

      TSmrdReconHit *reconHit =
          new ((*fSmrdReconHits)[fNSmrdReconHits++]) TSmrdReconHit;

      reconHit->Position = TLorentzVector(hit->GetPosition(), hit->GetTime());
      reconHit->PositionUncertainty =
          TLorentzVector(hit->GetUncertainty(), hit->GetTimeUncertainty());
      reconHit->Charge = hit->GetCharge();

      reconHit->dZ = (*h)->GetPosition().Z() - contrhit1->GetPosition().Z();
      reconHit->dT = contrhit1->GetTime() - contrhit2->GetTime();
      reconHit->ContribCharge[0] = contrhit1->GetCharge();
      reconHit->ContribCharge[1] = contrhit2->GetCharge();

      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) {
        reconHit->Wall = kTop;
      }
      if (isBottomWall) {
        reconHit->Wall = kBottom;
      }
      if (isRightWall) {
        reconHit->Wall = kRight;
      }
      if (isLeftWall) {
        reconHit->Wall = kLeft;
      }

      reconHit->Yoke = ND::TGeomInfo::SMRD().GetYokeRingNumber(hit);
      reconHit->Layer = ND::TGeomInfo::SMRD().GetSMRDRingNumber(hit);
      reconHit->Tower = ND::TGeomInfo::SMRD().GetSMRDTowerNumber(hit);
      reconHit->Counter = ND::TGeomInfo::SMRD().GetSMRDScintNumber(hit);

      // fill tfb and rmm info
      ND::TChannelId chan = contrhit1->GetChannelId();
      ND::TTFBChannelId tfbid(chan);

      reconHit->RMM = tfbid.GetRMM();
      reconHit->TFB = tfbid.GetTFB();

      reconHit->IsUsed = false;

      reconHit->IsInnerMatched = false;

      // should be always unmatched but let`s double check
      if (smrdReconHitsMatched) {
        if (std::find(smrdReconHitsMatched->begin(),
                      smrdReconHitsMatched->end(),
                      *h) != smrdReconHitsMatched->end()) {
          reconHit->IsInnerMatched = true;
        }
      }
    }  // end of hits loop
  }    // end of if(smrdReconHitsUnused)

  // smrdReconHitsUsed
  if (smrdReconHitsUsed) {
    for (ND::THitSelection::iterator h = smrdReconHitsUsed->begin();
         h != smrdReconHitsUsed->end(); h++) {
      ND::THandle<ND::THit> hit = *h;
      if (!hit) {
        continue;
      }

      // get hit contributors
      // only deal with "fully-calibrated" two-sided hits
      if ((*h)->GetContributorCount() != 2) {
        ND280Warn("warning: more or less than two smrd single hits contributed"
                  << " to this TReconHit; skipping this TReconHit");
        continue;
      }

      //******************get two single hits-contributors
      ND::THandle<ND::THit> contrhit1 = (*h)->GetContributor(0);
      ND::THandle<ND::THit> contrhit2 = (*h)->GetContributor(1);

      if (!contrhit1 || !contrhit2) {
        ND280Warn("warning: failed get constituents for SMRD TRecon hit; "
                  << "skip this TReconHit");
        continue;
      }

      TSmrdReconHit *reconHit =
          new ((*fSmrdReconHits)[fNSmrdReconHits++]) TSmrdReconHit;

      reconHit->Position = TLorentzVector(hit->GetPosition(), hit->GetTime());
      reconHit->PositionUncertainty =
          TLorentzVector(hit->GetUncertainty(), hit->GetTimeUncertainty());
      reconHit->Charge = hit->GetCharge();

      reconHit->dZ = (*h)->GetPosition().Z() - contrhit1->GetPosition().Z();
      reconHit->dT = contrhit1->GetTime() - contrhit2->GetTime();

      reconHit->ContribCharge[0] = contrhit1->GetCharge();
      reconHit->ContribCharge[1] = contrhit2->GetCharge();

      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) {
        reconHit->Wall = kTop;
      }
      if (isBottomWall) {
        reconHit->Wall = kBottom;
      }
      if (isRightWall) {
        reconHit->Wall = kRight;
      }
      if (isLeftWall) {
        reconHit->Wall = kLeft;
      }

      reconHit->Yoke = ND::TGeomInfo::SMRD().GetYokeRingNumber(hit);
      reconHit->Layer = ND::TGeomInfo::SMRD().GetSMRDRingNumber(hit);
      reconHit->Tower = ND::TGeomInfo::SMRD().GetSMRDTowerNumber(hit);
      reconHit->Counter = ND::TGeomInfo::SMRD().GetSMRDScintNumber(hit);

      // fill tfb and rmm info
      ND::TChannelId chan = contrhit1->GetChannelId();
      ND::TTFBChannelId tfbid(chan);

      reconHit->RMM = tfbid.GetRMM();
      reconHit->TFB = tfbid.GetTFB();

      reconHit->IsUsed = true;

      reconHit->IsInnerMatched = false;

      if (smrdReconHitsMatched) {
        if (std::find(smrdReconHitsMatched->begin(),
                      smrdReconHitsMatched->end(),
                      *h) != smrdReconHitsMatched->end()) {
          reconHit->IsInnerMatched = true;
        }
      }
    }  // end of hits loop
  }    // end of if(smrdReconHitsUsed)

  //----------------------------------------------------------------------------
  // Fill tracks information
  //----------------------------------------------------------------------------
  if (xyzReconTracks) {
    bool isMC = event.GetContext().IsMC();

    for (ND::TReconObjectContainer::iterator t = xyzReconTracks->begin();
         t != xyzReconTracks->end(); t++) {
      ND::THandle<ND::TReconTrack> track = *t;
      if (!track) {
        continue;
      }

      ND::THandle<ND::TTrackState> trackState = track->GetState();
      if (!trackState) {
        continue;
      }

      TSmrdIsoTrack *isoTrack =
          new ((*fSmrdIsoTracks)[fNSmrdIsoTracks++]) TSmrdIsoTrack;

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

      ND::THandle<ND::TTrackState> frontState =
          TrackingUtils::GetFirstState(*track);
      if (!frontState) {
        continue;
      }

      ND::THandle<ND::TTrackState> backState =
          TrackingUtils::GetLastState(*track);
      if (!backState) {
        continue;
      }

      ND::THandle<ND::THitSelection> hits = track->GetHits();

      if (!hits) {
        continue;
      }

      isoTrack->AlgorithmName = track->GetAlgorithmName();
      isoTrack->Chi2 = track->GetQuality();

      // some of all contributed hits charges
      isoTrack->EDeposit = -9999.0;
      ND::THandle<ND::TRealDatum> totalCharge =
          track->Get<ND::TRealDatum>("totalCharge");
      if (totalCharge) {
        isoTrack->EDeposit =
            track->Get<ND::TRealDatum>("totalCharge")->GetValue();
      }

      // hits averaged time
      isoTrack->avgtime = -9999.0;
      ND::THandle<ND::TRealDatum> smrdTimeRD =
          track->Get<ND::TRealDatum>("averageHitTime");
      if (smrdTimeRD) {
        isoTrack->avgtime =
            track->Get<ND::TRealDatum>("averageHitTime")->GetValue();
      }

      isoTrack->NHits = hits->size();
      isoTrack->NNodes = nodes.size();
      isoTrack->Status = track->CheckStatus(track->kSuccess);
      isoTrack->NDOF = track->GetNDOF();
      isoTrack->UniqueID = track->GetUniqueID();

      // Check the Kalman refitting status
      std::string name("RECPACK_KalmanFilter");
      if ((track->GetAlgorithmName().find(name) != std::string::npos) &&
          (track->CheckStatus(track->kSuccess))) {
        isoTrack->KalmanStatus = 1;
      } else {
        isoTrack->KalmanStatus = 0;
      }

      isoTrack->FrontPos = frontState->GetPosition();
      isoTrack->FrontPosVariance = frontState->GetPositionVariance();

      isoTrack->BackPos = backState->GetPosition();
      isoTrack->BackPosVariance = backState->GetPositionVariance();

      TVector3 backpos = backState->GetPosition().Vect();
      TVector3 frontpos = frontState->GetPosition().Vect();

      // isoTrack->Direction = (backpos - frontpos).Unit();
      isoTrack->Direction = trackState->GetDirection();
      isoTrack->DirectionVariance = trackState->GetDirectionVariance();

      // Calculate length of the track assuming straight line
      TVector3 diff = backpos - frontpos;
      isoTrack->Range = diff.Mag();

      // Angles...
      isoTrack->ThetaAngle = trackState->GetDirection().Theta();
      isoTrack->PhiAngle = trackState->GetDirection().Phi();

      // Set truth info
      int trueId = -999;
      TLorentzVector trueInitPos(-9999, -9999, -9999, -9999);
      TLorentzVector trueFinalPos(-9999, -9999, -9999, -9999);
      TLorentzVector trueInitMom(-9999, -9999, -9999, -9999);
      int truePDG = 0;
      int trueParentId = -999;
      double truePurity = -9999.0;
      double trueEff = -9999.0;

      if (isMC) {
        trueId = TrackTruthInfo::GetG4TrajIDHits(*hits, trueEff, truePurity);

        ND::THandle<ND::TG4TrajectoryContainer> trajectoryContainer =
            event.Get<ND::TG4TrajectoryContainer>("truth/G4Trajectories");
        ND::THandle<ND::TG4Trajectory> g4traj(0);
        if (trajectoryContainer)
          g4traj = trajectoryContainer->GetTrajectory(trueId);

        if (g4traj) {
          trueInitPos = g4traj->GetInitialPosition();
          trueFinalPos = g4traj->GetFinalPosition();
          trueInitMom = g4traj->GetInitialMomentum();
          truePDG = g4traj->GetPDGEncoding();
          trueParentId = g4traj->GetParentId();
        }
      }

      isoTrack->TrueInitPos = trueInitPos;

      gGeoManager->FindNode(trueInitPos.X(), trueInitPos.Y(), trueInitPos.Z());
      gGeoManager->GetCurrentNodeId();
      isoTrack->TrueInitDet =
          eventAnalysisEnums::PathToSubdetector(gGeoManager->GetPath());

      isoTrack->TrueFinalPos = trueFinalPos;

      gGeoManager->FindNode(trueFinalPos.X(), trueFinalPos.Y(),
                            trueFinalPos.Z());
      gGeoManager->GetCurrentNodeId();
      isoTrack->TrueFinalDet =
          eventAnalysisEnums::PathToSubdetector(gGeoManager->GetPath());

      isoTrack->TrueInitMom = trueInitMom;
      isoTrack->TrueId = trueId;
      isoTrack->TruePDG = truePDG;
      isoTrack->TrueParentId = trueParentId;
      // isoTrack->TrueParentPDG =;
      isoTrack->TrueHitPurity = truePurity;
      isoTrack->TrueHitEff = trueEff;

    }  // end of tracks loop
  }    // end of if(xyzTracks)
  return true;
}