LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Private Member Functions | Private Attributes | List of all members
shower::TrackShowerSeparationAlg Class Reference

#include "TrackShowerSeparationAlg.h"

Public Member Functions

 TrackShowerSeparationAlg (fhicl::ParameterSet const &pset)
 
void reconfigure (fhicl::ParameterSet const &pset)
 Read in configurable parameters from provided parameter set. More...
 
std::vector< art::Ptr< recob::Hit > > SelectShowerHits (int event, const std::vector< art::Ptr< recob::Hit > > &hits, const std::vector< art::Ptr< recob::Track > > &tracks, const std::vector< art::Ptr< recob::SpacePoint > > &spacePoints, const art::FindManyP< recob::Hit > &fmht, const art::FindManyP< recob::Track > &fmth, const art::FindManyP< recob::SpacePoint > &fmspt, const art::FindManyP< recob::Track > &fmtsp)
 

Private Member Functions

std::vector< int > InitialTrackLikeSegment (std::map< int, std::unique_ptr< ReconTrack > > &reconTracks)
 
TVector3 Gradient (const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir)
 
TVector3 Gradient (const art::Ptr< recob::Track > &track)
 
TVector3 Gradient (const std::vector< art::Ptr< recob::SpacePoint > > &spacePoints)
 
TVector3 ProjPoint (const TVector3 &point, const TVector3 &direction, const TVector3 &origin=TVector3(0, 0, 0))
 
TVector3 SpacePointPos (const art::Ptr< recob::SpacePoint > &spacePoint)
 Return 3D point of this space point. More...
 
double SpacePointsRMS (const std::vector< art::Ptr< recob::SpacePoint > > &spacePoints)
 

Private Attributes

int fDebug
 
double fConeAngle
 
double fCylinderRadius
 
double fTrackVertexCut
 
double fCylinderCut
 
double fShowerConeCut
 

Detailed Description

Definition at line 169 of file TrackShowerSeparationAlg.h.

Constructor & Destructor Documentation

shower::TrackShowerSeparationAlg::TrackShowerSeparationAlg ( fhicl::ParameterSet const &  pset)

Definition at line 14 of file TrackShowerSeparationAlg.cxx.

References reconfigure().

14  {
15  this->reconfigure(pset);
16 }
shower::TrackShowerSeparationAlg::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Read in configurable parameters from provided parameter set.
Definition: TrackShowerSeparationAlg.cxx:18

Member Function Documentation

TVector3 shower::TrackShowerSeparationAlg::Gradient ( const std::vector< TVector3 > &  points,
const std::unique_ptr< TVector3 > &  dir 
)
private

Definition at line 346 of file TrackShowerSeparationAlg.cxx.

Referenced by Gradient(), SelectShowerHits(), and SpacePointsRMS().

346  {
347 
348  int nhits = 0;
349  double sumx=0., sumy=0., sumz=0., sumx2=0., sumy2=0., sumxy=0., sumxz=0., sumyz=0.;
350  for (std::vector<TVector3>::const_iterator pointIt = points.begin(); pointIt != points.end(); ++pointIt) {
351  ++nhits;
352  sumx += pointIt->X();
353  sumy += pointIt->Y();
354  sumz += pointIt->Z();
355  sumx2 += pointIt->X() * pointIt->X();
356  sumy2 += pointIt->Y() * pointIt->Y();
357  sumxy += pointIt->X() * pointIt->Y();
358  sumxz += pointIt->X() * pointIt->Z();
359  sumyz += pointIt->Y() * pointIt->Z();
360  }
361 
362  double dydx = (nhits * sumxy - sumx * sumy) / (nhits * sumx2 - sumx * sumx);
363  double yint = (sumy * sumx2 - sumx * sumxy) / (nhits * sumx2 - sumx * sumx);
364  double dzdx = (nhits * sumxz - sumx * sumz) / (nhits * sumx2 - sumx * sumx);
365  double zint = (sumz * sumx2 - sumx * sumxz) / (nhits * sumx2 - sumx * sumx);
366  TVector2 directionXY = TVector2(1,dydx).Unit(), directionXZ = TVector2(1,dzdx).Unit();
367  TVector3 direction = TVector3(1,dydx,dzdx).Unit();
368  TVector3 intercept = TVector3(0,yint,zint);
369 
370  // Make sure the best fit direction is pointing correctly
371  if (dir and TMath::Abs(direction.Angle(*dir)) > TMath::Pi() / 2.)
372  direction *= -1;
373 
374  return direction;
375 
376 }
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:25
TVector3 shower::TrackShowerSeparationAlg::Gradient ( const art::Ptr< recob::Track > &  track)
private

