TrackActivityRemover.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file:   TrackActivityRemover.cxx
// \brief   Remove hits associated with track activity, including track
//          hits, hits surrounding all track bodies, and hits surrounding
//          the end points of stopping tracks (Michels).
// \author  Justin Vasel <justin.vasel@gmail.com>
// \date    2019-11-21
////////////////////////////////////////////////////////////////////////

// ART includes
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "fhiclcpp/ParameterSet.h"

// ROOT includes
#include "TVector3.h"

// NOvASOFT includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "RawData/RawDigit.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Supernova/BackgroundRejection/HitVetoMap.h"
#include "Supernova/BackgroundRejection/TrackActivityRemover.h"

// ............................................................................
sn::TrackActivityRemover::TrackActivityRemover(sn::HitVetoMap& vetoMap, fhicl::ParameterSet const & p, art::ServiceHandle<geo::Geometry>& geom, art::FindOneP<rb::Cluster>& clusterFromTrack):
fVetoMap(vetoMap),
fGeom(geom),
fClusterFromTrack(clusterFromTrack),
fMuonEndTimeCut(p.get<int>("MuonEndTimeCut")),
fMuonBodyTimeCut(p.get<int>("MuonBodyTimeCut")),
fMuonEndSphereRadius(p.get<int>("MuonEndSphereRadius")),
fMuonBodyCylinderRadius(p.get<int>("MuonBodyCylinderRadius")),
fFiducialVolume(p.get<fhicl::ParameterSet>("FiducialVolume")),
fFVMinX(fFiducialVolume.get<int>("MinX")),
fFVMaxX(fFiducialVolume.get<int>("MaxX")),
fFVMinY(fFiducialVolume.get<int>("MinY")),
fFVMaxY(fFiducialVolume.get<int>("MaxY")),
fFVMinZ(fFiducialVolume.get<int>("MinZ")),
fFVMaxZ(fFiducialVolume.get<int>("MaxZ")),
fNumRemovedMichel(0),
fNumRemovedBody(0),
fNumRemovedTrack(0)
{
}

// ............................................................................
sn::TrackActivityRemover::~TrackActivityRemover() {}


// ............................................................................
void sn::TrackActivityRemover::SetVetoRegion(int muonEndTimeCut, int muonBodyTimeCut, int muonEndSphereRadius, int muonBodyCylinderRadius)
{
  this->fMuonEndTimeCut         = muonEndTimeCut;
  this->fMuonBodyTimeCut        = muonBodyTimeCut;
  this->fMuonEndSphereRadius    = muonEndSphereRadius;
  this->fMuonBodyCylinderRadius = muonBodyCylinderRadius;

  return;
}


// ............................................................................
double sn::TrackActivityRemover::DistanceHitToPoint(const art::Ptr<rb::CellHit>& h, const TVector3& p)
{
  // Get h and track vectors in XYZ-space
  TVector3 r0;

  // Note: r0 is passed by reference here, and thus assigned at this line
  fGeom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(r0);

  // Define ZW coordinate system based on the view of the hit
  TVector3 w(0,0,0);
  if (h->View() == geo::kX) w.SetX(1);
  if (h->View() == geo::kY) w.SetY(1);

  // Transform from XYZ-space into ZW-space
  TVector3 w0(r0.Z(), r0.Dot(w), 0);
  TVector3 w1( p.Z(),  p.Dot(w), 0);

  return (w1-w0).Mag();
}


// ............................................................................
double sn::TrackActivityRemover::DistanceHitToTrackBody(const art::Ptr<rb::CellHit>& h, const art::Ptr<rb::Track>& t)
{
  // Get hit position in XYZ-space, transform into ZW coordinate system
  TVector3 r0;
  fGeom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(r0);

  // Define ZW coordinate system based on the view of the hit
  TVector3 w(0,0,0);
  if (h->View() == geo::kX) w.SetX(1);
  if (h->View() == geo::kY) w.SetY(1);

  // From here on out, the W coordinate maps onto either detector X or detector Y
  // depending on the view
  double wHitZ = r0.Z();
  double wHitW = r0.Dot(w);

  // Loop over track trajectory positions and compute min. distance to hit
  std::vector<TVector3> const & trackTrajPoints = t->Trajectory();
  int min_square_distance = -1;
  for (TVector3 trajPoint : trackTrajPoints) {
    double wTrajZ = trajPoint.Z();
    double wTrajW = trajPoint.Dot(w);
    double square_distance = std::pow((wTrajZ - wHitZ), 2) + std::pow((wTrajW - wHitW), 2);

    // Only update min_square_distance if it hasn't been set yet, or if the
    // square_distance is smaller than it.
    if (min_square_distance < 0 || square_distance < min_square_distance) {
      min_square_distance = square_distance;
    }

  }

  return std::sqrt(min_square_distance);
}


