All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackShowerSeparationAlg.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: TrackShowerSeparationAlg
3 // File: TrackShowerSeparationAlg.cxx
4 // Author: Mike Wallbank (m.wallbank@sheffield.ac.uk), November 2015
5 //
6 // Track/shower separation class.
7 // Provides methods for removing hits associated with track-like
8 // objects.
9 // To be run after track reconstruction, before shower reconstruction.
10 ////////////////////////////////////////////////////////////////////////
11 
12 #include "fhiclcpp/ParameterSet.h"
13 #include "larreco/RecoAlg/TrackShowerSeparationAlg.h"
14 
15 #include "TMathBase.h"
16 #include "TVector2.h"
17 
18 shower::TrackShowerSeparationAlg::TrackShowerSeparationAlg(fhicl::ParameterSet const& pset) {
19  this->reconfigure(pset);
20 }
21 
22 void shower::TrackShowerSeparationAlg::reconfigure(fhicl::ParameterSet const& pset) {
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 }
31 
32 std::vector<art::Ptr<recob::Hit> > shower::TrackShowerSeparationAlg::SelectShowerHits(int event,
33  const std::vector<art::Ptr<recob::Hit> >& hits,
34  const std::vector<art::Ptr<recob::Track> >& tracks,
35  const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints,
36  const art::FindManyP<recob::Hit>& fmht,
37  const art::FindManyP<recob::Track>& fmth,
38  const art::FindManyP<recob::SpacePoint>& fmspt,
39  const art::FindManyP<recob::Track>& fmtsp) const {
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 }
283 
284 std::vector<int> shower::TrackShowerSeparationAlg::InitialTrackLikeSegment(std::map<int,std::unique_ptr<ReconTrack> >& reconTracks) const {
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 }
349 
350 TVector3 shower::TrackShowerSeparationAlg::Gradient(const std::vector<TVector3>& points, const std::unique_ptr<TVector3>& dir) const {
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 }
381 
382 TVector3 shower::TrackShowerSeparationAlg::Gradient(const art::Ptr<recob::Track>& track) const {
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 }
394 
395 TVector3 shower::TrackShowerSeparationAlg::Gradient(const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints) const {
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 }
406 
407 TVector3 shower::TrackShowerSeparationAlg::ProjPoint(const TVector3& point, const TVector3& direction, const TVector3& origin) const {
408  return (point-origin).Dot(direction) * direction + origin;
409 }
410 
411 TVector3 shower::TrackShowerSeparationAlg::SpacePointPos(const art::Ptr<recob::SpacePoint>& spacePoint) const {
412  const double* xyz = spacePoint->XYZ();
413  return TVector3(xyz[0], xyz[1], xyz[2]);
414 }
415 
416 double shower::TrackShowerSeparationAlg::SpacePointsRMS(const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints) const {
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 }
433 
434 
435 
436 
437 
438 // // --------------------------- OLD (late 2015) -------------------------------
439 
440 
441 // shower::TrackShowerSeparationAlg::TrackShowerSeparationAlg(fhicl::ParameterSet const& pset) {
442 // this->reconfigure(pset);
443 // // tree
444 // // ftree = tfs->make<TTree>("tree","tree");
445 // // ftree->Branch("Run",&Run);
446 // // ftree->Branch("Event",&Event);
447 // // ftree->Branch("Distance",&Distance);
448 // // ftree->Branch("Angle",&Angle);
449 // // ftree->Branch("Length",&Length);
450 // // ftree->Branch("TrackID",&TrackID);
451 // // ftree->Branch("PDG",&pdg);
452 // // ftree->Branch("NSpacePoints",&NSpacePoints);
453 // // ftree->Branch("AvDistance",&AvDistance);
454 // }
455 
456 // void shower::TrackShowerSeparationAlg::reconfigure(fhicl::ParameterSet const& pset) {
457 // fConeAngle = pset.get<double>("ConeAngle");
458 
459 // // fAngleCut = pset.get<double>("AngleCut");
460 // // fDistanceCut = pset.get<double>("DistanceCut");
461 // // fVertexProximityCut = pset.get<double>("VertexProximityCut");
462 // // fTrackProximityCut = pset.get<double>("TrackProximityCut");
463 // // fAvTrackHitDistance = pset.get<double>("AvTrackHitDistance");
464 // }
465 
466 // void shower::TrackShowerSeparationAlg::IdentifyTracksFromEventCentre(const std::vector<art::Ptr<recob::Track> >& tracks,
467 // const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints,
468 // const art::FindManyP<recob::Track>& fmtsp) {
469 
470 // // Find the charge weighted centre of the entire event!
471 // // ---- except space points don't have charge (could look at associated hits, but cba!)
472 // TVector3 centre = TVector3(0,0,0);
473 // for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt)
474 // centre += TVector3((*spacePointIt)->XYZ()[0], (*spacePointIt)->XYZ()[1], (*spacePointIt)->XYZ()[2]);
475 // centre *= 1/(double)spacePoints.size();
476 
477 // TVector3 trackVertex, trackEnd, trackDirection;
478 
479 // // Look at all tracks
480 // for (std::vector<art::Ptr<recob::Track> >::const_iterator trackIt = tracks.begin(); trackIt != tracks.end(); ++trackIt) {
481 
482 // trackVertex = ( ((*trackIt)->Vertex()-centre).Mag() > ((*trackIt)->End()-centre).Mag() ) ? (*trackIt)->Vertex() : (*trackIt)->End();
483 // trackEnd = ( ((*trackIt)->Vertex()-centre).Mag() > ((*trackIt)->End()-centre).Mag() ) ? (*trackIt)->End() : (*trackIt)->Vertex();
484 // trackDirection = ( ((*trackIt)->Vertex()-centre).Mag() > ((*trackIt)->End()-centre).Mag() ) ? (*trackIt)->VertexDirection() : (-1)*(*trackIt)->EndDirection();
485 
486 // // Get the space points not associated with the current track
487 // std::vector<art::Ptr<recob::SpacePoint> > surroundingSpacePoints = GetSurroundingSpacePoints(spacePoints, fmtsp, (*trackIt)->ID());
488 
489 // // // TRUE -- use truth to find which particle this track is associated with --------
490 // // std::vector<art::Ptr<recob::Hit> > trackHits = fmht.at(trackIt->key());
491 // // int trueTrackID = FindTrueTrack(trackHits);
492 // // const simb::MCParticle* trueParticle = backtracker->TrackIDToParticle(trueTrackID);
493 // // pdg = trueParticle->PdgCode();
494 // // Length = (*trackIt)->Length();
495 // // TrackID = (*trackIt)->ID();
496 // // // -------------------------------------------------------------------------------
497 
498 // bool showerLike = IdentifyShowerLikeTrack(trackEnd, trackDirection, surroundingSpacePoints);
499 
500 // if (showerLike) {
501 // fShowerLikeIDs.push_back((*trackIt)->ID());
502 // continue;
503 // }
504 
505 // else
506 // fTrackLikeIDs.push_back((*trackIt)->ID());
507 
508 // }
509 
510 // return;
511 
512 // }
513 
514 // void shower::TrackShowerSeparationAlg::IdentifyTracksNearTracks(const std::vector<art::Ptr<recob::Track> >& tracks) {
515 
516 // std::vector<art::Ptr<recob::Track> > identifiedTracks;
517 
518 // for (std::vector<art::Ptr<recob::Track> >::const_iterator trackIt = tracks.begin(); trackIt != tracks.end(); ++trackIt)
519 // if (std::find(fTrackLikeIDs.begin(), fTrackLikeIDs.end(), (*trackIt)->ID()) != fTrackLikeIDs.end())
520 // identifiedTracks.push_back(*trackIt);
521 
522 // // Look through tracks
523 // bool allTracksRemoved = false;
524 // while (!allTracksRemoved) {
525 
526 // int tracksRemoved = 0;
527 
528 // // Look through all the tracks
529 // for (std::vector<art::Ptr<recob::Track> >::const_iterator trackIt = tracks.begin(); trackIt != tracks.end(); ++trackIt) {
530 
531 // // Don't consider tracks already tagged as tracks!
532 // if (std::find(fTrackLikeIDs.begin(), fTrackLikeIDs.end(), (*trackIt)->ID()) != fTrackLikeIDs.end())
533 // continue;
534 
535 // bool trackShouldBeRemoved = false;
536 
537 // // Find tracks which are close to previously identified tracks
538 // for (std::vector<art::Ptr<recob::Track> >::iterator identifiedTrackIt = identifiedTracks.begin(); identifiedTrackIt != identifiedTracks.end(); ++identifiedTrackIt) {
539 // if (std::find(fShowerLikeIDs.begin(), fShowerLikeIDs.end(), (*trackIt)->ID()) != fShowerLikeIDs.end())
540 // continue;
541 // if ( ( ((*trackIt)->Vertex() - (*identifiedTrackIt)->Vertex()).Mag() < fTrackProximityCut ) or
542 // ( ((*trackIt)->Vertex() - (*identifiedTrackIt)->End() ).Mag() < fTrackProximityCut ) or
543 // ( ((*trackIt)->End() - (*identifiedTrackIt)->Vertex()).Mag() < fTrackProximityCut ) or
544 // ( ((*trackIt)->End() - (*identifiedTrackIt)->End() ).Mag() < fTrackProximityCut ) )
545 // trackShouldBeRemoved = true;
546 // }
547 
548 // // Tag this track as a 'track-like' object
549 // if (trackShouldBeRemoved) {
550 // fTrackLikeIDs.push_back((*trackIt)->ID());
551 // identifiedTracks.push_back(*trackIt);
552 // ++tracksRemoved;
553 // }
554 
555 // }
556 
557 // // If there were no tracks removed then we'll call it a day
558 // if (tracksRemoved == 0)
559 // allTracksRemoved = true;
560 
561 // }
562 
563 // return;
564 
565 // }
566 
567 // void shower::TrackShowerSeparationAlg::IdentifyTracksNearVertex(const art::Ptr<recob::Vertex>& vertex,
568 // const std::vector<art::Ptr<recob::Track> >& tracks,
569 // const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints,
570 // const art::FindManyP<recob::Track>& fmtsp) {
571 
572 // double xyz[3];
573 // vertex->XYZ(xyz);
574 // TVector3 vertexPos = TVector3(xyz[0], xyz[1], xyz[2]);
575 
576 // // Look at each track
577 // for (std::vector<art::Ptr<recob::Track> >::const_iterator trackIt = tracks.begin(); trackIt != tracks.end(); ++trackIt) {
578 
579 // TVector3 end, direction;
580 
581 // // See if either end is close to the vertex
582 // if ( ((*trackIt)->Vertex() - vertexPos).Mag() < fVertexProximityCut or
583 // ((*trackIt)->End() - vertexPos).Mag() < fVertexProximityCut) {
584 // end = ((*trackIt)->Vertex() - vertexPos).Mag() < ((*trackIt)->End() - vertexPos).Mag() ? (*trackIt)->End() : (*trackIt)->VertexDirection();
585 // direction = ((*trackIt)->Vertex() - vertexPos).Mag() < ((*trackIt)->End() - vertexPos).Mag() ? (*trackIt)->VertexDirection() : (-1)*(*trackIt)->VertexDirection();
586 // }
587 
588 // else
589 // continue;
590 
591 // // Get the space points not associated with the current track
592 // std::vector<art::Ptr<recob::SpacePoint> > surroundingSpacePoints = GetSurroundingSpacePoints(spacePoints, fmtsp, (*trackIt)->ID());
593 
594 // // Make sure this track start isn't a shower start
595 // bool showerLike = IdentifyShowerLikeTrack(end, direction, surroundingSpacePoints);
596 
597 // if (showerLike) {
598 // fShowerLikeIDs.push_back((*trackIt)->ID());
599 // continue;
600 // }
601 
602 // // This track originates from near the interaction vertex and is not shower-like
603 // fTrackLikeIDs.push_back((*trackIt)->ID());
604 
605 // }
606 
607 // return;
608 
609 // }
610 
611 // bool shower::TrackShowerSeparationAlg::IdentifyShowerLikeTrack(const TVector3& end,
612 // const TVector3& direction,
613 // const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints) {
614 
615 // std::vector<art::Ptr<recob::SpacePoint> > spacePointsInCone = GetSpacePointsInCone(spacePoints, end, direction);
616 
617 // if (spacePointsInCone.size() < 2)
618 // return false;
619 
620 // // Get the average spread of these space points
621 // double spread = SpacePointSpread(spacePointsInCone);
622 
623 // if (spread < fAvTrackHitDistance)
624 // return false;
625 
626 // return true;
627 
628 // }
629 
630 // std::vector<art::Ptr<recob::Hit> > shower::TrackShowerSeparationAlg::FillHitsToCluster(const std::vector<art::Ptr<recob::Hit> >& initialHits,
631 // const art::FindManyP<recob::Track>& fmt) {
632 
633 // // Container to fill with shower-like hits
634 // std::vector<art::Ptr<recob::Hit> > hitsToCluster;
635 
636 // for (std::vector<art::Ptr<recob::Hit> >::const_iterator initialHit = initialHits.begin(); initialHit != initialHits.end(); ++initialHit) {
637 // std::vector<art::Ptr<recob::Track> > showerTracks = fmt.at(initialHit->key());
638 // if ( (showerTracks.size() and (std::find(fTrackLikeIDs.begin(), fTrackLikeIDs.end(), showerTracks.at(0)->ID()) == fTrackLikeIDs.end()))//hit on track-like track
639 // || !showerTracks.size() )//hit not on any track
640 // hitsToCluster.push_back(*initialHit);
641 // }
642 
643 // return hitsToCluster;
644 
645 // }
646 
647 // int shower::TrackShowerSeparationAlg::FindTrackID(const art::Ptr<recob::Hit>& hit) {
648 // double particleEnergy = 0;
649 // int likelyTrackID = 0;
650 // std::vector<sim::TrackIDE> trackIDs = backtracker->HitToTrackID(hit);
651 // for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
652 // if (trackIDs.at(idIt).energy > particleEnergy) {
653 // particleEnergy = trackIDs.at(idIt).energy;
654 // likelyTrackID = TMath::Abs(trackIDs.at(idIt).trackID);
655 // }
656 // }
657 // return likelyTrackID;
658 // }
659 
660 // int shower::TrackShowerSeparationAlg::FindTrueTrack(const std::vector<art::Ptr<recob::Hit> >& trackHits) {
661 // std::map<int,double> trackMap;
662 // for (std::vector<art::Ptr<recob::Hit> >::const_iterator trackHitIt = trackHits.begin(); trackHitIt != trackHits.end(); ++trackHitIt) {
663 // art::Ptr<recob::Hit> hit = *trackHitIt;
664 // int trackID = FindTrackID(hit);
665 // trackMap[trackID] += hit->Integral();
666 // }
667 // //return std::max_element(trackMap.begin(), trackMap.end(), [](const std::pair<int,double>& p1, const std::pair<int,double>& p2) {return p1.second < p2.second;} )->first;
668 // double highestCharge = 0;
669 // int clusterTrack = 0;
670 // for (std::map<int,double>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end(); ++trackIt)
671 // if (trackIt->second > highestCharge) {
672 // highestCharge = trackIt->second;
673 // clusterTrack = trackIt->first;
674 // }
675 // return clusterTrack;
676 // }
677 
678 // std::vector<art::Ptr<recob::SpacePoint> > shower::TrackShowerSeparationAlg::GetSurroundingSpacePoints(const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints,
679 // const art::FindManyP<recob::Track>& fmt,
680 // unsigned int trackID) {
681 
682 // // The space points to return
683 // std::vector<art::Ptr<recob::SpacePoint> > surroundingSpacePoints;
684 
685 // for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
686 
687 // bool spacePointIsInCurrentTrack = false;
688 
689 // std::vector<art::Ptr<recob::Track> > spacePointTracks = fmt.at(spacePointIt->key());
690 // for (std::vector<art::Ptr<recob::Track> >::iterator spacePointTrackIt = spacePointTracks.begin(); spacePointTrackIt != spacePointTracks.end(); ++spacePointTrackIt)
691 // if (spacePointTrackIt->key() == trackID) spacePointIsInCurrentTrack = true;
692 
693 // if (!spacePointIsInCurrentTrack)
694 // surroundingSpacePoints.push_back(*spacePointIt);
695 
696 // }
697 
698 // return surroundingSpacePoints;
699 
700 // }
701 
702 // std::vector<art::Ptr<recob::SpacePoint> > shower::TrackShowerSeparationAlg::GetSpacePointsInCone(const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints,
703 // const TVector3& trackEnd,
704 // const TVector3& trackDirection) {
705 
706 // // The space points in cone to return
707 // std::vector<art::Ptr<recob::SpacePoint> > spacePointsInCone;
708 
709 // TVector3 spacePointPos, spacePointProj;
710 // double displacementFromAxis, distanceFromTrackEnd, projDistanceFromTrackEnd, angleFromTrackEnd;
711 
712 // for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
713 
714 // // Get the properties of this space point
715 // const double* xyz = (*spacePointIt)->XYZ();
716 // spacePointPos = TVector3(xyz[0], xyz[1], xyz[2]);
717 // spacePointProj = ((spacePointPos-trackEnd).Dot(trackDirection))*trackDirection + trackEnd;
718 // displacementFromAxis = (spacePointProj - spacePointPos).Mag();
719 // distanceFromTrackEnd = (spacePointPos - trackEnd).Mag();
720 // projDistanceFromTrackEnd = (spacePointProj - trackEnd).Mag();
721 // angleFromTrackEnd = TMath::ASin(displacementFromAxis/distanceFromTrackEnd) * 180 / TMath::Pi();
722 
723 // if ( (projDistanceFromTrackEnd < fDistanceCut) and (angleFromTrackEnd < fAngleCut) and ((spacePointProj-trackEnd).Dot(trackDirection) > 0) )
724 // spacePointsInCone.push_back(*spacePointIt);
725 
726 // }
727 
728 // return spacePointsInCone;
729 
730 // }
731 
732 // std::vector<art::Ptr<recob::Hit> > shower::TrackShowerSeparationAlg::RemoveTrackHits(const std::vector<art::Ptr<recob::Hit> >& initialHits,
733 // const std::vector<art::Ptr<recob::Track> >& tracks,
734 // const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints,
735 // const std::vector<art::Ptr<recob::Vertex> >& vertices,
736 // const art::FindManyP<recob::Track>& fmth,
737 // const art::FindManyP<recob::Track>& fmtsp,
738 // const art::FindManyP<recob::Hit>& fmh,
739 // int event,
740 // int run) {
741 
742 // Event = event;
743 // Run = run;
744 
745 // if (spacePoints.size() == 0)
746 // return initialHits;
747 
748 // // Container for shower-like hits
749 // std::vector<art::Ptr<recob::Hit> > hitsToCluster;
750 
751 // fTrackLikeIDs.clear();
752 // fShowerLikeIDs.clear();
753 
754 // // Find the vertex furthest upstream (if it exists)
755 // art::Ptr<recob::Vertex> vertex;
756 // for (std::vector<art::Ptr<recob::Vertex> >::const_iterator vertexIt = vertices.begin(); vertexIt != vertices.end(); ++vertexIt) {
757 // double xyzNew[3], xyzOld[3];
758 // (*vertexIt)->XYZ(xyzNew);
759 // if (vertex.isNull())
760 // vertex = *vertexIt;
761 // else {
762 // vertex->XYZ(xyzOld);
763 // if (xyzNew[2] < xyzOld[2])
764 // vertex = *vertexIt;
765 // }
766 // }
767 
768 // // If we have a vertex then, for the time being, use it to remove tracks which are too close!
769 // if (!vertex.isNull())
770 // this->IdentifyTracksNearVertex(vertex, tracks, spacePoints, fmtsp);
771 
772 // // If no vertex, things are harder
773 // else
774 // this->IdentifyTracksFromEventCentre(tracks, spacePoints, fmtsp);
775 
776 // // Once we've identified some tracks, can look for others at the ends
777 // this->IdentifyTracksNearTracks(tracks);
778 
779 // hitsToCluster = FillHitsToCluster(initialHits, fmth);
780 
781 // return hitsToCluster;
782 
783 // }
784 
785 // std::vector<art::Ptr<recob::Hit> > shower::TrackShowerSeparationAlg::RemoveTrackHits(const std::vector<art::Ptr<recob::Hit> >& hits,
786 // const std::vector<art::Ptr<recob::PFParticle> > pfParticles,
787 // const art::FindManyP<recob::Cluster>& fmc,
788 // const art::FindManyP<recob::Hit>& fmh) {
789 
790 // std::vector<art::Ptr<recob::Hit> > showerHits;
791 
792 // // Use information from Pandora to identify shower-like hits
793 // for (std::vector<art::Ptr<recob::PFParticle> >::const_iterator pfParticleIt = pfParticles.begin(); pfParticleIt != pfParticles.end(); ++pfParticleIt) {
794 
795 // // See if this is a shower particle
796 // if ((*pfParticleIt)->PdgCode() == 11) {
797 // std::vector<art::Ptr<recob::Cluster> > clusters = fmc.at(pfParticleIt->key());
798 // for (std::vector<art::Ptr<recob::Cluster> >::iterator clusterIt = clusters.begin(); clusterIt != clusters.end(); ++clusterIt) {
799 // std::vector<art::Ptr<recob::Hit> > hits = fmh.at(clusterIt->key());
800 // for (std::vector<art::Ptr<recob::Hit> >::iterator hitIt = hits.begin(); hitIt != hits.end(); ++hitIt)
801 // showerHits.push_back(*hitIt);
802 // }
803 // }
804 
805 // }
806 
807 // return showerHits;
808 
809 // }
810 
811 // double shower::TrackShowerSeparationAlg::SpacePointSpread(const std::vector<art::Ptr<recob::SpacePoint> >& spacePoints) {
812 
813 // /// Finds the spread of a set of space points about their central axis
814 
815 // // Find the centre of the space points
816 // TVector3 centre = TVector3(0,0,0);
817 // for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt)
818 // centre += TVector3((*spacePointIt)->XYZ()[0], (*spacePointIt)->XYZ()[1], (*spacePointIt)->XYZ()[2]);
819 // centre *= 1/(double)spacePoints.size();
820 
821 // // Find the central axis of the space points
822 // TPrincipal* pca = new TPrincipal(3,"");
823 // for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt)
824 // pca->AddRow((*spacePointIt)->XYZ());
825 // pca->MakePrincipals();
826 // const TMatrixD* eigenvectors = pca->GetEigenVectors();
827 // TVector3 spacePointDir = TVector3((*eigenvectors)[0][0], (*eigenvectors)[1][0], (*eigenvectors)[2][0]);
828 // delete pca;
829 
830 // // See if the space points form something which may resemble a track (i.e. straight line) or a shower
831 // double avDistanceFromCentralAxis = 0;
832 // for (std::vector<art::Ptr<recob::SpacePoint> >::const_iterator spacePointIt = spacePoints.begin(); spacePointIt != spacePoints.end(); ++spacePointIt) {
833 // const double* xyz = (*spacePointIt)->XYZ();
834 // TVector3 spacePointPos = TVector3(xyz[0], xyz[1], xyz[2]);
835 // TVector3 projectionOntoCentralAxis = ((spacePointPos-centre).Dot(spacePointDir))*spacePointDir + centre;
836 // avDistanceFromCentralAxis += (spacePointPos - projectionOntoCentralAxis).Mag();
837 // }
838 // avDistanceFromCentralAxis /= spacePoints.size();
839 
840 // return avDistanceFromCentralAxis;
841 
842 // }
shower::ReconTrack::SetDirection
void SetDirection(TVector3 direction)
Definition: TrackShowerSeparationAlg.h:55
shower::TrackShowerSeparationAlg::InitialTrackLikeSegment
std::vector< int > InitialTrackLikeSegment(std::map< int, std::unique_ptr< ReconTrack > > &reconTracks) const
Definition: TrackShowerSeparationAlg.cxx:284
shower::ReconTrack::SetEnd
void SetEnd(TVector3 end)
Definition: TrackShowerSeparationAlg.h:52
tracks
ClusterModuleLabel join with tracks
Definition: trackfinderalgorithms.fcl:56
shower::TrackShowerSeparationAlg::reconfigure
void reconfigure(fhicl::ParameterSet const &pset)
Read in configurable parameters from provided parameter set.
Definition: TrackShowerSeparationAlg.cxx:22
shower::ReconTrack::SetLength
void SetLength(double length)
Definition: TrackShowerSeparationAlg.h:53
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::ReconTrack::SetSpacePoints
void SetSpacePoints(std::vector< art::Ptr< recob::SpacePoint > > spacePoints)
Definition: TrackShowerSeparationAlg.h:57
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::ReconTrack::SetVertex
void SetVertex(TVector3 vertex)
Definition: TrackShowerSeparationAlg.h:51
shower::TrackShowerSeparationAlg::SpacePointsRMS
double SpacePointsRMS(const std::vector< art::Ptr< recob::SpacePoint > > &spacePoints) const
Definition: TrackShowerSeparationAlg.cxx:416
shower::TrackShowerSeparationAlg::TrackShowerSeparationAlg
TrackShowerSeparationAlg(fhicl::ParameterSet const &pset)
Definition: TrackShowerSeparationAlg.cxx:18
shower::TrackShowerSeparationAlg::SelectShowerHits
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
Definition: TrackShowerSeparationAlg.cxx:32
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::ReconTrack::SetVertexDir
void SetVertexDir(TVector3 vertexDir)
Definition: TrackShowerSeparationAlg.h:54
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
shower::ReconTrack::SetHits
void SetHits(std::vector< art::Ptr< recob::Hit > > hits)
Definition: TrackShowerSeparationAlg.h:56
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
geo::origin
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
Generated by   doxygen 1.8.5