LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
cluster::ClusterMatchTQ Class Reference

#include "ClusterMatchTQ.h"

Public Member Functions

 ClusterMatchTQ (fhicl::ParameterSet const &pset)
 
std::vector< std::vector< unsigned int > > MatchedClusters (const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Cluster >> &clusterlist, const art::FindManyP< recob::Hit > &fm) const
 

Private Attributes

double fKSCut
 
bool fEnableU
 
bool fEnableV
 
bool fEnableZ
 

Detailed Description

Definition at line 30 of file ClusterMatchTQ.h.

Constructor & Destructor Documentation

cluster::ClusterMatchTQ::ClusterMatchTQ ( fhicl::ParameterSet const &  pset)

Definition at line 40 of file ClusterMatchTQ.cxx.

References fhicl::ParameterSet::get().

41  {
42  fKSCut = pset.get<double>("KSCut");
43  fEnableU = pset.get<bool>("EnableU");
44  fEnableV = pset.get<bool>("EnableV");
45  fEnableZ = pset.get<bool>("EnableZ");
46  }

Member Function Documentation

std::vector< std::vector< unsigned int > > cluster::ClusterMatchTQ::MatchedClusters ( const detinfo::DetectorPropertiesData detProp,
const std::vector< art::Ptr< recob::Cluster >> &  clusterlist,
const art::FindManyP< recob::Hit > &  fm 
) const

Definition at line 49 of file ClusterMatchTQ.cxx.

References util::begin(), bin, c1, c2, detinfo::DetectorPropertiesData::ConvertTicksToX(), util::end(), detinfo::DetectorPropertiesData::GetXTicksOffset(), geo::kU, geo::kV, geo::kZ, geo::GeometryCore::Nplanes(), detinfo::DetectorPropertiesData::NumberTimeSamples(), recob::tracking::Plane::Plane(), util::size(), SortByWire(), t1, tmp, geo::GeometryCore::WirePitch(), and x1.

Referenced by trkf::CosmicTracker::produce().

