From e96d9eaf3392035b80bc4f49cbf784f5d51be39e Mon Sep 17 00:00:00 2001 From: GNiendorf Date: Mon, 20 Mar 2023 13:28:04 -0400 Subject: [PATCH] remove tracklet and pixeltracklet files --- SDL/Constants.cu | 2 +- SDL/Constants.cuh | 1 + SDL/Kernels.cuh | 2 - SDL/PixelTracklet.cu | 143 ------- SDL/PixelTracklet.cuh | 545 --------------------------- SDL/PixelTriplet.cu | 1 - SDL/PixelTriplet.cuh | 3 - SDL/Quintuplet.cuh | 1 - SDL/TrackCandidate.cuh | 2 - SDL/Tracklet.cu | 835 ----------------------------------------- SDL/Tracklet.cuh | 156 -------- SDL/Triplet.cuh | 1 - 12 files changed, 2 insertions(+), 1690 deletions(-) delete mode 100644 SDL/PixelTracklet.cu delete mode 100644 SDL/PixelTracklet.cuh delete mode 100644 SDL/Tracklet.cu delete mode 100644 SDL/Tracklet.cuh diff --git a/SDL/Constants.cu b/SDL/Constants.cu index a5016227..70323eae 100644 --- a/SDL/Constants.cu +++ b/SDL/Constants.cu @@ -16,4 +16,4 @@ CUDA_CONST_VAR const float SDL::ptCut = 0.8; CUDA_CONST_VAR const float SDL::deltaZLum = 15.0; CUDA_CONST_VAR const float SDL::pixelPSZpitch = 0.15; CUDA_CONST_VAR const float SDL::strip2SZpitch = 5.0; - +CUDA_CONST_VAR const float SDL::pt_betaMax = 7.0f; diff --git a/SDL/Constants.cuh b/SDL/Constants.cuh index 1b35d62f..b926a4b7 100644 --- a/SDL/Constants.cuh +++ b/SDL/Constants.cuh @@ -72,5 +72,6 @@ namespace SDL extern CUDA_CONST_VAR const float deltaZLum; extern CUDA_CONST_VAR const float pixelPSZpitch; extern CUDA_CONST_VAR const float strip2SZpitch; + extern CUDA_CONST_VAR const float pt_betaMax; } #endif diff --git a/SDL/Kernels.cuh b/SDL/Kernels.cuh index be2dba34..5702d0cd 100644 --- a/SDL/Kernels.cuh +++ b/SDL/Kernels.cuh @@ -18,8 +18,6 @@ #include "Hit.cuh" #include "MiniDoublet.cuh" #include "Segment.cuh" -#include "Tracklet.cuh" -#include "PixelTracklet.cuh" #include "Triplet.cuh" #include "TrackCandidate.cuh" #include "Quintuplet.cuh" diff --git a/SDL/PixelTracklet.cu b/SDL/PixelTracklet.cu deleted file mode 100644 index be25f635..00000000 --- a/SDL/PixelTracklet.cu +++ /dev/null @@ -1,143 +0,0 @@ -# include "PixelTracklet.cuh" - -#ifdef __CUDACC__ -#define CUDA_CONST_VAR __device__ -#endif - -//#ifdef CACHE_ALLOC -#include "allocate.h" -//#endif - -void SDL::createPixelTrackletsInExplicitMemory(struct pixelTracklets& pixelTrackletsInGPU, unsigned int maxPixelTracklets,cudaStream_t stream) -{ -#ifdef CACHE_ALLOC -// cudaStream_t stream = 0; - int dev; - cudaGetDevice(&dev); - - pixelTrackletsInGPU.segmentIndices = (unsigned int*)cms::cuda::allocate_device(dev, maxPixelTracklets * sizeof(unsigned int) * 2,stream); - pixelTrackletsInGPU.lowerModuleIndices = (unsigned int*)cms::cuda::allocate_device(dev, maxPixelTracklets * sizeof(unsigned int) * 2,stream);//split up to avoid runtime error of exceeding max byte allocation at a time - pixelTrackletsInGPU.nPixelTracklets = (unsigned int*)cms::cuda::allocate_device(dev, sizeof(unsigned int),stream); - pixelTrackletsInGPU.zOut = (float*)cms::cuda::allocate_device(dev, maxPixelTracklets * sizeof(float) * 4,stream); - pixelTrackletsInGPU.betaIn = (float*)cms::cuda::allocate_device(dev, maxPixelTracklets * sizeof(float) * 3,stream); - -#else - cudaMalloc(&pixelTrackletsInGPU.segmentIndices, 2 * maxPixelTracklets * sizeof(unsigned int)); - cudaMalloc(&pixelTrackletsInGPU.lowerModuleIndices, 2 * maxPixelTracklets * sizeof(unsigned int)); - cudaMalloc(&pixelTrackletsInGPU.nPixelTracklets, sizeof(unsigned int)); - cudaMalloc(&pixelTrackletsInGPU.zOut, maxPixelTracklets *4* sizeof(float)); - cudaMalloc(&pixelTrackletsInGPU.betaIn, maxPixelTracklets *3* sizeof(float)); -#endif - pixelTrackletsInGPU.rtOut = pixelTrackletsInGPU.zOut + maxPixelTracklets; - pixelTrackletsInGPU.deltaPhiPos = pixelTrackletsInGPU.zOut + maxPixelTracklets * 2; - pixelTrackletsInGPU.deltaPhi = pixelTrackletsInGPU.zOut + maxPixelTracklets * 3; - pixelTrackletsInGPU.betaOut = pixelTrackletsInGPU.betaIn + maxPixelTracklets; - pixelTrackletsInGPU.pt_beta = pixelTrackletsInGPU.betaIn + maxPixelTracklets * 2; - - cudaMemsetAsync(pixelTrackletsInGPU.nPixelTracklets, 0, sizeof(unsigned int),stream); -} - -#ifdef CUT_VALUE_DEBUG -__device__ void SDL::addPixelTrackletToMemory(struct pixelTracklets& pixelTrackletsInGPU, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int innerInnerLowerModuleIndex, unsigned int innerOuterLowerModuleIndex, unsigned int outerInnerLowerModuleIndex, unsigned int outerOuterLowerModuleIndex, 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, unsigned int pixelTrackletIndex) -#else -__device__ void SDL::addPixelTrackletToMemory(struct pixelTracklets& pixelTrackletsInGPU, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int innerInnerLowerModuleIndex, unsigned int innerOuterLowerModuleIndex, unsigned int outerInnerLowerModuleIndex, unsigned int outerOuterLowerModuleIndex, float& zOut, float& rtOut, float& deltaPhiPos, float& deltaPhi, float& betaIn, float& betaOut, float pt_beta, unsigned int pixelTrackletIndex) -#endif -{ - pixelTrackletsInGPU.segmentIndices[2 * pixelTrackletIndex] = innerSegmentIndex; - pixelTrackletsInGPU.segmentIndices[2 * pixelTrackletIndex + 1] = outerSegmentIndex; - pixelTrackletsInGPU.lowerModuleIndices[2 * pixelTrackletIndex] = outerInnerLowerModuleIndex; - pixelTrackletsInGPU.lowerModuleIndices[2 * pixelTrackletIndex + 1] = outerOuterLowerModuleIndex; - - pixelTrackletsInGPU.zOut[pixelTrackletIndex] = zOut; - pixelTrackletsInGPU.rtOut[pixelTrackletIndex] = rtOut; - pixelTrackletsInGPU.deltaPhiPos[pixelTrackletIndex] = deltaPhiPos; - pixelTrackletsInGPU.deltaPhi[pixelTrackletIndex] = deltaPhi; - - pixelTrackletsInGPU.betaIn[pixelTrackletIndex] = betaIn; - pixelTrackletsInGPU.betaOut[pixelTrackletIndex] = betaOut; - pixelTrackletsInGPU.pt_beta[pixelTrackletIndex] = pt_beta; - -#ifdef CUT_VALUE_DEBUG - pixelTrackletsInGPU.zLo[pixelTrackletIndex] = zLo; - pixelTrackletsInGPU.zHi[pixelTrackletIndex] = zHi; - pixelTrackletsInGPU.rtLo[pixelTrackletIndex] = rtLo; - pixelTrackletsInGPU.rtHi[pixelTrackletIndex] = rtHi; - pixelTrackletsInGPU.zLoPointed[pixelTrackletIndex] = zLoPointed; - pixelTrackletsInGPU.zHiPointed[pixelTrackletIndex] = zHiPointed; - pixelTrackletsInGPU.sdlCut[pixelTrackletIndex] = sdlCut; - pixelTrackletsInGPU.betaInCut[pixelTrackletIndex] = betaInCut; - pixelTrackletsInGPU.betaOutCut[pixelTrackletIndex] = betaOutCut; - pixelTrackletsInGPU.deltaBetaCut[pixelTrackletIndex] = deltaBetaCut; - pixelTrackletsInGPU.kZ[pixelTrackletIndex] = kZ; -#endif - -} - -void SDL::pixelTracklets::freeMemoryCache() -{ - int dev; - cudaGetDevice(&dev); - cms::cuda::free_device(dev,segmentIndices); - cms::cuda::free_device(dev,lowerModuleIndices); - cms::cuda::free_device(dev,zOut); - cms::cuda::free_device(dev,betaIn); - cms::cuda::free_device(dev,nPixelTracklets); -} - -void SDL::pixelTracklets::freeMemory() -{ - cudaFree(segmentIndices); - cudaFree(lowerModuleIndices); - cudaFree(nPixelTracklets); - cudaFree(zOut); - cudaFree(betaIn); -#ifdef CUT_VALUE_DEBUG - cudaFree(zLo); - cudaFree(zHi); - cudaFree(rtLo); - cudaFree(rtHi); - cudaFree(zLoPointed); - cudaFree(zHiPointed); - cudaFree(sdlCut); - cudaFree(betaInCut); - cudaFree(betaOutCut); - cudaFree(deltaBetaCut); - cudaFree(kZ); -#endif -} - -SDL::pixelTracklets::pixelTracklets() -{ - segmentIndices = nullptr; - lowerModuleIndices = nullptr; - nPixelTracklets = nullptr; - zOut = nullptr; - rtOut = nullptr; - - deltaPhiPos = nullptr; - deltaPhi = nullptr; - betaIn = nullptr; - betaOut = nullptr; - pt_beta = nullptr; -#ifdef CUT_VALUE_DEBUG - zLo = nullptr; - zHi = nullptr; - rtLo = nullptr; - rtHi = nullptr; - zLoPointed = nullptr; - zHiPointed = nullptr; - sdlCut = nullptr; - betaInCut = nullptr; - betaOutCut = nullptr; - deltaBetaCut = nullptr; - kZ = nullptr; -#endif - -} - -SDL::pixelTracklets::~pixelTracklets() -{ - -} - diff --git a/SDL/PixelTracklet.cuh b/SDL/PixelTracklet.cuh deleted file mode 100644 index 52589a23..00000000 --- a/SDL/PixelTracklet.cuh +++ /dev/null @@ -1,545 +0,0 @@ -#ifndef Pixel_Tracklet_cuh -#define Pixel_Tracklet_cuh - -#ifdef __CUDACC__ -#define CUDA_HOSTDEV __host__ __device__ -#define CUDA_DEV __device__ -#define CUDA_CONST_VAR __device__ -#else -#define CUDA_HOSTDEV -#define CUDA_DEV -#define CUDA_CONST_VAR -#endif - -#include "Constants.cuh" -#include "EndcapGeometry.cuh" -#include "TiltedGeometry.h" -#include "Segment.cuh" -#include "MiniDoublet.cuh" -#include "Module.cuh" -#include "Hit.cuh" -#include "PrintUtil.h" -#include "Tracklet.cuh" - - -namespace SDL -{ - struct pixelTracklets - { - unsigned int* segmentIndices; - unsigned int* lowerModuleIndices; - unsigned int* nPixelTracklets; - float* zOut; - float* rtOut; - - float* deltaPhiPos; - float* deltaPhi; - float* betaIn; - float* betaOut; - float* pt_beta; - -#ifdef CUT_VALUE_DEBUG - //debug variables - float* zLo; - float* zHi; - float* zLoPointed; - float* zHiPointed; - float* sdlCut; - float* betaInCut; - float* betaOutCut; - float* deltaBetaCut; - float* rtLo; - float* rtHi; - float* kZ; -#endif - - pixelTracklets(); - ~pixelTracklets(); - void freeMemory(); - void freeMemoryCache(); - }; - - void createPixelTrackletsInExplicitMemory(struct pixelTracklets& pixelTrackletsInGPU, unsigned int maxPixelTracklets,cudaStream_t stream); - -#ifdef CUT_VALUE_DEBUG - CUDA_DEV void addPixelTrackletToMemory(struct pixelTracklets& pixelTrackletsInGPU, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int innerInnerLowerModuleIndex, unsigned int innerOuterLowerModuleIndex, unsigned int outerInnerLowerModuleIndex, unsigned int outerOuterLowerModuleIndex, 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, unsigned int pixelTrackletIndex); - -#else - CUDA_DEV void addPixelTrackletToMemory(struct pixelTracklets& pixelTrackletsInGPU, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int innerInnerLowerModuleIndex, unsigned int innerOuterLowerModuleIndex, unsigned int outerInnerLowerModuleIndex, unsigned int outerOuterLowerModuleIndex, float& zOut, float& rtOut, float& deltaPhiPos, float& deltaPhi, float& betaIn, float& betaOut, float pt_beta, unsigned int pixelTrackletIndex); -#endif - -CUDA_DEV bool inline runTrackletDefaultAlgoPPBB(struct modules& modulesInGPU, struct objectRanges& rangesInGPU, struct miniDoublets& mdsInGPU ,struct 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_InUp = mdsInGPU.anchorZ[secondMDIndex]; - z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - - 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; - //float& z_InOut = z_InUp; - - pass = pass and (fabsf(deltaPhi(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 = fmaxf(ptCut, ptSLo - 10.0f*fmaxf(ptErr, 0.005f*ptSLo)); - ptSLo = fminf(10.0f, ptSLo); - - - float alpha1GeV_OutLo = asinf(fminf(rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - //float rtRatio_OutLoInLo = rt_OutLo / rt_InLo; // Outer segment beginning rt divided by inner segment beginning rt; - const float rtRatio_OutLoInOut = rt_OutLo / rt_InOut; // Outer segment beginning rt divided by inner segment beginning rt; - - float dzDrtScale = tanf(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 = sqrtf(ptIn * ptIn + pz * pz) / ptIn; - const float drt_OutLo_InUp = (rt_OutLo - rt_InUp); - const float r3_InUp = sqrtf(z_InUp * z_InUp + rt_InUp * rt_InUp); - - float drt_InSeg = rt_InOut - rt_InLo; - - const float sdlThetaMulsF = 0.015f * sqrtf(0.1f + 0.2f * (rt_OutLo - rt_InUp) / 50.f) * sqrtf(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 = sqrtf(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 + sqrtf(sdlMuls * sdlMuls + sdlPVoff * sdlPVoff); - - dPhiPos = deltaPhi(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 diffX = x_OutLo - x_InLo; - float diffY = y_OutLo - y_InLo; - - - dPhi = deltaPhi(midPointX, midPointY, diffX, diffY); - - pass = pass and (fabsf(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 = deltaPhi(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_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 = -deltaPhi(px, py, tl_axis_x, tl_axis_y); - float betaInRHmin = betaIn; - float betaInRHmax = betaIn; - - betaOut = -alpha_OutUp + deltaPhi(x_OutUp, y_OutUp, tl_axis_x, tl_axis_y); - - float betaOutRHmin = betaOut; - float betaOutRHmax = betaOut; - - if(isEC_lastLayer) - { - alpha_OutUp_highEdge = deltaPhi(mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], mdsInGPU.anchorHighEdgeX[fourthMDIndex] - x_OutLo, mdsInGPU.anchorHighEdgeY[fourthMDIndex] - y_OutLo); - alpha_OutUp_lowEdge = deltaPhi(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 + deltaPhi(mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], tl_axis_highEdge_x, tl_axis_highEdge_y); - betaOutRHmax = -alpha_OutUp_lowEdge + deltaPhi(mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], tl_axis_lowEdge_x, tl_axis_lowEdge_y); - } - - //beta computation - float drt_tl_axis = sqrtf(tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - //float drt_tl_lowEdge = sqrtf(tl_axis_lowEdge_x * tl_axis_lowEdge_x + tl_axis_lowEdge_y * tl_axis_lowEdge_y); - //float drt_tl_highEdge = sqrtf(tl_axis_highEdge_x * tl_axis_highEdge_x + tl_axis_highEdge_y * tl_axis_highEdge_y); - - //innerOuterAnchor - innerInnerAnchor - const float rt_InSeg = sqrtf((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 = sqrtf((x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo)); - float sdOut_d = rt_OutUp - rt_OutLo; - - //const float diffDr = fabsf(rt_InSeg - sdOut_dr) / fabsf(rt_InSeg + sdOut_dr); - - runDeltaBetaIterations(betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn); - - const float betaInMMSF = (fabsf(betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / fabsf(betaInRHmin + betaInRHmax)) : 0.; //mean value of min,max is the old betaIn - const float betaOutMMSF = (fabsf(betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / fabsf(betaOutRHmin + betaOutRHmax)) : 0.; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / fminf(fabsf(pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - const float alphaInAbsReg = fmaxf(fabsf(alpha_InLo), asinf(fminf(rt_InUp * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = fmaxf(fabsf(alpha_OutLo), asinf(fminf(rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : fabsf(alphaInAbsReg*deltaZLum / z_InUp); - const float dBetaOutLum = lOut < 11 ? 0.0f : fabsf(alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - - const float sinDPhi = sinf(dPhi); - const float dBetaRIn2 = 0; // TODO-RH - - float dBetaROut = 0; - if(isEC_lastLayer) - { - dBetaROut = (sqrtf(mdsInGPU.anchorHighEdgeX[fourthMDIndex] * mdsInGPU.anchorHighEdgeX[fourthMDIndex] + mdsInGPU.anchorHighEdgeY[fourthMDIndex] * mdsInGPU.anchorHighEdgeY[fourthMDIndex]) - sqrtf(mdsInGPU.anchorLowEdgeX[fourthMDIndex] * mdsInGPU.anchorLowEdgeX[fourthMDIndex] + mdsInGPU.anchorLowEdgeY[fourthMDIndex] * mdsInGPU.anchorLowEdgeY[fourthMDIndex])) * sinDPhi / drt_tl_axis; - } - - const float dBetaROut2 = dBetaROut * dBetaROut; - - betaOutCut = asinf(fminf(drt_tl_axis*k2Rinv1GeVf / ptCut, sinAlphaMax)) //FIXME: need faster version - + (0.02f / sdOut_d) + sqrtf(dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and (fabsf(betaOut) < betaOutCut); - if(not pass) return pass; - - //const float pt_betaOut = drt_tl_axis * k2Rinv1GeVf / sin(betaOut); - const float dBetaRes = 0.02f / fminf(sdOut_d, drt_InSeg); - const float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 - + 0.25f * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax)) * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax))); - float dBeta = betaIn - betaOut; - deltaBetaCut = sqrtf(dBetaCut2); - - pass = pass and (dBeta * dBeta <= dBetaCut2); - - return pass; -} - -CUDA_DEV bool inline runTrackletDefaultAlgoPPEE(struct modules& modulesInGPU, struct objectRanges& rangesInGPU, struct miniDoublets& mdsInGPU ,struct 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_InUp = mdsInGPU.anchorZ[secondMDIndex]; - z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - - 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 = fmaxf(ptCut, ptSLo - 10.0f*fmaxf(ptErr, 0.005f*ptSLo)); - ptSLo = fminf(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 = asinf(fminf(rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - const float dzDrtScale = tanf(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 = 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 = 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 = 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 = fabsf(z_OutLo - z_InUp); - const float coshEta = hypotf(ptIn, pz) / ptIn; - const float multDzDr = dzOutInAbs*coshEta/(coshEta*coshEta - 1.f); - const float r3_InUp = sqrtf(z_InUp * z_InUp + rt_InUp * rt_InUp); - const float sdlThetaMulsF = 0.015f * sqrtf(0.1f + 0.2f * (rt_OutLo - rt_InUp) / 50.f) * sqrtf(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 = sqrtf(drtErr); - const float drtDzIn = fabsf(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 = asinf(fminf(rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - const float sdlPVoff = 0.1f / rt_OutLo; - sdlCut = alpha1GeV_OutLo + sqrtf(sdlMuls * sdlMuls + sdlPVoff * sdlPVoff); - - deltaPhiPos = deltaPhi(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 diffX = x_OutLo - x_InLo; - float diffY = y_OutLo - y_InLo; - - dPhi = deltaPhi(midPointX, midPointY, diffX, diffY); - - // Cut #5: deltaPhiChange - pass = pass and (fabsf(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 = deltaPhi(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_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 = -deltaPhi(px, py, tl_axis_x, tl_axis_y); - float betaInRHmin = betaIn; - float betaInRHmax = betaIn; - - betaOut = -alpha_OutUp + deltaPhi(x_OutUp, y_OutUp, tl_axis_x, tl_axis_y); - float betaOutRHmin = betaOut; - float betaOutRHmax = betaOut; - - if(isEC_lastLayer) - { - - alpha_OutUp_highEdge = deltaPhi(mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], mdsInGPU.anchorHighEdgeX[fourthMDIndex] - x_OutLo, mdsInGPU.anchorHighEdgeY[fourthMDIndex] - y_OutLo); - alpha_OutUp_lowEdge = deltaPhi(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 + deltaPhi(mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], tl_axis_highEdge_x, tl_axis_highEdge_y); - betaOutRHmax = -alpha_OutUp_lowEdge + deltaPhi(mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], tl_axis_lowEdge_x, tl_axis_lowEdge_y); - } - - //beta computation - float drt_tl_axis = sqrtf(tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - //float drt_tl_lowEdge = sqrtf(tl_axis_lowEdge_x * tl_axis_lowEdge_x + tl_axis_lowEdge_y * tl_axis_lowEdge_y); - //float drt_tl_highEdge = sqrtf(tl_axis_highEdge_x * tl_axis_highEdge_x + tl_axis_highEdge_y * tl_axis_highEdge_y); -//no betaIn cut for the pixels - const float rt_InSeg = sqrtf((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 = sqrtf((x_OutUp - x_OutLo) * (x_OutUp - x_OutLo) + (y_OutUp - y_OutLo) * (y_OutUp - y_OutLo)); - float sdOut_d = rt_OutUp - rt_OutLo; - - //const float diffDr = fabsf(rt_InSeg - sdOut_dr) / fabsf(rt_InSeg + sdOut_dr); - - runDeltaBetaIterations(betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn); - - const float betaInMMSF = (fabsf(betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / fabsf(betaInRHmin + betaInRHmax)) : 0.; //mean value of min,max is the old betaIn - const float betaOutMMSF = (fabsf(betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / fabsf(betaOutRHmin + betaOutRHmax)) : 0.; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / fminf(fabsf(pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - - const float alphaInAbsReg = fmaxf(fabsf(alpha_InLo), asinf(fminf(rt_InUp * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = fmaxf(fabsf(alpha_OutLo), asinf(fminf(rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : fabsf(alphaInAbsReg*deltaZLum / z_InUp); - const float dBetaOutLum = lOut < 11 ? 0.0f : fabsf(alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - - const float sinDPhi = sinf(dPhi); - const float dBetaRIn2 = 0; // TODO-RH - - float dBetaROut = 0; - if(isEC_lastLayer) - { - dBetaROut = (sqrtf(mdsInGPU.anchorHighEdgeX[fourthMDIndex] * mdsInGPU.anchorHighEdgeX[fourthMDIndex] + mdsInGPU.anchorHighEdgeY[fourthMDIndex] * mdsInGPU.anchorHighEdgeY[fourthMDIndex]) - sqrtf(mdsInGPU.anchorLowEdgeX[fourthMDIndex] * mdsInGPU.anchorLowEdgeX[fourthMDIndex] + mdsInGPU.anchorLowEdgeY[fourthMDIndex] * mdsInGPU.anchorLowEdgeY[fourthMDIndex])) * sinDPhi / drt_tl_axis; - } - - const float dBetaROut2 = dBetaROut * dBetaROut; - - betaOutCut = asinf(fminf(drt_tl_axis*k2Rinv1GeVf / ptCut, sinAlphaMax)) //FIXME: need faster version - + (0.02f / sdOut_d) + sqrtf(dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and (fabsf(betaOut) < betaOutCut); - if(not pass) return pass; - - // const float pt_betaOut = drt_tl_axis * k2Rinv1GeVf / sin(betaOut); - float drt_InSeg = rt_InUp - rt_InLo; - - const float dBetaRes = 0.02f / fminf(sdOut_d, drt_InSeg); - const float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 - + 0.25f * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax)) * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax))); - float dBeta = betaIn - betaOut; - deltaBetaCut = sqrtf(dBetaCut2); - - pass = pass and (dBeta * dBeta <= dBetaCut2); - return pass; -} - -CUDA_DEV bool inline runPixelTrackletDefaultAlgo(struct modules& modulesInGPU, struct objectRanges& rangesInGPU, struct miniDoublets& mdsInGPU, struct 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 runTrackletDefaultAlgoPPBB(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 runTrackletDefaultAlgoPPBB(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 runTrackletDefaultAlgoPPEE(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 diff --git a/SDL/PixelTriplet.cu b/SDL/PixelTriplet.cu index 9d4c4b11..2e44ebf9 100644 --- a/SDL/PixelTriplet.cu +++ b/SDL/PixelTriplet.cu @@ -3,7 +3,6 @@ #define CUDA_CONST_VAR __device__ #endif # include "PixelTriplet.cuh" -# include "PixelTracklet.cuh" #include "allocate.h" #include "Kernels.cuh" diff --git a/SDL/PixelTriplet.cuh b/SDL/PixelTriplet.cuh index cb952840..1a858cdd 100644 --- a/SDL/PixelTriplet.cuh +++ b/SDL/PixelTriplet.cuh @@ -12,7 +12,6 @@ #include "Constants.cuh" #include "EndcapGeometry.cuh" #include "TiltedGeometry.h" -#include "Tracklet.cuh" #include "Triplet.cuh" #include "Segment.cuh" #include "MiniDoublet.cuh" @@ -20,7 +19,6 @@ #include "Hit.cuh" #include "PrintUtil.h" #include "Quintuplet.cuh" -#include "PixelTracklet.cuh" namespace SDL { @@ -131,7 +129,6 @@ namespace SDL #include "Hit.cuh" #include "PrintUtil.h" #include "Triplet.cuh" -#include "Tracklet.cuh" #include "Quintuplet.cuh" #include "PixelTriplet.cuh" diff --git a/SDL/Quintuplet.cuh b/SDL/Quintuplet.cuh index ee663d3b..250d20fd 100644 --- a/SDL/Quintuplet.cuh +++ b/SDL/Quintuplet.cuh @@ -20,7 +20,6 @@ #include "Hit.cuh" #include "PrintUtil.h" #include "Triplet.cuh" -#include "Tracklet.cuh" namespace SDL { diff --git a/SDL/TrackCandidate.cuh b/SDL/TrackCandidate.cuh index 470e9bb7..50d8de96 100644 --- a/SDL/TrackCandidate.cuh +++ b/SDL/TrackCandidate.cuh @@ -15,8 +15,6 @@ #include "EndcapGeometry.cuh" #include "TiltedGeometry.h" #include "Triplet.cuh" -#include "Tracklet.cuh" -#include "PixelTracklet.cuh" #include "Segment.cuh" #include "MiniDoublet.cuh" #include "PixelTriplet.cuh" diff --git a/SDL/Tracklet.cu b/SDL/Tracklet.cu deleted file mode 100644 index e6e76bc7..00000000 --- a/SDL/Tracklet.cu +++ /dev/null @@ -1,835 +0,0 @@ -#ifdef __CUDACC__ -#define CUDA_CONST_VAR __device__ -#endif -#include "Tracklet.cuh" - -#include "allocate.h" - -CUDA_CONST_VAR float SDL::pt_betaMax = 7.0f; - - -void SDL::createTrackletsInExplicitMemory(struct tracklets& trackletsInGPU, unsigned int maxTracklets, uint16_t nLowerModules,cudaStream_t stream) -{ - - unsigned int nMemoryLocations = maxTracklets * nLowerModules; -#ifdef CACHE_ALLOC -// cudaStream_t stream =0; - int dev; - cudaGetDevice(&dev); - trackletsInGPU.segmentIndices = (unsigned int*)cms::cuda::allocate_device(dev,nMemoryLocations * sizeof(unsigned int) * 2,stream); - trackletsInGPU.lowerModuleIndices = (unsigned int*)cms::cuda::allocate_device(dev,nMemoryLocations * sizeof(unsigned int) * 4,stream);//split up to avoid runtime error of exceeding max byte allocation at a time - trackletsInGPU.nTracklets = (unsigned int*)cms::cuda::allocate_device(dev,nLowerModules * sizeof(unsigned int),stream); - trackletsInGPU.zOut = (float*)cms::cuda::allocate_device(dev,nMemoryLocations * sizeof(float) * 4,stream); - trackletsInGPU.betaIn = (float*)cms::cuda::allocate_device(dev,nMemoryLocations * sizeof(float) * 3,stream); -#else - cudaMalloc(&trackletsInGPU.segmentIndices, 2 * nMemoryLocations * sizeof(unsigned int)); - cudaMalloc(&trackletsInGPU.lowerModuleIndices, 4 * nMemoryLocations * sizeof(unsigned int)); - cudaMalloc(&trackletsInGPU.nTracklets,nLowerModules * sizeof(unsigned int)); - cudaMalloc(&trackletsInGPU.zOut, nMemoryLocations *4* sizeof(float)); - cudaMalloc(&trackletsInGPU.betaIn, nMemoryLocations *3* sizeof(float)); -#endif - cudaMemsetAsync(trackletsInGPU.nTracklets,0,nLowerModules*sizeof(unsigned int),stream); - trackletsInGPU.rtOut = trackletsInGPU.zOut + nMemoryLocations; - trackletsInGPU.deltaPhiPos = trackletsInGPU.zOut + nMemoryLocations * 2; - trackletsInGPU.deltaPhi = trackletsInGPU.zOut + nMemoryLocations * 3; - trackletsInGPU.betaOut = trackletsInGPU.betaIn + nMemoryLocations; - trackletsInGPU.pt_beta = trackletsInGPU.betaIn + nMemoryLocations * 2; -} - -#ifdef CUT_VALUE_DEBUG -__device__ void SDL::addTrackletToMemory(struct tracklets& trackletsInGPU, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int innerInnerLowerModuleIndex, unsigned int innerOuterLowerModuleIndex, unsigned int outerInnerLowerModuleIndex, unsigned int outerOuterLowerModuleIndex, 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, unsigned int trackletIndex) -#else -__device__ void SDL::addTrackletToMemory(struct tracklets& trackletsInGPU, unsigned int innerSegmentIndex, unsigned int outerSegmentIndex, unsigned int innerInnerLowerModuleIndex, unsigned int innerOuterLowerModuleIndex, unsigned int outerInnerLowerModuleIndex, unsigned int outerOuterLowerModuleIndex, float& zOut, float& rtOut, float& deltaPhiPos, float& deltaPhi, float& betaIn, float& betaOut, float pt_beta, unsigned int trackletIndex) -#endif -{ - trackletsInGPU.segmentIndices[trackletIndex * 2] = innerSegmentIndex; - trackletsInGPU.segmentIndices[trackletIndex * 2 + 1] = outerSegmentIndex; - trackletsInGPU.lowerModuleIndices[trackletIndex * 4] = innerInnerLowerModuleIndex; - trackletsInGPU.lowerModuleIndices[trackletIndex * 4 + 1] = innerOuterLowerModuleIndex; - trackletsInGPU.lowerModuleIndices[trackletIndex * 4 + 2] = outerInnerLowerModuleIndex; - trackletsInGPU.lowerModuleIndices[trackletIndex * 4 + 3] = outerOuterLowerModuleIndex; - - trackletsInGPU.zOut[trackletIndex] = zOut; - trackletsInGPU.rtOut[trackletIndex] = rtOut; - trackletsInGPU.deltaPhiPos[trackletIndex] = deltaPhiPos; - trackletsInGPU.deltaPhi[trackletIndex] = deltaPhi; - - trackletsInGPU.betaIn[trackletIndex] = betaIn; - trackletsInGPU.betaOut[trackletIndex] = betaOut; - trackletsInGPU.pt_beta[trackletIndex] = pt_beta; - -#ifdef CUT_VALUE_DEBUG - trackletsInGPU.zLo[trackletIndex] = zLo; - trackletsInGPU.zHi[trackletIndex] = zHi; - trackletsInGPU.rtLo[trackletIndex] = rtLo; - trackletsInGPU.rtHi[trackletIndex] = rtHi; - trackletsInGPU.zLoPointed[trackletIndex] = zLoPointed; - trackletsInGPU.zHiPointed[trackletIndex] = zHiPointed; - trackletsInGPU.sdlCut[trackletIndex] = sdlCut; - trackletsInGPU.betaInCut[trackletIndex] = betaInCut; - trackletsInGPU.betaOutCut[trackletIndex] = betaOutCut; - trackletsInGPU.deltaBetaCut[trackletIndex] = deltaBetaCut; - trackletsInGPU.kZ[trackletIndex] = kZ; -#endif -} - -SDL::tracklets::tracklets() -{ - segmentIndices = nullptr; - lowerModuleIndices = nullptr; - zOut = nullptr; - rtOut = nullptr; - - deltaPhiPos = nullptr; - deltaPhi = nullptr; - betaIn = nullptr; - betaOut = nullptr; - pt_beta = nullptr; -#ifdef CUT_VALUE_DEBUG - zLo = nullptr; - zHi = nullptr; - rtLo = nullptr; - rtHi = nullptr; - zLoPointed = nullptr; - zHiPointed = nullptr; - sdlCut = nullptr; - betaInCut = nullptr; - betaOutCut = nullptr; - deltaBetaCut = nullptr; - kZ = nullptr; -#endif -} - -SDL::tracklets::~tracklets() -{ -} - -void SDL::tracklets::freeMemoryCache() -{ - int dev; - cudaGetDevice(&dev); - cms::cuda::free_device(dev,segmentIndices); - cms::cuda::free_device(dev,lowerModuleIndices); - cms::cuda::free_device(dev,zOut); - cms::cuda::free_device(dev,betaIn); - cms::cuda::free_device(dev,nTracklets); -} -void SDL::tracklets::freeMemory() -{ - cudaFree(segmentIndices); - cudaFree(lowerModuleIndices); - cudaFree(nTracklets); - cudaFree(zOut); - cudaFree(betaIn); -#ifdef CUT_VALUE_DEBUG - cudaFree(zLo); - cudaFree(zHi); - cudaFree(rtLo); - cudaFree(rtHi); - cudaFree(zLoPointed); - cudaFree(zHiPointed); - cudaFree(sdlCut); - cudaFree(betaInCut); - cudaFree(betaOutCut); - cudaFree(deltaBetaCut); - cudaFree(kZ); -#endif -} - - -__device__ bool SDL::runTrackletDefaultAlgoBBBB(struct modules& modulesInGPU, struct miniDoublets& mdsInGPU, struct segments& segmentsInGPU, uint16_t& innerInnerLowerModuleIndex, uint16_t& innerOuterLowerModuleIndex, 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& zOut, float& rtOut, float& deltaPhiPos, float& dPhi, float& betaIn, float& - betaOut, float& pt_beta, float& zLo, float& zHi, float& zLoPointed, float& zHiPointed, float& sdlCut, float& betaInCut, float& betaOutCut, float& deltaBetaCut) -{ - bool pass = true; - - bool isPS_InLo = (modulesInGPU.moduleType[innerInnerLowerModuleIndex] == SDL::PS); - bool isPS_OutLo = (modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS); - - float rt_InLo = mdsInGPU.anchorRt[firstMDIndex]; - float rt_InOut = mdsInGPU.anchorRt[secondMDIndex]; - float rt_OutLo = mdsInGPU.anchorRt[thirdMDIndex]; - - float z_InLo = mdsInGPU.anchorZ[firstMDIndex]; - float z_InOut = mdsInGPU.anchorZ[secondMDIndex]; - float z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - - float alpha1GeV_OutLo = asinf(fminf(rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - - float rtRatio_OutLoInLo = rt_OutLo / rt_InLo; // Outer segment beginning rt divided by inner segment beginning rt; - float dzDrtScale = tanf(alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly - float zpitch_InLo = (isPS_InLo ? pixelPSZpitch : strip2SZpitch); - float zpitch_OutLo = (isPS_OutLo ? pixelPSZpitch : strip2SZpitch); - - zHi = z_InLo + (z_InLo + deltaZLum) * (rtRatio_OutLoInLo - 1.f) * (z_InLo < 0.f ? 1.f : dzDrtScale) + (zpitch_InLo + zpitch_OutLo); - zLo = z_InLo + (z_InLo - deltaZLum) * (rtRatio_OutLoInLo - 1.f) * (z_InLo > 0.f ? 1.f : dzDrtScale) - (zpitch_InLo + zpitch_OutLo); - - - //Cut 1 - z compatibility - zOut = z_OutLo; - rtOut = rt_OutLo; - pass = pass and ((z_OutLo >= zLo) & (z_OutLo <= zHi)); - if(not pass) return pass; - - float drt_OutLo_InLo = (rt_OutLo - rt_InLo); - float r3_InLo = sqrtf(z_InLo * z_InLo + rt_InLo * rt_InLo); - float drt_InSeg = rt_InOut - rt_InLo; - float dz_InSeg = z_InOut - z_InLo; - float dr3_InSeg = sqrtf(rt_InOut * rt_InOut + z_InOut * z_InOut) - sqrtf(rt_InLo * rt_InLo + z_InLo * z_InLo); - - float coshEta = dr3_InSeg/drt_InSeg; - float dzErr = (zpitch_InLo + zpitch_OutLo) * (zpitch_InLo + zpitch_OutLo) * 2.f; - - float sdlThetaMulsF = 0.015f * sqrtf(0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f) * sqrtf(r3_InLo / rt_InLo); - float sdlMuls = sdlThetaMulsF * 3.f / ptCut * 4.f; // will need a better guess than x4? - dzErr += sdlMuls * sdlMuls * drt_OutLo_InLo * drt_OutLo_InLo / 3.f * coshEta * coshEta; //sloppy - dzErr = sqrtf(dzErr); - - // Constructing upper and lower bound - const float dzMean = dz_InSeg / drt_InSeg * drt_OutLo_InLo; - const float zWindow = dzErr / drt_InSeg * drt_OutLo_InLo + (zpitch_InLo + zpitch_OutLo); //FIXME for ptCut lower than ~0.8 need to add curv path correction - zLoPointed = z_InLo + dzMean * (z_InLo > 0.f ? 1.f : dzDrtScale) - zWindow; - zHiPointed = z_InLo + dzMean * (z_InLo < 0.f ? 1.f : dzDrtScale) + zWindow; - - // Cut #2: Pointed Z (Inner segment two MD points to outer segment inner MD) - pass = pass and ((z_OutLo >= zLoPointed) & (z_OutLo <= zHiPointed)); - if(not pass) return pass; - - float sdlPVoff = 0.1f/rt_OutLo; - sdlCut = alpha1GeV_OutLo + sqrtf(sdlMuls * sdlMuls + sdlPVoff * sdlPVoff); - - deltaPhiPos = deltaPhi(mdsInGPU.anchorX[secondMDIndex], mdsInGPU.anchorY[secondMDIndex], mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex]); - // Cut #3: FIXME:deltaPhiPos can be tighter - pass = pass and (fabsf(deltaPhiPos) <= sdlCut); - if(not pass) return pass; - - float midPointX = 0.5f*(mdsInGPU.anchorX[firstMDIndex] + mdsInGPU.anchorX[thirdMDIndex]); - float midPointY = 0.5f* (mdsInGPU.anchorY[firstMDIndex] + mdsInGPU.anchorY[thirdMDIndex]); - float midPointZ = 0.5f*(mdsInGPU.anchorZ[firstMDIndex] + mdsInGPU.anchorZ[thirdMDIndex]); - float diffX = mdsInGPU.anchorX[thirdMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - float diffY = mdsInGPU.anchorY[thirdMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - float diffZ = mdsInGPU.anchorZ[thirdMDIndex] - mdsInGPU.anchorZ[firstMDIndex]; - - dPhi = deltaPhi(midPointX, midPointY, diffX, diffY); - - // Cut #4: deltaPhiChange - pass = pass and (fabsf(dPhi) <= sdlCut); - //lots of array accesses below. Cut here! - if(not pass) return pass; - - // First obtaining the raw betaIn and betaOut values without any correction and just purely based on the mini-doublet hit positions - - 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 = deltaPhi(mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex], mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex], mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]); - - alpha_OutUp_highEdge = alpha_OutUp; - alpha_OutUp_lowEdge = alpha_OutUp; - - float tl_axis_x = mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - float tl_axis_y = mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - float tl_axis_z = mdsInGPU.anchorZ[fourthMDIndex] - mdsInGPU.anchorZ[firstMDIndex]; - 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 = alpha_InLo - deltaPhi(mdsInGPU.anchorX[firstMDIndex], mdsInGPU.anchorY[firstMDIndex], tl_axis_x, tl_axis_y); - - float betaInRHmin = betaIn; - float betaInRHmax = betaIn; - betaOut = -alpha_OutUp + deltaPhi(mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex], tl_axis_x, tl_axis_y); - - float betaOutRHmin = betaOut; - float betaOutRHmax = betaOut; - - if(isEC_lastLayer) - { - alpha_OutUp_highEdge = deltaPhi(mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], mdsInGPU.anchorHighEdgeX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]); - alpha_OutUp_lowEdge = deltaPhi(mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], mdsInGPU.anchorLowEdgeX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]); - - tl_axis_highEdge_x = mdsInGPU.anchorHighEdgeX[fourthMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - tl_axis_highEdge_y = mdsInGPU.anchorHighEdgeY[fourthMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - tl_axis_lowEdge_x = mdsInGPU.anchorLowEdgeX[fourthMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - tl_axis_lowEdge_y = mdsInGPU.anchorLowEdgeY[fourthMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - - - betaOutRHmin = -alpha_OutUp_highEdge + deltaPhi(mdsInGPU.anchorHighEdgeX[fourthMDIndex], mdsInGPU.anchorHighEdgeY[fourthMDIndex], tl_axis_highEdge_x, tl_axis_highEdge_y); - betaOutRHmax = -alpha_OutUp_lowEdge + deltaPhi(mdsInGPU.anchorLowEdgeX[fourthMDIndex], mdsInGPU.anchorLowEdgeY[fourthMDIndex], tl_axis_lowEdge_x, tl_axis_lowEdge_y); - } - - //beta computation - float drt_tl_axis = sqrtf(tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - float drt_tl_lowEdge = sqrtf(tl_axis_lowEdge_x * tl_axis_lowEdge_x + tl_axis_lowEdge_y * tl_axis_lowEdge_y); - float drt_tl_highEdge = sqrtf(tl_axis_highEdge_x * tl_axis_highEdge_x + tl_axis_highEdge_y * tl_axis_highEdge_y); - - float corrF = 1.f; - //innerOuterAnchor - innerInnerAnchor - const float rt_InSeg = sqrtf((mdsInGPU.anchorX[secondMDIndex] - mdsInGPU.anchorX[firstMDIndex]) * (mdsInGPU.anchorX[secondMDIndex] - mdsInGPU.anchorX[firstMDIndex]) + (mdsInGPU.anchorY[secondMDIndex] - mdsInGPU.anchorY[firstMDIndex]) * (mdsInGPU.anchorY[secondMDIndex] - mdsInGPU.anchorY[firstMDIndex])); - betaInCut = asinf(fminf((-rt_InSeg * corrF + drt_tl_axis) * k2Rinv1GeVf / ptCut, sinAlphaMax)) + (0.02f / drt_InSeg); - - //Cut #5: first beta cut - pass = pass and (fabsf(betaInRHmin) < betaInCut); - if(not pass) return pass; - - float betaAv = 0.5f * (betaIn + betaOut); - pt_beta = drt_tl_axis * k2Rinv1GeVf/sinf(betaAv); - int lIn = 5; - int lOut = isEC_lastLayer ? 11 : 5; - float sdOut_dr = sqrtf((mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex]) * (mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex]) + (mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]) * (mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex])); - float sdOut_d = mdsInGPU.anchorRt[fourthMDIndex] - mdsInGPU.anchorRt[thirdMDIndex]; - - const float diffDr = fabsf(rt_InSeg - sdOut_dr) / fabsf(rt_InSeg + sdOut_dr); - - runDeltaBetaIterations(betaIn, betaOut, betaAv, pt_beta, rt_InSeg, sdOut_dr, drt_tl_axis, lIn); - - const float betaInMMSF = (fabsf(betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / fabsf(betaInRHmin + betaInRHmax)) : 0.f; //mean value of min,max is the old betaIn - const float betaOutMMSF = (fabsf(betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / fabsf(betaOutRHmin + betaOutRHmax)) : 0.f; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / fminf(fabsf(pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - - - const float alphaInAbsReg = fmaxf(fabsf(alpha_InLo), asinf(fminf(rt_InLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = fmaxf(fabs(alpha_OutLo), asinf(fminf(rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : fabsf(alphaInAbsReg*deltaZLum / z_InLo); - const float dBetaOutLum = lOut < 11 ? 0.0f : fabsf(alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - const float sinDPhi = sinf(dPhi); - - const float dBetaRIn2 = 0; // TODO-RH - // const float dBetaROut2 = 0; // TODO-RH - float dBetaROut = 0; - if(isEC_lastLayer) - { - dBetaROut = (sqrtf(mdsInGPU.anchorHighEdgeX[fourthMDIndex] * mdsInGPU.anchorHighEdgeX[fourthMDIndex] + mdsInGPU.anchorHighEdgeY[fourthMDIndex] * mdsInGPU.anchorHighEdgeY[fourthMDIndex]) - sqrtf(mdsInGPU.anchorLowEdgeX[fourthMDIndex] * mdsInGPU.anchorLowEdgeX[fourthMDIndex] + mdsInGPU.anchorLowEdgeY[fourthMDIndex] * mdsInGPU.anchorLowEdgeY[fourthMDIndex])) * sinDPhi / drt_tl_axis; - } - - const float dBetaROut2 = dBetaROut * dBetaROut; - - betaOutCut = asinf(fminf(drt_tl_axis*k2Rinv1GeVf / ptCut, sinAlphaMax)) //FIXME: need faster version - + (0.02f / sdOut_d) + sqrtf(dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and ((fabsf(betaOut) < betaOutCut)); - if(not pass) return pass; - - float pt_betaIn = drt_tl_axis * k2Rinv1GeVf/sinf(betaIn); - float pt_betaOut = drt_tl_axis * k2Rinv1GeVf / sinf(betaOut); - float dBetaRes = 0.02f/fminf(sdOut_d,drt_InSeg); - float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 - + 0.25f * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax)) * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax))); - - float dBeta = betaIn - betaOut; - deltaBetaCut = sqrtf(dBetaCut2); - pass = pass and (dBeta * dBeta <= dBetaCut2); - - return pass; -} - -__device__ bool SDL::runTrackletDefaultAlgoBBEE(struct modules& modulesInGPU, struct miniDoublets& mdsInGPU, struct segments& segmentsInGPU, uint16_t& innerInnerLowerModuleIndex, uint16_t& innerOuterLowerModuleIndex, 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& zOut, float& rtOut, 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) -{ - bool pass = true; - bool isPS_InLo = (modulesInGPU.moduleType[innerInnerLowerModuleIndex] == SDL::PS); - bool isPS_OutLo = (modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS); - - float rt_InLo = mdsInGPU.anchorRt[firstMDIndex]; - float rt_InOut = mdsInGPU.anchorRt[secondMDIndex]; - float rt_OutLo = mdsInGPU.anchorRt[thirdMDIndex]; - - float z_InLo = mdsInGPU.anchorZ[firstMDIndex]; - float z_InOut = mdsInGPU.anchorZ[secondMDIndex]; - float z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - - float alpha1GeV_OutLo = asinf(fminf(rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - - float rtRatio_OutLoInLo = rt_OutLo / rt_InLo; // Outer segment beginning rt divided by inner segment beginning rt; - float dzDrtScale = tanf(alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly - float zpitch_InLo = (isPS_InLo ? pixelPSZpitch : strip2SZpitch); - float zpitch_OutLo = (isPS_OutLo ? pixelPSZpitch : strip2SZpitch); - float zGeom = zpitch_InLo + zpitch_OutLo; - - zLo = z_InLo + (z_InLo - deltaZLum) * (rtRatio_OutLoInLo - 1.f) * (z_InLo > 0.f ? 1.f : dzDrtScale) - zGeom; - - // Cut #0: Preliminary (Only here in endcap case) - pass = pass and (z_InLo * z_OutLo > 0); - if(not pass) return pass; - - float dLum = copysignf(deltaZLum, z_InLo); - bool isOutSgInnerMDPS = modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS; - float rtGeom1 = isOutSgInnerMDPS ? pixelPSZpitch : strip2SZpitch; - float zGeom1 = copysignf(zGeom,z_InLo); - rtLo = rt_InLo * (1.f + (z_OutLo - z_InLo - zGeom1) / (z_InLo + zGeom1 + dLum) / dzDrtScale) - rtGeom1; //slope correction only on the lower end - zOut = z_OutLo; - rtOut = rt_OutLo; - - //Cut #1: rt condition - pass = pass and (rtOut >= rtLo); - if(not pass) return pass; - - float zInForHi = z_InLo - zGeom1 - dLum; - if(zInForHi * z_InLo < 0) - { - zInForHi = copysignf(0.1f,z_InLo); - } - rtHi = rt_InLo * (1.f + (z_OutLo - z_InLo + zGeom1) / zInForHi) + rtGeom1; - - //Cut #2: rt condition - pass = pass and ((rt_OutLo >= rtLo) & (rt_OutLo <= rtHi)); - if(not pass) return pass; - - float rIn = sqrtf(z_InLo * z_InLo + rt_InLo * rt_InLo); - const float drtSDIn = rt_InOut - rt_InLo; - const float dzSDIn = z_InOut - z_InLo; - const float dr3SDIn = sqrtf(rt_InOut * rt_InOut + z_InOut * z_InOut) - sqrtf(rt_InLo * rt_InLo + z_InLo * z_InLo); - - const float coshEta = dr3SDIn / drtSDIn; //direction estimate - const float dzOutInAbs = fabsf(z_OutLo - z_InLo); - const float multDzDr = dzOutInAbs * coshEta / (coshEta * coshEta - 1.f); - const float zGeom1_another = pixelPSZpitch; //What's this? - kZ = (z_OutLo - z_InLo) / dzSDIn; - float drtErr = zGeom1_another * zGeom1_another * drtSDIn * drtSDIn / dzSDIn / dzSDIn * (1.f - 2.f * kZ + 2.f * kZ * kZ); //Notes:122316 - const float sdlThetaMulsF = 0.015f * sqrtf(0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f) * sqrtf(rIn / rt_InLo); - const float sdlMuls = sdlThetaMulsF * 3.f / ptCut * 4.f; //will need a better guess than x4? - drtErr += sdlMuls * sdlMuls * multDzDr * multDzDr / 3.f * coshEta * coshEta; //sloppy: relative muls is 1/3 of total muls - drtErr = sqrtf(drtErr); - const float drtMean = drtSDIn * dzOutInAbs / fabsf(dzSDIn); // - const float rtWindow = drtErr + rtGeom1; - const float rtLo_another = rt_InLo + drtMean / dzDrtScale - rtWindow; - const float rtHi_another = rt_InLo + drtMean + rtWindow; - - //Cut #3: rt-z pointed - pass = pass and ((kZ >= 0) & (rtOut >= rtLo) & (rtOut <= rtHi)); - if(not pass) return pass; - - const float sdlPVoff = 0.1f / rt_OutLo; - sdlCut = alpha1GeV_OutLo + sqrtf(sdlMuls * sdlMuls + sdlPVoff*sdlPVoff); - - - deltaPhiPos = deltaPhi(mdsInGPU.anchorX[secondMDIndex], mdsInGPU.anchorY[secondMDIndex], mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex]); - - - //Cut #4: deltaPhiPos can be tighter - pass = pass and (fabsf(deltaPhiPos) <= sdlCut); - if(not pass) return pass; - - float midPointX = 0.5f*(mdsInGPU.anchorX[firstMDIndex] + mdsInGPU.anchorX[thirdMDIndex]); - float midPointY = 0.5f* (mdsInGPU.anchorY[firstMDIndex] + mdsInGPU.anchorY[thirdMDIndex]); - float midPointZ = 0.5f*(mdsInGPU.anchorZ[firstMDIndex] + mdsInGPU.anchorZ[thirdMDIndex]); - float diffX = mdsInGPU.anchorX[thirdMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - float diffY = mdsInGPU.anchorY[thirdMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - float diffZ = mdsInGPU.anchorZ[thirdMDIndex] - mdsInGPU.anchorZ[firstMDIndex]; - - dPhi = deltaPhi(midPointX, midPointY, diffX, diffY); - // Cut #5: deltaPhiChange - pass = pass and (fabsf(dPhi) <= sdlCut); - if(not pass) return pass; - - float sdIn_alpha = __H2F(segmentsInGPU.dPhiChanges[innerSegmentIndex]); - float sdIn_alpha_min = __H2F(segmentsInGPU.dPhiChangeMins[innerSegmentIndex]); - float sdIn_alpha_max = __H2F(segmentsInGPU.dPhiChangeMaxs[innerSegmentIndex]); - float sdOut_alpha = sdIn_alpha; //weird - - float sdOut_alphaOut = deltaPhi(mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex], mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex], mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]); - - float sdOut_alphaOut_min = phi_mpi_pi(__H2F(segmentsInGPU.dPhiChangeMins[outerSegmentIndex]) - __H2F(segmentsInGPU.dPhiMins[outerSegmentIndex])); - float sdOut_alphaOut_max = phi_mpi_pi(__H2F(segmentsInGPU.dPhiChangeMaxs[outerSegmentIndex]) - __H2F(segmentsInGPU.dPhiMaxs[outerSegmentIndex])); - - float tl_axis_x = mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - float tl_axis_y = mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - float tl_axis_z = mdsInGPU.anchorZ[fourthMDIndex] - mdsInGPU.anchorZ[firstMDIndex]; - - betaIn = sdIn_alpha - deltaPhi(mdsInGPU.anchorX[firstMDIndex], mdsInGPU.anchorY[firstMDIndex], tl_axis_x, tl_axis_y); - - float betaInRHmin = betaIn; - float betaInRHmax = betaIn; - betaOut = -sdOut_alphaOut + deltaPhi(mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex], tl_axis_x, tl_axis_y); - - float betaOutRHmin = betaOut; - float betaOutRHmax = betaOut; - - bool isEC_secondLayer = (modulesInGPU.subdets[innerOuterLowerModuleIndex] == SDL::Endcap) and (modulesInGPU.moduleType[innerOuterLowerModuleIndex] == SDL::TwoS); - - if(isEC_secondLayer) - { - betaInRHmin = betaIn - sdIn_alpha_min + sdIn_alpha; - betaInRHmax = betaIn - sdIn_alpha_max + sdIn_alpha; - } - - betaOutRHmin = betaOut - sdOut_alphaOut_min + sdOut_alphaOut; - betaOutRHmax = betaOut - sdOut_alphaOut_max + sdOut_alphaOut; - - float swapTemp; - if(fabsf(betaOutRHmin) > fabsf(betaOutRHmax)) - { - swapTemp = betaOutRHmin; - betaOutRHmin = betaOutRHmax; - betaOutRHmax = swapTemp; - } - - if(fabsf(betaInRHmin) > fabsf(betaInRHmax)) - { - swapTemp = betaInRHmin; - betaInRHmin = betaInRHmax; - betaInRHmax = swapTemp; - } - - float sdIn_dr = sqrtf((mdsInGPU.anchorX[secondMDIndex] - mdsInGPU.anchorX[firstMDIndex]) * (mdsInGPU.anchorX[secondMDIndex] - mdsInGPU.anchorX[firstMDIndex]) + (mdsInGPU.anchorY[secondMDIndex] - mdsInGPU.anchorY[firstMDIndex]) * (mdsInGPU.anchorY[secondMDIndex] - mdsInGPU.anchorY[firstMDIndex])); - float sdIn_d = rt_InOut - rt_InLo; - - float dr = sqrtf(tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - const float corrF = 1.f; - betaInCut = asinf(fminf((-sdIn_dr * corrF + dr) * k2Rinv1GeVf / ptCut, sinAlphaMax)) + (0.02f / sdIn_d); - - //Cut #6: first beta cut - pass = pass and (fabsf(betaInRHmin) < betaInCut); - if(not pass) return pass; - - float betaAv = 0.5f * (betaIn + betaOut); - pt_beta = dr * k2Rinv1GeVf / sinf(betaAv); - - float lIn = 5; - float lOut = 11; - - float sdOut_dr = sqrtf((mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex]) * (mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex]) + (mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]) * (mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex])); - float sdOut_d = mdsInGPU.anchorRt[fourthMDIndex] - mdsInGPU.anchorRt[thirdMDIndex]; - - runDeltaBetaIterations(betaIn, betaOut, betaAv, pt_beta, sdIn_dr, sdOut_dr, dr, lIn); - - const float betaInMMSF = (fabsf(betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / fabsf(betaInRHmin + betaInRHmax)) : 0.; //mean value of min,max is the old betaIn - const float betaOutMMSF = (fabsf(betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / fabsf(betaOutRHmin + betaOutRHmax)) : 0.; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / fminf(fabsf(pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - - const float alphaInAbsReg = fmaxf(fabsf(sdIn_alpha), asinf(fminf(rt_InLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = fmaxf(fabsf(sdOut_alpha), asinf(fminf(rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : fabsf(alphaInAbsReg*deltaZLum / z_InLo); - const float dBetaOutLum = lOut < 11 ? 0.0f : fabsf(alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - const float sinDPhi = sinf(dPhi); - - const float dBetaRIn2 = 0; // TODO-RH - // const float dBetaROut2 = 0; // TODO-RH - float dBetaROut = 0; - if(modulesInGPU.moduleType[outerOuterLowerModuleIndex] == SDL::TwoS) - { - dBetaROut = (sqrtf(mdsInGPU.anchorHighEdgeX[fourthMDIndex] * mdsInGPU.anchorHighEdgeX[fourthMDIndex] + mdsInGPU.anchorHighEdgeY[fourthMDIndex] * mdsInGPU.anchorHighEdgeY[fourthMDIndex]) - sqrtf(mdsInGPU.anchorLowEdgeX[fourthMDIndex] * mdsInGPU.anchorLowEdgeX[fourthMDIndex] + mdsInGPU.anchorLowEdgeY[fourthMDIndex] * mdsInGPU.anchorLowEdgeY[fourthMDIndex])) * sinDPhi / dr; - } - - const float dBetaROut2 = dBetaROut * dBetaROut; - betaOutCut = asinf(fminf(dr*k2Rinv1GeVf / ptCut, sinAlphaMax)) //FIXME: need faster version - + (0.02f / sdOut_d) + sqrtf(dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and (fabsf(betaOut) < betaOutCut); - if(not pass) return pass; - - float pt_betaIn = dr * k2Rinv1GeVf/sinf(betaIn); - float pt_betaOut = dr * k2Rinv1GeVf / sinf(betaOut); - float dBetaRes = 0.02f/fminf(sdOut_d,sdIn_d); - float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 - + 0.25f * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax)) * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax))); - float dBeta = betaIn - betaOut; - deltaBetaCut = sqrtf(dBetaCut2); - //Cut #7: Cut on dBet - pass = pass and (dBeta * dBeta <= dBetaCut2); - - return pass; -} - -__device__ bool SDL::runTrackletDefaultAlgoEEEE(struct modules& modulesInGPU, struct miniDoublets& mdsInGPU, struct segments& segmentsInGPU, uint16_t& innerInnerLowerModuleIndex, uint16_t& innerOuterLowerModuleIndex, 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& zOut, float& rtOut, 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) -{ - bool pass = true; - - bool isPS_InLo = (modulesInGPU.moduleType[innerInnerLowerModuleIndex] == SDL::PS); - bool isPS_OutLo = (modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS); - - float rt_InLo = mdsInGPU.anchorRt[firstMDIndex]; - float rt_InOut = mdsInGPU.anchorRt[secondMDIndex]; - float rt_OutLo = mdsInGPU.anchorRt[thirdMDIndex]; - - float z_InLo = mdsInGPU.anchorZ[firstMDIndex]; - float z_InOut = mdsInGPU.anchorZ[secondMDIndex]; - float z_OutLo = mdsInGPU.anchorZ[thirdMDIndex]; - - float alpha1GeV_OutLo = asinf(fminf(rt_OutLo * k2Rinv1GeVf / ptCut, sinAlphaMax)); - - float rtRatio_OutLoInLo = rt_OutLo / rt_InLo; // Outer segment beginning rt divided by inner segment beginning rt; - float dzDrtScale = tanf(alpha1GeV_OutLo) / alpha1GeV_OutLo; // The track can bend in r-z plane slightly - float zpitch_InLo = (isPS_InLo ? pixelPSZpitch : strip2SZpitch); - float zpitch_OutLo = (isPS_OutLo ? pixelPSZpitch : strip2SZpitch); - float zGeom = zpitch_InLo + zpitch_OutLo; - - zLo = z_InLo + (z_InLo - deltaZLum) * (rtRatio_OutLoInLo - 1.f) * (z_InLo > 0.f ? 1.f : dzDrtScale) - zGeom; //slope-correction only on outer end - - // Cut #0: Preliminary (Only here in endcap case) - pass = pass and ((z_InLo * z_OutLo) > 0); - if(not pass) return pass; - - float dLum = copysignf(deltaZLum, z_InLo); - bool isOutSgInnerMDPS = modulesInGPU.moduleType[outerInnerLowerModuleIndex] == SDL::PS; - bool isInSgInnerMDPS = modulesInGPU.moduleType[innerInnerLowerModuleIndex] == SDL::PS; - - float rtGeom = (isInSgInnerMDPS and isOutSgInnerMDPS) ? 2.f * pixelPSZpitch : (isInSgInnerMDPS or isOutSgInnerMDPS) ? pixelPSZpitch + strip2SZpitch : 2.f * strip2SZpitch; - - float zGeom1 = copysignf(zGeom,z_InLo); - float dz = z_OutLo - z_InLo; - rtLo = rt_InLo * (1.f + dz / (z_InLo + dLum) / dzDrtScale) - rtGeom; //slope correction only on the lower end - - zOut = z_OutLo; - rtOut = rt_OutLo; - - //Cut #1: rt condition - - rtHi = rt_InLo * (1.f + dz / (z_InLo - dLum)) + rtGeom; - - pass = pass and ((rtOut >= rtLo) & (rtOut <= rtHi)); - if(not pass) return pass; - - bool isInSgOuterMDPS = modulesInGPU.moduleType[innerOuterLowerModuleIndex] == SDL::PS; - - float drOutIn = rtOut - rt_InLo; - const float drtSDIn = rt_InOut - rt_InLo; - const float dzSDIn = z_InOut - z_InLo; - const float dr3SDIn = sqrtf(rt_InOut * rt_InOut + z_InOut * z_InOut) - sqrtf(rt_InLo * rt_InLo + z_InLo * z_InLo); - float coshEta = dr3SDIn / drtSDIn; //direction estimate - float dzOutInAbs = fabsf(z_OutLo - z_InLo); - float multDzDr = dzOutInAbs * coshEta / (coshEta * coshEta - 1.f); - - kZ = (z_OutLo - z_InLo) / dzSDIn; - float sdlThetaMulsF = 0.015f * sqrtf(0.1f + 0.2f * (rt_OutLo - rt_InLo) / 50.f); - - float sdlMuls = sdlThetaMulsF * 3.f / ptCut * 4.f; //will need a better guess than x4? - - float drtErr = sqrtf(pixelPSZpitch * pixelPSZpitch * 2.f / (dzSDIn * dzSDIn) * (dzOutInAbs * dzOutInAbs) + sdlMuls * sdlMuls * multDzDr * multDzDr / 3.f * coshEta * coshEta); - - float drtMean = drtSDIn * dzOutInAbs/fabsf(dzSDIn); - float rtWindow = drtErr + rtGeom; - float rtLo_point = rt_InLo + drtMean / dzDrtScale - rtWindow; - float rtHi_point = rt_InLo + drtMean + rtWindow; - - // Cut #3: rt-z pointed - // https://github.com/slava77/cms-tkph2-ntuple/blob/superDoubletLinked-91X-noMock/doubletAnalysis.C#L3765 - - if (isInSgInnerMDPS and isInSgOuterMDPS) // If both PS then we can point - { - pass = pass and (kZ >= 0 and rtOut >= rtLo_point and rtOut <= rtHi_point); - if(not pass) return pass; - } - - float sdlPVoff = 0.1f/rtOut; - sdlCut = alpha1GeV_OutLo + sqrtf(sdlMuls * sdlMuls + sdlPVoff * sdlPVoff); - - deltaPhiPos = deltaPhi(mdsInGPU.anchorX[secondMDIndex], mdsInGPU.anchorY[secondMDIndex], mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex]); - - pass = pass and (fabsf(deltaPhiPos) <= sdlCut); - if(not pass) return pass; - - float midPointX = 0.5f*(mdsInGPU.anchorX[firstMDIndex] + mdsInGPU.anchorX[thirdMDIndex]); - float midPointY = 0.5f* (mdsInGPU.anchorY[firstMDIndex] + mdsInGPU.anchorY[thirdMDIndex]); - float midPointZ = 0.5f*(mdsInGPU.anchorZ[firstMDIndex] + mdsInGPU.anchorZ[thirdMDIndex]); - float diffX = mdsInGPU.anchorX[thirdMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - float diffY = mdsInGPU.anchorY[thirdMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - float diffZ = mdsInGPU.anchorZ[thirdMDIndex] - mdsInGPU.anchorZ[firstMDIndex]; - - dPhi = deltaPhi(midPointX, midPointY, diffX, diffY); - - // Cut #5: deltaPhiChange - pass = pass and ((fabsf(dPhi) <= sdlCut)); - if(not pass) return pass; - - float sdIn_alpha = __H2F(segmentsInGPU.dPhiChanges[innerSegmentIndex]); - float sdOut_alpha = sdIn_alpha; //weird - float sdOut_dPhiPos = deltaPhi(mdsInGPU.anchorX[thirdMDIndex], mdsInGPU.anchorY[thirdMDIndex], mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex]); - - float sdOut_dPhiChange = __H2F(segmentsInGPU.dPhiChanges[outerSegmentIndex]); - float sdOut_dPhiChange_min = __H2F(segmentsInGPU.dPhiChangeMins[outerSegmentIndex]); - float sdOut_dPhiChange_max = __H2F(segmentsInGPU.dPhiChangeMaxs[outerSegmentIndex]); - - float sdOut_alphaOutRHmin = phi_mpi_pi(sdOut_dPhiChange_min - sdOut_dPhiPos); - float sdOut_alphaOutRHmax = phi_mpi_pi(sdOut_dPhiChange_max - sdOut_dPhiPos); - float sdOut_alphaOut = phi_mpi_pi(sdOut_dPhiChange - sdOut_dPhiPos); - - float tl_axis_x = mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[firstMDIndex]; - float tl_axis_y = mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[firstMDIndex]; - float tl_axis_z = mdsInGPU.anchorZ[fourthMDIndex] - mdsInGPU.anchorZ[firstMDIndex]; - - betaIn = sdIn_alpha - deltaPhi(mdsInGPU.anchorX[firstMDIndex], mdsInGPU.anchorY[firstMDIndex], tl_axis_x, tl_axis_y); - - float sdIn_alphaRHmin = __H2F(segmentsInGPU.dPhiChangeMins[innerSegmentIndex]); - float sdIn_alphaRHmax = __H2F(segmentsInGPU.dPhiChangeMaxs[innerSegmentIndex]); - float betaInRHmin = betaIn + sdIn_alphaRHmin - sdIn_alpha; - float betaInRHmax = betaIn + sdIn_alphaRHmax - sdIn_alpha; - - betaOut = -sdOut_alphaOut + deltaPhi(mdsInGPU.anchorX[fourthMDIndex], mdsInGPU.anchorY[fourthMDIndex], tl_axis_x, tl_axis_y); - - float betaOutRHmin = betaOut - sdOut_alphaOutRHmin + sdOut_alphaOut; - float betaOutRHmax = betaOut - sdOut_alphaOutRHmax + sdOut_alphaOut; - - float swapTemp; - if(fabsf(betaOutRHmin) > fabsf(betaOutRHmax)) - { - swapTemp = betaOutRHmin; - betaOutRHmin = betaOutRHmax; - betaOutRHmax = swapTemp; - } - - if(fabsf(betaInRHmin) > fabsf(betaInRHmax)) - { - swapTemp = betaInRHmin; - betaInRHmin = betaInRHmax; - betaInRHmax = swapTemp; - } - float sdIn_dr = sqrtf((mdsInGPU.anchorX[secondMDIndex] - mdsInGPU.anchorX[firstMDIndex]) * (mdsInGPU.anchorX[secondMDIndex] - mdsInGPU.anchorX[firstMDIndex]) + (mdsInGPU.anchorY[secondMDIndex] - mdsInGPU.anchorY[firstMDIndex]) * (mdsInGPU.anchorY[secondMDIndex] - mdsInGPU.anchorY[firstMDIndex])); - float sdIn_d = rt_InOut - rt_InLo; - - float dr = sqrtf(tl_axis_x * tl_axis_x + tl_axis_y * tl_axis_y); - const float corrF = 1.f; - betaInCut = asinf(fminf((-sdIn_dr * corrF + dr) * k2Rinv1GeVf / ptCut, sinAlphaMax)) + (0.02f / sdIn_d); - - //Cut #6: first beta cut - pass = pass and (fabsf(betaInRHmin) < betaInCut); - if(not pass) return pass; - - float betaAv = 0.5f * (betaIn + betaOut); - pt_beta = dr * k2Rinv1GeVf / sinf(betaAv); - - - int lIn= 11; //endcap - int lOut = 13; //endcap - - float sdOut_dr = sqrtf((mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex]) * (mdsInGPU.anchorX[fourthMDIndex] - mdsInGPU.anchorX[thirdMDIndex]) + (mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex]) * (mdsInGPU.anchorY[fourthMDIndex] - mdsInGPU.anchorY[thirdMDIndex])); - float sdOut_d = mdsInGPU.anchorRt[fourthMDIndex] - mdsInGPU.anchorRt[thirdMDIndex]; - - float diffDr = fabsf(sdIn_dr - sdOut_dr)/fabs(sdIn_dr + sdOut_dr); - - runDeltaBetaIterations(betaIn, betaOut, betaAv, pt_beta, sdIn_dr, sdOut_dr, dr, lIn); - - const float betaInMMSF = (fabsf(betaInRHmin + betaInRHmax) > 0) ? (2.f * betaIn / fabsf(betaInRHmin + betaInRHmax)) : 0.; //mean value of min,max is the old betaIn - const float betaOutMMSF = (fabsf(betaOutRHmin + betaOutRHmax) > 0) ? (2.f * betaOut / fabsf(betaOutRHmin + betaOutRHmax)) : 0.; - betaInRHmin *= betaInMMSF; - betaInRHmax *= betaInMMSF; - betaOutRHmin *= betaOutMMSF; - betaOutRHmax *= betaOutMMSF; - - const float dBetaMuls = sdlThetaMulsF * 4.f / fminf(fabsf(pt_beta), SDL::pt_betaMax); //need to confirm the range-out value of 7 GeV - - const float alphaInAbsReg = fmaxf(fabsf(sdIn_alpha), asinf(fminf(rt_InLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float alphaOutAbsReg = fmaxf(fabsf(sdOut_alpha), asinf(fminf(rt_OutLo * k2Rinv1GeVf / 3.0f, sinAlphaMax))); - const float dBetaInLum = lIn < 11 ? 0.0f : fabsf(alphaInAbsReg*deltaZLum / z_InLo); - const float dBetaOutLum = lOut < 11 ? 0.0f : fabsf(alphaOutAbsReg*deltaZLum / z_OutLo); - const float dBetaLum2 = (dBetaInLum + dBetaOutLum) * (dBetaInLum + dBetaOutLum); - const float sinDPhi = sinf(dPhi); - - const float dBetaRIn2 = 0; // TODO-RH - // const float dBetaROut2 = 0; // TODO-RH - float dBetaROut2 = 0;//TODO-RH - betaOutCut = asinf(fminf(dr*k2Rinv1GeVf / ptCut, sinAlphaMax)) //FIXME: need faster version - + (0.02f / sdOut_d) + sqrtf(dBetaLum2 + dBetaMuls*dBetaMuls); - - //Cut #6: The real beta cut - pass = pass and (fabsf(betaOut) < betaOutCut); - if(not pass) return pass; - - float pt_betaIn = dr * k2Rinv1GeVf/sinf(betaIn); - float pt_betaOut = dr * k2Rinv1GeVf / sinf(betaOut); - float dBetaRes = 0.02f/fminf(sdOut_d,sdIn_d); - float dBetaCut2 = (dBetaRes*dBetaRes * 2.0f + dBetaMuls * dBetaMuls + dBetaLum2 + dBetaRIn2 + dBetaROut2 - + 0.25f * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax)) * (fabsf(betaInRHmin - betaInRHmax) + fabsf(betaOutRHmin - betaOutRHmax))); - float dBeta = betaIn - betaOut; - //Cut #7: Cut on dBeta - deltaBetaCut = sqrtf(dBetaCut2); - - pass = pass and (dBeta * dBeta <= dBetaCut2); - - return pass; -} - -__device__ void SDL::runDeltaBetaIterations(float& betaIn, float& betaOut, float& betaAv, float & pt_beta, float sdIn_dr, float sdOut_dr, float dr, float lIn) -{ - if (lIn == 0) - { - betaOut += copysign(asinf(fminf(sdOut_dr * k2Rinv1GeVf / fabsf(pt_beta), sinAlphaMax)), betaOut); - return; - } - - if (betaIn * betaOut > 0.f and (fabsf(pt_beta) < 4.f * SDL::pt_betaMax or (lIn >= 11 and fabsf(pt_beta) < 8.f * SDL::pt_betaMax))) //and the pt_beta is well-defined; less strict for endcap-endcap - { - - const float betaInUpd = betaIn + copysignf(asinf(fminf(sdIn_dr * k2Rinv1GeVf / fabsf(pt_beta), sinAlphaMax)), betaIn); //FIXME: need a faster version - const float betaOutUpd = betaOut + copysignf(asinf(fminf(sdOut_dr * k2Rinv1GeVf / fabsf(pt_beta), sinAlphaMax)), betaOut); //FIXME: need a faster version - betaAv = 0.5f * (betaInUpd + betaOutUpd); - - //1st update - //pt_beta = dr * k2Rinv1GeVf / sinf(betaAv); //get a better pt estimate - const float pt_beta_inv = 1.f/fabsf(dr * k2Rinv1GeVf / sinf(betaAv)); //get a better pt estimate - - betaIn += copysignf(asinf(fminf(sdIn_dr * k2Rinv1GeVf *pt_beta_inv, sinAlphaMax)), betaIn); //FIXME: need a faster version - betaOut += copysignf(asinf(fminf(sdOut_dr * k2Rinv1GeVf *pt_beta_inv, sinAlphaMax)), betaOut); //FIXME: need a faster version - //update the av and pt - betaAv = 0.5f * (betaIn + betaOut); - //2nd update - pt_beta = dr * k2Rinv1GeVf / sinf(betaAv); //get a better pt estimate - } - else if (lIn < 11 && fabsf(betaOut) < 0.2f * fabsf(betaIn) && fabsf(pt_beta) < 12.f * SDL::pt_betaMax) //use betaIn sign as ref - { - - const float pt_betaIn = dr * k2Rinv1GeVf / sinf(betaIn); - - const float betaInUpd = betaIn + copysignf(asinf(fminf(sdIn_dr * k2Rinv1GeVf / fabsf(pt_betaIn), sinAlphaMax)), betaIn); //FIXME: need a faster version - const float betaOutUpd = betaOut + copysignf(asinf(fminf(sdOut_dr * k2Rinv1GeVf / fabsf(pt_betaIn), sinAlphaMax)), betaIn); //FIXME: need a faster version - betaAv = (fabsf(betaOut) > 0.2f * fabsf(betaIn)) ? (0.5f * (betaInUpd + betaOutUpd)) : betaInUpd; - - //1st update - pt_beta = dr * k2Rinv1GeVf / sin(betaAv); //get a better pt estimate - betaIn += copysignf(asinf(fminf(sdIn_dr * k2Rinv1GeVf / fabsf(pt_beta), sinAlphaMax)), betaIn); //FIXME: need a faster version - betaOut += copysignf(asinf(fminf(sdOut_dr * k2Rinv1GeVf / fabsf(pt_beta), sinAlphaMax)), betaIn); //FIXME: need a faster version - //update the av and pt - betaAv = 0.5f * (betaIn + betaOut); - //2nd update - pt_beta = dr * k2Rinv1GeVf / sin(betaAv); //get a better pt estimate - - } -} - -void SDL::printTracklet(struct SDL::tracklets& trackletsInGPU, struct SDL::segments& segmentsInGPU, struct SDL::miniDoublets& mdsInGPU, struct SDL::hits& hitsInGPU, struct SDL::modules& modulesInGPU, unsigned int trackletIndex) -{ - unsigned int innerSegmentIndex = trackletsInGPU.segmentIndices[trackletIndex * 2]; - unsigned int outerSegmentIndex = trackletsInGPU.segmentIndices[trackletIndex * 2 + 1]; - - std::cout<