10 #include <tracking/trackFindingCDC/findlets/minimal/TrackCombiner.h>
12 #include <tracking/trackFindingCDC/eventdata/tracks/CDCTrack.h>
13 #include <tracking/trackFindingCDC/eventdata/segments/CDCSegment3D.h>
14 #include <tracking/trackFindingCDC/eventdata/hits/CDCRLWireHit.h>
15 #include <tracking/trackFindingCDC/eventdata/hits/CDCWireHit.h>
17 #include <tracking/trackFindingCDC/numerics/Index.h>
19 #include <tracking/trackFindingCDC/utilities/WeightedRelation.h>
20 #include <tracking/trackFindingCDC/utilities/Functional.h>
21 #include <tracking/trackFindingCDC/utilities/Range.h>
22 #include <tracking/trackFindingCDC/utilities/StringManipulation.h>
24 #include <framework/core/ModuleParamList.templateDetails.h>
30 using namespace TrackFindingCDC;
33 std::array<int, ISuperLayerUtil::c_N> getNHitsByISuperLayer(
const CDCTrack& track)
35 std::array<int, ISuperLayerUtil::c_N> nHitsBySLayer = {0};
36 for (
const CDCRecoHit3D& hit : track) {
37 nHitsBySLayer[hit.getISuperLayer()]++;
42 CDCTrack condense(
const TrackFindingCDC::Path<const CDCSegment3D>& segmentPath)
45 for (
const CDCSegment3D* segment : segmentPath) {
46 for (
const CDCRecoHit3D& recoHit3D : *segment) {
47 result.push_back(recoHit3D);
50 if (segmentPath.size()) {
51 result.setStartTrajectory3D(segmentPath.front()->getTrajectory3D());
56 std::vector<CDCRLWireHit> getRLWireHitSegment(
const CDCSegment3D& segment3D)
58 std::vector<CDCRLWireHit> result;
59 result.reserve(segment3D.size());
60 for (
const CDCRecoHit3D& recoHit3D : segment3D) {
62 result.push_back(recoHit3D.getRLWireHit());
67 std::vector<std::pair<std::pair<Index, Index>, CDCRLWireHit>> getCommonRLWireHits(
const std::vector<CDCRLWireHit>& rlWireHits1,
68 const std::vector<CDCRLWireHit>& rlWireHits2)
70 std::vector<std::pair<std::pair<Index, Index>, CDCRLWireHit>> result;
72 for (
const CDCRLWireHit& rlWireHit : rlWireHits1) {
74 auto itRLWireHits2 = std::find(rlWireHits2.begin(), rlWireHits2.end(), rlWireHit);
75 if (itRLWireHits2 == rlWireHits2.end())
continue;
76 Index index2 = std::distance(rlWireHits2.begin(), itRLWireHits2);
78 result.push_back({{index1, index2}, rlWireHit});
88 return "Combines two sets of tracks to one final set by merging tracks that have large overlaps.";
92 const std::string& prefix)
94 moduleParamList->
addParameter(prefixed(prefix,
"identifyCommonSegments"),
96 "Activate the identification of common segments",
101 const std::vector<CDCTrack>& secondInputTracks,
102 std::vector<CDCTrack>& outputTracks)
104 using InTracks = std::array<const CDCTrack*, 2>;
107 std::map<const CDCWireHit*, InTracks> inTracksByWireHit;
108 for (
const CDCTrack& track : inputTracks) {
110 const CDCWireHit& wireHit = recoHit3D.getWireHit();
112 if (inTracksByWireHit.count(&wireHit) == 0) inTracksByWireHit[&wireHit] = {
nullptr,
nullptr};
113 inTracksByWireHit[&wireHit][0] = &track;
120 for (
const CDCTrack& track : secondInputTracks) {
122 const CDCWireHit& wireHit = recoHit3D.getWireHit();
124 if (inTracksByWireHit.count(&wireHit) == 0) inTracksByWireHit[&wireHit] = {
nullptr,
nullptr};
125 inTracksByWireHit[&wireHit][1] = &track;
133 std::vector<std::pair<std::pair<int, size_t>,
const CDCTrack*> > rankedTracks;
134 for (
const CDCTrack& track : inputTracks) {
135 std::array<int, ISuperLayerUtil::c_N> nHitsBySLayer = getNHitsByISuperLayer(track);
136 int nSL = std::count_if(nHitsBySLayer.begin(), nHitsBySLayer.end(),
Id() > 0);
137 rankedTracks.push_back({{nSL, track.size()}, &track});
139 for (
const CDCTrack& track : secondInputTracks) {
140 std::array<int, ISuperLayerUtil::c_N> nHitsBySLayer = getNHitsByISuperLayer(track);
141 int nSL = std::count_if(nHitsBySLayer.begin(), nHitsBySLayer.end(),
Id() > 0);
142 rankedTracks.push_back({{nSL, track.size()}, &track});
147 std::deque<CDCSegment3D> segments;
150 std::vector<WeightedRelation<const CDCSegment3D>> segmentRelations;
153 std::map<const CDCSegment3D*, InTracks> inTracksBySegment;
156 for (
const std::pair<std::pair<int, size_t>,
const CDCTrack*> rankedTrack : rankedTracks) {
158 const CDCTrack* track = rankedTrack.second;
159 std::vector<std::pair<InTracks, CDCSegment3D> > segmentsInTrack;
163 ISuperLayer lastISuperLayer = -1;
164 InTracks lastTracksForHit{{
nullptr,
nullptr}};
167 std::array<const CDCTrack*, 2> tracksForWireHit = inTracksByWireHit[&recoHit3D.getWireHit()];
168 if (not currentSegment or iSuperLayer != lastISuperLayer or tracksForWireHit != lastTracksForHit) {
170 segmentsInTrack.push_back({tracksForWireHit,
CDCSegment3D()});
171 currentSegment = &segmentsInTrack.back().second;
172 currentSegment->setTrajectory3D(track->getStartTrajectory3D());
174 currentSegment->push_back(recoHit3D);
175 lastISuperLayer = iSuperLayer;
176 lastTracksForHit = tracksForWireHit;
180 auto absorbsSegment = [](std::pair<InTracks, CDCSegment3D>& segment1,
181 std::pair<InTracks, CDCSegment3D>& segment2) {
183 if (segment1.second.getISuperLayer() != segment2.second.getISuperLayer())
return false;
184 if (segment1.second.size() < 3) {
185 segment1.first = segment2.first;
188 if (segment1.second.size() < 3 or segment2.second.size() < 3 or
189 segment1.first == segment2.first) {
190 segment1.second.insert(segment1.second.end(),
191 segment2.second.begin(),
192 segment2.second.end());
197 auto itSegmentToDrop = std::unique(segmentsInTrack.begin(), segmentsInTrack.end(), absorbsSegment);
198 segmentsInTrack.erase(itSegmentToDrop, segmentsInTrack.end());
201 for (std::pair<InTracks, CDCSegment3D>& segment : segmentsInTrack) {
202 nHits += segment.second.size();
204 B2ASSERT(
"Expected the number of hits to be the same", nHits == track->size());
208 for (std::pair<InTracks, CDCSegment3D>& segmentInTrack : segmentsInTrack) {
209 const InTracks& inTracks = segmentInTrack.first;
211 segment->setCellWeight(segment.size());
212 segments.push_back(std::move(segment));
213 inTracksBySegment[&segments.back()] = inTracks;
214 if (lastSegment !=
nullptr) {
215 segmentRelations.push_back({lastSegment, 0, &segments.back()});
217 lastSegment = &segments.back();
222 std::sort(segmentRelations.begin(), segmentRelations.end());
226 std::vector<std::pair<const CDCSegment3D*, const CDCSegment3D*>> segmentContainmentRelation;
227 for (
auto itSegment1 = segments.begin(); itSegment1 != segments.end(); ++itSegment1) {
229 std::vector<CDCRLWireHit> rlWireHits1 = getRLWireHitSegment(segment1);
231 InTracks inTracks1 = inTracksBySegment[&segment1];
232 for (
auto itSegment2 = itSegment1 + 1; itSegment2 != segments.end(); ++itSegment2) {
236 if (&segment1 == &segment2)
continue;
239 if (iSL1 != iSL2)
continue;
241 InTracks inTracks2 = inTracksBySegment[&segment2];
242 if (inTracks1 != inTracks2)
continue;
245 std::vector<CDCRLWireHit> rlWireHits2 = getRLWireHitSegment(segment2);
246 std::vector<std::pair<std::pair<int, int>,
CDCRLWireHit>> commonRLWireHits =
247 getCommonRLWireHits(rlWireHits1, rlWireHits2);
251 bool contains = commonRLWireHits.size() > 2 and
252 commonRLWireHits.size() > rlWireHits1.size() * 0.8 and
253 commonRLWireHits.size() > rlWireHits2.size() * 0.8;
254 if (not contains)
continue;
257 double n = commonRLWireHits.size();
263 for (
const auto& commonRLWireHit : commonRLWireHits) {
264 const std::pair<int, int>& indices = commonRLWireHit.first;
265 prod11 += indices.first * indices.first;
266 prod12 += indices.first * indices.second;
267 prod22 += indices.second * indices.second;
269 sum1 += indices.first;
270 sum2 += indices.second;
272 double var1 = (prod11 - sum1 * sum1 / n);
273 double var2 = (prod22 - sum2 * sum2 / n);
274 double cov12 = (prod12 - sum1 * sum2 / n);
275 double cor = cov12 / std::sqrt(var1 * var2);
276 if (not(cor > 0.75))
continue;
278 segmentContainmentRelation.push_back({&segment1, &segment2});
283 std::vector<WeightedRelation<const CDCSegment3D>> additionalSegmentRelations;
284 for (std::pair<const CDCSegment3D*, const CDCSegment3D*>& rel : segmentContainmentRelation) {
286 asRange(std::equal_range(segmentRelations.begin(), segmentRelations.end(), rel.first))) {
287 additionalSegmentRelations.push_back({rel.second, 0, rel1.getTo()});
291 std::equal_range(segmentRelations.begin(), segmentRelations.end(), rel.second))) {
292 additionalSegmentRelations.push_back({rel.first, 0, rel2.getFrom()});
295 segmentRelations.insert(segmentRelations.end(),
296 additionalSegmentRelations.begin(),
297 additionalSegmentRelations.end());
298 std::sort(segmentRelations.begin(), segmentRelations.end());
303 std::vector<const CDCSegment3D*> segmentPtrs = as_pointers<const CDCSegment3D>(segments);
306 std::vector<TrackFindingCDC::Path<const CDCSegment3D>> segmentPaths;
310 outputTracks.clear();
311 for (
const std::vector<const CDCSegment3D*>& segmentPath : segmentPaths) {
313 if (segmentPath.size() < 2)
continue;
314 outputTracks.push_back(condense(segmentPath));