Public Member Functions | Private Attributes | List of all members
cluster::ClusterMatchTQ Class Reference

#include <ClusterMatchTQ.h>

Public Member Functions

 ClusterMatchTQ (fhicl::ParameterSet const &pset)
 
std::vector< std::vector< unsigned int > > MatchedClusters (const detinfo::DetectorClocksData &clockdata, 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 31 of file ClusterMatchTQ.h.

Constructor & Destructor Documentation

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

Definition at line 42 of file ClusterMatchTQ.cxx.

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

Member Function Documentation

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

Definition at line 52 of file ClusterMatchTQ.cxx.

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

Member Data Documentation

bool cluster::ClusterMatchTQ::fEnableU
private

Definition at line 43 of file ClusterMatchTQ.h.

bool cluster::ClusterMatchTQ::fEnableV
private

Definition at line 44 of file ClusterMatchTQ.h.

bool cluster::ClusterMatchTQ::fEnableZ
private

Definition at line 45 of file ClusterMatchTQ.h.

double cluster::ClusterMatchTQ::fKSCut
private

Definition at line 42 of file ClusterMatchTQ.h.


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