1 #include <tracking/trackFindingCDC/processing/TrackQualityTools.h>
2 #include <tracking/trackFindingCDC/eventdata/tracks/CDCTrack.h>
3 #include <tracking/trackFindingCDC/eventdata/hits/CDCWireHit.h>
4 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory3D.h>
5 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory2D.h>
7 #include <tracking/trackFindingCDC/topology/CDCWire.h>
8 #include <tracking/trackFindingCDC/utilities/ReversedRange.h>
13 using namespace TrackFindingCDC;
19 if (recoHit.getStereoKind() == EStereoKind::c_Axial) {
22 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
34 if (arcLength2DOfApogee < 0) {
35 arcLength2DOfApogee += 2 * TMath::Pi() * radius;
42 double currentArcLength2D = recoHit.getArcLength2D();
43 if (currentArcLength2D < 0) B2INFO(
"Below 0");
44 if (currentArcLength2D > arcLength2DOfApogee) {
45 splittedCDCTrack.push_back(recoHit);
50 if (splittedCDCTrack.size() > 0) {
51 tracks.push_back(splittedCDCTrack);
57 if (track.empty())
return;
62 const Vector3D startPosition(0, 0, 0);
70 unsigned int numberOfPositiveHits = 0;
72 const double currentArcLength = trajectory2D.
calcArcLength2D(recoHit.getRecoPos2D());
73 if (currentArcLength > 0) {
74 numberOfPositiveHits++;
77 const bool reverseTrajectory = 2 * numberOfPositiveHits < track.size();
79 if (reverseTrajectory) {
86 Vector2D recoPos2D = recoHit.getRLWireHit().reconstruct3D(trajectory2D).xy();
89 if (arcLength2D < 0) {
90 arcLength2D += arcLength2DPeriod;
92 recoHit.setArcLength2D(arcLength2D);
93 recoHit.getWireHit().getAutomatonCell().unsetAssignedFlag();
94 recoHit.getWireHit().getAutomatonCell().setTakenFlag();
98 track.sortByArcLength2D();
101 Vector3D frontPosition = track.front().getRLWireHit().reconstruct3D(trajectory2D);
102 double arcLengthOffset = trajectory3D.
setLocalOrigin(frontPosition);
103 track.setStartTrajectory3D(trajectory3D);
105 recoHit.shiftArcLength2D(-arcLengthOffset);
109 Vector3D backPosition = track.back().getRLWireHit().reconstruct3D(trajectory2D);
110 double backArcLength2D = trajectory3D.
setLocalOrigin(backPosition);
111 if (backArcLength2D < 0) {
114 track.setEndTrajectory3D(trajectory3D);
120 const CDCTrajectory2D& trajectory2D = track.getStartTrajectory3D().getTrajectory2D();
124 if (trajectory2D.
isCurler(outerCylindricalRFactor)) {
130 double arcLength2DOfExitWithFactor = trajectory2D.
calcArcLength2D(outerExitWithFactor);
131 if (arcLength2DOfExitWithFactor < 0) {
132 arcLength2DOfExitWithFactor += 2 * TMath::Pi() * radius;
134 bool removeAfterThis =
false;
137 if (removeAfterThis) {
138 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
142 double currentArcLength2D = recoHit.getArcLength2D();
143 if (currentArcLength2D < 0) {
144 currentArcLength2D += 2 * TMath::Pi() * radius;
147 if (currentArcLength2D > arcLength2DOfExitWithFactor) {
148 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
149 removeAfterThis =
true;
156 ILayer lastLayer = -1;
159 std::vector<std::vector<const CDCRecoHit3D*>> trackletList;
160 trackletList.reserve(3);
162 std::vector<const CDCRecoHit3D*>* currentTracklet =
nullptr;
165 if (currentTracklet ==
nullptr) {
166 trackletList.emplace_back();
167 currentTracklet = &(trackletList.back());
170 const ILayer currentLayer = recoHit.getWire().getICLayer();
171 const Vector2D& currentPosition = recoHit.getRecoPos2D();
172 if (lastLayer != -1) {
173 const ILayer delta = currentLayer - lastLayer;
174 const double distance = (currentPosition - lastWirePosition).norm();
175 if (abs(delta) > 4 or distance > 50) {
176 trackletList.emplace_back();
177 currentTracklet = &(trackletList.back());
181 lastWirePosition = currentPosition;
182 lastLayer = currentLayer;
184 currentTracklet->push_back(&recoHit);
187 if (trackletList.size() > 1) {
188 for (
const std::vector<const CDCRecoHit3D*>& tracklet : trackletList) {
189 if (tracklet.size() < 5) {
191 recoHit->getWireHit().getAutomatonCell().setAssignedFlag();
204 if (std::isnan(radius)) {
208 std::vector<std::vector<const CDCRecoHit3D*>> trackletList;
210 double lastArcLength2D = NAN;
212 std::vector<const CDCRecoHit3D*>* currentTracklet =
nullptr;
215 if (currentTracklet ==
nullptr) {
216 trackletList.emplace_back();
217 currentTracklet = &(trackletList.back());
220 const double currentArcLength2D = recoHit.getArcLength2D();
221 if (not std::isnan(lastArcLength2D)) {
222 const double delta = (currentArcLength2D - lastArcLength2D);
223 if (std::fabs(delta) > m_maximumArcLength2DDistance) {
224 trackletList.emplace_back();
225 currentTracklet = &(trackletList.back());
229 lastArcLength2D = currentArcLength2D;
231 currentTracklet->push_back(&recoHit);
234 if (trackletList.size() > 1) {
236 while (trackletList.size() > 0) {
237 if (trackletList.back().size() > 4) {
241 for (
const CDCRecoHit3D* recoHit : trackletList.back()) {
242 recoHit->getWireHit().getAutomatonCell().setAssignedFlag();
245 trackletList.pop_back();
248 std::reverse(trackletList.begin(), trackletList.end());
250 while (trackletList.size() > 0) {
251 if (trackletList.back().size() > 4) {
255 for (
const CDCRecoHit3D* recoHit : trackletList.back()) {
256 recoHit->getWireHit().getAutomatonCell().setAssignedFlag();
259 trackletList.pop_back();
266 const bool deleteTrack = track.size() < minimalHits;
270 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
277 double lastAngle = NAN;
278 bool removeAfterThis =
false;
280 for (
const CDCRecoHit3D& recoHit : reversedRange(track)) {
281 if (removeAfterThis) {
282 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
286 const double currentAngle = recoHit.getRecoPos2D().phi();
287 if (not std::isnan(lastAngle)) {
288 const double delta = currentAngle - lastAngle;
289 const double normalizedDelta = std::min(TVector2::Phi_0_2pi(delta), TVector2::Phi_0_2pi(-delta));
290 if (fabs(normalizedDelta) > maximalAngle) {
291 removeAfterThis =
true;
292 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
296 lastAngle = currentAngle;
303 ISuperLayer lastLayer = -1;
304 bool deleteTrack =
true;
307 const ISuperLayer currentLayer = recoHit.getISuperLayer();
308 if (lastLayer != -1 and lastLayer != currentLayer) {
313 lastLayer = currentLayer;
318 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
334 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
345 if (std::isnan(radius)) {
349 double lastArcLength2D = NAN;
350 bool removeAfterThis =
false;
353 if (removeAfterThis) {
354 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
358 const double currentArcLength2D = recoHit.getArcLength2D();
359 if (not std::isnan(lastArcLength2D)) {
360 const double delta = (currentArcLength2D - lastArcLength2D) / radius;
361 if (fabs(delta) > m_maximumArcLength2DDistance) {
362 removeAfterThis =
true;
363 recoHit.getWireHit().getAutomatonCell().setAssignedFlag();
368 lastArcLength2D = currentArcLength2D;