Extensive Code Duplication

[GNiendorf]

To simplify the code, we should remove the extensive use of code duplication (copy-paste).

Here is a gist showing some of the current duplication, automatically generated with PMD. Not all of these are meaningful duplications, but it catches many obvious repeats:

https://gist.github.com/GNiendorf/fd4fe952f38149b5bfd1269a7039370e

This was generated with the following:

1. Download PMD

wget https://github.com/pmd/pmd/releases/download/pmd_releases%2F7.5.0/pmd-dist-7.5.0-bin.zip

2. Extract the PMD archive

unzip pmd-dist-7.5.0-bin.zip -d pmd-dist

3. Navigate to the PMD bin directory

cd pmd-dist/pmd-bin-7.5.0/bin

4. Create a list of .cc and .h files in specific directories

find
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/interface/alpaka
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/interface
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/src/alpaka
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/src
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/standalone/bin
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/standalone/code/core
CMSSW_14_1_0_pre3/src/RecoTracker/LSTCore/standalone/efficiency/src
-type f ( -name ".cc" -o -name ".h" ) > file_list.txt

5. Run CPD to generate the duplication report

./pmd cpd --minimum-tokens 100 --file-list file_list.txt --language cpp --format text > duplication_report.txt

[GNiendorf]

Some examples below from the report:

Constants repeated in Segment.h

cmssw/RecoTracker/LSTCore/src/alpaka/Segment.h

Lines 205 to 213 in ff8446c

	static constexpr float miniDeltaTilted[3] = {0.26f, 0.26f, 0.26f};
	static constexpr float miniDeltaFlat[6] = {0.26f, 0.16f, 0.16f, 0.18f, 0.18f, 0.18f};
	static constexpr float miniDeltaLooseTilted[3] = {0.4f, 0.4f, 0.4f};
	static constexpr float miniDeltaEndcap[5][15] = {
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f}};

cmssw/RecoTracker/LSTCore/src/alpaka/Segment.h

Lines 235 to 243 in ff8446c

	static constexpr float miniDeltaTilted[3] = {0.26f, 0.26f, 0.26f};
	static constexpr float miniDeltaFlat[6] = {0.26f, 0.16f, 0.16f, 0.18f, 0.18f, 0.18f};
	static constexpr float miniDeltaLooseTilted[3] = {0.4f, 0.4f, 0.4f};
	static constexpr float miniDeltaEndcap[5][15] = {
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f},
	{0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.4f, 0.18f, /*10*/ 0.18f, 0.18f, 0.18f, 0.18f, 0.18f}};

[GNiendorf]

Duplicate functions in Quintuplet.h and PixelQuintuplet.h

cmssw/RecoTracker/LSTCore/src/alpaka/Quintuplet.h