// ............................................................................
bool sn::TrackActivityRemover::PointIsContained(const TVector3& p)
{
  return (p.X() > fFVMinX) && (p.X() < fFVMaxX) &&
         (p.Y() > fFVMinY) && (p.Y() < fFVMaxY) &&
         (p.Z() > fFVMinZ) && (p.Z() < fFVMaxZ);
}


// ............................................................................
sn::TrackType sn::TrackActivityRemover::GetTrackType(const art::Ptr<rb::Track>& track)
{
  if      (!PointIsContained(track->Start()) &&  PointIsContained(track->Stop())) return sn::TrackType::kStopping;
  else if (!PointIsContained(track->Start()) && !PointIsContained(track->Stop())) return sn::TrackType::kThroughGoing;
  else if ( PointIsContained(track->Start()) &&  PointIsContained(track->Stop())) return sn::TrackType::kContained;
  else if ( PointIsContained(track->Start()) && !PointIsContained(track->Stop())) return sn::TrackType::kExiting;
  else                                                                            return sn::TrackType::kUnknown;
}


// ............................................................................
void sn::TrackActivityRemover::remove(std::vector<art::Ptr<rb::CellHit>>& hits, art::Handle<std::vector<rb::Track>>& hdlTracks)
{
  std::vector<art::Ptr<rb::Track>> tracks;
  art::fill_ptr_vector(tracks, hdlTracks);
  rb::SortTracksByTime(tracks);
  rb::SortByTime(hits);

  for (size_t trackId=0; trackId<tracks.size(); ++trackId) {
    art::Ptr<rb::Track> track = tracks.at(trackId);

    art::PtrVector<rb::CellHit> trackHits = this->fClusterFromTrack.at(trackId)->AllCells();
    for (art::Ptr<rb::CellHit> trackHit : trackHits) {
      this->fVetoMap.AddHit(*trackHit);
      ++this->fNumRemovedTrack;
    }

    sn::TrackType trackType = this->GetTrackType(track);

    double trackMinTNS = track->MinTNS();
    double trackMaxTNS = track->MaxTNS();

    /* Find nearby hits */
    unsigned int idxStartCellHit = 0;
    for (size_t idx=idxStartCellHit; idx<hits.size(); ++idx) {
      art::Ptr<rb::CellHit> hit = hits.at(idx);
      double hitTNS = hit->TNS();

      // Case: Hit is earlier than track
      if (hitTNS < trackMinTNS) {
        idxStartCellHit = idx;
        continue;
      }

      // Case: Hit is part of the track, and has already been removed
      if (std::find(trackHits.begin(), trackHits.end(), hit) != trackHits.end()) {
        idxStartCellHit = idx;
        continue;
      }

      bool kPassedHitTrackEndTime   = hitTNS <= (trackMaxTNS + fMuonEndTimeCut);
      bool kPassedHitTrackBodyTime  = hitTNS <= (trackMaxTNS + fMuonBodyTimeCut);
      bool kPassedHitTrackEndSpace  = this->DistanceHitToPoint(hit, track->Stop()) <= fMuonEndSphereRadius;
      bool kPassedHitTrackBodySpace = this->DistanceHitToTrackBody(hit, track) <= fMuonBodyCylinderRadius;

      // Case: Hit is outside of the temporal region of the track
      if (!kPassedHitTrackEndTime && !kPassedHitTrackBodyTime) {
        break;
      }

      // Case: Hit is outside of the spatial region of the track
      if (!kPassedHitTrackEndSpace && !kPassedHitTrackBodySpace) {
        continue;
      }

      // Case: Hit is in the Michel veto region
      if (kPassedHitTrackEndTime && kPassedHitTrackEndSpace && trackType == sn::TrackType::kStopping) {
        this->fVetoMap.AddHit(*hit);
        ++this->fNumRemovedMichel;
        continue;
      }

      // Case: Hit is within the track body veto region
      if (kPassedHitTrackBodyTime && kPassedHitTrackBodySpace) {
        this->fVetoMap.AddHit(*hit);
        ++this->fNumRemovedBody;
      }

    }
  }

  return;
}