CosmicTrackSelection.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 /// \brief Cosmic track selection code
3 /// \author Brian Rebel - brebel@fnal.gov
4 /// \date
5 //////////////////////////////////////////////////////////////////////////
6 #include <vector>
7 #include <initializer_list>
8 
9 #include "TVector3.h"
10 
12 
14 #include "RecoBase/Track.h"
15 #include "RecoBase/Cluster.h"
17 
18 namespace trk{
19 
20  //----------------------------------------------------------------------------
22  {
24  this->reconfigure(pset);
25 
26  return;
27  }
28 
29  //----------------------------------------------------------------------------
31  {
32  }
33 
34  //----------------------------------------------------------------------------
36  {
37  fMaxDeltaEndPlane = pset.get< double >("MaxDeltaEndPlane", 3 );
38  fMaxPlaneAsymmetry = pset.get< double >("MaxPlaneAsymmetry", 0.1 );
39  fMaxCellsPerPlane = pset.get< double >("MaxCellsPerPlane", 6 );
40  fBadStepSizeLimit = pset.get< double >("BadStepSizeLimit", 3. );
41  fMaxCosTheta = pset.get< double >("MaxCosTheta", 0.95 );
42  fTrkLengthCut = pset.get< double >("TrkLengthCut", 50. );
43  fMaxVtxDistFromEdge = pset.get< double >("MaxVtxDistFromEdge", 10. );
44  fMaxEndDistFromEdge = pset.get< double >("MaxEndDistFromEdge", 10. );
45 
46  return;
47  }
48 
49  //----------------------------------------------------------------------------
50  bool CosmicTrackSelection::Stopper(TVector3 const& finalPoint)
51  {
52  // call it a stopper if the track ends at least 50 cm from an edge of
53  // the detector
54  // distEdge will be negative if the point is outside of the detector
55 
56  return (this->DistToEdge(finalPoint) > 50.);
57  }
58 
59  //----------------------------------------------------------------------------
60  // tracks that have huge changes in the step size from one trajectory
61  // point to another are probably poorly reconstructed, and should be removed
62  // from the analysis. Usually the biggest steps are at either the start
63  // or end of the track in cases of really bad reconstruction
65  std::vector<float> & steps,
66  float & median)
67  {
68 
69  // first find the median step size
70  size_t size = track.NTrajectoryPoints() - 1;
71 
72  // track with only two distances between trajectory points,
73  // probably not very useful for the analysis so return false
74  if(size < 3) return false;
75 
76  std::vector<float> stepsize(size, 0.);
77  float largestStep = -1.;
78  for(size_t t = 1; t < size+1; ++t){
79  TVector3 const& tp1 = track.TrajectoryPoint(t-1);
80  TVector3 const& tp2 = track.TrajectoryPoint(t);
81  stepsize[t-1] = (tp1 - tp2).Mag();
82  largestStep = std::max(stepsize[t-1], largestStep);
83  }
84 
85  // sort the vector of step sizes and find the median
86  std::sort(stepsize.begin(), stepsize.end());
87 
88  median = 0.;
89  if(size%2 == 0)
90  median = 0.5*(stepsize[size/2] + stepsize[(size/2) - 1]);
91  else
92  median = stepsize[size/2];
93 
94  steps.swap(stepsize);
95 
96  // the median represents the value of the step size
97  // distribution where half of the distribution is
98  // below that value and half is above it. If the
99  // largest step is more than twice the median, then
100  // it is very far outside the distribution and this
101  // track was not well reconstructed
102  return (largestStep < median*fBadStepSizeLimit);
103  }
104 
105  //----------------------------------------------------------------------------
107  {
108  std::vector<float> steps;
109  float median = 0.;
110 
111  return (std::abs(track.Stop().Z() - track.Start().Z()) >= fTrkLengthCut &&
112  -track.Dir().Y() < fMaxCosTheta*track.Dir().Mag() &&
113  this->GoodSteps(track, steps, median) );
114  }
115 
116  //----------------------------------------------------------------------------
118  {
119  return (ri.stop &&
120  ri.goodVtx &&
125  }
126 
127  //----------------------------------------------------------------------------
129  rb::Cluster const& slice)
130  {
131 
132  //Track must have at least two X cells and 2 Y cells:
133  if(track.NXCell() < 2 || track.NYCell() < 2) {
134  return false;
135  }
136 
137  //Check that direction vector is well-defined:
138  const TVector3 dr = track.Dir().Unit();
139  const double dx = std::abs(dr.X());
140  const double dy = std::abs(dr.Y());
141  const double dz = std::abs(dr.Z());
142 
143  //No undefined directions:
144  if(std::isnan(dx) ||
145  std::isnan(dy) ||
146  std::isnan(dz)) return false;
147 
148  //At least one direction must be non-zero:
149  if(dx == 0 &&
150  dy == 0 &&
151  dz == 0) {
152  return false;
153  }
154 
155  // Additional quality cuts. These remove tracks that often have
156  // badly-reconstructed W positions.
157 
158  // Cut very steep tracks two ways. These can be misreconstructed, or
159  // they can simply make it difficult for the W estimation to do a good
160  // job. OK to use the rb::Track::Dir() method here because the initial
161  // direction of the track is what we want to cut on
162  if(std::abs(track.Start().Z() - track.Stop().Z() ) <= 70 ||
163  std::abs(track.Dir().Z()) < 0.2) return false;
164 
165  // Cut out actual reconstruction failures, where the track doesn't use
166  // a significant amount of the slice hits. If either view is not
167  // complete enough discard the whole thing
168  if(track.NXCell() < 0.8*slice.NXCell() ||
169  track.NYCell() < 0.8*slice.NYCell()) return false;
170 
171  // Check that the average number of cells per plane is similar in
172  // each view. If the value is very high in either view, or larger in
173  // one view than the other then you might have some form of FEB flash
174  // or cosmic ray brem that will mess up the calibration
175  std::set<unsigned short> xPlaneList;
176  std::set<unsigned short> yPlaneList;
177  for(size_t c = 0; c < track.NCell(); ++c){
178  if (track.Cell(c)->View() == geo::kX) xPlaneList.insert(track.Cell(c)->Plane());
179  else if(track.Cell(c)->View() == geo::kY) yPlaneList.insert(track.Cell(c)->Plane());
180  }// end loop over cells in track
181  double xCells = track.NXCell();
182  double yCells = track.NYCell();
183  double xPlanes = 1. * xPlaneList.size();
184  double yPlanes = 1. * yPlaneList.size();
185  double cellsPerPlaneX = (xPlanes > 0) ? xCells/xPlanes : 999.;
186  double cellsPerPlaneY = (yPlanes > 0) ? yCells/yPlanes : 999.;
187  if( std::max(cellsPerPlaneX, cellsPerPlaneY) > fMaxCellsPerPlane ) return false;
188 
189  // remove tracks where there is an asymmetry between the number
190  // of planes hit in each view. One class of asymmetry is where
191  // there are extra planes off either end of the track in one view.
192  // Since there is no information in the other view for those planes
193  // you can't be sure you got the length/containment correct
194  double xStart = track.MaxPlane(geo::kX);
195  double yStart = track.MaxPlane(geo::kY);
196  double xEnd = track.MinPlane(geo::kX);
197  double yEnd = track.MinPlane(geo::kY);
198  if( std::abs(xStart - yStart) > fMaxDeltaEndPlane ||
199  std::abs(xEnd - yEnd ) > fMaxDeltaEndPlane ||
200  std::abs(xPlanes - yPlanes) > fMaxPlaneAsymmetry*(xPlanes + yPlanes) ) return false;
201 
202 
203  // Check that the track doesn't take one or two big steps
204  // compared to the rest, that is usually a reconstruction failure.
205  std::vector<float> steps;
206  float median = 0.;
207  if( !this->GoodSteps(track, steps, median) ) return false;
208 
209  // Remove tracks whose vertex position is well away from the edge of the
210  // detector. Muon catcher is not counted as part of the detector proper in
211  // ND so have to check for that as well.
212  if( std::abs(this->DistToEdge(track.Start())) > fMaxVtxDistFromEdge)
213  return false;
214 
215 
216  // Remove tracks whose end position is more than fMaxEndDistFromEdge
217  // outside of the detector The minus sign is because LiveGeometry returns a
218  // negative value if the point is outside of the active detector. A track
219  // can stop anywhere inside the detector and still be OK.
220  if( this->DistToEdge(track.Stop()) < -fMaxEndDistFromEdge)
221  return false;
222 
223  return true;
224  }
225 
226  //----------------------------------------------------------------------------
227  float CosmicTrackSelection::DistToEdge(TVector3 const& point)
228  {
229  return std::min(fLiveGeom->DistToEdgeXY(point), fLiveGeom->DistToEdgeZ(point));
230  }
231 
232 } // namespace trk
T max(const caf::Proxy< T > &a, T b)
size_t NTrajectoryPoints() const
Definition: Track.h:83
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void reconfigure(fhicl::ParameterSet const &pset)
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
unsigned short Plane() const
Definition: CellHit.h:39
double DistToEdgeXY(TVector3 vertex)
geo::View_t View() const
Definition: CellHit.h:41
double fMaxVtxDistFromEdge
maximum distance of the track vertex from the edge of the detector
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
double fMaxCellsPerPlane
to remove shower-y events
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
bool Stopper(TVector3 const &finalPoint)
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
double fTrkLengthCut
Minimum length of track to accept (cm)
Track finder for cosmic rays.
double dy[NP][NC]
double dx[NP][NC]
double DistToEdgeZ(TVector3 vertex)
double dz[NP][NC]
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
double fMaxCosTheta
reject tracks coming from larger cos(theta)
double median(TH1D *h)
Definition: absCal.cxx:524
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
double fMaxDeltaEndPlane
biggest allowed difference in end plane number
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
bool IsCalibTrack(rb::Track const &track, rb::Cluster const &slice)
float Mag() const
double fMaxPlaneAsymmetry
biggest allowed asymmetry in number of planes in each view
bool GoodReco(rb::Track const &track)
double fMaxEndDistFromEdge
maximum distance of the track end from the edge of the detector
float DistToEdge(TVector3 const &point)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
bool CleanSample(CosmicTrackInfo const &ri)
T min(const caf::Proxy< T > &a, T b)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
bool GoodSteps(rb::Track const &track, std::vector< float > &steps, float &median)