Definition at line 378 of file TrackShowerSeparationAlg.cxx.

References dir, Gradient(), recob::Track::LocationAtPoint(), recob::Track::NumberTrajectoryPoints(), and recob::Track::VertexDirection().

378  {
379 
380  std::vector<TVector3> points;
381  std::unique_ptr<TVector3> dir;
382 
383  for (unsigned int traj = 0; traj < track->NumberTrajectoryPoints(); ++traj)
384  points.push_back(track->LocationAtPoint<TVector3>(traj));
385  dir = std::make_unique<TVector3>(track->VertexDirection<TVector3>());
386 
387  return Gradient(points, dir);
388 
389 }
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Access to track position at different points.
Definition: Track.h:129
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:105
recob::Track::VertexDirection
Vector_t VertexDirection() const
Access to track direction at different points.
Definition: Track.h:135
shower::TrackShowerSeparationAlg::Gradient
TVector3 Gradient(const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir)
Definition: TrackShowerSeparationAlg.cxx:346
dir
TDirectory * dir
Definition: macro.C:5
TVector3 shower::TrackShowerSeparationAlg::Gradient ( const std::vector< art::Ptr< recob::SpacePoint > > &  spacePoints)
private

Definition at line 391 of file TrackShowerSeparationAlg.cxx.

References dir, Gradient(), SpacePointPos(), and lar::dump::vector().

391  {
392 
393  std::vector<TVector3> points;
394  std::unique_ptr<TVector3> dir;
395 
396  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt)
397  points.push_back(SpacePointPos(*spacePointIt));
398 
399  return Gradient(points, dir);
400 
401 }
shower::TrackShowerSeparationAlg::SpacePointPos
TVector3 SpacePointPos(const art::Ptr< recob::SpacePoint > &spacePoint)
Return 3D point of this space point.
Definition: TrackShowerSeparationAlg.cxx:407
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
shower::TrackShowerSeparationAlg::Gradient
TVector3 Gradient(const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir)
Definition: TrackShowerSeparationAlg.cxx:346
dir
TDirectory * dir
Definition: macro.C:5
std::vector< int > shower::TrackShowerSeparationAlg::InitialTrackLikeSegment ( std::map< int, std::unique_ptr< ReconTrack > > &  reconTracks)
private

Definition at line 280 of file TrackShowerSeparationAlg.cxx.

References fConeAngle.