53  {
54  std::vector<std::vector<unsigned int>> matchedclusters;
55 
56  // get services
58  int nplanes = geom->Nplanes();
59  int nts = detProp.NumberTimeSamples();
60 
61  std::vector<std::vector<TH1D>> signals(nplanes);
62 
63  std::vector<std::vector<unsigned int>> Cls(nplanes);
64  std::vector<std::vector<CluLen>> clulens(nplanes);
65 
66  for (size_t iclu = 0; iclu < clusterlist.size(); ++iclu) {
67 
68  float wire_pitch = geom->WirePitch(clusterlist[iclu]->Plane());
69 
70  float w0 = clusterlist[iclu]->StartWire();
71  float w1 = clusterlist[iclu]->EndWire();
72  float t0 = clusterlist[iclu]->StartTick();
73  float t1 = clusterlist[iclu]->EndTick();
74 
75  CluLen clulen;
76  clulen.index = iclu;
77 
78  auto const x0 = detProp.ConvertTicksToX(t0,
79  clusterlist[iclu]->Plane().Plane,
80  clusterlist[iclu]->Plane().TPC,
81  clusterlist[iclu]->Plane().Cryostat);
82  auto const x1 = detProp.ConvertTicksToX(t1,
83  clusterlist[iclu]->Plane().Plane,
84  clusterlist[iclu]->Plane().TPC,
85  clusterlist[iclu]->Plane().Cryostat);
86  clulen.length = std::hypot((w0 - w1) * wire_pitch, x0 - x1);
87 
88  switch (clusterlist[iclu]->View()) {
89  case geo::kU:
90  if (fEnableU) clulens[0].push_back(clulen);
91  break;
92  case geo::kV:
93  if (fEnableV) clulens[1].push_back(clulen);
94  break;
95  case geo::kZ:
96  if (fEnableZ) clulens[2].push_back(clulen);
97  break;
98  default: break;
99  }
100  }
101 
102  //sort clusters based on 2D length
103  for (size_t i = 0; i < clulens.size(); ++i) {
104  std::sort(clulens[i].begin(), clulens[i].end(), SortByLength);
105  for (size_t j = 0; j < clulens[i].size(); ++j) {
106  Cls[i].push_back(clulens[i][j].index);
107  }
108  }
109 
110  for (int i = 0; i < nplanes; ++i) {
111  for (size_t ic = 0; ic < Cls[i].size(); ++ic) {
112  TH1D sig(
113  Form("sig_%d_%d", i, int(ic)), Form("sig_%d_%d", i, int(ic)), nts + 100, -100, nts);
114  TH1D sigint(
115  Form("sigint_%d_%d", i, int(ic)), Form("sigint_%d_%d", i, int(ic)), nts + 100, -100, nts);
116  std::vector<art::Ptr<recob::Hit>> hitlist = fm.at(Cls[i][ic]);
117  std::sort(hitlist.begin(), hitlist.end(), SortByWire);
118 
119  for (auto theHit = hitlist.begin(); theHit != hitlist.end(); theHit++) {
120 
121  double time = (*theHit)->PeakTime();
122  time -= detProp.GetXTicksOffset(
123  (*theHit)->WireID().Plane, (*theHit)->WireID().TPC, (*theHit)->WireID().Cryostat);
124 
125  double charge = (*theHit)->Integral();
126  int bin = sig.FindBin(time);
127  sig.SetBinContent(bin, sig.GetBinContent(bin) + charge);
128  for (int j = bin; j <= sig.GetNbinsX(); ++j) {
129  sigint.SetBinContent(j, sigint.GetBinContent(j) + charge);
130  }
131  }
132  if (sigint.Integral()) sigint.Scale(1. / sigint.GetBinContent(sigint.GetNbinsX()));
133  signals[i].push_back(sigint);
134  }
135  }
136 
137  //matching clusters between different views
138  std::vector<int> matched(clusterlist.size());
139  for (size_t i = 0; i < clusterlist.size(); ++i)
140  matched[i] = 0;
141 
142  for (int i = 0; i < nplanes - 1; ++i) {
143  for (int j = i + 1; j < nplanes; ++j) {
144  for (size_t c1 = 0; c1 < Cls[i].size(); ++c1) {
145  for (size_t c2 = 0; c2 < Cls[j].size(); ++c2) {
146 
147  // check if both are the same view
148  if (clusterlist[Cls[i][c1]]->View() == clusterlist[Cls[j][c2]]->View()) continue;
149  // check if both are in the same cryostat and tpc
150  if (clusterlist[Cls[i][c1]]->Plane().Cryostat !=
151  clusterlist[Cls[j][c2]]->Plane().Cryostat)
152  continue;
153  if (clusterlist[Cls[i][c1]]->Plane().TPC != clusterlist[Cls[j][c2]]->Plane().TPC)
154  continue;
155  // check if both are already in the matched list
156  if (matched[Cls[i][c1]] == 1 && matched[Cls[j][c2]] == 1) continue;
157  // KS test between two views in time
158  double ks = 0;
159  if (signals[i][c1].Integral() && signals[j][c2].Integral())
160  ks = signals[i][c1].KolmogorovTest(&signals[j][c2]);
161  else {
162  mf::LogWarning("ClusterMatchTQ")
163  << "One of the two clusters appears to be empty: " << clusterlist[Cls[i][c1]]->ID()
164  << " " << clusterlist[Cls[j][c2]]->ID() << " " << i << " " << j << " " << c1 << " "
165  << c2 << " " << signals[i][c1].Integral() << " " << signals[j][c2].Integral();
166  }
167  //hks->Fill(ks);
168  int imatch = -1; //track candidate index
169  int iadd = -1; //cluster index to be inserted
170  if (ks > fKSCut) { //pass KS test
171  // check both clusters with all matched clusters
172  // if one is already matched,
173  // check if need to add the other to the same track candidate
174  for (size_t l = 0; l < matchedclusters.size(); ++l) {
175  for (size_t m = 0; m < matchedclusters[l].size(); ++m) {
176  if (matchedclusters[l][m] == Cls[i][c1]) {
177  imatch = l; //track candidate
178  iadd = j; //consider the other cluster
179  }
180  else if (matchedclusters[l][m] == Cls[j][c2]) {
181  imatch = l; //track candidate
182  iadd = i; //consider the other cluster
183  }
184  }
185  }
186  if (imatch >= 0) {
187  if (iadd == i) {
188  bool matchview = false;
189  // check if one matched cluster has the same view
190  for (size_t ii = 0; ii < matchedclusters[imatch].size(); ++ii) {
191  if (clusterlist[matchedclusters[imatch][ii]]->View() ==
192  clusterlist[Cls[i][c1]]->View()) {
193  matchview = true;
194  //replace if the new cluster has more hits
195  if (fm.at(Cls[i][c1]).size() > fm.at(matchedclusters[imatch][ii]).size()) {
196  matched[matchedclusters[imatch][ii]] = 0;
197  matchedclusters[imatch][ii] = Cls[i][c1];
198  matched[Cls[i][c1]] = 1;
199  }
200  }
201  }
202  if (!matchview) { //not matched view found, just add
203  matchedclusters[imatch].push_back(Cls[i][c1]);
204  matched[Cls[i][c1]] = 1;
205  }
206  }
207  else {
208  bool matchview = false;
209  for (size_t jj = 0; jj < matchedclusters[imatch].size(); ++jj) {
210  if (clusterlist[matchedclusters[imatch][jj]]->View() ==
211  clusterlist[Cls[j][c2]]->View()) {
212  matchview = true;
213  //replace if it has more hits
214  if (fm.at(Cls[j][c2]).size() > fm.at(matchedclusters[imatch][jj]).size()) {
215  matched[matchedclusters[imatch][jj]] = 0;
216  matchedclusters[imatch][jj] = Cls[j][c2];
217  matched[Cls[j][c2]] = 1;
218  }
219  }
220  }
221  if (!matchview) {
222  matchedclusters[imatch].push_back(Cls[j][c2]);
223  matched[Cls[j][c2]] = 1;
224  }
225  }
226  }
227  else {
228  std::vector<unsigned int> tmp;
229  tmp.push_back(Cls[i][c1]);
230  tmp.push_back(Cls[j][c2]);
231  matchedclusters.push_back(tmp);
232  matched[Cls[i][c1]] = 1;
233  matched[Cls[j][c2]] = 1;
234  }
235  } //pass KS test
236  } //c2
237  } //c1
238  } //j
239  } //i
240 
241  for (size_t i = 0; i < matchedclusters.size(); ++i) {
242  if (matchedclusters[i].size())
243  mf::LogVerbatim("ClusterMatchTQ") << "Cluster group " << i << ":";
244  for (size_t j = 0; j < matchedclusters[i].size(); ++j) {
245  mf::LogVerbatim("ClusterMatchTQ") << matchedclusters[i][j];
246  }
247  }
248 
249  return matchedclusters;
250  }
code to link reconstructed objects back to the MC truth information
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
TTree * t1
Definition: plottest35.C:26
Planes which measure V.
Definition: geo_types.h:136
double GetXTicksOffset(int p, int t, int c) const
Float_t x1[n_points_granero]
Definition: compare.C:5
Planes which measure Z direction.
Definition: geo_types.h:138
Float_t tmp
Definition: plot.C:35
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
Planes which measure U.
Definition: geo_types.h:135
bool SortByWire(art::Ptr< recob::Hit > const &h1, art::Ptr< recob::Hit > const &h2)
TCanvas * c1
Definition: plotHisto.C:7
TCanvas * c2
Definition: plot_hist.C:75
float bin[41]
Definition: plottest35.C:14
double ConvertTicksToX(double ticks, int p, int t, int c) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
recob::tracking::Plane Plane
Definition: TrackState.h:17
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.

Member Data Documentation

bool cluster::ClusterMatchTQ::fEnableU
private

Definition at line 41 of file ClusterMatchTQ.h.

bool cluster::ClusterMatchTQ::fEnableV
private

Definition at line 42 of file ClusterMatchTQ.h.

bool cluster::ClusterMatchTQ::fEnableZ
private

Definition at line 43 of file ClusterMatchTQ.h.

double cluster::ClusterMatchTQ::fKSCut
private

Definition at line 40 of file ClusterMatchTQ.h.


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