TrackActivityRemover.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file: TrackActivityRemover.cxx
3 // \brief Remove hits associated with track activity, including track
4 // hits, hits surrounding all track bodies, and hits surrounding
5 // the end points of stopping tracks (Michels).
6 // \author Justin Vasel <justin.vasel@gmail.com>
7 // \date 2019-11-21
8 ////////////////////////////////////////////////////////////////////////
9 
10 // ART includes
16 #include "fhiclcpp/ParameterSet.h"
17 
18 // ROOT includes
19 #include "TVector3.h"
20 
21 // NOvASOFT includes
22 #include "Geometry/Geometry.h"
25 #include "RawData/RawDigit.h"
26 #include "RecoBase/CellHit.h"
27 #include "RecoBase/Cluster.h"
28 #include "RecoBase/Track.h"
31 
32 // ............................................................................
34 fVetoMap(vetoMap),
35 fGeom(geom),
36 fClusterFromTrack(clusterFromTrack),
37 fMuonEndTimeCut(p.get<int>("MuonEndTimeCut")),
38 fMuonBodyTimeCut(p.get<int>("MuonBodyTimeCut")),
39 fMuonEndSphereRadius(p.get<int>("MuonEndSphereRadius")),
40 fMuonBodyCylinderRadius(p.get<int>("MuonBodyCylinderRadius")),
41 fFiducialVolume(p.get<fhicl::ParameterSet>("FiducialVolume")),
42 fFVMinX(fFiducialVolume.get<int>("MinX")),
43 fFVMaxX(fFiducialVolume.get<int>("MaxX")),
44 fFVMinY(fFiducialVolume.get<int>("MinY")),
45 fFVMaxY(fFiducialVolume.get<int>("MaxY")),
46 fFVMinZ(fFiducialVolume.get<int>("MinZ")),
47 fFVMaxZ(fFiducialVolume.get<int>("MaxZ")),
48 fNumRemovedMichel(0),
49 fNumRemovedBody(0),
50 fNumRemovedTrack(0)
51 {
52 }
53 
54 // ............................................................................
56 
57 
58 // ............................................................................
59 void sn::TrackActivityRemover::SetVetoRegion(int muonEndTimeCut, int muonBodyTimeCut, int muonEndSphereRadius, int muonBodyCylinderRadius)
60 {
61  this->fMuonEndTimeCut = muonEndTimeCut;
62  this->fMuonBodyTimeCut = muonBodyTimeCut;
63  this->fMuonEndSphereRadius = muonEndSphereRadius;
64  this->fMuonBodyCylinderRadius = muonBodyCylinderRadius;
65 
66  return;
67 }
68 
69 
70 // ............................................................................
72 {
73  // Get h and track vectors in XYZ-space
74  TVector3 r0;
75 
76  // Note: r0 is passed by reference here, and thus assigned at this line
77  fGeom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(r0);
78 
79  // Define ZW coordinate system based on the view of the hit
80  TVector3 w(0,0,0);
81  if (h->View() == geo::kX) w.SetX(1);
82  if (h->View() == geo::kY) w.SetY(1);
83 
84  // Transform from XYZ-space into ZW-space
85  TVector3 w0(r0.Z(), r0.Dot(w), 0);
86  TVector3 w1( p.Z(), p.Dot(w), 0);
87 
88  return (w1-w0).Mag();
89 }
90 
91 
92 // ............................................................................
94 {
95  // Get hit position in XYZ-space, transform into ZW coordinate system
96  TVector3 r0;
97  fGeom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(r0);
98 
99  // Define ZW coordinate system based on the view of the hit
100  TVector3 w(0,0,0);
101  if (h->View() == geo::kX) w.SetX(1);
102  if (h->View() == geo::kY) w.SetY(1);
103 
104  // From here on out, the W coordinate maps onto either detector X or detector Y
105  // depending on the view
106  double wHitZ = r0.Z();
107  double wHitW = r0.Dot(w);
108 
109  // Loop over track trajectory positions and compute min. distance to hit
110  std::vector<TVector3> const & trackTrajPoints = t->Trajectory();
111  int min_square_distance = -1;
112  for (TVector3 trajPoint : trackTrajPoints) {
113  double wTrajZ = trajPoint.Z();
114  double wTrajW = trajPoint.Dot(w);
115  double square_distance = std::pow((wTrajZ - wHitZ), 2) + std::pow((wTrajW - wHitW), 2);
116 
117  // Only update min_square_distance if it hasn't been set yet, or if the
118  // square_distance is smaller than it.
119  if (min_square_distance < 0 || square_distance < min_square_distance) {
120  min_square_distance = square_distance;
121  }
122 
123  }
124 
125  return std::sqrt(min_square_distance);
126 }
127 
128 
129 // ............................................................................
131 {
132  return (p.X() > fFVMinX) && (p.X() < fFVMaxX) &&
133  (p.Y() > fFVMinY) && (p.Y() < fFVMaxY) &&
134  (p.Z() > fFVMinZ) && (p.Z() < fFVMaxZ);
135 }
136 
137 
138 // ............................................................................
140 {
141  if (!PointIsContained(track->Start()) && PointIsContained(track->Stop())) return sn::TrackType::kStopping;
142  else if (!PointIsContained(track->Start()) && !PointIsContained(track->Stop())) return sn::TrackType::kThroughGoing;
143  else if ( PointIsContained(track->Start()) && PointIsContained(track->Stop())) return sn::TrackType::kContained;
144  else if ( PointIsContained(track->Start()) && !PointIsContained(track->Stop())) return sn::TrackType::kExiting;
145  else return sn::TrackType::kUnknown;
146 }
147 
148 
149 // ............................................................................
150 void sn::TrackActivityRemover::remove(std::vector<art::Ptr<rb::CellHit>>& hits, art::Handle<std::vector<rb::Track>>& hdlTracks)
151 {
152  std::vector<art::Ptr<rb::Track>> tracks;
153  art::fill_ptr_vector(tracks, hdlTracks);
154  rb::SortTracksByTime(tracks);
156 
157  for (size_t trackId=0; trackId<tracks.size(); ++trackId) {
158  art::Ptr<rb::Track> track = tracks.at(trackId);
159 
160  art::PtrVector<rb::CellHit> trackHits = this->fClusterFromTrack.at(trackId)->AllCells();
161  for (art::Ptr<rb::CellHit> trackHit : trackHits) {
162  this->fVetoMap.AddHit(*trackHit);
163  ++this->fNumRemovedTrack;
164  }
165 
166  sn::TrackType trackType = this->GetTrackType(track);
167 
168  double trackMinTNS = track->MinTNS();
169  double trackMaxTNS = track->MaxTNS();
170 
171  /* Find nearby hits */
172  unsigned int idxStartCellHit = 0;
173  for (size_t idx=idxStartCellHit; idx<hits.size(); ++idx) {
175  double hitTNS = hit->TNS();
176 
177  // Case: Hit is earlier than track
178  if (hitTNS < trackMinTNS) {
179  idxStartCellHit = idx;
180  continue;
181  }
182 
183  // Case: Hit is part of the track, and has already been removed
184  if (std::find(trackHits.begin(), trackHits.end(), hit) != trackHits.end()) {
185  idxStartCellHit = idx;
186  continue;
187  }
188 
189  bool kPassedHitTrackEndTime = hitTNS <= (trackMaxTNS + fMuonEndTimeCut);
190  bool kPassedHitTrackBodyTime = hitTNS <= (trackMaxTNS + fMuonBodyTimeCut);
191  bool kPassedHitTrackEndSpace = this->DistanceHitToPoint(hit, track->Stop()) <= fMuonEndSphereRadius;
192  bool kPassedHitTrackBodySpace = this->DistanceHitToTrackBody(hit, track) <= fMuonBodyCylinderRadius;
193 
194  // Case: Hit is outside of the temporal region of the track
195  if (!kPassedHitTrackEndTime && !kPassedHitTrackBodyTime) {
196  break;
197  }
198 
199  // Case: Hit is outside of the spatial region of the track
200  if (!kPassedHitTrackEndSpace && !kPassedHitTrackBodySpace) {
201  continue;
202  }
203 
204  // Case: Hit is in the Michel veto region
205  if (kPassedHitTrackEndTime && kPassedHitTrackEndSpace && trackType == sn::TrackType::kStopping) {
206  this->fVetoMap.AddHit(*hit);
207  ++this->fNumRemovedMichel;
208  continue;
209  }
210 
211  // Case: Hit is within the track body veto region
212  if (kPassedHitTrackBodyTime && kPassedHitTrackBodySpace) {
213  this->fVetoMap.AddHit(*hit);
214  ++this->fNumRemovedBody;
215  }
216 
217  }
218  }
219 
220  return;
221 }
float TNS() const
Definition: CellHit.h:46
void SetVetoRegion(int muonEndTimeCut, int muonBodyTimeCut, int muonEndSphereRadius, int muonBodyCylinderRadius)
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
bool PointIsContained(const TVector3 &p)
Track enters the detector and stops inside.
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
TrackActivityRemover(sn::HitVetoMap &vetoMap, fhicl::ParameterSet const &p, art::ServiceHandle< geo::Geometry > &geom, art::FindOneP< rb::Cluster > &clusterFromTrack)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::FindOneP< rb::Cluster > & fClusterFromTrack
int fMuonBodyTimeCut
Min. time difference between hit and track to reject [ns].
unsigned short Cell() const
Definition: CellHit.h:40
double MinTNS() const
Definition: Cluster.cxx:482
void SortTracksByTime(std::vector< art::Ptr< rb::Track >> &t)
Sort t in time order (earliest to latest).
Definition: Track.cxx:539
Track starts inside the detector and exits.
void hits()
Definition: readHits.C:15
double DistanceHitToTrackBody(const art::Ptr< rb::CellHit > &h, const art::Ptr< rb::Track > &t)
const Cut kContained([](const caf::SRProxy *sr){return(sr->slc.boxmin.Z() > 100 && sr->slc.boxmax.Z()< 800);})
2 di-block era contained-type event selector
Definition: Cuts.h:42
void AddHit(rawdata::RawDigit h)
Definition: HitVetoMap.cxx:27
sn::TrackType GetTrackType(const art::Ptr< rb::Track > &track)
int fMuonEndSphereRadius
Min. distance of hit from track end to reject [cm].
art::ServiceHandle< geo::Geometry > & fGeom
double MaxTNS() const
Definition: Cluster.cxx:528
int fMuonBodyCylinderRadius
Min. distance of hit from track to reject [cm].
Definition: event.h:1
int fMuonEndTimeCut
Min. time difference between hit and track end to reject [ns].
double DistanceHitToPoint(const art::Ptr< rb::CellHit > &h, const TVector3 &p)
void geom(int which=0)
Definition: geom.C:163
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
const Cut kThroughGoing([](const caf::SRProxy *sr){return(sr->slc.boxmin.Z()< 100 && sr->slc.boxmax.Z() > 800);})
2 di-block era through-going-type event selector
Definition: Cuts.h:45
void remove(std::vector< art::Ptr< rb::CellHit >> &hits, art::Handle< std::vector< rb::Track >> &hdlTracks)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Float_t w
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)