280  {
281 
282  // Consider the cones for each track
283  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
284 
285  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator otherTrackIt = reconTracks.begin(); otherTrackIt != reconTracks.end(); ++otherTrackIt) {
286 
287  if (trackIt->first == otherTrackIt->first)
288  continue;
289  if ((otherTrackIt->second->Vertex() - trackIt->second->Vertex()).Angle(trackIt->second->VertexDirection()) < fConeAngle * TMath::Pi() / 180 or
290  (otherTrackIt->second->End() - trackIt->second->Vertex()).Angle(trackIt->second->VertexDirection()) < fConeAngle * TMath::Pi() / 180) {
291  trackIt->second->AddForwardTrack(otherTrackIt->first);
292  otherTrackIt->second->AddShowerTrack(trackIt->first);
293  }
294  if ((otherTrackIt->second->Vertex() - trackIt->second->Vertex()).Angle(-1*trackIt->second->VertexDirection()) < fConeAngle * TMath::Pi() / 180 or
295  (otherTrackIt->second->End() - trackIt->second->Vertex()).Angle(-1*trackIt->second->VertexDirection()) < fConeAngle * TMath::Pi() / 180)
296  trackIt->second->AddBackwardTrack(otherTrackIt->first);
297 
298  }
299 
300  }
301 
302  // Determine if any of these tracks are actually shower tracks
303  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
304 
305  if (!trackIt->second->ShowerTrackCandidate())
306  continue;
307 
308  std::cout << "Track " << trackIt->first << " is a candidate, with shower cone tracks:" << std::endl;
309  const std::vector<int>& showerConeTracks = trackIt->second->ForwardConeTracks();
310  for (std::vector<int>::const_iterator showerConeTrackIt = showerConeTracks.begin(); showerConeTrackIt != showerConeTracks.end(); ++showerConeTrackIt)
311  std::cout << " " << *showerConeTrackIt << std::endl;
312 
313  bool isBestCandidate = true;
314  const std::vector<int>& showerTracks = trackIt->second->ShowerTracks();
315  for (std::vector<int>::const_iterator showerTrackIt = showerTracks.begin(); showerTrackIt != showerTracks.end(); ++showerTrackIt) {
316  if (!reconTracks[*showerTrackIt]->ShowerTrackCandidate())
317  continue;
318  if (std::find(showerConeTracks.begin(), showerConeTracks.end(), *showerTrackIt) == showerConeTracks.end())
319  continue;
320  if (reconTracks[*showerTrackIt]->IsShowerTrack())
321  continue;
322  if (trackIt->second->TrackConeSize() < reconTracks[*showerTrackIt]->TrackConeSize())
323  isBestCandidate = false;
324  }
325 
326  if (isBestCandidate)
327  trackIt->second->MakeShowerTrack();
328 
329  }
330 
331  // Determine which tracks are shower cones
332  std::vector<int> showerTracks;
333  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
334  if (trackIt->second->IsShowerTrack()) {
335  showerTracks.push_back(trackIt->first);
336  const std::vector<int>& coneTracks = trackIt->second->ForwardConeTracks();
337  for (std::vector<int>::const_iterator coneTrackIt = coneTracks.begin(); coneTrackIt != coneTracks.end(); ++coneTrackIt)
338  reconTracks[*coneTrackIt]->MakeShowerCone();
339  }
340  }
341 
342  return showerTracks;
343 
344 }
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:25
TVector3 shower::TrackShowerSeparationAlg::ProjPoint ( const TVector3 &  point,
const TVector3 &  direction,
const TVector3 &  origin = TVector3(0,0,0) 
)
private

Projects the 3D point given along the line defined by the specified direction Coordinate system is assumed to be centred at the origin unless a difference point is specified

Definition at line 403 of file TrackShowerSeparationAlg.cxx.

References geo::origin().

Referenced by SelectShowerHits(), and SpacePointsRMS().

403  {
404  return (point-origin).Dot(direction) * direction + origin;
405 }
geo::origin
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:230
void shower::TrackShowerSeparationAlg::reconfigure ( fhicl::ParameterSet const &  pset)

Read in configurable parameters from provided parameter set.

Definition at line 18 of file TrackShowerSeparationAlg.cxx.

References fConeAngle, fCylinderCut, fCylinderRadius, fDebug, fShowerConeCut, fTrackVertexCut, and fhicl::ParameterSet::get().

Referenced by TrackShowerSeparationAlg().

18  {
19  fConeAngle = pset.get<double>("ConeAngle");
20  fCylinderRadius = pset.get<double>("CylinderRadius");
21  fTrackVertexCut = pset.get<double>("TrackVertexCut");
22  fCylinderCut = pset.get<double>("CylinderCut");
23  fShowerConeCut = pset.get<double>("ShowerConeCut");
24 
25  fDebug = pset.get<int>("Debug",0);
26 }
std::vector< art::Ptr< recob::Hit > > shower::TrackShowerSeparationAlg::SelectShowerHits ( int  event,
const std::vector< art::Ptr< recob::Hit > > &  hits,
const std::vector< art::Ptr< recob::Track > > &  tracks,
const std::vector< art::Ptr< recob::SpacePoint > > &  spacePoints,
const art::FindManyP< recob::Hit > &  fmht,
const art::FindManyP< recob::Track > &  fmth,
const art::FindManyP< recob::SpacePoint > &  fmspt,
const art::FindManyP< recob::Track > &  fmtsp 
)

Definition at line 28 of file TrackShowerSeparationAlg.cxx.

