All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) const
 

Private Member Functions

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

Private Attributes

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

Detailed Description

Definition at line 160 of file TrackShowerSeparationAlg.h.

Constructor & Destructor Documentation

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

Definition at line 18 of file TrackShowerSeparationAlg.cxx.

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

Member Function Documentation

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

Definition at line 350 of file TrackShowerSeparationAlg.cxx.

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

Definition at line 382 of file TrackShowerSeparationAlg.cxx.

382  {
383 
384  std::vector<TVector3> points;
385  std::unique_ptr<TVector3> dir;
386 
387  for (unsigned int traj = 0; traj < track->NumberTrajectoryPoints(); ++traj)
388  points.push_back(track->LocationAtPoint<TVector3>(traj));
389  dir = std::make_unique<TVector3>(track->VertexDirection<TVector3>());
390 
391  return Gradient(points, dir);
392 
393 }
track
process_name use argoneut_mc_hitfinder track
Definition: standard_reco.fcl:36
dropbox.dir
tuple dir
Definition: dropbox.py:28
shower::TrackShowerSeparationAlg::Gradient
TVector3 Gradient(const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir) const
Definition: TrackShowerSeparationAlg.cxx:350
TVector3 shower::TrackShowerSeparationAlg::Gradient ( const std::vector< art::Ptr< recob::SpacePoint > > &  spacePoints) const
private

Definition at line 395 of file TrackShowerSeparationAlg.cxx.

395  {
396 
397  std::vector<TVector3> points;
398  std::unique_ptr<TVector3> dir;
399 
400  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt)
401  points.push_back(SpacePointPos(*spacePointIt));
402 
403  return Gradient(points, dir);
404 
405 }
shower::TrackShowerSeparationAlg::SpacePointPos
TVector3 SpacePointPos(const art::Ptr< recob::SpacePoint > &spacePoint) const
Return 3D point of this space point.
Definition: TrackShowerSeparationAlg.cxx:411
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
dropbox.dir
tuple dir
Definition: dropbox.py:28
shower::TrackShowerSeparationAlg::Gradient
TVector3 Gradient(const std::vector< TVector3 > &points, const std::unique_ptr< TVector3 > &dir) const
Definition: TrackShowerSeparationAlg.cxx:350
std::vector< int > shower::TrackShowerSeparationAlg::InitialTrackLikeSegment ( std::map< int, std::unique_ptr< ReconTrack > > &  reconTracks) const
private

Definition at line 284 of file TrackShowerSeparationAlg.cxx.

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

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

Read in configurable parameters from provided parameter set.

Definition at line 22 of file TrackShowerSeparationAlg.cxx.

22  {
23  fConeAngle = pset.get<double>("ConeAngle");
24  fCylinderRadius = pset.get<double>("CylinderRadius");
25  fTrackVertexCut = pset.get<double>("TrackVertexCut");
26  fCylinderCut = pset.get<double>("CylinderCut");
27  fShowerConeCut = pset.get<double>("ShowerConeCut");
28 
29  fDebug = pset.get<int>("Debug",0);
30 }
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 
) const

Definition at line 32 of file TrackShowerSeparationAlg.cxx.

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

Return 3D point of this space point.

Definition at line 411 of file TrackShowerSeparationAlg.cxx.

411  {
412  const double* xyz = spacePoint->XYZ();
413  return TVector3(xyz[0], xyz[1], xyz[2]);
414 }
double shower::TrackShowerSeparationAlg::SpacePointsRMS ( const std::vector< art::Ptr< recob::SpacePoint > > &  spacePoints) const
private

Definition at line 416 of file TrackShowerSeparationAlg.cxx.

416  {
417 
418  TVector3 point = SpacePointPos(spacePoints.at(0));
419  TVector3 direction = Gradient(spacePoints);
420 
421  std::vector<double> distances;
422  for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
423  TVector3 pos = SpacePointPos(*spacePointIt);
424  TVector3 proj = ProjPoint(pos, direction, point);
425  distances.push_back((pos-proj).Mag());
426  }
427 
428  double rms = TMath::RMS(distances.begin(), distances.end());
429 
430  return rms;
431 
432 }
shower::TrackShowerSeparationAlg::ProjPoint
TVector3 ProjPoint(const TVector3 &point, const TVector3 &direction, const TVector3 &origin=TVector3(0, 0, 0)) const
Definition: TrackShowerSeparationAlg.cxx:407
shower::TrackShowerSeparationAlg::SpacePointPos
TVector3 SpacePointPos(const art::Ptr< recob::SpacePoint > &spacePoint) const
Return 3D point of this space point.
Definition: TrackShowerSeparationAlg.cxx:411
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) const
Definition: TrackShowerSeparationAlg.cxx:350

Member Data Documentation

double shower::TrackShowerSeparationAlg::fConeAngle
private

Definition at line 205 of file TrackShowerSeparationAlg.h.

double shower::TrackShowerSeparationAlg::fCylinderCut
private

Definition at line 210 of file TrackShowerSeparationAlg.h.

double shower::TrackShowerSeparationAlg::fCylinderRadius
private

Definition at line 206 of file TrackShowerSeparationAlg.h.

int shower::TrackShowerSeparationAlg::fDebug
private

Definition at line 202 of file TrackShowerSeparationAlg.h.

double shower::TrackShowerSeparationAlg::fShowerConeCut
private

Definition at line 211 of file TrackShowerSeparationAlg.h.

double shower::TrackShowerSeparationAlg::fTrackVertexCut
private

Definition at line 209 of file TrackShowerSeparationAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5