LArSoft  v10_04_05
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 41 of file ClusterMatchTQ.cxx.

References fhicl::ParameterSet::get().

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

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 50 of file ClusterMatchTQ.cxx.

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

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

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

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: