Main algorithm.
100{
101 using InTracks = std::array<const CDCTrack*, 2>;
102
103
104 std::map<const CDCWireHit*, InTracks> inTracksByWireHit;
105 for (const CDCTrack& track : inputTracks) {
106 for (const CDCRecoHit3D& recoHit3D : track) {
107 const CDCWireHit& wireHit = recoHit3D.getWireHit();
108 if (inTracksByWireHit.count(&wireHit) == 0) inTracksByWireHit[&wireHit] = {nullptr, nullptr};
109 inTracksByWireHit[&wireHit][0] = &track;
110
113 }
114 }
115
116 for (const CDCTrack& track : secondInputTracks) {
117 for (const CDCRecoHit3D& recoHit3D : track) {
118 const CDCWireHit& wireHit = recoHit3D.getWireHit();
119 if (inTracksByWireHit.count(&wireHit) == 0) inTracksByWireHit[&wireHit] = {nullptr, nullptr};
120 inTracksByWireHit[&wireHit][1] = &track;
121
124 }
125 }
126
127
128 std::vector<std::pair<std::pair<int, size_t>, const CDCTrack*> > rankedTracks;
129 for (const CDCTrack& track : inputTracks) {
130 std::array<int, ISuperLayerUtil::c_N> nHitsBySLayer = getNHitsByISuperLayer(track);
131 int nSL = std::count_if(nHitsBySLayer.begin(), nHitsBySLayer.end(), Id() > 0);
132 rankedTracks.push_back({{nSL, track.size()}, &track});
133 }
134 for (const CDCTrack& track : secondInputTracks) {
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});
138 }
139 std::sort(rankedTracks.begin(), rankedTracks.end(), GreaterOf<First>());
140
141
142 std::deque<CDCSegment3D> segments;
143
144
145 std::vector<WeightedRelation<const CDCSegment3D>> segmentRelations;
146
147
148 std::map<const CDCSegment3D*, InTracks> inTracksBySegment;
149
150
151 for (const std::pair<std::pair<int, size_t>, const CDCTrack*>& rankedTrack : rankedTracks) {
152
153 const CDCTrack* track = rankedTrack.second;
154 std::vector<std::pair<InTracks, CDCSegment3D> > segmentsInTrack;
155
156
157 CDCSegment3D* currentSegment = nullptr;
158 ISuperLayer lastISuperLayer = -1;
159 InTracks lastTracksForHit{{nullptr, nullptr}};
160 for (const CDCRecoHit3D& recoHit3D : *track) {
161 ISuperLayer iSuperLayer = recoHit3D.getISuperLayer();
162 std::array<const CDCTrack*, 2> tracksForWireHit = inTracksByWireHit[&recoHit3D.getWireHit()];
163 if (not currentSegment or iSuperLayer != lastISuperLayer or tracksForWireHit != lastTracksForHit) {
164
165 segmentsInTrack.push_back({tracksForWireHit, CDCSegment3D()});
166 currentSegment = &segmentsInTrack.back().second;
167 currentSegment->setTrajectory3D(track->getStartTrajectory3D());
168 }
169 currentSegment->push_back(recoHit3D);
170 lastISuperLayer = iSuperLayer;
171 lastTracksForHit = tracksForWireHit;
172 }
173
174
175 auto absorbsSegment = [](std::pair<InTracks, CDCSegment3D>& segment1,
176 std::pair<InTracks, CDCSegment3D>& segment2) {
177
178 if (segment1.second.getISuperLayer() != segment2.second.getISuperLayer()) return false;
179 if (segment1.second.size() < 3) {
180 segment1.first = segment2.first;
181 }
182
183 if (segment1.second.size() < 3 or segment2.second.size() < 3 or
184 segment1.first == segment2.first) {
185 segment1.second.insert(segment1.second.end(),
186 segment2.second.begin(),
187 segment2.second.end());
188 return true;
189 }
190 return false;
191 };
192 auto itSegmentToDrop = std::unique(segmentsInTrack.begin(), segmentsInTrack.end(), absorbsSegment);
193 segmentsInTrack.erase(itSegmentToDrop, segmentsInTrack.end());
194
195 size_t nHits = 0;
196 for (std::pair<InTracks, CDCSegment3D>& segment : segmentsInTrack) {
197 nHits += segment.second.size();
198 }
199 B2ASSERT("Expected the number of hits to be the same", nHits == track->size());
200
201
202 const CDCSegment3D* lastSegment = nullptr;
203 for (std::pair<InTracks, CDCSegment3D>& segmentInTrack : segmentsInTrack) {
204 const InTracks& inTracks = segmentInTrack.first;
205 const CDCSegment3D& segment = segmentInTrack.second;
206 segment->setCellWeight(segment.size());
207 segments.push_back(std::move(segment));
208 inTracksBySegment[&segments.back()] = inTracks;
209 if (lastSegment != nullptr) {
210 segmentRelations.push_back({lastSegment, 0, &segments.back()});
211 }
212 lastSegment = &segments.back();
213 }
214 }
215
216
217 std::sort(segmentRelations.begin(), segmentRelations.end());
218
219
221 std::vector<std::pair<const CDCSegment3D*, const CDCSegment3D*>> segmentContainmentRelation;
222 for (auto itSegment1 = segments.begin(); itSegment1 != segments.end(); ++itSegment1) {
223 const CDCSegment3D& segment1 = *itSegment1;
224 std::vector<CDCRLWireHit> rlWireHits1 = getRLWireHitSegment(segment1);
226 InTracks inTracks1 = inTracksBySegment[&segment1];
227 for (auto itSegment2 = itSegment1 + 1; itSegment2 != segments.end(); ++itSegment2) {
228 const CDCSegment3D& segment2 = *itSegment2;
229
230
231 if (&segment1 == &segment2) continue;
232
234 if (iSL1 != iSL2) continue;
235
236 InTracks inTracks2 = inTracksBySegment[&segment2];
237 if (inTracks1 != inTracks2) continue;
238
239
240 std::vector<CDCRLWireHit> rlWireHits2 = getRLWireHitSegment(segment2);
241 std::vector<std::pair<std::pair<int, int>, CDCRLWireHit>> commonRLWireHits =
242 getCommonRLWireHits(rlWireHits1, rlWireHits2);
243
244
245
246 bool contains = commonRLWireHits.size() > 2 and
247 commonRLWireHits.size() > rlWireHits1.size() * 0.8 and
248 commonRLWireHits.size() > rlWireHits2.size() * 0.8;
249 if (not contains) continue;
250
251
252 double n = commonRLWireHits.size();
253 double prod11 = 0;
254 double prod12 = 0;
255 double prod22 = 0;
256 double sum1 = 0;
257 double sum2 = 0;
258 for (const auto& commonRLWireHit : commonRLWireHits) {
259 const std::pair<int, int>& indices = commonRLWireHit.first;
260 prod11 += indices.first * indices.first;
261 prod12 += indices.first * indices.second;
262 prod22 += indices.second * indices.second;
263
264 sum1 += indices.first;
265 sum2 += indices.second;
266 }
267 double var1 = (prod11 - sum1 * sum1 / n);
268 double var2 = (prod22 - sum2 * sum2 / n);
269 double cov12 = (prod12 - sum1 * sum2 / n);
270 double cor = cov12 / std::sqrt(var1 * var2);
271 if (not(cor > 0.75)) continue;
272
273 segmentContainmentRelation.push_back({&segment1, &segment2});
274 }
275 }
276
277
278 std::vector<WeightedRelation<const CDCSegment3D>> additionalSegmentRelations;
279 for (std::pair<const CDCSegment3D*, const CDCSegment3D*>& rel : segmentContainmentRelation) {
280 for (WeightedRelation<const CDCSegment3D> rel1 :
281 asRange(std::equal_range(segmentRelations.begin(), segmentRelations.end(), rel.first))) {
282 additionalSegmentRelations.push_back({rel.second, 0, rel1.getTo()});
283 }
284
285 for (WeightedRelation<const CDCSegment3D> rel2 : asRange(
286 std::equal_range(segmentRelations.begin(), segmentRelations.end(), rel.second))) {
287 additionalSegmentRelations.push_back({rel.first, 0, rel2.getFrom()});
288 }
289 }
290 segmentRelations.insert(segmentRelations.end(),
291 additionalSegmentRelations.begin(),
292 additionalSegmentRelations.end());
293 std::sort(segmentRelations.begin(), segmentRelations.end());
294 }
295
296
297
298 std::vector<const CDCSegment3D*> segmentPtrs = as_pointers<const CDCSegment3D>(segments);
299
300
301 std::vector<TrackFindingCDC::Path<const CDCSegment3D>> segmentPaths;
303
304
305 outputTracks.clear();
306 for (const std::vector<const CDCSegment3D*>& segmentPath : segmentPaths) {
307
308 if (segmentPath.size() < 2) continue;
309 outputTracks.push_back(condense(segmentPath));
310 }
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332}
void unsetMaskedFlag()
Resets the masked flag to false.
void unsetTemporaryFlags()
Resets the assigned, start and cycle marker flag.
ISuperLayer getISuperLayer() const
Returns the common super layer id of all stored tracking hits.
MultipassCellularPathFinder< const CDCSegment3D > m_cellularPathFinder
Instance of the cellular automaton path finder.
bool m_param_identifyCommonSegments
Parameter : Activate the identification of common segments.