Lines 1015 to 1099 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeSigmasForRegression(TAcc const& acc,
	lst::Modules const& modulesInGPU,
	const uint16_t* lowerModuleIndices,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	unsigned int nPoints = 5,
	bool anchorHits = true) {
	/*
	Bool anchorHits required to deal with a weird edge case wherein
	the hits ultimately used in the regression are anchor hits, but the
	lower modules need not all be Pixel Modules (in case of PS). Similarly,
	when we compute the chi squared for the non-anchor hits, the "partner module"
	need not always be a PS strip module, but all non-anchor hits sit on strip
	modules.
	*/
	ModuleType moduleType;
	short moduleSubdet, moduleSide;
	float inv1 = kWidthPS / kWidth2S;
	float inv2 = kPixelPSZpitch / kWidth2S;
	float inv3 = kStripPSZpitch / kWidth2S;
	for (size_t i = 0; i < nPoints; i++) {
	moduleType = modulesInGPU.moduleType[lowerModuleIndices[i]];
	moduleSubdet = modulesInGPU.subdets[lowerModuleIndices[i]];
	moduleSide = modulesInGPU.sides[lowerModuleIndices[i]];
	const float& drdz = modulesInGPU.drdzs[lowerModuleIndices[i]];
	slopes[i] = modulesInGPU.dxdys[lowerModuleIndices[i]];
	//category 1 - barrel PS flat
	if (moduleSubdet == Barrel and moduleType == PS and moduleSide == Center) {
	delta1[i] = inv1;
	delta2[i] = inv1;
	slopes[i] = -999.f;
	isFlat[i] = true;
	}
	//category 2 - barrel 2S
	else if (moduleSubdet == Barrel and moduleType == TwoS) {
	delta1[i] = 1.f;
	delta2[i] = 1.f;
	slopes[i] = -999.f;
	isFlat[i] = true;
	}
	//category 3 - barrel PS tilted
	else if (moduleSubdet == Barrel and moduleType == PS and moduleSide != Center) {
	delta1[i] = inv1;
	isFlat[i] = false;
	if (anchorHits) {
	delta2[i] = (inv2 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
	} else {
	delta2[i] = (inv3 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
	}
	}
	//category 4 - endcap PS
	else if (moduleSubdet == Endcap and moduleType == PS) {
	delta1[i] = inv1;
	isFlat[i] = false;
	/*
	despite the type of the module layer of the lower module index,
	all anchor hits are on the pixel side and all non-anchor hits are
	on the strip side!
	*/
	if (anchorHits) {
	delta2[i] = inv2;
	} else {
	delta2[i] = inv3;
	}
	}
	//category 5 - endcap 2S
	else if (moduleSubdet == Endcap and moduleType == TwoS) {
	delta1[i] = 1.f;
	delta2[i] = 500.f * inv1;
	isFlat[i] = false;
	} else {
	#ifdef WARNINGS
	printf("ERROR!!!!! I SHOULDN'T BE HERE!!!! subdet = %d, type = %d, side = %d\n",
	moduleSubdet,
	moduleType,
	moduleSide);
	#endif
	}
	}
	};

cmssw/RecoTracker/LSTCore/src/alpaka/PixelQuintuplet.h

Lines 430 to 514 in ff8446c

	template <typename TAcc>
	ALPAKA_FN_ACC ALPAKA_FN_INLINE void computeSigmasForRegression_pT5(TAcc const& acc,
	lst::Modules const& modulesInGPU,
	const uint16_t* lowerModuleIndices,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	unsigned int nPoints = 5,
	bool anchorHits = true) {
	/*
	bool anchorHits required to deal with a weird edge case wherein
	the hits ultimately used in the regression are anchor hits, but the
	lower modules need not all be Pixel Modules (in case of PS). Similarly,
	when we compute the chi squared for the non-anchor hits, the "partner module"
	need not always be a PS strip module, but all non-anchor hits sit on strip
	modules.
	*/
	ModuleType moduleType;
	short moduleSubdet, moduleSide;
	float inv1 = kWidthPS / kWidth2S;
	float inv2 = kPixelPSZpitch / kWidth2S;
	float inv3 = kStripPSZpitch / kWidth2S;
	for (size_t i = 0; i < nPoints; i++) {
	moduleType = modulesInGPU.moduleType[lowerModuleIndices[i]];
	moduleSubdet = modulesInGPU.subdets[lowerModuleIndices[i]];
	moduleSide = modulesInGPU.sides[lowerModuleIndices[i]];
	const float& drdz = modulesInGPU.drdzs[lowerModuleIndices[i]];
	slopes[i] = modulesInGPU.dxdys[lowerModuleIndices[i]];
	//category 1 - barrel PS flat
	if (moduleSubdet == Barrel and moduleType == PS and moduleSide == Center) {
	delta1[i] = inv1;
	delta2[i] = inv1;
	slopes[i] = -999.f;
	isFlat[i] = true;
	}
	//category 2 - barrel 2S
	else if (moduleSubdet == Barrel and moduleType == TwoS) {
	delta1[i] = 1.f;
	delta2[i] = 1.f;
	slopes[i] = -999.f;
	isFlat[i] = true;
	}
	//category 3 - barrel PS tilted
	else if (moduleSubdet == Barrel and moduleType == PS and moduleSide != Center) {
	delta1[i] = inv1;
	isFlat[i] = false;
	if (anchorHits) {
	delta2[i] = (inv2 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
	} else {
	delta2[i] = (inv3 * drdz / alpaka::math::sqrt(acc, 1 + drdz * drdz));
	}
	}
	//category 4 - endcap PS
	else if (moduleSubdet == Endcap and moduleType == PS) {
	delta1[i] = inv1;
	isFlat[i] = false;
	/*
	despite the type of the module layer of the lower module index,
	all anchor hits are on the pixel side and all non-anchor hits are
	on the strip side!
	*/
	if (anchorHits) {
	delta2[i] = inv2;
	} else {
	delta2[i] = inv3;
	}
	}
	//category 5 - endcap 2S
	else if (moduleSubdet == Endcap and moduleType == TwoS) {
	delta1[i] = 1.f;
	delta2[i] = 500.f * inv1;
	isFlat[i] = false;
	}
	#ifdef WARNINGS
	else {
	printf("ERROR!!!!! I SHOULDN'T BE HERE!!!! subdet = %d, type = %d, side = %d\n",
	moduleSubdet,
	moduleType,
	moduleSide);
	}
	#endif
	}
	};

[GNiendorf]

Triple duplication of chi-squared function. @YonsiG do you know if there is supposed to be any difference between these three, or can they be merged?

cmssw/RecoTracker/LSTCore/src/alpaka/PixelQuintuplet.h

Lines 385 to 428 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquaredpT5(TAcc const& acc,
	unsigned int nPoints,
	float* xs,
	float* ys,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	float g,
	float f,
	float radius) {
	/*
	Given values of (g, f, radius) and a set of points (and its uncertainties) compute chi squared
	*/
	float c = g * g + f * f - radius * radius;
	float chiSquared = 0.f;
	float absArctanSlope, angleM, xPrime, yPrime, sigma2;
	for (size_t i = 0; i < nPoints; i++) {
	absArctanSlope = ((slopes[i] != lst::lst_INF) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
	: 0.5f * float(M_PI));
	if (xs[i] > 0 and ys[i] > 0) {
	angleM = 0.5f * float(M_PI) - absArctanSlope;
	} else if (xs[i] < 0 and ys[i] > 0) {
	angleM = absArctanSlope + 0.5f * float(M_PI);
	} else if (xs[i] < 0 and ys[i] < 0) {
	angleM = -(absArctanSlope + 0.5f * float(M_PI));
	} else if (xs[i] > 0 and ys[i] < 0) {
	angleM = -(0.5f * float(M_PI) - absArctanSlope);
	} else {
	angleM = 0;
	}
	if (not isFlat[i]) {
	xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
	yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
	} else {
	xPrime = xs[i];
	yPrime = ys[i];
	}
	sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
	chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
	(xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / (sigma2);
	}
	return chiSquared;
	};

cmssw/RecoTracker/LSTCore/src/alpaka/Quintuplet.h

Lines 1200 to 1243 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquared(TAcc const& acc,
	unsigned int nPoints,
	float* xs,
	float* ys,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	float g,
	float f,
	float radius) {
	// given values of (g, f, radius) and a set of points (and its uncertainties)
	// compute chi squared
	float c = g * g + f * f - radius * radius;
	float chiSquared = 0.f;
	float absArctanSlope, angleM, xPrime, yPrime, sigma2;
	for (size_t i = 0; i < nPoints; i++) {
	absArctanSlope = ((slopes[i] != lst::lst_INF) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
	: 0.5f * float(M_PI));
	if (xs[i] > 0 and ys[i] > 0) {
	angleM = 0.5f * float(M_PI) - absArctanSlope;
	} else if (xs[i] < 0 and ys[i] > 0) {
	angleM = absArctanSlope + 0.5f * float(M_PI);
	} else if (xs[i] < 0 and ys[i] < 0) {
	angleM = -(absArctanSlope + 0.5f * float(M_PI));
	} else if (xs[i] > 0 and ys[i] < 0) {
	angleM = -(0.5f * float(M_PI) - absArctanSlope);
	} else {
	angleM = 0;
	}
	if (not isFlat[i]) {
	xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
	yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
	} else {
	xPrime = xs[i];
	yPrime = ys[i];
	}
	sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
	chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
	(xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / sigma2;
	}
	return chiSquared;
	};

cmssw/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h

Lines 325 to 368 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquaredpT3(TAcc const& acc,
	unsigned int nPoints,
	float* xs,
	float* ys,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	float g,
	float f,
	float radius) {
	//given values of (g, f, radius) and a set of points (and its uncertainties)
	//compute chi squared
	float c = g * g + f * f - radius * radius;
	float chiSquared = 0.f;
	float absArctanSlope, angleM, xPrime, yPrime, sigma2;
	for (size_t i = 0; i < nPoints; i++) {
	absArctanSlope = ((slopes[i] != lst::lst_INF) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
	: 0.5f * float(M_PI));
	if (xs[i] > 0 and ys[i] > 0) {
	angleM = 0.5f * float(M_PI) - absArctanSlope;
	} else if (xs[i] < 0 and ys[i] > 0) {
	angleM = absArctanSlope + 0.5f * float(M_PI);
	} else if (xs[i] < 0 and ys[i] < 0) {
	angleM = -(absArctanSlope + 0.5f * float(M_PI));
	} else if (xs[i] > 0 and ys[i] < 0) {
	angleM = -(0.5f * float(M_PI) - absArctanSlope);
	} else {
	angleM = 0;
	}
	if (not isFlat[i]) {
	xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
	yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
	} else {
	xPrime = xs[i];
	yPrime = ys[i];
	}
	sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
	chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
	(xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / sigma2;
	}
	return chiSquared;
	};

[GNiendorf]

runDeltaBetaIterations duplication.

cmssw/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h

Lines 1074 to 1178 in ff8446c

	template <typename TAcc>
	ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterationspT3(TAcc const& acc,
	float& betaIn,
	float& betaOut,
	float betaAv,
	float& pt_beta,
	float sdIn_dr,
	float sdOut_dr,
	float dr,
	float lIn) {
	if (lIn == 0) {
	betaOut += alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaOut);
	return;
	}
	if (betaIn * betaOut > 0.f and
	(alpaka::math::abs(acc, pt_beta) < 4.f * lst::kPt_betaMax or
	(lIn >= 11 and alpaka::math::abs(acc, pt_beta) <
	8.f * lst::kPt_betaMax))) //and the pt_beta is well-defined; less strict for endcap-endcap
	{
	const float betaInUpd =
	betaIn + alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdIn_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	const float betaOutUpd =
	betaOut + alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaOut); //FIXME: need a faster version
	betaAv = 0.5f * (betaInUpd + betaOutUpd);
	//1st update
	const float pt_beta_inv =
	1.f / alpaka::math::abs(acc, dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv)); //get a better pt estimate
	betaIn += alpaka::math::copysign(
	acc,
	alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * lst::k2Rinv1GeVf * pt_beta_inv, lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	betaOut += alpaka::math::copysign(
	acc,
	alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * lst::k2Rinv1GeVf * pt_beta_inv, lst::kSinAlphaMax)),
	betaOut); //FIXME: need a faster version
	//update the av and pt
	betaAv = 0.5f * (betaIn + betaOut);
	//2nd update
	pt_beta = dr * lst::k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
	} else if (lIn < 11 && alpaka::math::abs(acc, betaOut) < 0.2f * alpaka::math::abs(acc, betaIn) &&
	alpaka::math::abs(acc, pt_beta) < 12.f * lst::kPt_betaMax) //use betaIn sign as ref
	{
	const float pt_betaIn = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaIn);
	const float betaInUpd =
	betaIn + alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdIn_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	const float betaOutUpd =
	betaOut +
	alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	betaAv = (alpaka::math::abs(acc, betaOut) > 0.2f * alpaka::math::abs(acc, betaIn))
	? (0.5f * (betaInUpd + betaOutUpd))
	: betaInUpd;
	//1st update
	pt_beta = dr * lst::k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
	betaIn += alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(acc, sdIn_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	betaOut += alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	//update the av and pt
	betaAv = 0.5f * (betaIn + betaOut);
	//2nd update
	pt_beta = dr * lst::k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
	}
	};

cmssw/RecoTracker/LSTCore/src/alpaka/Quintuplet.h

Lines 1245 to 1349 in ff8446c

	template <typename TAcc>
	ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterationsT5(TAcc const& acc,
	float& betaIn,
	float& betaOut,
	float betaAv,
	float& pt_beta,
	float sdIn_dr,
	float sdOut_dr,
	float dr,
	float lIn) {
	if (lIn == 0) {
	betaOut += alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaOut);
	return;
	}
	if (betaIn * betaOut > 0.f and
	(alpaka::math::abs(acc, pt_beta) < 4.f * lst::kPt_betaMax or
	(lIn >= 11 and alpaka::math::abs(acc, pt_beta) <
	8.f * lst::kPt_betaMax))) //and the pt_beta is well-defined; less strict for endcap-endcap
	{
	const float betaInUpd =
	betaIn + alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdIn_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	const float betaOutUpd =
	betaOut + alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaOut); //FIXME: need a faster version
	betaAv = 0.5f * (betaInUpd + betaOutUpd);
	//1st update
	const float pt_beta_inv =
	1.f / alpaka::math::abs(acc, dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaAv)); //get a better pt estimate
	betaIn += alpaka::math::copysign(
	acc,
	alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * lst::k2Rinv1GeVf * pt_beta_inv, lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	betaOut += alpaka::math::copysign(
	acc,
	alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * lst::k2Rinv1GeVf * pt_beta_inv, lst::kSinAlphaMax)),
	betaOut); //FIXME: need a faster version
	//update the av and pt
	betaAv = 0.5f * (betaIn + betaOut);
	//2nd update
	pt_beta = dr * lst::k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
	} else if (lIn < 11 && alpaka::math::abs(acc, betaOut) < 0.2f * alpaka::math::abs(acc, betaIn) &&
	alpaka::math::abs(acc, pt_beta) < 12.f * lst::kPt_betaMax) //use betaIn sign as ref
	{
	const float pt_betaIn = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaIn);
	const float betaInUpd =
	betaIn + alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdIn_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	const float betaOutUpd =
	betaOut +
	alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(
	acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	betaAv = (alpaka::math::abs(acc, betaOut) > 0.2f * alpaka::math::abs(acc, betaIn))
	? (0.5f * (betaInUpd + betaOutUpd))
	: betaInUpd;
	//1st update
	pt_beta = dr * lst::k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
	betaIn += alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(acc, sdIn_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	betaOut += alpaka::math::copysign(
	acc,
	alpaka::math::asin(
	acc,
	alpaka::math::min(acc, sdOut_dr * lst::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), lst::kSinAlphaMax)),
	betaIn); //FIXME: need a faster version
	//update the av and pt
	betaAv = 0.5f * (betaIn + betaOut);
	//2nd update
	pt_beta = dr * lst::k2Rinv1GeVf / alpaka::math::sin(acc, betaAv); //get a better pt estimate
	}
	};