References fConeAngle, fCylinderCut, fCylinderRadius, fDebug, for(), fTrackVertexCut, Gradient(), hits(), proj, ProjPoint(), SpacePointPos(), track, and lar::dump::vector().

Referenced by cluster::BlurredClustering::produce().

35  {
36 
37  // Ok, here we are again
38  // Playing the game in which no one wins, but everyone loses
39  // Trying to separate tracks from showers
40  // Ode to track/shower separation
41  // M Wallbank, Oct 2016
42 
43  // Showers are comprised of two topologically different parts:
44  // -- the 'shower track' (the initial part of the shower)
45  // -- the 'shower cone' (the part of the shower after the initial track when showering happens)
46 
47  // -
48  // --
49  // - -- --
50  // -- --- ---
51  // -- ---- --
52  // ----------- ----- - -- ---
53  // --- ---
54  // -- ---
55  // {=========}{==============================}
56  // shower track shower cone
57 
58  std::map<int,std::unique_ptr<ReconTrack> > reconTracks;
59  for (std::vector<art::Ptr<recob::Track> >::const_iterator trackIt = tracks.begin(); trackIt != tracks.end(); ++trackIt) {
60  std::unique_ptr<ReconTrack> track = std::make_unique<ReconTrack>(trackIt->key());
61  track->SetVertex((*trackIt)->Vertex<TVector3>());
62  track->SetEnd((*trackIt)->End<TVector3>());
63  track->SetVertexDir((*trackIt)->VertexDirection<TVector3>());
64  track->SetLength((*trackIt)->Length());
65  track->SetDirection(Gradient(*trackIt));
66  track->SetHits(fmht.at(trackIt->key()));
67  track->SetSpacePoints(fmspt.at(trackIt->key()));
68  const std::vector<art::Ptr<recob::SpacePoint> > spsss = fmspt.at(trackIt->key());
69  // std::cout << "Track " << trackIt->key() << " has " << spsss.size() << " space points and " << (*trackIt)->NumberTrajectoryPoints() << " traj points" << std::endl;
70  // if (trackIt->key() == 5)
71  // for (unsigned int i = 0; i < (*trackIt)->NumberTrajectoryPoints(); ++i)
72  // std::cout << "Traj point " << i << " has position (" << (*trackIt)->LocationAtPoint(i).X() << ", " << (*trackIt)->LocationAtPoint(i).Y() << ", " << (*trackIt)->LocationAtPoint(i).Z() << ")" << std::endl;
73  reconTracks[trackIt->key()] = std::move(track);
74  }
75 
76  // std::vector<int> showerLikeTracks, trackLikeTracks;
77  // std::vector<int> showerTracks = InitialTrackLikeSegment(reconTracks);
78 
79  // Consider the space point cylinder situation
80  double avCylinderSpacePoints = 0;
81  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
82  // Get the 3D properties of the track
83  TVector3 point = trackIt->second->Vertex();
84  TVector3 direction = trackIt->second->Direction();
85  // if (trackIt->second->Vertex().X() > 250 and trackIt->second->Vertex().X() < 252 and
86  // trackIt->second->Vertex().Y() > -440 and trackIt->second->Vertex().Y() < -430 and
87  // trackIt->second->Vertex().Z() > 1080 and trackIt->second->Vertex().Z() < 1090)
88  // std::cout << "Track " << trackIt->first << " begins at the most upstream vertex" << std::endl;
89  // if (trackIt->second->Vertex().X() > 254 and trackIt->second->Vertex().X() < 258 and
90  // trackIt->second->Vertex().Y() > -440 and trackIt->second->Vertex().Y() < -420 and
91  // trackIt->second->Vertex().Z() > 1115 and trackIt->second->Vertex().Z() < 1130)
92  // std::cout << "Track " << trackIt->first << " begins at the supposed vertex" << std::endl;
93  // if (trackIt->second->End().X() > 254 and trackIt->second->End().X() < 258 and
94  // trackIt->second->End().Y() > -440 and trackIt->second->End().Y() < -420 and
95  // trackIt->second->End().Z() > 1115 and trackIt->second->End().Z() < 1130)
96  // std::cout << "Track " << trackIt->first << " ends at the supposed vertex" << std::endl;
97  // std::cout << "Track " << trackIt->first << " has vertex (" << trackIt->second->Vertex().X() << ", " << trackIt->second->Vertex().Y() << ", " << trackIt->second->Vertex().Z() << ") and end (" << trackIt->second->End().X() << ", " << trackIt->second->End().Y() << ", " << trackIt->second->End().Z() << "), with vertex direction (" << trackIt->second->VertexDirection().X() << ", " << trackIt->second->VertexDirection().Y() << ", " << trackIt->second->VertexDirection().Z() << ")" << std::endl;
98  // Count space points in the volume around the track
99  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
100  const std::vector<art::Ptr<recob::Track> > spTracks = fmtsp.at(spacePointIt->key());
101  if (find_if(spTracks.begin(), spTracks.end(), [&trackIt](const art::Ptr<recob::Track>& t){ return (int)t.key() == trackIt->first; }) != spTracks.end())
102  continue;
103  // Get the properties of this space point
104  TVector3 pos = SpacePointPos(*spacePointIt);
105  TVector3 proj = ProjPoint(pos, direction, point);
106  if ((pos-proj).Mag() < fCylinderRadius)
107  trackIt->second->AddCylinderSpacePoint(spacePointIt->key());
108  // if ((pos-proj).Mag() < fCylinderRadius and
109  // (pos-point)*direction > 0 and
110  // (pos-point)*direction < trackIt->second->Length())
111  // std::cout << "Space point " << spacePointIt->key() << " (" << pos.X() << ", " << pos.Y() << ", " << pos.Z() << ") in cylinder around track " << trackIt->first << " (assocatied with track " << spTracks.at(0).key() << "); point is (" << point.X() << ", " << point.Y() << ", " << point.Z() << "), proj end (" << ((trackIt->second->Length()*trackIt->second->VertexDirection())+point).X() << ", " << ((trackIt->second->Length()*trackIt->second->VertexDirection())+point).Y() << ", " << ((trackIt->second->Length()*trackIt->second->VertexDirection())+point).Z() << ")" << std::endl;
112  }
113  avCylinderSpacePoints += trackIt->second->CylinderSpacePointRatio();
114  }
115  avCylinderSpacePoints /= (double)reconTracks.size();
116 
117  if (fDebug > 1) {
118  std::cout << std::endl << "Cylinder space point ratios:" << std::endl;
119  for (std::map<int,std::unique_ptr<ReconTrack> >::const_iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt)
120  std::cout << " Track " << trackIt->first << " has cylinder space point ratio " << trackIt->second->CylinderSpacePointRatio() << " (event average " << avCylinderSpacePoints << ")" << std::endl;
121  }
122 
123  // Identify tracks
124  if (fDebug > 0)
125  std::cout << std::endl << "Identifying tracks:" << std::endl;
126  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt)
127  if (trackIt->second->CylinderSpacePointRatio() / avCylinderSpacePoints < fCylinderCut) {
128  if (fDebug > 0)
129  std::cout << " Making track " << trackIt->first << " a track (Type I)" << std::endl;
130  trackIt->second->MakeTrack();
131  }
132 
133  // for (std::map<int,std::unique_ptr<ReconTrack> >::const_iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
134  // if (!trackIt->second->IsTrack())
135  // continue;
136  // std::cout << "Track " << trackIt->first << " (tagged as track) is this close to other tracks:" << std::endl;
137  // for (std::map<int,std::unique_ptr<ReconTrack> >::const_iterator otherTrackIt = reconTracks.begin(); otherTrackIt != reconTracks.end(); ++otherTrackIt) {
138  // if (trackIt->first == otherTrackIt->first)
139  // continue;
140  // std::cout << " Other track " << otherTrackIt->first << " from track " << trackIt->first << " vertex: " << (otherTrackIt->second->Vertex()-trackIt->second->Vertex()).Mag() << " (vertex) and " << (otherTrackIt->second->End()-trackIt->second->Vertex()).Mag() << " (end); end: " << (otherTrackIt->second->Vertex()-trackIt->second->End()).Mag() << " (vertex) and " << (otherTrackIt->second->End()-trackIt->second->End()).Mag() << " (end)" << std::endl;
141  // }
142  // }
143 
144  // Identify further tracks
145  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
146  if (trackIt->second->IsTrack())
147  continue;
148  for (std::map<int,std::unique_ptr<ReconTrack> >::const_iterator otherTrackIt = reconTracks.begin(); otherTrackIt != reconTracks.end(); ++otherTrackIt) {
149  if (trackIt->first == otherTrackIt->first or !otherTrackIt->second->IsTrack())
150  continue;
151  if ((trackIt->second->Vertex() - otherTrackIt->second->Vertex()).Mag() < fTrackVertexCut or
152  (trackIt->second->Vertex() - otherTrackIt->second->End()).Mag() < fTrackVertexCut or
153  (trackIt->second->End() - otherTrackIt->second->Vertex()).Mag() < fTrackVertexCut or
154  (trackIt->second->End() - otherTrackIt->second->End()).Mag() < fTrackVertexCut) {
155  if (fDebug > 0)
156  std::cout << " Making track " << trackIt->first << " a track (Type II)" << std::endl;
157  trackIt->second->MakeTrack();
158  }
159  }
160  }
161 
162  // Consider removing false tracks by looking at their closest approach to any other track
163 
164  // Consider the space point cone situation
165  std::vector<art::Ptr<recob::SpacePoint> > showerSpacePoints;
166  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
167  bool showerSpacePoint = true;
168  const std::vector<art::Ptr<recob::Track> > spacePointTracks = fmtsp.at(spacePointIt->key());
169  for (std::vector<art::Ptr<recob::Track> >::const_iterator trackIt = spacePointTracks.begin(); trackIt != spacePointTracks.end(); ++trackIt)
170  if (reconTracks[trackIt->key()]->IsTrack())
171  showerSpacePoint = false;
172  if (showerSpacePoint)
173  showerSpacePoints.push_back(*spacePointIt);
174  }
175 
176  // Identify tracks which slipped through and shower tracks
177  // For the moment, until the track tagging gets better at least, don't try to identify tracks from this
178  double avConeSize = 0;
179  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
180  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = showerSpacePoints.begin(); spacePointIt != showerSpacePoints.end(); ++spacePointIt) {
181  bool associatedSpacePoint = false;
182  const std::vector<art::Ptr<recob::Track> > spTracks = fmtsp.at(spacePointIt->key());
183  for (std::vector<art::Ptr<recob::Track> >::const_iterator spTrackIt = spTracks.begin(); spTrackIt != spTracks.end(); ++spTrackIt)
184  if ((int)spTrackIt->key() == trackIt->first)
185  associatedSpacePoint = true;
186  if (associatedSpacePoint)
187  continue;
188  if ((SpacePointPos(*spacePointIt) - trackIt->second->Vertex()).Angle(trackIt->second->Direction()) < fConeAngle * TMath::Pi() / 180) {
189  trackIt->second->AddForwardSpacePoint(spacePointIt->key());
190  trackIt->second->AddForwardTrack(spTracks.at(0).key());
191  }
192  if ((SpacePointPos(*spacePointIt) - trackIt->second->Vertex()).Angle(-1*trackIt->second->Direction()) < fConeAngle * TMath::Pi() / 180) {
193  trackIt->second->AddBackwardSpacePoint(spacePointIt->key());
194  trackIt->second->AddBackwardTrack(spTracks.at(0).key());
195  }
196  }
197  avConeSize += trackIt->second->ConeSize();
198  }
199  avConeSize /= (double)reconTracks.size();
200  if (fDebug > 0)
201  std::cout << std::endl << "Identifying showers:" << std::endl;
202  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
203  // if (trackIt->second->ForwardSpacePoints() == 0)
204  // trackIt->second->MakeTrack();
205  // double distanceFromAverage = (trackIt->second->ConeSize() - avConeSize) / TMath::Abs(avConeSize);
206  // if (fDebug > 1)
207  // std::cout << " Track " << trackIt->first << " has cone size " << trackIt->second->ConeSize() << " (forward " << trackIt->second->ForwardSpacePoints() << ") and tracks " << trackIt->second->TrackConeSize() << " (average " << avConeSize << ", distance from average " << distanceFromAverage << ")" << std::endl;
208  // if (distanceFromAverage > fShowerConeCut) {
209  // trackIt->second->MakeShower();
210  // if (fDebug > 0)
211  // std::cout << " Making track " << trackIt->first << " a shower (Type I)" << std::endl;
212  // }
213  if (fDebug > 1)
214  std::cout << " Track " << trackIt->first << " has space point cone size " << trackIt->second->ConeSize() << " and track cone size " << trackIt->second->TrackConeSize() << std::endl;
215  if (TMath::Abs(trackIt->second->ConeSize()) > 30 and TMath::Abs(trackIt->second->TrackConeSize()) > 3) {
216  trackIt->second->MakeShower();
217  if (fDebug > 0)
218  std::cout << " Making track " << trackIt->first << " a shower (Type I)" << std::endl;
219  if (trackIt->second->ConeSize() < 0)
220  trackIt->second->FlipTrack();
221  }
222  }
223 
224  // Look for shower cones
225  std::cout << std::endl << " Shower cones:" << std::endl;
226  for (std::map<int,std::unique_ptr<ReconTrack> >::const_iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
227  if (trackIt->second->IsShower()) {
228  if (fDebug > 1)
229  std::cout << " Track " << trackIt->first << std::endl;
230  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator otherTrackIt = reconTracks.begin(); otherTrackIt != reconTracks.end(); ++otherTrackIt) {
231  if (trackIt->first == otherTrackIt->first or !otherTrackIt->second->IsUndetermined())
232  continue;
233  if ((otherTrackIt->second->Vertex()-trackIt->second->Vertex()).Angle(trackIt->second->Direction()) < fConeAngle * TMath::Pi() / 180 or
234  (otherTrackIt->second->End()-trackIt->second->Vertex()).Angle(trackIt->second->Direction()) < fConeAngle * TMath::Pi() / 180) {
235  std::cout << " " << otherTrackIt->first << std::endl;
236  otherTrackIt->second->MakeShower();
237  if (fDebug > 0)
238  std::cout << " Making track " << otherTrackIt->first << " a shower (Type II)" << std::endl;
239  }
240  }
241  }
242  }
243 
244  // Look at remaining undetermined tracks
245  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt) {
246  if (trackIt->second->IsUndetermined()) {
247  // std::cout << "Remaining undetermined track " << trackIt->first << std::endl;
248  // std::cout << "Length is " << trackIt->second->Length() << " and rms of space points is " << SpacePointsRMS(trackIt->second->SpacePoints()) << std::endl;
249  trackIt->second->MakeShower();
250  }
251  }
252 
253  std::cout << std::endl << "Event " << event << " track shower separation:" << std::endl;
254  std::cout << "Shower initial tracks are:" << std::endl;
255  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt)
256  if (trackIt->second->IsShowerTrack())
257  std::cout << " " << trackIt->first << std::endl;
258 
259  std::cout << "Track tracks are:" << std::endl;
260  for (std::map<int,std::unique_ptr<ReconTrack> >::iterator trackIt = reconTracks.begin(); trackIt != reconTracks.end(); ++trackIt)
261  if (trackIt->second->IsTrack())
262  std::cout << " " << trackIt->first << std::endl;
263 
264  // Shower hits -- select all hits which aren't associated with a determined track
265  std::vector<art::Ptr<recob::Hit> > showerHits;
266  for (std::vector<art::Ptr<recob::Hit> >::const_iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt) {
267  bool showerHit = true;
268  const std::vector<art::Ptr<recob::Track> > hitTracks = fmth.at(hitIt->key());
269  for (std::vector<art::Ptr<recob::Track> >::const_iterator hitTrackIt = hitTracks.begin(); hitTrackIt != hitTracks.end(); ++hitTrackIt)
270  if (reconTracks[hitTrackIt->key()]->IsTrack())
271  showerHit = false;
272  if (showerHit)
273  showerHits.push_back(*hitIt);
274  }
275 
276  return showerHits;
277 
278 }
shower::TrackShowerSeparationAlg::ProjPoint
TVector3 ProjPoint(const TVector3 &point, const TVector3 &direction, const TVector3 &origin=TVector3(0, 0, 0))
Definition: TrackShowerSeparationAlg.cxx:403
shower::TrackShowerSeparationAlg::SpacePointPos
TVector3 SpacePointPos(const art::Ptr< recob::SpacePoint > &spacePoint)
Return 3D point of this space point.
Definition: TrackShowerSeparationAlg.cxx:407
for
for(int i=0;i< 401;i++)
Definition: TG43_relative_dose.C:23
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
shower::TrackShowerSeparationAlg::Gradient
TVector3 Gradient(const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir)
Definition: TrackShowerSeparationAlg.cxx:346
proj
Float_t proj
Definition: plot.C:34
track
Float_t track
Definition: plot.C:34
TVector3 shower::TrackShowerSeparationAlg::SpacePointPos ( const art::Ptr< recob::SpacePoint > &  spacePoint)
private

