15 #include "TPrincipal.h" 21 fTree =
tfs->make<TTree>(
"MatchingVariables",
"MatchingVariables");
36 TVector2
const& centre,
37 TVector2
const& direction,
45 std::map<double, TVector2> hitProjection;
48 for (
auto&
hit : cluster) {
50 hitProjection[direction * pos] = pos;
54 start = hitProjection.begin()->second.Proj(direction) + centre;
55 end = hitProjection.rbegin()->second.Proj(direction) + centre;
61 TVector2
const& centre,
62 TVector2
const& start1,
64 TVector2
const& start2,
65 TVector2
const& end2)
const 70 double clusterOverlap = 0;
73 double s1 = (start1 - centre) * direction;
74 double e1 = (end1 - centre) * direction;
75 double s2 = (start2 - centre) * direction;
76 double e2 = (end2 - centre) * direction;
80 std::cout <<
"s1>e1: " << s1 <<
" and " << e1 << std::endl;
86 std::cout <<
"s1>e1: " << s1 <<
" and " << e1 << std::endl;
93 if ((e1 > s2) && (e2 > s1)) clusterOverlap = std::min((e1 - s2), (e2 - s1));
95 return clusterOverlap;
99 TVector2
const& centre1,
100 TVector2
const& direction2,
101 TVector2
const& centre2)
const 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);
113 double crossingDistance = std::min((centre1 - crossing).Mod(), (centre2 - crossing).Mod());
115 return crossingDistance;
124 double minDistance = 99999.;
127 for (
auto const& hit1 : cluster1) {
128 for (
auto const& hit2 : cluster2) {
133 double distance = (pos1 - pos2).Mod();
135 if (distance < minDistance) minDistance = distance;
143 TVector2
const& start1,
144 TVector2
const& end1,
145 TVector2
const& centre2,
146 TVector2
const& start2,
147 TVector2
const& end2)
const 153 TVector2 parallel = (centre2 - centre1).Unit();
154 TVector2 perpendicular = parallel.Rotate(TMath::Pi() / 2);
157 double s1 = (start1 - centre1) * perpendicular;
158 double e1 = (end1 - centre1) * perpendicular;
159 double s2 = (start2 - centre2) * perpendicular;
160 double e2 = (end2 - centre2) * perpendicular;
163 double projectionStart = std::max(TMath::Abs(s1), TMath::Abs(s2));
164 double projectionEnd = std::max(TMath::Abs(e1), TMath::Abs(e2));
166 double projectionWidth = projectionStart + projectionEnd;
168 return projectionWidth;
197 std::vector<unsigned int> mergedClusters;
199 std::vector<art::PtrVector<recob::Hit>> oldClusters = planeClusters;
202 std::sort(oldClusters.begin(),
205 return a.
size() > b.size();
209 unsigned int nclusters = 0;
210 for (
auto&
cluster : oldClusters)
214 bool mergedAllClusters =
false;
215 while (!mergedAllClusters) {
221 for (
unsigned int initCluster = 0; initCluster < oldClusters.size(); ++initCluster) {
223 std::find(mergedClusters.begin(), mergedClusters.end(), initCluster) !=
224 mergedClusters.end())
226 cluster = oldClusters.
at(initCluster);
227 mergedClusters.push_back(initCluster);
232 bool mergedAllToThisCluster =
false;
233 while (!mergedAllToThisCluster) {
237 for (
unsigned int trialCluster = 0; trialCluster < oldClusters.size(); ++trialCluster) {
240 std::find(mergedClusters.begin(), mergedClusters.end(), trialCluster) !=
241 mergedClusters.end())
245 TPrincipal *pca1 =
new TPrincipal(2,
""), *pca2 =
new TPrincipal(2,
"");
250 TVector2 chargePoint1 = TVector2(0, 0), chargePoint2 = TVector2(0, 0);
251 double totalCharge1 = 0, totalCharge2 = 0;
253 for (
auto& hit1 : cluster) {
258 chargePoint1 += hit1->Integral() * pos;
259 totalCharge1 += hit1->Integral();
261 for (
auto& hit2 : oldClusters.at(trialCluster)) {
266 chargePoint2 += hit2->Integral() * pos;
267 totalCharge2 += hit2->Integral();
270 pca1->MakePrincipals();
271 pca2->MakePrincipals();
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;
286 double length1 = (end1 - start1).Mod();
287 double length2 = (end2 - start2).Mod();
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;
298 angle, crossingDistance, projectedWidth, separation, TMath::Max(length1, length2))) {
300 for (
auto&
hit : oldClusters.at(trialCluster))
301 cluster.push_back(
hit);
303 mergedClusters.push_back(trialCluster);
312 if (nadded == 0) mergedAllToThisCluster =
true;
316 clusters.push_back(cluster);
317 if (mergedClusters.size() == nclusters) mergedAllClusters =
true;
320 return clusters.size();
324 double const crossingDistance,
325 double const projectedWidth,
326 double const separation,
327 double const longLength)
const 332 bool passCrossingDistanceAngle =
false;
333 if (crossingDistance < (-2 + (5 / (1 * TMath::Abs(angle)) - 0))) passCrossingDistanceAngle =
true;
335 bool passSeparationAngle =
false;
336 if (separation < (200 * TMath::Abs(angle) + 40)) passSeparationAngle =
true;
338 bool passProjectedWidth =
false;
340 passProjectedWidth =
true;
342 return passCrossingDistanceAngle and passSeparationAngle and passProjectedWidth;
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].
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.
void FindClusterEndPoints(art::PtrVector< recob::Hit > const &cluster, TVector2 const ¢re, TVector2 const &direction, TVector2 &start, TVector2 &end) const
art::ServiceHandle< art::TFileService const > tfs
CryostatID_t Cryostat
Index of cryostat.
WireID_t Wire
Index of the wire within its plane.
geo::WireID const & WireID() const
Initial tdc tick for hit.
Cluster finding and building.
void reconfigure(fhicl::ParameterSet const &p)
double FindClusterOverlap(TVector2 const &direction, TVector2 const ¢re, 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.
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Signal from induction planes.
T get(std::string const &key) const
double fProjWidthThreshold
reference at(size_type n)
PlaneID_t Plane
Index of the plane within its TPC.
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.
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 ¢re1, TVector2 const &direction2, TVector2 const ¢re2) 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.
double FindProjectedWidth(TVector2 const ¢re1, TVector2 const &start1, TVector2 const &end1, TVector2 const ¢re2, TVector2 const &start2, TVector2 const &end2) const
unsigned int fMinMergeClusterSize
double fMaxMergeSeparation