bug fixes PLS IGNORE THIS ENTIRE BRANCH IM JUST TOO LAZY TO FORK

Pr/PrVeloUT/src/PrVeloUtOpt.cpp

@@ -100,9 +100,10 @@ PrVeloUtOutput PrVeloUtOpt::operator()(const std::vector<TrackObj>& veloTracks,
 	std::array<simd_vector<uint32_t>, 4> gatheredOffests;
 
 	for (size_t i = 0; i < 4; ++i) {
-		auto [h, o] = GatherHits(hits, filteredStates, layerStates[i], 1.4f);
+		auto [h, o] = GatherHitsSimd(hits, filteredStates, layerStates[i], 1.4f);
 		gatheredHits[i] = std::move(h);
 		gatheredOffests[i] = std::move(o);
+		//auto [_ignore1, _ingore2] = GatherHitsSimd(hits, filteredStates, layerStates[i], 1.4f);
 	}
 
 	auto trackCandidates = FindLines(gatheredHits, gatheredOffests, filteredStates);
@@ -114,6 +115,7 @@ PrVeloUtOutput PrVeloUtOpt::operator()(const std::vector<TrackObj>& veloTracks,
 	//}
 	auto assembledTracks = AssembleTracks(bestCandidates, trackCandidates);
 
+
 	// Ugly hacky code to produce acceptable output.
 	// TODO: refactor to proper format and function.
 	const auto numTracks = assembledTracks.trackIndex.size();
@@ -252,11 +254,69 @@ simd_vector<float> GetHitTolerances(const States& states, const ExtrapolatedStat
 	return result;
 }
 
-
+//*
 std::pair<GatheredHits, simd_vector<uint32_t>> PrVeloUtOpt::GatherHits(const UtHitSpacePartition& hits, const States& states, const ExtrapolatedStates& layerStates, float tolerance) {
 	const UtHits& hitsContainer = hits.GetPartitionedHits();
 	const size_t numLayerStates = layerStates.size();
 
+	size_t writeIndex = 0;
+
+	size_t sumBinSizes;
+	StateBinLocations binLocations = GetBinLocations(hits, layerStates, sumBinSizes);
+
+	GatheredHits gatheredHits;
+	simd_vector<uint32_t> offsets;
+	gatheredHits.resize(sumBinSizes);
+	offsets.resize(states.size());
+
+	for (size_t i = 0; i < numLayerStates; ++i) {
+		const auto index = binLocations.index[i];
+		auto count = binLocations.count[i];
+		float trackX = layerStates.x[i];
+		float trackY = layerStates.y[i];
+		float trackTx = states.tx[i];
+		float trackTy = states.ty[i];
+		float trackZ = layerStates.z;
+
+		constexpr float adaptiveToleranceConstant = 83.f; // This is a magic number obtained from old UT. Roughly 1/500*distToMomentum*minPt, I don't know what it means or if it's needed.
+		float adaptiveTolerance = adaptiveToleranceConstant * std::sqrt(trackTx * trackTx + trackTy * trackTy);
+		adaptiveTolerance = std::clamp(adaptiveTolerance, 0.2f, 16.f); // Limits from old UT as well.
+
+		while (count) {
+			auto totalIndex = index + count - 1;
+
+			const float hitZ0 = hitsContainer.z0[totalIndex];
+			const float diffZ = hitZ0 - trackZ;
+			const float trackAccurateX = trackX + diffZ * trackTx;
+			const float trackAccurateY = trackY + diffZ * trackTy;
+
+			const float hitX0 = hitsContainer.x0[totalIndex];
+			const float hitDxdy = hitsContainer.dxdy[totalIndex];
+			const float hitAccurateX = hitX0 + hitDxdy * trackAccurateY;
+
+			if (std::abs(hitAccurateX - trackAccurateX) < adaptiveTolerance) {
+				gatheredHits.x[writeIndex] = hitAccurateX;
+				gatheredHits.z[writeIndex] = hitZ0;
+				gatheredHits.index[writeIndex] = totalIndex;
+				++writeIndex;
+			}
+
+			--count;
+		}
+
+		offsets[i] = writeIndex;
+	}
+
+	gatheredHits.resize(writeIndex);
+	
+	return { std::move(gatheredHits), std::move(offsets) };
+}
+
+//*/
+std::pair<GatheredHits, simd_vector<uint32_t>> PrVeloUtOpt::GatherHitsSimd(const UtHitSpacePartition& hits, const States& states, const ExtrapolatedStates& layerStates, float tolerance) {
+	const UtHits& hitsContainer = hits.GetPartitionedHits();
+	const size_t numLayerStates = layerStates.size();
+
 	size_t sumBinSizes;
 	StateBinLocations binLocations = GetBinLocations(hits, layerStates, sumBinSizes);
 	simd_vector<float> tolerances = GetHitTolerances(states, layerStates, tolerance);
@@ -281,7 +341,7 @@ std::pair<GatheredHits, simd_vector<uint32_t>> PrVeloUtOpt::GatherHits(const UtH
 
 		for (int i = 0; i < count; i += SimdWidth) {
 			auto totalIndex = index + i;
-			uint32_v totalIndexSimd = uint32_v(index) + uint32_v::IndexesFromZero();
+			uint32_v totalIndexSimd = uint32_v(totalIndex) + uint32_v::IndexesFromZero();
 
 			float_v hitX0;
 			float_v hitZ0;
@@ -295,7 +355,7 @@ std::pair<GatheredHits, simd_vector<uint32_t>> PrVeloUtOpt::GatherHits(const UtH
 			const float_v trackAccurateY = trackY + diffZ * trackTy;
 			const float_v hitAccurateX = hitX0 + hitDxdy * trackAccurateY;
 
-			const float_mask_v mask = Vc::abs(hitAccurateX - trackAccurateX) < float_v(adaptiveTolerance);
+			const float_mask_v mask = Vc::abs(hitAccurateX - trackAccurateX) < float_v(adaptiveTolerance) && float_mask_v(totalIndexSimd < uint32_v(index + count));
 			int maski = mask.toInt();
 #ifdef _MSC_VER
 			int popcount = __popcnt(maski);
@@ -318,9 +378,10 @@ std::pair<GatheredHits, simd_vector<uint32_t>> PrVeloUtOpt::GatherHits(const UtH
 	}
 
 	gatheredHits.resize(writeIndex);
-
+	
 	return { std::move(gatheredHits), std::move(offsets) };
 }
+//*/
 
 
 
@@ -359,7 +420,7 @@ uint32_t FindLinesForTrack(const std::array<GatheredHits, 4>& gatheredHits,
 	const size_t writeIndexStart = writeIndex;
 	uint32_t minHitsFound = 3;
 	uint32_t passedCount = 0;
-	constexpr float maxDiff = 0.4f;
+	constexpr float maxDiff = 0.5f;
 
 	for (size_t ibase1 = beginIndices[layerConfig.base1]; ibase1 < endIndices[layerConfig.base1]; ++ibase1) {
 		const float xbase1 = gatheredHits[layerConfig.base1].x[ibase1];
@@ -372,8 +433,10 @@ uint32_t FindLinesForTrack(const std::array<GatheredHits, 4>& gatheredHits,
 			const float dxdz = (xbase1 - xbase2) / (zbase1 - zbase2);
 			const float x0 = xbase2 - dxdz * zbase2;
 
+			//std::cout << xbase1 << " " << xbase2 << " -> ";
 			auto [minDiffExtra1, idxExtra1] = ScanExtraHits(gatheredHits[layerConfig.extra1], beginIndices[layerConfig.extra1], endIndices[layerConfig.extra1], x0, dxdz, maxDiff);
 			auto [minDiffExtra2, idxExtra2] = ScanExtraHits(gatheredHits[layerConfig.extra2], beginIndices[layerConfig.extra2], endIndices[layerConfig.extra2], x0, dxdz, maxDiff);
+			//std::cout << minDiffExtra1 << " " << minDiffExtra2 << "   ";
 
 			uint32_t numHitsFound = 2 + (minDiffExtra1 < maxDiff) + (minDiffExtra2 < maxDiff);
 			if (numHitsFound == 4 && minHitsFound == 3) {
@@ -407,6 +470,7 @@ uint32_t FindLinesForTrack(const std::array<GatheredHits, 4>& gatheredHits,
 		}
 	}
 
+	//std::cout << "((" << passedCount << "))" << std::endl;
 	return passedCount;
 }
 
@@ -425,8 +489,8 @@ inline std::array<int, 2> GetLineCombinatorics(const std::array<uint32_t, 4>& be
 		endIndices[3] - beginIndices[3],
 	};
 
-	const int combinatorics1 = counts[0] * counts[2] * (counts[1] + counts[3]);
-	const int combinatorics2 = counts[1] * counts[3] * (counts[0] + counts[2]);
+	const int combinatorics1 = counts[lineSearchConfigs[0].base1] * counts[lineSearchConfigs[0].base2] * (counts[lineSearchConfigs[0].extra1] + counts[lineSearchConfigs[0].extra2]);
+	const int combinatorics2 = counts[lineSearchConfigs[1].base1] * counts[lineSearchConfigs[1].base2] * (counts[lineSearchConfigs[1].extra1] + counts[lineSearchConfigs[1].extra2]);
 
 	return { combinatorics1, combinatorics2 };
 }
@@ -434,7 +498,7 @@ inline std::array<int, 2> GetLineCombinatorics(const std::array<uint32_t, 4>& be
 
 inline std::pair<const SearchLayers&, const SearchLayers&> GetSearchConfig(int c1, int c2) {
 	using Ret = std::pair<const SearchLayers&, const SearchLayers&>;
-	return c1 > c2 ? Ret{ lineSearchConfigs[0], lineSearchConfigs[1] } : Ret{ lineSearchConfigs[0], lineSearchConfigs[1] };
+	return c1 > c2 ? Ret{ lineSearchConfigs[0], lineSearchConfigs[1] } : Ret{ lineSearchConfigs[1], lineSearchConfigs[0] };
 }
 
 
@@ -455,7 +519,7 @@ UtTrackCandidates FindLines(const std::array<GatheredHits, 4>& gatheredHits, con
 
 		const auto& [combinatorics1, combinatorics2] = GetLineCombinatorics(hitIndices, currentIndices);
 		const auto& [config1, config2] = GetSearchConfig(combinatorics1, combinatorics2);
-		const auto& [maxComb, minComb] = std::minmax(combinatorics1, combinatorics2);
+		const auto& [minComb, maxComb] = std::minmax(combinatorics1, combinatorics2);
 
 		if (trackCandidates.size() < writeIndex + maxComb) {
 			trackCandidates.resize(writeIndex + maxComb);
@@ -642,7 +706,7 @@ FittedTrackProperties FitTrackCandidates(const UtTrackCandidates& candidates) {
 
 		const FittedTrackParameters fitParams = {
 			fitSystem.b[0],
-			fitSystem.b[1],
+			fitSystem.b[1] * zUnitScaling,
 			fitSystem.b[2],
 		};
 

Pr/PrVeloUT/src/PrVeloUtOpt.h

@@ -32,6 +32,7 @@ private:
 	static auto FilterStates(const std::vector<TrackObj>& veloTracks, const States& veloEndStates) -> std::tuple<std::vector<TrackObj>, States>;
 	static StateBinLocations GetBinLocations(const UtHitSpacePartition& hits, const ExtrapolatedStates& layerStates, size_t& totalNumBins);
 	static std::pair<GatheredHits, simd_vector<uint32_t>> GatherHits(const UtHitSpacePartition& hits, const States& states, const ExtrapolatedStates& layerStates, float tolerance);
+	static std::pair<GatheredHits, simd_vector<uint32_t>> GatherHitsSimd(const UtHitSpacePartition& hits, const States& states, const ExtrapolatedStates& layerStates, float tolerance);
 };
 
 

Pr/PrVeloUT/src/UtHitSpacePartition.h

@@ -10,7 +10,7 @@ namespace opt {
 
 
 constexpr size_t numRows = 14;
-constexpr size_t numBins = 128;
+constexpr size_t numBins = 64;
 constexpr size_t numLayers = 4;
 constexpr float detectorBottom = -1386.f / 2.0f;
 constexpr float detectorTop = +1386.f / 2.0f;