[GNiendorf]

See PR #97

cmssw/RecoTracker/LSTCore/src/alpaka/PixelQuintuplet.h

Lines 852 to 878 in ff8446c

	// calculation is detailed documented here https://indico.cern.ch/event/1185895/contributions/4982756/attachments/2526561/4345805/helix%20pT3%20summarize.pdf
	float diffr, diffz;
	float p = alpaka::math::sqrt(acc, Px * Px + Py * Py + Pz * Pz);
	float rou = a / p;
	if (moduleSubdet == lst::Endcap) {
	float s = (zsi - z1) * p / Pz;
	float x = x1 + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
	float y = y1 + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
	diffr = alpaka::math::abs(acc, rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
	}
	if (moduleSubdet == lst::Barrel) {
	float paraA = r1 * r1 + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y1 * Px - x1 * Py) / a - rtsi * rtsi;
	float paraB = 2 * (x1 * Px + y1 * Py) / a;
	float paraC = 2 * (y1 * Px - x1 * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
	float A = paraB * paraB + paraC * paraC;
	float B = 2 * paraA * paraB;
	float C = paraA * paraA - paraC * paraC;
	float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
	float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
	float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z1;
	float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z1;
	float diffz1 = alpaka::math::abs(acc, solz1 - zsi) * 100;
	float diffz2 = alpaka::math::abs(acc, solz2 - zsi) * 100;
	diffz = alpaka::math::min(acc, diffz1, diffz2);
	}

cmssw/RecoTracker/LSTCore/src/alpaka/Quintuplet.h

Lines 532 to 565 in ff8446c

	// calculation is copied from PixelTriplet.cc lst::computePT3RZChiSquared
	float diffr = 0, diffz = 0;
	float rou = a / p;
	// for endcap
	float s = (zsi - z_init) * p / Pz;
	float x = x_init + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
	float y = y_init + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
	diffr = (rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
	// for barrel
	if (layeri <= 6) {
	float paraA =
	rt_init * rt_init + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y_init * Px - x_init * Py) / a - rtsi * rtsi;
	float paraB = 2 * (x_init * Px + y_init * Py) / a;
	float paraC = 2 * (y_init * Px - x_init * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
	float A = paraB * paraB + paraC * paraC;
	float B = 2 * paraA * paraB;
	float C = paraA * paraA - paraC * paraC;
	float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
	float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
	float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z_init;
	float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z_init;
	float diffz1 = (solz1 - zsi) * 100;
	float diffz2 = (solz2 - zsi) * 100;
	// Alpaka : Needs to be moved over
	if (alpaka::math::isnan(acc, diffz1))
	diffz = diffz2;
	else if (alpaka::math::isnan(acc, diffz2))
	diffz = diffz1;
	else {
	diffz = (alpaka::math::abs(acc, diffz1) < alpaka::math::abs(acc, diffz2)) ? diffz1 : diffz2;
	}
	}

cmssw/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h

Lines 726 to 754 in ff8446c

	// calculation is detailed documented here https://indico.cern.ch/event/1185895/contributions/4982756/attachments/2526561/4345805/helix%20pT3%20summarize.pdf
	float diffr, diffz;
	float p = alpaka::math::sqrt(acc, Px * Px + Py * Py + Pz * Pz);
	float rou = a / p;
	if (moduleSubdet == lst::Endcap) {
	float s = (zsi - z1) * p / Pz;
	float x = x1 + Px / a * alpaka::math::sin(acc, rou * s) - Py / a * (1 - alpaka::math::cos(acc, rou * s));
	float y = y1 + Py / a * alpaka::math::sin(acc, rou * s) + Px / a * (1 - alpaka::math::cos(acc, rou * s));
	diffr = alpaka::math::abs(acc, rtsi - alpaka::math::sqrt(acc, x * x + y * y)) * 100;
	residual = diffr;
	}
	if (moduleSubdet == lst::Barrel) {
	float paraA = r1 * r1 + 2 * (Px * Px + Py * Py) / (a * a) + 2 * (y1 * Px - x1 * Py) / a - rtsi * rtsi;
	float paraB = 2 * (x1 * Px + y1 * Py) / a;
	float paraC = 2 * (y1 * Px - x1 * Py) / a + 2 * (Px * Px + Py * Py) / (a * a);
	float A = paraB * paraB + paraC * paraC;
	float B = 2 * paraA * paraB;
	float C = paraA * paraA - paraC * paraC;
	float sol1 = (-B + alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
	float sol2 = (-B - alpaka::math::sqrt(acc, B * B - 4 * A * C)) / (2 * A);
	float solz1 = alpaka::math::asin(acc, sol1) / rou * Pz / p + z1;
	float solz2 = alpaka::math::asin(acc, sol2) / rou * Pz / p + z1;
	float diffz1 = alpaka::math::abs(acc, solz1 - zsi) * 100;
	float diffz2 = alpaka::math::abs(acc, solz2 - zsi) * 100;
	diffz = alpaka::math::min(acc, diffz1, diffz2);
	residual = diffz;
	}

[GNiendorf]

PR #98 addresses some of the above, with a 60% speedup on GPU for MD creation.

[YonsiG]

Triple duplication of chi-squared function. @YonsiG do you know if there is supposed to be any difference between these three, or can they be merged?

cmssw/RecoTracker/LSTCore/src/alpaka/PixelQuintuplet.h

Lines 385 to 428 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquaredpT5(TAcc const& acc,
	unsigned int nPoints,
	float* xs,
	float* ys,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	float g,
	float f,
	float radius) {
	/*
	Given values of (g, f, radius) and a set of points (and its uncertainties) compute chi squared
	*/
	float c = g * g + f * f - radius * radius;
	float chiSquared = 0.f;
	float absArctanSlope, angleM, xPrime, yPrime, sigma2;
	for (size_t i = 0; i < nPoints; i++) {
	absArctanSlope = ((slopes[i] != lst::lst_INF) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
	: 0.5f * float(M_PI));
	if (xs[i] > 0 and ys[i] > 0) {
	angleM = 0.5f * float(M_PI) - absArctanSlope;
	} else if (xs[i] < 0 and ys[i] > 0) {
	angleM = absArctanSlope + 0.5f * float(M_PI);
	} else if (xs[i] < 0 and ys[i] < 0) {
	angleM = -(absArctanSlope + 0.5f * float(M_PI));
	} else if (xs[i] > 0 and ys[i] < 0) {
	angleM = -(0.5f * float(M_PI) - absArctanSlope);
	} else {
	angleM = 0;
	}
	if (not isFlat[i]) {
	xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
	yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
	} else {
	xPrime = xs[i];
	yPrime = ys[i];
	}
	sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
	chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
	(xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / (sigma2);
	}
	return chiSquared;
	};

cmssw/RecoTracker/LSTCore/src/alpaka/Quintuplet.h

Lines 1200 to 1243 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquared(TAcc const& acc,
	unsigned int nPoints,
	float* xs,
	float* ys,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	float g,
	float f,
	float radius) {
	// given values of (g, f, radius) and a set of points (and its uncertainties)
	// compute chi squared
	float c = g * g + f * f - radius * radius;
	float chiSquared = 0.f;
	float absArctanSlope, angleM, xPrime, yPrime, sigma2;
	for (size_t i = 0; i < nPoints; i++) {
	absArctanSlope = ((slopes[i] != lst::lst_INF) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
	: 0.5f * float(M_PI));
	if (xs[i] > 0 and ys[i] > 0) {
	angleM = 0.5f * float(M_PI) - absArctanSlope;
	} else if (xs[i] < 0 and ys[i] > 0) {
	angleM = absArctanSlope + 0.5f * float(M_PI);
	} else if (xs[i] < 0 and ys[i] < 0) {
	angleM = -(absArctanSlope + 0.5f * float(M_PI));
	} else if (xs[i] > 0 and ys[i] < 0) {
	angleM = -(0.5f * float(M_PI) - absArctanSlope);
	} else {
	angleM = 0;
	}
	if (not isFlat[i]) {
	xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
	yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
	} else {
	xPrime = xs[i];
	yPrime = ys[i];
	}
	sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
	chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
	(xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / sigma2;
	}
	return chiSquared;
	};

cmssw/RecoTracker/LSTCore/src/alpaka/PixelTriplet.h

Lines 325 to 368 in ff8446c

	ALPAKA_FN_ACC ALPAKA_FN_INLINE float computeChiSquaredpT3(TAcc const& acc,
	unsigned int nPoints,
	float* xs,
	float* ys,
	float* delta1,
	float* delta2,
	float* slopes,
	bool* isFlat,
	float g,
	float f,
	float radius) {
	//given values of (g, f, radius) and a set of points (and its uncertainties)
	//compute chi squared
	float c = g * g + f * f - radius * radius;
	float chiSquared = 0.f;
	float absArctanSlope, angleM, xPrime, yPrime, sigma2;
	for (size_t i = 0; i < nPoints; i++) {
	absArctanSlope = ((slopes[i] != lst::lst_INF) ? alpaka::math::abs(acc, alpaka::math::atan(acc, slopes[i]))
	: 0.5f * float(M_PI));
	if (xs[i] > 0 and ys[i] > 0) {
	angleM = 0.5f * float(M_PI) - absArctanSlope;
	} else if (xs[i] < 0 and ys[i] > 0) {
	angleM = absArctanSlope + 0.5f * float(M_PI);
	} else if (xs[i] < 0 and ys[i] < 0) {
	angleM = -(absArctanSlope + 0.5f * float(M_PI));
	} else if (xs[i] > 0 and ys[i] < 0) {
	angleM = -(0.5f * float(M_PI) - absArctanSlope);
	} else {
	angleM = 0;
	}
	if (not isFlat[i]) {
	xPrime = xs[i] * alpaka::math::cos(acc, angleM) + ys[i] * alpaka::math::sin(acc, angleM);
	yPrime = ys[i] * alpaka::math::cos(acc, angleM) - xs[i] * alpaka::math::sin(acc, angleM);
	} else {
	xPrime = xs[i];
	yPrime = ys[i];
	}
	sigma2 = 4 * ((xPrime * delta1[i]) * (xPrime * delta1[i]) + (yPrime * delta2[i]) * (yPrime * delta2[i]));
	chiSquared += (xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) *
	(xs[i] * xs[i] + ys[i] * ys[i] - 2 * g * xs[i] - 2 * f * ys[i] + c) / sigma2;
	}
	return chiSquared;
	};

Hi Gavin, for these chi-square function, they are doing chi2 based on x-y plane and radius calculations. I think those can be merged, just the input number of MDs need to be changed.

[GNiendorf]

PR #124 fixes an issue caused by code duplication.