LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MergeClusterAlg.cxx
Go to the documentation of this file.
1 // Merge Cluster algorithm
3 //
4 // Runs on the output of previous clustering algorithms to merge
5 // clusters together which lie on a straight line and are within
6 // some separation threshold.
7 // Runs recursively over all clusters, including new ones formed
8 // in the algorithm.
9 //
10 // M Wallbank (m.wallbank@sheffield.ac.uk), July 2015
12 
14 
15 #include "TPrincipal.h"
16 #include "TTree.h"
17 
19 {
20  this->reconfigure(pset);
21  fTree = tfs->make<TTree>("MatchingVariables", "MatchingVariables");
22  fTree->Branch("Angle", &fAngle);
23  fTree->Branch("Eigenvalue", &fEigenvalue);
24  fTree->Branch("Cluster1Size", &fCluster1Size);
25  fTree->Branch("Cluster2Size", &fCluster2Size);
26  fTree->Branch("Length1", &fLength1);
27  fTree->Branch("Length2", &fLength2);
28  fTree->Branch("Separation", &fSeparation);
29  fTree->Branch("CrossingDistance", &fCrossingDistance);
30  fTree->Branch("ProjectedWidth", &fProjectedWidth);
31  fTree->Branch("Overlap", &fOverlap);
32  fTree->Branch("TrueMerge", &fTrueMerge);
33 }
34 
36  TVector2 const& centre,
37  TVector2 const& direction,
38  TVector2& start,
39  TVector2& end) const
40 {
41 
43 
44  TVector2 pos;
45  std::map<double, TVector2> hitProjection;
46 
47  // Project all hits onto line to determine end points
48  for (auto& hit : cluster) {
49  pos = HitCoordinates(hit) - centre;
50  hitProjection[direction * pos] = pos;
51  }
52 
53  // Project end points onto line which passes through centre of cluster
54  start = hitProjection.begin()->second.Proj(direction) + centre;
55  end = hitProjection.rbegin()->second.Proj(direction) + centre;
56 
57  return;
58 }
59 
60 double cluster::MergeClusterAlg::FindClusterOverlap(TVector2 const& direction,
61  TVector2 const& centre,
62  TVector2 const& start1,
63  TVector2 const& end1,
64  TVector2 const& start2,
65  TVector2 const& end2) const
66 {
67 
69 
70  double clusterOverlap = 0;
71 
72  // Project onto the average direction through both clusters
73  double s1 = (start1 - centre) * direction;
74  double e1 = (end1 - centre) * direction;
75  double s2 = (start2 - centre) * direction;
76  double e2 = (end2 - centre) * direction;
77 
78  // Make sure end > start
79  if (s1 > e1) {
80  std::cout << "s1>e1: " << s1 << " and " << e1 << std::endl;
81  double tmp = e1;
82  e1 = s1;
83  s1 = tmp;
84  }
85  if (s2 > e2) {
86  std::cout << "s1>e1: " << s1 << " and " << e1 << std::endl;
87  double tmp = e2;
88  e2 = s2;
89  s2 = tmp;
90  }
91 
92  // Find the overlap of the clusters on the centre line
93  if ((e1 > s2) && (e2 > s1)) clusterOverlap = std::min((e1 - s2), (e2 - s1));
94 
95  return clusterOverlap;
96 }
97 
98 double cluster::MergeClusterAlg::FindCrossingDistance(TVector2 const& direction1,
99  TVector2 const& centre1,
100  TVector2 const& direction2,
101  TVector2 const& centre2) const
102 {
103 
105 
106  // Find intersection point of two lines drawn through the centre of the clusters
107  double dcross = (direction1.X() * direction2.Y()) - (direction1.Y() * direction2.X());
108  TVector2 p = centre2 - centre1;
109  double pcrossd = (p.X() * direction2.Y()) - (p.Y() * direction2.X());
110  TVector2 crossing = centre1 + ((pcrossd / dcross) * direction1);
111 
112  // Get distance from this point to the clusters
113  double crossingDistance = std::min((centre1 - crossing).Mod(), (centre2 - crossing).Mod());
114 
115  return crossingDistance;
116 }
117 
119  art::PtrVector<recob::Hit> const& cluster2) const
120 {
121 
123 
124  double minDistance = 99999.;
125 
126  // Loop over the two clusters to find the smallest distance
127  for (auto const& hit1 : cluster1) {
128  for (auto const& hit2 : cluster2) {
129 
130  TVector2 pos1 = HitCoordinates(hit1);
131  TVector2 pos2 = HitCoordinates(hit2);
132 
133  double distance = (pos1 - pos2).Mod();
134 
135  if (distance < minDistance) minDistance = distance;
136  }
137  }
138 
139  return minDistance;
140 }
141 
142 double cluster::MergeClusterAlg::FindProjectedWidth(TVector2 const& centre1,
143  TVector2 const& start1,
144  TVector2 const& end1,
145  TVector2 const& centre2,
146  TVector2 const& start2,
147  TVector2 const& end2) const
148 {
149 
151 
152  // Get the line running through the centre of the two clusters
153  TVector2 parallel = (centre2 - centre1).Unit();
154  TVector2 perpendicular = parallel.Rotate(TMath::Pi() / 2);
155 
156  // Project the cluster vector onto this perpendicular line
157  double s1 = (start1 - centre1) * perpendicular;
158  double e1 = (end1 - centre1) * perpendicular;
159  double s2 = (start2 - centre2) * perpendicular;
160  double e2 = (end2 - centre2) * perpendicular;
161 
162  // Find the width in each direction
163  double projectionStart = std::max(TMath::Abs(s1), TMath::Abs(s2));
164  double projectionEnd = std::max(TMath::Abs(e1), TMath::Abs(e2));
165 
166  double projectionWidth = projectionStart + projectionEnd;
167 
168  return projectionWidth;
169 }
170 
172 {
174 
175  auto const wireCenter = fGeom->WireIDToWireGeo(wireID).GetCenter();
176  geo::PlaneID const planeID{wireID.Cryostat, wireID.TPC % 2, wireID.Plane};
177 
178  if (fGeom->SignalType(wireID) == geo::kInduction) {
179  return fGeom->WireCoordinate(wireCenter, planeID);
180  }
181  return wireID.Wire + ((wireID.TPC / 2) * fGeom->Nwires(planeID));
182 }
183 
185 {
187 
188  return TVector2(GlobalWire(hit->WireID()), hit->PeakTime());
189 }
190 
192  std::vector<art::PtrVector<recob::Hit>> const& planeClusters,
193  std::vector<art::PtrVector<recob::Hit>>& clusters) const
194 {
196 
197  std::vector<unsigned int> mergedClusters;
198 
199  std::vector<art::PtrVector<recob::Hit>> oldClusters = planeClusters;
200 
201  // Sort the clusters by size
202  std::sort(oldClusters.begin(),
203  oldClusters.end(),
205  return a.size() > b.size();
206  });
207 
208  // Find the numbers of clusters above size threshold
209  unsigned int nclusters = 0;
210  for (auto& cluster : oldClusters)
211  if (cluster.size() >= fMinMergeClusterSize) ++nclusters;
212 
213  // Until all clusters are merged, create new clusters
214  bool mergedAllClusters = false;
215  while (!mergedAllClusters) {
216 
217  // New cluster
219 
220  // Put the largest unmerged cluster in this new cluster
221  for (unsigned int initCluster = 0; initCluster < oldClusters.size(); ++initCluster) {
222  if (oldClusters.at(initCluster).size() < fMinMergeClusterSize or
223  std::find(mergedClusters.begin(), mergedClusters.end(), initCluster) !=
224  mergedClusters.end())
225  continue;
226  cluster = oldClusters.at(initCluster);
227  mergedClusters.push_back(initCluster);
228  break;
229  }
230 
231  // Merge all aligned clusters to this
232  bool mergedAllToThisCluster = false;
233  while (!mergedAllToThisCluster) {
234 
235  // Look at all clusters and merge
236  int nadded = 0;
237  for (unsigned int trialCluster = 0; trialCluster < oldClusters.size(); ++trialCluster) {
238 
239  if (oldClusters.at(trialCluster).size() < fMinMergeClusterSize or
240  std::find(mergedClusters.begin(), mergedClusters.end(), trialCluster) !=
241  mergedClusters.end())
242  continue;
243 
244  // Calculate the PCA for each
245  TPrincipal *pca1 = new TPrincipal(2, ""), *pca2 = new TPrincipal(2, "");
246  double hits[2];
247  TVector2 pos;
248 
249  // Cluster centre
250  TVector2 chargePoint1 = TVector2(0, 0), chargePoint2 = TVector2(0, 0);
251  double totalCharge1 = 0, totalCharge2 = 0;
252 
253  for (auto& hit1 : cluster) {
254  pos = HitCoordinates(hit1);
255  hits[0] = pos.X();
256  hits[1] = pos.Y();
257  pca1->AddRow(hits);
258  chargePoint1 += hit1->Integral() * pos;
259  totalCharge1 += hit1->Integral();
260  }
261  for (auto& hit2 : oldClusters.at(trialCluster)) {
262  pos = HitCoordinates(hit2);
263  hits[0] = pos.X();
264  hits[1] = pos.Y();
265  pca2->AddRow(hits);
266  chargePoint2 += hit2->Integral() * pos;
267  totalCharge2 += hit2->Integral();
268  }
269 
270  pca1->MakePrincipals();
271  pca2->MakePrincipals();
272 
273  // Properties of these clusters
274  TVector2 direction1 =
275  TVector2((*pca1->GetEigenVectors())[0][0], (*pca1->GetEigenVectors())[1][0]).Unit();
276  TVector2 direction2 =
277  TVector2((*pca2->GetEigenVectors())[0][0], (*pca2->GetEigenVectors())[1][0]).Unit();
278  TVector2 direction = ((direction1 + direction2) / 2).Unit();
279  TVector2 centre1 = chargePoint1 / totalCharge1;
280  TVector2 centre2 = chargePoint2 / totalCharge2;
281  TVector2 centre = (centre1 + centre2) / 2;
282  TVector2 start1, end1;
283  TVector2 start2, end2;
284  FindClusterEndPoints(cluster, centre1, direction1, start1, end1);
285  FindClusterEndPoints(oldClusters.at(trialCluster), centre2, direction2, start2, end2);
286  double length1 = (end1 - start1).Mod();
287  double length2 = (end2 - start2).Mod();
288 
289  // Properties of the pair of clusters
290  double crossingDistance = FindCrossingDistance(direction1, centre1, direction2, centre2);
291  double projectedWidth = FindProjectedWidth(centre1, start1, end1, centre2, start2, end2);
292  double angle = direction1.DeltaPhi(direction2);
293  if (angle > 1.57) angle = 3.14159 - angle;
294  double separation = FindMinSeparation(cluster, oldClusters.at(trialCluster));
295 
296  if (separation > fMaxMergeSeparation) continue;
297  if (PassCuts(
298  angle, crossingDistance, projectedWidth, separation, TMath::Max(length1, length2))) {
299 
300  for (auto& hit : oldClusters.at(trialCluster))
301  cluster.push_back(hit);
302 
303  mergedClusters.push_back(trialCluster);
304  ++nadded;
305  }
306 
307  delete pca1;
308  delete pca2;
309 
310  } // loop over clusters to add
311 
312  if (nadded == 0) mergedAllToThisCluster = true;
313 
314  } // while loop
315 
316  clusters.push_back(cluster);
317  if (mergedClusters.size() == nclusters) mergedAllClusters = true;
318  }
319 
320  return clusters.size();
321 }
322 
323 bool cluster::MergeClusterAlg::PassCuts(double const angle,
324  double const crossingDistance,
325  double const projectedWidth,
326  double const separation,
327  double const longLength) const
328 {
329 
331 
332  bool passCrossingDistanceAngle = false;
333  if (crossingDistance < (-2 + (5 / (1 * TMath::Abs(angle)) - 0))) passCrossingDistanceAngle = true;
334 
335  bool passSeparationAngle = false;
336  if (separation < (200 * TMath::Abs(angle) + 40)) passSeparationAngle = true;
337 
338  bool passProjectedWidth = false;
339  if (((double)projectedWidth / (double)longLength) < fProjWidthThreshold)
340  passProjectedWidth = true;
341 
342  return passCrossingDistanceAngle and passSeparationAngle and passProjectedWidth;
343 }
344 
346 {
347  fMinMergeClusterSize = p.get<int>("MinMergeClusterSize");
348  fMaxMergeSeparation = p.get<double>("MaxMergeSeparation");
349  fProjWidthThreshold = p.get<double>("ProjWidthThreshold");
350 }
MergeClusterAlg(fhicl::ParameterSet const &pset)
Length_t WireCoordinate(Point_t const &pos, PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Point_t const & GetCenter() const
Returns the world coordinate of the center of the wire [cm].
Definition: WireGeo.h:221
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
int MergeClusters(std::vector< art::PtrVector< recob::Hit >> const &planeClusters, std::vector< art::PtrVector< recob::Hit >> &clusters) const
art::ServiceHandle< geo::Geometry const > fGeom
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
void FindClusterEndPoints(art::PtrVector< recob::Hit > const &cluster, TVector2 const &centre, TVector2 const &direction, TVector2 &start, TVector2 &end) const
art::ServiceHandle< art::TFileService const > tfs
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
Float_t tmp
Definition: plot.C:35
Cluster finding and building.
void reconfigure(fhicl::ParameterSet const &p)
double FindClusterOverlap(TVector2 const &direction, TVector2 const &centre, TVector2 const &start1, TVector2 const &end1, TVector2 const &start2, TVector2 const &end2) const
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void hits()
Definition: readHits.C:15
Signal from induction planes.
Definition: geo_types.h:151
T get(std::string const &key) const
Definition: ParameterSet.h:314
reference at(size_type n)
Definition: PtrVector.h:359
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
size_type size() const
Definition: PtrVector.h:302
Detector simulation of raw signals on wires.
TVector2 HitCoordinates(art::Ptr< recob::Hit > const &hit) const
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:220
SigType_t SignalType(PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
double FindMinSeparation(art::PtrVector< recob::Hit > const &cluster1, art::PtrVector< recob::Hit > const &cluster2) const
double FindCrossingDistance(TVector2 const &direction1, TVector2 const &centre1, TVector2 const &direction2, TVector2 const &centre2) const
unsigned int Nwires(PlaneID const &planeid) const
Returns the total number of wires in the specified plane.
double GlobalWire(geo::WireID const &wireID) const
bool PassCuts(double angle, double crossingDistance, double projectedWidth, double separation, double longLength) const
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
double FindProjectedWidth(TVector2 const &centre1, TVector2 const &start1, TVector2 const &end1, TVector2 const &centre2, TVector2 const &start2, TVector2 const &end2) const
unsigned int fMinMergeClusterSize