21 #include "RStarTree/RStarBoundingBox.h" 58 Visitor() : count(0), vResult(), sResult(), ContinueVisiting(true){};
61 vResult.push_back(leaf->leaf);
62 sResult.insert(leaf->leaf);
81 c[0] = (m_bound.edges[0].second + m_bound.edges[0].first) / 2.0;
82 c[1] = (m_bound.edges[1].second + m_bound.edges[1].first) / 2.0;
83 d[0] = (m_bound.edges[0].second - m_bound.edges[0].first) / 2.0;
84 d[1] = (m_bound.edges[1].second - m_bound.edges[1].first) / 2.0;
89 return m_bound.overlaps(node->bound);
95 C[0] = (leaf->bound.edges[0].second + leaf->bound.edges[0].first) / 2.0;
96 C[1] = (leaf->bound.edges[1].second + leaf->bound.edges[1].first) / 2.0;
97 D[0] = (leaf->bound.edges[0].second - leaf->bound.edges[0].first) / 2.0;
98 D[1] = (leaf->bound.edges[1].second - leaf->bound.edges[1].first) / 2.0;
100 for (
int i = 0; i < 2; ++i) {
103 t += ((c[i] - C[i]) * (c[i] - C[i])) / ((d[i] + D[i]) * (d[i] + D[i]));
114 c[0] = (m_bound.edges[0].second - m_bound.edges[0].first) / 2.0;
115 c[1] = (m_bound.edges[1].second - m_bound.edges[1].first) / 2.0;
143 std::vector<unsigned int>& badWireSum)
144 : fBound(b), fEps(), fMaxWidth(maxWidth), fWireDist(wireDist), fBadWireSum(badWireSum)
153 c.edges[0].first = c.edges[0].second = (b.edges[0].first + b.edges[0].second) / 2.0;
154 c.edges[1].first = c.edges[1].second = (b.edges[1].first + b.edges[1].second) / 2.0;
163 double bCenter0 = center(b).edges[0].first;
164 double bCenter1 = center(b).edges[1].first;
165 double tCenter0 = center().edges[0].first;
166 double tCenter1 = center().edges[1].first;
168 double bWidth =
std::abs(b.edges[1].second - b.edges[1].first);
169 double tWidth =
std::abs(fBound.edges[1].second - fBound.edges[1].first);
171 unsigned int wire1 = (
unsigned int)(tCenter0 / fWireDist + 0.5);
172 unsigned int wire2 = (
unsigned int)(bCenter0 / fWireDist + 0.5);
175 if (wire1 < fBadWireSum.size()) wire1 = fBadWireSum.size();
176 if (wire2 < fBadWireSum.size()) wire2 = fBadWireSum.size();
180 unsigned int wirestobridge =
lar::util::absDiff(fBadWireSum[wire1], fBadWireSum[wire2]);
181 double cmtobridge = wirestobridge * fWireDist;
183 double sim =
std::abs(tCenter0 - bCenter0) - cmtobridge;
186 if (
std::abs(tCenter0 - bCenter0) > 1
e-10) {
187 cmtobridge *=
std::abs((tCenter1 - bCenter1) / (tCenter0 - bCenter0));
189 double sim2 =
std::abs(tCenter1 - bCenter1) - cmtobridge;
193 double WFactor = (exp(4.6 * ((tWidth * tWidth) + (bWidth * bWidth)))) * k;
195 if (WFactor < 1.0) WFactor = 1.0;
196 if (WFactor > 6.25) WFactor = 6.25;
199 return (((sim) / (fEps[0] * fEps[0])) + ((sim2) / (fEps[1] * fEps[1] * (WFactor))) <= 1);
205 for (
int i = 0; i < 2; ++i) {
208 if (b.edges[i].first > c.edges[i].second) {
210 n.edges[i].first = n.edges[i].second = b.edges[i].first;
212 else if (b.edges[0].second < c.edges[0].first) {
214 n.edges[i].first = n.edges[i].second = b.edges[i].second;
218 n.edges[i].first = n.edges[i].second = c.edges[i].first;
222 n.edges[1].first -= fMaxWidth / 2.0;
223 n.edges[1].second += fMaxWidth / 2.0;
229 if (fBound.overlaps(node->bound))
return true;
233 return isNear(nearestPoint(node->bound));
235 bool operator()(
const RTree::Leaf*
const leaf)
const {
return isNear(leaf->bound); }
248 fEps = p.
get<
double>(
"eps");
249 fEps2 = p.
get<
double>(
"epstwo");
250 fMinPts = p.
get<
int>(
"minPts");
251 fClusterMethod = p.
get<
int>(
"Method");
252 fDistanceMetric = p.
get<
int>(
"Metric");
259 std::set<uint32_t> badChannels,
260 const std::vector<geo::WireID>& wireids)
262 if (wireids.size() && wireids.size() != allhits.size()) {
263 throw cet::exception(
"DBScanAlg") <<
"allhits size = " << allhits.size()
264 <<
" wireids size = " << wireids.size() <<
" do not match\n";
268 fpointId_to_clusterId.clear();
277 fBadChannels = badChannels;
281 fRTree.Remove(RTree::AcceptAny(), RTree::RemoveLeaf());
292 for (
auto const& plane : wireReadoutGeom.Iterate<
geo::PlaneGeo>(tpcid))
293 fWirePitch.push_back(plane.WirePitch());
296 if (fClusterMethod) {
297 auto const nchannels = wireReadoutGeom.Nchannels();
298 fBadWireSum.resize(nchannels);
299 unsigned int count = 0;
300 for (
unsigned int i = 0; i < fBadWireSum.size(); ++i) {
301 count += fBadChannels.count(i);
302 fBadWireSum[i] = count;
309 for (
unsigned int j = 0; j < allhits.size(); ++j) {
311 std::vector<double> p(dims);
316 p[0] = (allhits[j]->
WireID().Wire) * fWirePitch[allhits[j]->
WireID().Plane];
318 p[0] = (wireids[j].Wire) * fWirePitch[allhits[j]->
WireID().Plane];
319 p[1] = allhits[j]->PeakTime() * tickToDist;
320 p[2] = 2. * allhits[j]->RMS() * tickToDist;
323 if (p[2] > fMaxWidth) fMaxWidth = p[2];
327 if (fClusterMethod) {
329 dbsPoint pp(p[0], p[1], 0.0, p[2] / 2.0);
330 fRTree.Insert(j, pp.
bounds());
336 fpointId_to_clusterId.resize(fps.size(),
kNO_CLUSTER);
337 fnoise.resize(fps.size(),
false);
338 fvisited.resize(fps.size(),
false);
340 if (fClusterMethod) {
342 mf::LogInfo(
"DBscan") <<
"InitScan: hits RTree loaded with " << visitor.
count <<
" items.";
344 mf::LogInfo(
"DBscan") <<
"InitScan: hits vector size is " << fps.size();
362 double wire_dist = fWirePitch[0];
365 (
unsigned int)(v1[0] / wire_dist + 0.5);
366 unsigned int wire2 = (
unsigned int)(v2[0] / wire_dist + 0.5);
367 int wirestobridge = 0;
370 unsigned int wire = wire1;
375 for (
unsigned int i = wire1; i < wire2; i++) {
376 if (fBadChannels.find(i) != fBadChannels.end()) wirestobridge++;
379 double cmtobridge = wirestobridge * wire_dist;
381 return ((
std::abs(v2[0] - v1[0]) - cmtobridge) *
382 (
std::abs(v2[0] - v1[0]) - cmtobridge));
387 const std::vector<double> v2)
395 double wire_dist = fWirePitch[0];
398 (
unsigned int)(v1[0] / wire_dist + 0.5);
399 unsigned int wire2 = (
unsigned int)(v2[0] / wire_dist + 0.5);
400 int wirestobridge = 0;
403 unsigned int wire = wire1;
408 for (
unsigned int i = wire1; i < wire2; i++) {
409 if (fBadChannels.find(i) != fBadChannels.end()) wirestobridge++;
412 double cmtobridge = wirestobridge * wire_dist;
415 cmtobridge *=
std::abs((v2[1] - v1[1]) / (v2[0] - v1[0]));
420 return ((
std::abs(v2[1] - v1[1]) - cmtobridge) *
421 (
std::abs(v2[1] - v1[1]) - cmtobridge));
426 const std::vector<double> v2)
437 double WFactor = (exp(4.6 * ((v1[2] * v1[2]) + (v2[2] * v2[2])))) * k;
458 std::vector<unsigned int> ne;
460 for (
int unsigned j = 0; j < fsim.size(); j++) {
462 (((fsim[pid][j]) / (threshold * threshold)) +
463 ((fsim2[pid][j]) / (threshold2 * threshold2 * (fsim3[pid][j])))) < 1) {
474 int size = fps.size();
475 fsim.resize(size, std::vector<double>(size));
476 for (
int i = 0; i <
size; i++) {
477 for (
int j = i + 1; j <
size; j++) {
478 fsim[j][i] = fsim[i][j] = getSimilarity(fps[i], fps[j]);
486 int size = fps.size();
487 fsim2.resize(size, std::vector<double>(size));
488 for (
int i = 0; i <
size; i++) {
489 for (
int j = i + 1; j <
size; j++) {
490 fsim2[j][i] = fsim2[i][j] = getSimilarity2(fps[i], fps[j]);
498 int size = fps.size();
499 fsim3.resize(size, std::vector<double>(size));
501 for (
int i = 0; i <
size; i++) {
502 for (
int j = i + 1; j <
size; j++) {
503 fsim3[j][i] = fsim3[i][j] = getWidthFactor(fps[i], fps[j]);
514 switch (fClusterMethod) {
515 case 2:
return run_dbscan_cluster();
516 case 1:
return run_FN_cluster();
519 computeSimilarity2();
520 computeWidthFactor();
521 return run_FN_naive_cluster();
532 unsigned int cid = 0;
534 for (
size_t pid = 0;
pid < fps.size();
pid++) {
537 if (ExpandCluster(
pid, cid)) { cid++; }
544 fclusters.resize(cid);
545 for (
size_t y = 0;
y < fpointId_to_clusterId.size(); ++
y) {
554 unsigned int c = fpointId_to_clusterId[
y];
556 mf::LogWarning(
"DBscan") <<
"Point in cluster " << c <<
" when only " << cid
557 <<
" clusters wer found [0-" << cid - 1 <<
"]";
559 fclusters[c].push_back(
y);
562 mf::LogInfo(
"DBscan") <<
"DWM (R*-tree): Found " << cid <<
" clusters...";
563 for (
unsigned int c = 0; c < cid; ++c) {
565 <<
"Cluster " << c <<
":\t" << fclusters[c].size();
568 <<
"...and " << noise <<
" noise points.";
599 std::vector<unsigned int>& v = visitor.
vResult;
603 v.erase(std::remove(v.begin(), v.end(), point), v.end());
612 std::set<unsigned int>
seeds = RegionQuery(point);
615 if (seeds.size() < fMinPts) {
621 fpointId_to_clusterId[point] = clusterID;
623 fpointId_to_clusterId[*itr] = clusterID;
626 while (!seeds.empty()) {
627 unsigned int currentP = *(seeds.begin());
628 std::set<unsigned int> result = RegionQuery(currentP);
630 if (result.size() >= fMinPts) {
632 unsigned int resultP = *itr;
634 if (fpointId_to_clusterId[resultP] ==
kNO_CLUSTER ||
636 if (fpointId_to_clusterId[resultP] ==
kNO_CLUSTER) { seeds.insert(resultP); }
637 fpointId_to_clusterId[resultP] = clusterID;
641 seeds.erase(currentP);
656 unsigned int cid = 0;
658 for (
size_t pid = 0;
pid < fps.size();
pid++) {
660 if (!fvisited[
pid]) {
662 fvisited[pid] =
true;
664 std::vector<unsigned int> ne = RegionQuery_vector(pid);
667 if (ne.size() < fMinPts) { fnoise[pid] =
true; }
671 std::vector<unsigned int> c;
674 fpointId_to_clusterId[pid] = cid;
676 for (
size_t i = 0; i < ne.size(); ++i) {
677 unsigned int nPid = ne[i];
680 if (!fvisited[nPid]) {
681 fvisited[nPid] =
true;
683 std::vector<unsigned int> ne1 = RegionQuery_vector(nPid);
685 if (ne1.size() >= fMinPts) {
689 for (
size_t i = 0; i < ne1.size(); ++i) {
691 ne.push_back(ne1[i]);
699 fpointId_to_clusterId[nPid] = cid;
703 fclusters.push_back(c);
713 for (
size_t y = 0;
y < fpointId_to_clusterId.size(); ++
y) {
716 mf::LogInfo(
"DBscan") <<
"FindNeighbors (R*-tree): Found " << cid <<
" clusters...";
717 for (
unsigned int c = 0; c < cid; ++c) {
719 <<
"Cluster " << c <<
":\t" << fclusters[c].size();
722 <<
"...and " << noise <<
" noise points.";
733 unsigned int cid = 0;
735 for (
size_t pid = 0;
pid < fps.size(); ++
pid) {
737 if (!fvisited[
pid]) {
739 fvisited[pid] =
true;
741 std::vector<unsigned int> ne = findNeighbors(pid, fEps, fEps2);
744 if (ne.size() < fMinPts) { fnoise[pid] =
true; }
748 std::vector<unsigned int> c;
751 fpointId_to_clusterId[pid] = cid;
753 for (
size_t i = 0; i < ne.size(); ++i) {
754 unsigned int nPid = ne[i];
757 if (!fvisited[nPid]) {
758 fvisited[nPid] =
true;
760 std::vector<unsigned int> ne1 = findNeighbors(nPid, fEps, fEps2);
762 if (ne1.size() >= fMinPts) {
766 for (
unsigned int i = 0; i < ne1.size(); i++) {
768 ne.push_back(ne1[i]);
776 fpointId_to_clusterId[nPid] = cid;
780 fclusters.push_back(c);
789 for (
size_t y = 0;
y < fpointId_to_clusterId.size(); ++
y) {
792 mf::LogInfo(
"DBscan") <<
"FindNeighbors (naive): Found " << cid <<
" clusters...";
793 for (
unsigned int c = 0; c < cid; ++c) {
795 <<
"Cluster " << c <<
":\t" << fclusters[c].size() <<
" points";
798 <<
"...and " << noise <<
" noise points.";
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
double getSimilarity(const std::vector< double > v1, const std::vector< double > v2)
const unsigned int kNO_CLUSTER
Functions to help with numbers.
Utilities related to art service access.
void computeWidthFactor()
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
bool operator()(const RTree::Leaf *const leaf) const
Declaration of signal hit object.
AcceptFindNeighbors(const BoundingBox &b, double eps, double eps2, double maxWidth, double wireDist, std::vector< unsigned int > &badWireSum)
double Temperature() const
In kelvin.
std::vector< unsigned int > & fBadWireSum
static const BoundingBox EmptyBoundingBox
constexpr auto abs(T v)
Returns the absolute value of the argument.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Cluster finding and building.
void run_dbscan_cluster()
std::set< unsigned int > RegionQuery(unsigned int point)
void run_FN_naive_cluster()
bool isNear(const BoundingBox &b) const
double Efield(unsigned int planegap=0) const
kV/cm
bool operator()(const RTree::Leaf *const leaf) const
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
const unsigned int kNOISE_CLUSTER
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
RTree::BoundingBox BoundingBox
std::set< unsigned int > sResult
const BoundingBox & fBound
void operator()(const RTree::Leaf *const leaf)
BoundingBox bounds() const
T get(std::string const &key) const
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
const BoundingBox & m_bound
BoundingBox center(const BoundingBox &b) const
bool operator()(const RTree::Node *const node) const
std::vector< unsigned int > vResult
double DriftVelocity(double efield=0., double temperature=0.) const
cm/us
The data type to uniquely identify a TPC.
DBScanAlg(fhicl::ParameterSet const &pset)
double getWidthFactor(const std::vector< double > v1, const std::vector< double > v2)
bool operator()(const RTree::Node *const node) const
void computeSimilarity2()
std::vector< unsigned int > findNeighbors(unsigned int pid, double threshold, double threshold2)
AcceptEllipse(const BoundingBox &b, double r1, double r2)
std::vector< TrajPoint > seeds
std::vector< unsigned int > RegionQuery_vector(unsigned int point)
void InitScan(const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &allhits, std::set< uint32_t > badChannels, const std::vector< geo::WireID > &wireids=std::vector< geo::WireID >())
bool ExpandCluster(unsigned int point, unsigned int clusterID)
Contains all timing reference information for the detector.
constexpr auto absDiff(A const &a, B const &b)
Returns the absolute value of the difference between two values.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double getSimilarity2(const std::vector< double > v1, const std::vector< double > v2)
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
BoundingBox center() const
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
const bool ContinueVisiting
BoundingBox nearestPoint(const BoundingBox &b) const