diff --git a/SDL/PixelTriplet.h b/SDL/PixelTriplet.h index 98ad0208..3b399f43 100644 --- a/SDL/PixelTriplet.h +++ b/SDL/PixelTriplet.h @@ -2325,512 +2325,5 @@ namespace SDL } // end i_pLS } }; - - template - ALPAKA_FN_ACC ALPAKA_FN_INLINE void runDeltaBetaIterationspT5(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 += SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), SDL::sinAlphaMax)), betaOut); - return; - } - - if (betaIn * betaOut > 0.f and (alpaka::math::abs(acc, pt_beta) < 4.f * SDL::pt_betaMax or (lIn >= 11 and alpaka::math::abs(acc, pt_beta) < 8.f * SDL::pt_betaMax))) //and the pt_beta is well-defined; less strict for endcap-endcap - { - const float betaInUpd = betaIn + SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), SDL::sinAlphaMax)), betaIn); //FIXME: need a faster version - const float betaOutUpd = betaOut + SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), SDL::sinAlphaMax)), 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 += SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * SDL::k2Rinv1GeVf *pt_beta_inv, SDL::sinAlphaMax)), betaIn); //FIXME: need a faster version - betaOut += SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * SDL::k2Rinv1GeVf *pt_beta_inv, SDL::sinAlphaMax)), betaOut); //FIXME: need a faster version - //update the av and pt - betaAv = 0.5f * (betaIn + betaOut); - //2nd update - pt_beta = dr * SDL::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 * SDL::pt_betaMax) //use betaIn sign as ref - { - const float pt_betaIn = dr * k2Rinv1GeVf / alpaka::math::sin(acc, betaIn); - - const float betaInUpd = betaIn + SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), SDL::sinAlphaMax)), betaIn); //FIXME: need a faster version - const float betaOutUpd = betaOut + SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_betaIn), SDL::sinAlphaMax)), 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 * SDL::k2Rinv1GeVf /alpaka::math::sin(acc, betaAv); //get a better pt estimate - betaIn += SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdIn_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), SDL::sinAlphaMax)), betaIn); //FIXME: need a faster version - betaOut += SDL::copysignf(alpaka::math::asin(acc, alpaka::math::min(acc, sdOut_dr * SDL::k2Rinv1GeVf / alpaka::math::abs(acc, pt_beta), SDL::sinAlphaMax)), betaIn); //FIXME: need a faster version - //update the av and pt - betaAv = 0.5f * (betaIn + betaOut); - //2nd update - pt_beta = dr * SDL::k2Rinv1GeVf /alpaka::math::sin(acc, betaAv); //get a better pt estimate - } - }; - - ALPAKA_FN_ACC ALPAKA_FN_INLINE bool checkIntervalOverlappT5(const float& firstMin, const float& firstMax, const float& secondMin, const float& secondMax) - { - return ((firstMin <= secondMin) & (secondMin < firstMax)) | ((secondMin < firstMin) & (firstMin < secondMax)); - }; - - template - ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runpT5DefaultAlgoPPBB(TAcc const & acc, struct SDL::modules& modulesInGPU, struct SDL::objectRanges& rangesInGPU, struct SDL::miniDoublets& mdsInGPU, struct SDL::segments& segmentsInGPU, uint16_t& pixelModuleIndex, uint16_t& outerInnerLowerModuleIndex, uint16_t& outerOuterLowerModuleIndex, unsigned int& innerSegmentIndex, unsigned int& outerSegmentIndex, unsigned int& firstMDIndex, unsigned int& secondMDIndex, unsigned int thirdMDIndex, unsigned int& fourthMDIndex, float& z_OutLo, float& rt_OutLo, float& dPhiPos, float& dPhi, float& betaIn, float& betaOut, float& pt_beta, float& zLo, float& zHi, float& zLoPointed, float& zHiPointed, float& sdlCut, float& betaOutCut, float& deltaBetaCut) // pixel to BB and BE segments - { - bool pass = true; - - bool isPS_OutLo = (modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS); - - float rt_InLo = mdsInGPU.anchorRt[firstMDIndex]; - float rt_InUp = mdsInGPU.anchorRt[secondMDIndex]; - rt_OutLo = mdsInGPU.anchorRt[thirdMDIndex]; - float rt_OutUp = mdsInGPU.anchorRt[fourthMDIndex]; - - float z_InLo = mdsInGPU.anchorZ[firstMDIndex]; - float z_InUp = mdsInGPU.anchorZ[secondMDIndex]; - z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - float z_OutUp = mdsInGPU.anchorZ[fourthMDIndex]; - - float x_InLo = mdsInGPU.anchorX[firstMDIndex]; - float x_InUp = mdsInGPU.anchorX[secondMDIndex]; - float x_OutLo = mdsInGPU.anchorX[thirdMDIndex]; - float x_OutUp = mdsInGPU.anchorX[fourthMDIndex]; - - float y_InLo = mdsInGPU.anchorY[firstMDIndex]; - float y_InUp = mdsInGPU.anchorY[secondMDIndex]; - float y_OutLo = mdsInGPU.anchorY[thirdMDIndex]; - float y_OutUp = mdsInGPU.anchorY[fourthMDIndex]; - - float& rt_InOut = rt_InUp; - - pass = pass and (alpaka::math::abs(acc, SDL::deltaPhi(acc, x_InUp, y_InUp, x_OutLo, y_OutLo)) <= 0.5f * float(M_PI)); - if(not pass) return pass; - - unsigned int pixelSegmentArrayIndex = innerSegmentIndex - rangesInGPU.segmentModuleIndices[pixelModuleIndex]; - float ptIn = segmentsInGPU.ptIn[pixelSegmentArrayIndex]; - float ptSLo = ptIn; - float px = segmentsInGPU.px[pixelSegmentArrayIndex]; - float py = segmentsInGPU.py[pixelSegmentArrayIndex]; - float pz = segmentsInGPU.pz[pixelSegmentArrayIndex]; - float ptErr = segmentsInGPU.ptErr[pixelSegmentArrayIndex]; - float etaErr = segmentsInGPU.etaErr[pixelSegmentArrayIndex]; - ptSLo = alpaka::math::max(acc, ptCut, ptSLo - 10.0f*alpaka::math::max(acc, ptErr, 0.005f*ptSLo)); - ptSLo = alpaka::math::min(acc, 10.0f, ptSLo); - - float alpha1GeV_OutLo = alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - const float rtRatio_OutLoInOut = rt_OutLo / rt_InOut; // Outer segment beginning rt divided by inner segment beginning rt; - - float dzDrtScale = alpaka::math::tan(acc, alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly - const float zpitch_InLo = 0.05f; - const float zpitch_InOut = 0.05f; - float zpitch_OutLo = (isPS_OutLo ? pixelPSZpitch : strip2SZpitch); - float zGeom = zpitch_InLo + zpitch_OutLo; - zHi = z_InUp + (z_InUp + deltaZLum) * (rtRatio_OutLoInOut - 1.f) * (z_InUp < 0.f ? 1.f : dzDrtScale) + (zpitch_InOut + zpitch_OutLo); - zLo = z_InUp + (z_InUp - deltaZLum) * (rtRatio_OutLoInOut - 1.f) * (z_InUp > 0.f ? 1.f : dzDrtScale) - (zpitch_InOut + zpitch_OutLo); //slope-correction only on outer end - - pass = pass and ((z_OutLo >= zLo) & (z_OutLo <= zHi)); - if(not pass) return pass; - - const float coshEta = alpaka::math::sqrt(acc, ptIn * ptIn + pz * pz) / ptIn; - const float drt_OutLo_InUp = (rt_OutLo - rt_InUp); - const float r3_InUp = alpaka::math::sqrt(acc, z_InUp * z_InUp + rt_InUp * rt_InUp); - float drt_InSeg = rt_InOut - rt_InLo; - - const float sdlThetaMulsF = 0.015f * alpaka::math::sqrt(acc, 0.1f + 0.2f * (rt_OutLo - rt_InUp) / 50.f) * alpaka::math::sqrt(acc, r3_InUp / rt_InUp); - const float sdlMuls = sdlThetaMulsF * 3.f / ptCut * 4.f; // will need a better guess than x4? - - float dzErr = drt_OutLo_InUp*etaErr*coshEta; //FIXME: check with the calc in the endcap - dzErr *= dzErr; - dzErr += 0.03f*0.03f; // pixel size x2. ... random for now - dzErr *= 9.f; //3 sigma - dzErr += sdlMuls*sdlMuls*drt_OutLo_InUp*drt_OutLo_InUp/3.f*coshEta*coshEta;//sloppy - dzErr += zGeom*zGeom; - dzErr = alpaka::math::sqrt(acc, dzErr); - - const float dzDrIn = pz / ptIn; - const float zWindow = dzErr / drt_InSeg * drt_OutLo_InUp + zGeom; - const float dzMean = dzDrIn * drt_OutLo_InUp * (1.f + drt_OutLo_InUp * drt_OutLo_InUp * 4 * k2Rinv1GeVf * k2Rinv1GeVf / ptIn / ptIn / 24.f); // with curved path correction - // Constructing upper and lower bound - zLoPointed = z_InUp + dzMean - zWindow; - zHiPointed = z_InUp + dzMean + zWindow; - - pass = pass and ((z_OutLo >= zLoPointed) & (z_OutLo <= zHiPointed)); - if(not pass) return pass; - - const float sdlPVoff = 0.1f / rt_OutLo; - sdlCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, sdlMuls * sdlMuls + sdlPVoff * sdlPVoff); - - dPhiPos = SDL::deltaPhi(acc, x_InUp, y_InUp, x_OutUp, y_OutUp); - - //no dphipos cut - float midPointX = 0.5f * (x_InLo + x_OutLo); - float midPointY = 0.5f * (y_InLo + y_OutLo); - float midPointZ = 0.5f * (z_InLo + z_OutLo); - - float diffX = x_OutLo - x_InLo; - float diffY = y_OutLo - y_InLo; - float diffZ = z_OutLo - z_InLo; - - dPhi = SDL::deltaPhi(acc, midPointX, midPointY, diffX, diffY); - - pass = pass and (alpaka::math::abs(acc, dPhi) <= sdlCut); - if(not pass) return pass; - - //lots of array accesses below this... - - float alpha_InLo = __H2F(segmentsInGPU.dPhiChanges[innerSegmentIndex]); - float alpha_OutLo = __H2F(segmentsInGPU.dPhiChanges[outerSegmentIndex]); - - bool isEC_lastLayer = modulesInGPU.subdets[outerOuterLowerModuleIndex] == SDL::Endcap and modulesInGPU.moduleType[outerOuterLowerModuleIndex] == SDL::TwoS; - - float alpha_OutUp,alpha_OutUp_highEdge,alpha_OutUp_lowEdge; - alpha_OutUp = SDL::deltaPhi(acc, x_OutUp, y_OutUp, x_OutUp - x_OutLo, y_OutUp - y_OutLo); - - alpha_OutUp_highEdge = alpha_OutUp; - alpha_OutUp_lowEdge = alpha_OutUp; - - float tl_axis_x = x_OutUp - x_InUp; - float tl_axis_y = y_OutUp - y_InUp; - float tl_axis_z = z_OutUp - z_InUp; - - float tl_axis_highEdge_x = tl_axis_x; - float tl_axis_highEdge_y = tl_axis_y; - - float tl_axis_lowEdge_x = tl_axis_x; - float tl_axis_lowEdge_y = tl_axis_y; - - betaIn = -SDL::deltaPhi(acc, px, py, tl_axis_x, tl_axis_y); - float betaInRHmin = betaIn; - float betaInRHmax = betaIn; - - betaOut = -alpha_OutUp + SDL::deltaPhi(acc, x_OutUp, y_OutUp, tl_axis_x, tl_axis_y); - - float betaOutRHmin = betaOut; - float betaOutRHmax = betaOut; - - if(isEC_lastLayer) - { - alpha_OutUp_highEdge = SDL::deltaPhi(acc, mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], mdsInGPU.anchorHighEdgeX[fourthMDIndex] - x_OutLo, mdsInGPU.anchorHighEdgeY[fourthMDIndex] - y_OutLo); - alpha_OutUp_lowEdge = SDL::deltaPhi(acc, mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], mdsInGPU.anchorLowEdgeX[fourthMDIndex] - x_OutLo, mdsInGPU.anchorLowEdgeY[fourthMDIndex] - y_OutLo); - - tl_axis_highEdge_x = mdsInGPU.anchorHighEdgeX[fourthMDIndex] - x_InUp; - tl_axis_highEdge_y = mdsInGPU.anchorHighEdgeY[fourthMDIndex] - y_InUp; - tl_axis_lowEdge_x = mdsInGPU.anchorLowEdgeX[fourthMDIndex] - x_InUp; - tl_axis_lowEdge_y = mdsInGPU.anchorLowEdgeY[fourthMDIndex] - y_InUp; - - betaOutRHmin = -alpha_OutUp_highEdge + SDL::deltaPhi(acc, mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], tl_axis_highEdge_x, tl_axis_highEdge_y); - betaOutRHmax = -alpha_OutUp_lowEdge + SDL::deltaPhi(acc, mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], tl_axis_lowEdge_x, tl_axis_lowEdge_y); - } - - //beta computation - float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - - //innerOuterAnchor - innerInnerAnchor - const float rt_InSeg = alpaka::math::sqrt(acc, (x_InUp - x_InLo) * (x_InUp - x_InLo) + (y_InUp - y_InLo) * (y_InUp - y_InLo)); - - //no betaIn cut for the pixels - float betaAv = 0.5f * (betaIn + betaOut); - pt_beta = ptIn; - - int lIn = 0; - int lOut = isEC_lastLayer ? 11 : 5; - float sdOut_dr = alpaka::math::sqrt(acc, (x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo)); - float sdOut_d = rt_OutUp - rt_OutLo; - - runDeltaBetaIterationspT5(acc, betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn); - - const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax)) : 0.; //mean value of min,max is the old betaIn - const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax)) : 0.; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / alpaka::math::min(acc, alpaka::math::abs(acc, pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - const float alphaInAbsReg = alpaka::math::max(acc, alpaka::math::abs(acc, alpha_InLo), alpaka::math::asin(acc, alpaka::math::min(acc, rt_InUp * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = alpaka::math::max(acc, alpaka::math::abs(acc, alpha_OutLo), alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg*deltaZLum / z_InUp); - const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - - const float sinDPhi = alpaka::math::sin(acc, dPhi); - const float dBetaRIn2 = 0; // TODO-RH - - float dBetaROut = 0; - if(isEC_lastLayer) - { - dBetaROut = (alpaka::math::sqrt(acc, mdsInGPU.anchorHighEdgeX[fourthMDIndex] * mdsInGPU.anchorHighEdgeX[fourthMDIndex] + mdsInGPU.anchorHighEdgeY[fourthMDIndex] * mdsInGPU.anchorHighEdgeY[fourthMDIndex]) - alpaka::math::sqrt(acc, mdsInGPU.anchorLowEdgeX[fourthMDIndex] * mdsInGPU.anchorLowEdgeX[fourthMDIndex] + mdsInGPU.anchorLowEdgeY[fourthMDIndex] * mdsInGPU.anchorLowEdgeY[fourthMDIndex])) * sinDPhi / drt_tl_axis; - } - - const float dBetaROut2 = dBetaROut * dBetaROut; - - //FIXME: need faster version - betaOutCut = alpaka::math::asin(acc, alpaka::math::min(acc, drt_tl_axis*k2Rinv1GeVf / ptCut, sinAlphaMax)) + (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and (alpaka::math::abs(acc, betaOut) < betaOutCut); - if(not pass) return pass; - - const float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, drt_InSeg); - const float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 + 0.25f * (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) * (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax))); - float dBeta = betaIn - betaOut; - deltaBetaCut = alpaka::math::sqrt(acc, dBetaCut2); - - pass = pass and (dBeta * dBeta <= dBetaCut2); - - return pass; - }; - - template - ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runpT5DefaultAlgoPPEE(TAcc const & acc, struct SDL::modules& modulesInGPU, struct SDL::objectRanges& rangesInGPU, struct SDL::miniDoublets& mdsInGPU ,struct SDL::segments& segmentsInGPU, uint16_t& pixelModuleIndex, uint16_t& outerInnerLowerModuleIndex, uint16_t& outerOuterLowerModuleIndex, unsigned int& innerSegmentIndex, unsigned int& outerSegmentIndex, unsigned int& firstMDIndex, unsigned int& secondMDIndex, unsigned int& thirdMDIndex, unsigned int& fourthMDIndex, float& z_OutLo, float& rt_OutLo, float& deltaPhiPos, float& dPhi, float& betaIn, float& betaOut, float& pt_beta, float& zLo, float& rtLo, float& rtHi, float& sdlCut, float& betaInCut, float& betaOutCut, float& deltaBetaCut, float& kZ) // pixel to EE segments - { - bool pass = true; - bool isPS_OutLo = (modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS); - - float z_InLo = mdsInGPU.anchorZ[firstMDIndex]; - float z_InUp = mdsInGPU.anchorZ[secondMDIndex]; - z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - float z_OutUp = mdsInGPU.anchorZ[fourthMDIndex]; - - pass = pass and (z_InUp * z_OutLo > 0); - if(not pass) return pass; - - float rt_InLo = mdsInGPU.anchorRt[firstMDIndex]; - float rt_InUp = mdsInGPU.anchorRt[secondMDIndex]; - rt_OutLo = mdsInGPU.anchorRt[thirdMDIndex]; - float rt_OutUp = mdsInGPU.anchorRt[fourthMDIndex]; - - float x_InLo = mdsInGPU.anchorX[firstMDIndex]; - float x_InUp = mdsInGPU.anchorX[secondMDIndex]; - float x_OutLo = mdsInGPU.anchorX[thirdMDIndex]; - float x_OutUp = mdsInGPU.anchorX[fourthMDIndex]; - - float y_InLo = mdsInGPU.anchorY[firstMDIndex]; - float y_InUp = mdsInGPU.anchorY[secondMDIndex]; - float y_OutLo = mdsInGPU.anchorY[thirdMDIndex]; - float y_OutUp = mdsInGPU.anchorY[fourthMDIndex]; - - unsigned int pixelSegmentArrayIndex = innerSegmentIndex - rangesInGPU.segmentModuleIndices[pixelModuleIndex]; - - float ptIn = segmentsInGPU.ptIn[pixelSegmentArrayIndex]; - float ptSLo = ptIn; - float px = segmentsInGPU.px[pixelSegmentArrayIndex]; - float py = segmentsInGPU.py[pixelSegmentArrayIndex]; - float pz = segmentsInGPU.pz[pixelSegmentArrayIndex]; - float ptErr = segmentsInGPU.ptErr[pixelSegmentArrayIndex]; - float etaErr = segmentsInGPU.etaErr[pixelSegmentArrayIndex]; - - ptSLo = alpaka::math::max(acc, ptCut, ptSLo - 10.0f*alpaka::math::max(acc, ptErr, 0.005f*ptSLo)); - ptSLo = alpaka::math::min(acc, 10.0f, ptSLo); - - float rtOut_o_rtIn = rt_OutLo/rt_InUp; - const float zpitch_InLo = 0.05f; - float zpitch_OutLo = (isPS_OutLo ? pixelPSZpitch : strip2SZpitch); - float zGeom = zpitch_InLo + zpitch_OutLo; - - const float sdlSlope = alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - const float dzDrtScale = alpaka::math::tan(acc, sdlSlope) / sdlSlope;//FIXME: need approximate value - zLo = z_InUp + (z_InUp - deltaZLum) * (rtOut_o_rtIn - 1.f) * (z_InUp > 0.f ? 1.f : dzDrtScale) - zGeom; //slope-correction only on outer end - - const float dLum = SDL::copysignf(deltaZLum, z_InUp); - bool isOutSgInnerMDPS = modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS; - - const float rtGeom1 = isOutSgInnerMDPS ? pixelPSZpitch : strip2SZpitch;//FIXME: make this chosen by configuration for lay11,12 full PS - const float zGeom1 = SDL::copysignf(zGeom, z_InUp); //used in B-E region - rtLo = rt_InUp * (1.f + (z_OutLo- z_InUp - zGeom1) / (z_InUp + zGeom1 + dLum) / dzDrtScale) - rtGeom1; //slope correction only on the lower end - - float zInForHi = z_InUp - zGeom1 - dLum; - if (zInForHi * z_InUp < 0) - zInForHi = SDL::copysignf(0.1f, z_InUp); - rtHi = rt_InUp * (1.f + (z_OutLo - z_InUp + zGeom1) / zInForHi) + rtGeom1; - - // Cut #2: rt condition - pass = pass and ((rt_OutLo >= rtLo) & (rt_OutLo <= rtHi)); - if(not pass) return pass; - - const float dzOutInAbs = alpaka::math::abs(acc, z_OutLo - z_InUp); - const float coshEta = alpaka::math::sqrt(acc, ptIn*ptIn + pz*pz) / ptIn; - const float multDzDr = dzOutInAbs*coshEta/(coshEta*coshEta - 1.f); - const float r3_InUp = alpaka::math::sqrt(acc, z_InUp * z_InUp + rt_InUp * rt_InUp); - const float sdlThetaMulsF = 0.015f * alpaka::math::sqrt(acc, 0.1f + 0.2f * (rt_OutLo - rt_InUp) / 50.f) * alpaka::math::sqrt(acc, r3_InUp / rt_InUp); - const float sdlMuls = sdlThetaMulsF * 3.f / ptCut * 4.f; // will need a better guess than x4? - - float drtErr = etaErr*multDzDr; - drtErr *= drtErr; - drtErr += 0.03f*0.03f; // pixel size x2. ... random for now - drtErr *= 9.f; //3 sigma - drtErr += sdlMuls*sdlMuls*multDzDr*multDzDr/3.f*coshEta*coshEta;//sloppy: relative muls is 1/3 of total muls - drtErr = alpaka::math::sqrt(acc, drtErr); - const float drtDzIn = alpaka::math::abs(acc, ptIn / pz);//all tracks are out-going in endcaps? - - const float drt_OutLo_InUp = (rt_OutLo - rt_InUp); // drOutIn - - const float rtWindow = drtErr + rtGeom1; - const float drtMean = drtDzIn * dzOutInAbs * (1.f - drt_OutLo_InUp * drt_OutLo_InUp * 4 * k2Rinv1GeVf * k2Rinv1GeVf / ptIn / ptIn / 24.f); // with curved path correction - const float rtLo_point = rt_InUp + drtMean - rtWindow; - const float rtHi_point = rt_InUp + drtMean + rtWindow; - - // Cut #3: rt-z pointed - pass = pass and ((rt_OutLo >= rtLo_point) & (rt_OutLo <= rtHi_point)); - if(not pass) return pass; - - const float alpha1GeV_OutLo = alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - const float sdlPVoff = 0.1f / rt_OutLo; - sdlCut = alpha1GeV_OutLo + alpaka::math::sqrt(acc, sdlMuls * sdlMuls + sdlPVoff * sdlPVoff); - - deltaPhiPos = SDL::deltaPhi(acc, x_InUp, y_InUp, x_OutUp, y_OutUp); - - float midPointX = 0.5f * (x_InLo + x_OutLo); - float midPointY = 0.5f * (y_InLo + y_OutLo); - float midPointZ = 0.5f * (z_InLo + z_OutLo); - - float diffX = x_OutLo - x_InLo; - float diffY = y_OutLo - y_InLo; - float diffZ = z_OutLo - z_InLo; - - dPhi = SDL::deltaPhi(acc, midPointX, midPointY, diffX, diffY); - - // Cut #5: deltaPhiChange - pass = pass and (alpaka::math::abs(acc, dPhi) <= sdlCut); - if(not pass) return pass; - - float alpha_InLo = __H2F(segmentsInGPU.dPhiChanges[innerSegmentIndex]); - float alpha_OutLo = __H2F(segmentsInGPU.dPhiChanges[outerSegmentIndex]); - - bool isEC_lastLayer = modulesInGPU.subdets[outerOuterLowerModuleIndex] == SDL::Endcap and modulesInGPU.moduleType[outerOuterLowerModuleIndex] == SDL::TwoS; - - float alpha_OutUp,alpha_OutUp_highEdge,alpha_OutUp_lowEdge; - - alpha_OutUp = SDL::deltaPhi(acc, x_OutUp, y_OutUp, x_OutUp - x_OutLo, y_OutUp - y_OutLo); - alpha_OutUp_highEdge = alpha_OutUp; - alpha_OutUp_lowEdge = alpha_OutUp; - - float tl_axis_x = x_OutUp - x_InUp; - float tl_axis_y = y_OutUp - y_InUp; - float tl_axis_z = z_OutUp - z_InUp; - - float tl_axis_highEdge_x = tl_axis_x; - float tl_axis_highEdge_y = tl_axis_y; - - float tl_axis_lowEdge_x = tl_axis_x; - float tl_axis_lowEdge_y = tl_axis_y; - - betaIn = -SDL::deltaPhi(acc, px, py, tl_axis_x, tl_axis_y); - float betaInRHmin = betaIn; - float betaInRHmax = betaIn; - - betaOut = -alpha_OutUp + SDL::deltaPhi(acc, x_OutUp, y_OutUp, tl_axis_x, tl_axis_y); - float betaOutRHmin = betaOut; - float betaOutRHmax = betaOut; - - if(isEC_lastLayer) - { - alpha_OutUp_highEdge = SDL::deltaPhi(acc, mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], mdsInGPU.anchorHighEdgeX[fourthMDIndex] - x_OutLo, mdsInGPU.anchorHighEdgeY[fourthMDIndex] - y_OutLo); - alpha_OutUp_lowEdge = SDL::deltaPhi(acc, mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], mdsInGPU.anchorLowEdgeX[fourthMDIndex] - x_OutLo, mdsInGPU.anchorLowEdgeY[fourthMDIndex] - y_OutLo); - - tl_axis_highEdge_x = mdsInGPU.anchorHighEdgeX[fourthMDIndex] - x_InUp; - tl_axis_highEdge_y = mdsInGPU.anchorHighEdgeY[fourthMDIndex] - y_InUp; - tl_axis_lowEdge_x = mdsInGPU.anchorLowEdgeX[fourthMDIndex] - x_InUp; - tl_axis_lowEdge_y = mdsInGPU.anchorLowEdgeY[fourthMDIndex] - y_InUp; - - betaOutRHmin = -alpha_OutUp_highEdge + SDL::deltaPhi(acc, mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], tl_axis_highEdge_x, tl_axis_highEdge_y); - betaOutRHmax = -alpha_OutUp_lowEdge + SDL::deltaPhi(acc, mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], tl_axis_lowEdge_x, tl_axis_lowEdge_y); - } - - //beta computation - float drt_tl_axis = alpaka::math::sqrt(acc, tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - //no betaIn cut for the pixels - const float rt_InSeg = alpaka::math::sqrt(acc, (x_InUp - x_InLo) * (x_InUp - x_InLo) + (y_InUp - y_InLo) * (y_InUp - y_InLo)); - - float betaAv = 0.5f * (betaIn + betaOut); - pt_beta = ptIn; - - int lIn = 0; - int lOut = isEC_lastLayer ? 11 : 5; - float sdOut_dr = alpaka::math::sqrt(acc, (x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo)); - float sdOut_d = rt_OutUp - rt_OutLo; - - runDeltaBetaIterationspT5(acc, betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn); - - const float betaInMMSF = (alpaka::math::abs(acc, betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / alpaka::math::abs(acc, betaInRHmin + betaInRHmax)) : 0.; //mean value of min,max is the old betaIn - const float betaOutMMSF = (alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / alpaka::math::abs(acc, betaOutRHmin + betaOutRHmax)) : 0.; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / alpaka::math::min(acc, alpaka::math::abs(acc, pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - - const float alphaInAbsReg = alpaka::math::max(acc, alpaka::math::abs(acc, alpha_InLo), alpaka::math::asin(acc, alpaka::math::min(acc, rt_InUp * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = alpaka::math::max(acc, alpaka::math::abs(acc, alpha_OutLo), alpaka::math::asin(acc, alpaka::math::min(acc, rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : alpaka::math::abs(acc, alphaInAbsReg*deltaZLum / z_InUp); - const float dBetaOutLum = lOut < 11 ? 0.0f : alpaka::math::abs(acc, alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - - const float sinDPhi = alpaka::math::sin(acc, dPhi); - const float dBetaRIn2 = 0; // TODO-RH - - float dBetaROut = 0; - if(isEC_lastLayer) - { - dBetaROut = (alpaka::math::sqrt(acc, mdsInGPU.anchorHighEdgeX[fourthMDIndex] * mdsInGPU.anchorHighEdgeX[fourthMDIndex] + mdsInGPU.anchorHighEdgeY[fourthMDIndex] * mdsInGPU.anchorHighEdgeY[fourthMDIndex]) - alpaka::math::sqrt(acc, mdsInGPU.anchorLowEdgeX[fourthMDIndex] * mdsInGPU.anchorLowEdgeX[fourthMDIndex] + mdsInGPU.anchorLowEdgeY[fourthMDIndex] * mdsInGPU.anchorLowEdgeY[fourthMDIndex])) * sinDPhi / drt_tl_axis; - } - - const float dBetaROut2 = dBetaROut * dBetaROut; - - betaOutCut = alpaka::math::asin(acc, alpaka::math::min(acc, drt_tl_axis*k2Rinv1GeVf / ptCut, sinAlphaMax)) //FIXME: need faster version - + (0.02f / sdOut_d) + alpaka::math::sqrt(acc, dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and (alpaka::math::abs(acc, betaOut) < betaOutCut); - if(not pass) return pass; - - float drt_InSeg = rt_InUp - rt_InLo; - - const float dBetaRes = 0.02f / alpaka::math::min(acc, sdOut_d, drt_InSeg); - const float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 + 0.25f * (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax)) * (alpaka::math::abs(acc, betaInRHmin - betaInRHmax) + alpaka::math::abs(acc, betaOutRHmin - betaOutRHmax))); - float dBeta = betaIn - betaOut; - deltaBetaCut = alpaka::math::sqrt(acc, dBetaCut2); - - pass = pass and (dBeta * dBeta <= dBetaCut2); - return pass; - }; - - template - ALPAKA_FN_ACC ALPAKA_FN_INLINE bool runpT5DefaultAlgo(TAcc const & acc, struct SDL::modules& modulesInGPU, struct SDL::objectRanges& rangesInGPU, struct SDL::miniDoublets& mdsInGPU, struct SDL::segments& segmentsInGPU, uint16_t& pixelLowerModuleIndex, uint16_t& outerInnerLowerModuleIndex, uint16_t& outerOuterLowerModuleIndex, unsigned int& innerSegmentIndex, unsigned int& outerSegmentIndex, float& zOut, float& rtOut, float& deltaPhiPos, float& deltaPhi, float& betaIn, float& betaOut, float& pt_beta, float& zLo, float& zHi, float& rtLo, float& rtHi, float& zLoPointed, float& zHiPointed, float& sdlCut, float& betaInCut, float& betaOutCut, float& deltaBetaCut, float& kZ) - { - zLo = -999; - zHi = -999; - rtLo = -999; - rtHi = -999; - zLoPointed = -999; - zHiPointed = -999; - kZ = -999; - betaInCut = -999; - - short outerInnerLowerModuleSubdet = modulesInGPU.subdets[outerInnerLowerModuleIndex]; - short outerOuterLowerModuleSubdet = modulesInGPU.subdets[outerOuterLowerModuleIndex]; - - unsigned int firstMDIndex = segmentsInGPU.mdIndices[2 * innerSegmentIndex]; - - unsigned int secondMDIndex = segmentsInGPU.mdIndices[2 * innerSegmentIndex + 1]; - unsigned int thirdMDIndex = segmentsInGPU.mdIndices[2 * outerSegmentIndex]; - unsigned int fourthMDIndex = segmentsInGPU.mdIndices[2 * outerSegmentIndex + 1]; - - if(outerInnerLowerModuleSubdet == SDL::Barrel and outerOuterLowerModuleSubdet == SDL::Barrel) - { - return runpT5DefaultAlgoPPBB(acc, modulesInGPU, rangesInGPU, mdsInGPU, segmentsInGPU, pixelLowerModuleIndex, outerInnerLowerModuleIndex, outerOuterLowerModuleIndex, innerSegmentIndex, outerSegmentIndex, firstMDIndex, secondMDIndex, thirdMDIndex, fourthMDIndex, zOut,rtOut,deltaPhiPos,deltaPhi,betaIn,betaOut,pt_beta, zLo, zHi, zLoPointed, zHiPointed, sdlCut, betaOutCut, deltaBetaCut); - } - else if(outerInnerLowerModuleSubdet == SDL::Barrel and outerOuterLowerModuleSubdet == SDL::Endcap) - { - return runpT5DefaultAlgoPPBB(acc, modulesInGPU, rangesInGPU, mdsInGPU, segmentsInGPU, pixelLowerModuleIndex, outerInnerLowerModuleIndex, outerOuterLowerModuleIndex, innerSegmentIndex, outerSegmentIndex, firstMDIndex, secondMDIndex, thirdMDIndex, fourthMDIndex,zOut,rtOut,deltaPhiPos,deltaPhi,betaIn,betaOut,pt_beta, zLo, zHi, zLoPointed, zHiPointed, sdlCut, betaOutCut, deltaBetaCut); - } - else if(outerInnerLowerModuleSubdet == SDL::Endcap and outerOuterLowerModuleSubdet == SDL::Endcap) - { - return runpT5DefaultAlgoPPEE(acc, modulesInGPU, rangesInGPU, mdsInGPU, segmentsInGPU, pixelLowerModuleIndex, outerInnerLowerModuleIndex, outerOuterLowerModuleIndex, innerSegmentIndex, outerSegmentIndex, firstMDIndex, secondMDIndex, thirdMDIndex, fourthMDIndex, zOut,rtOut,deltaPhiPos,deltaPhi,betaIn,betaOut,pt_beta, zLo, rtLo, rtHi, sdlCut, betaInCut, betaOutCut, deltaBetaCut, kZ); - } - return false; - }; } #endif \ No newline at end of file