Return 3D point of this space point.

Definition at line 407 of file TrackShowerSeparationAlg.cxx.

References recob::SpacePoint::XYZ().

Referenced by Gradient(), SelectShowerHits(), and SpacePointsRMS().

407  {
408  const double* xyz = spacePoint->XYZ();
409  return TVector3(xyz[0], xyz[1], xyz[2]);
410 }
recob::SpacePoint::XYZ
const Double32_t * XYZ() const
Definition: SpacePoint.h:65
double shower::TrackShowerSeparationAlg::SpacePointsRMS ( const std::vector< art::Ptr< recob::SpacePoint > > &  spacePoints)
private

Definition at line 412 of file TrackShowerSeparationAlg.cxx.

References Gradient(), proj, ProjPoint(), SpacePointPos(), and lar::dump::vector().

412  {
413 
414  TVector3 point = SpacePointPos(spacePoints.at(0));
415  TVector3 direction = Gradient(spacePoints);
416 
417  std::vector<double> distances;
418  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
419  TVector3 pos = SpacePointPos(*spacePointIt);
420  TVector3 proj = ProjPoint(pos, direction, point);
421  distances.push_back((pos-proj).Mag());
422  }
423 
424  double rms = TMath::RMS(distances.begin(), distances.end());
425 
426  return rms;
427 
428 }
shower::TrackShowerSeparationAlg::ProjPoint
TVector3 ProjPoint(const TVector3 &point, const TVector3 &direction, const TVector3 &origin=TVector3(0, 0, 0))
Definition: TrackShowerSeparationAlg.cxx:403
shower::TrackShowerSeparationAlg::SpacePointPos
TVector3 SpacePointPos(const art::Ptr< recob::SpacePoint > &spacePoint)
Return 3D point of this space point.
Definition: TrackShowerSeparationAlg.cxx:407
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
shower::TrackShowerSeparationAlg::Gradient
TVector3 Gradient(const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir)
Definition: TrackShowerSeparationAlg.cxx:346
proj
Float_t proj
Definition: plot.C:34

Member Data Documentation

double shower::TrackShowerSeparationAlg::fConeAngle
private

Definition at line 214 of file TrackShowerSeparationAlg.h.

Referenced by InitialTrackLikeSegment(), reconfigure(), and SelectShowerHits().

double shower::TrackShowerSeparationAlg::fCylinderCut
private

Definition at line 219 of file TrackShowerSeparationAlg.h.

Referenced by reconfigure(), and SelectShowerHits().

double shower::TrackShowerSeparationAlg::fCylinderRadius
private

Definition at line 215 of file TrackShowerSeparationAlg.h.

Referenced by reconfigure(), and SelectShowerHits().

int shower::TrackShowerSeparationAlg::fDebug
private

Definition at line 211 of file TrackShowerSeparationAlg.h.

Referenced by reconfigure(), and SelectShowerHits().

double shower::TrackShowerSeparationAlg::fShowerConeCut
private

Definition at line 220 of file TrackShowerSeparationAlg.h.

Referenced by reconfigure().

double shower::TrackShowerSeparationAlg::fTrackVertexCut
private

Definition at line 218 of file TrackShowerSeparationAlg.h.

Referenced by reconfigure(), and SelectShowerHits().

The documentation for this class was generated from the following files: