From b10131f78b90e31ba8aaa5a8528b028a3449430a Mon Sep 17 00:00:00 2001 From: Vincenzo Innocente Date: Tue, 31 Jul 2018 12:03:15 +0200 Subject: [PATCH] Heterogeneous ClusterTPAssociation (#105) Implement a heterogeneous Cluster-to-TrackingParticle associator running on the GPU. --- .../PixelTriplets/plugins/gpuPixelDoublets.h | 161 ++++++++++++++++++ .../RecoTrack/python/TrackValidation_cff.py | 15 +- 2 files changed, 174 insertions(+), 2 deletions(-) create mode 100644 RecoPixelVertexing/PixelTriplets/plugins/gpuPixelDoublets.h diff --git a/RecoPixelVertexing/PixelTriplets/plugins/gpuPixelDoublets.h b/RecoPixelVertexing/PixelTriplets/plugins/gpuPixelDoublets.h new file mode 100644 index 0000000000000..f09ae6aba5efb --- /dev/null +++ b/RecoPixelVertexing/PixelTriplets/plugins/gpuPixelDoublets.h @@ -0,0 +1,161 @@ +#ifndef RecoLocalTracker_SiPixelRecHits_plugins_gpuPixelDoublets_h +#define RecoLocalTracker_SiPixelRecHits_plugins_gpuPixelDouplets_h + +#include +#include +#include +#include +#include +#include + +#include "DataFormats/Math/interface/approx_atan2.h" +#include "RecoLocalTracker/SiPixelRecHits/plugins/siPixelRecHitsHeterogeneousProduct.h" + +namespace gpuPixelDoublets { + + __device__ + std::pair + findPhiLimits(int16_t phiMe, int16_t * iphi, uint16_t * index, uint16_t size, int16_t iphicut) { + + assert(iphicut>0); + + // find extreemes in top + int16_t minPhi = phiMe-iphicut; + int16_t maxPhi = phiMe+iphicut; + + // std::cout << "\n phi min/max " << phiMe << ' ' << minPhi << ' ' << maxPhi << std::endl; + + // guess and adjust + auto findLimit = [&](int16_t mPhi) { + int jm = float(0.5f*size)*(1.f+float(mPhi)/float(std::numeric_limits::max())); + // std::cout << "jm for " << mPhi << ' ' << jm << std::endl; + jm = std::min(size-1,std::max(0,jm)); + bool notDone=true; + while(jm>0 && mPhiiphi[index[++jm]]){} + jm = std::min(size-1,std::max(0,jm)); + return jm; + }; + + auto jmin = findLimit(minPhi); + auto jmax = findLimit(maxPhi); + + + /* + std::cout << "j min/max " << jmin << ' ' << jmax << std::endl; + std::cout << "found min/max " << iphi[index[jmin]] << ' ' << iphi[index[jmax]] << std::endl; + std::cout << "found min/max +1 " << iphi[index[jmin+1]] << ' ' << iphi[index[jmax+1]] << std::endl; + std::cout << "found min/max -1 " << iphi[index[jmin-1]] << ' ' << iphi[index[jmax-1]] << std::endl; + */ + + return std::make_pair(jmin,jmax); + } + + + __global__ + void getDoubletsFromSorted(int16_t * iphi, uint16_t * index, uint32_t * offsets, float phiCut) { + auto iphicut = phi2short(phiCut); + auto i = blockIdx.x*blockDim.x + threadIdx.x; + if (i>=offsets[9]) { + // get rid of last layer + return; + } + + assert(0==offsets[0]); + int top = (i>offsets[5]) ? 5: 0; + while (i>=offsets[++top]){}; + assert(top<10); + auto bottom = top-1; + if (bottom == 3 or bottom == 6) { + // do not have UP... (9 we got rid already) + return; + } + assert(i >= offsets[bottom]); + assert(i < offsets[top]); + + if (index[i]>= (offsets[top]-offsets[bottom])) { + printf("index problem: %d %d %d %d %d\n",i, offsets[top], offsets[bottom], offsets[top]-offsets[bottom], index[i]); + return; + } + + assert(index[i]::max()); + + auto jLimits = findPhiLimits(phiMe, iphi+offsets[top],index+offsets[top],size,iphicut); + + auto slidingWindow = [&](uint16_t mysize, uint16_t mymin,uint16_t mymax) { + auto topPhi = iphi+offsets[top]; + uint16_t imax = std::numeric_limits::max(); + uint16_t offset = (mymin>mymax) ? imax-(mysize-1) : 0; + int n=0; + for (uint16_t i = mymin+offset; i!=mymax; i++) { + assert(i<=imax); + uint16_t k = (i>mymax) ? i-offset : i; + assert(k=mymin || k2*iphicut && int16_t(phiMe-topPhi[k])>2*iphicut) + printf("deltaPhi problem: %d %d %d %d, deltas %d:%d cut %d\n",i,k,phiMe,topPhi[k],int16_t(topPhi[k]-phiMe),int16_t(phiMe-topPhi[k]),iphicut); + n++; + } + int tot = (mymin>mymax) ? (mysize-mymin)+mymax : mymax-mymin; + assert(n==tot); + }; + + slidingWindow(size,jLimits.first,jLimits.second); + } + + template + __device__ + void doubletsFromHisto(int16_t const * iphi, Hist const * hist, uint32_t const * offsets, float phiCut) { + auto iphicut = phi2short(phiCut); + auto i = blockIdx.x*blockDim.x + threadIdx.x; + if (i>=offsets[9]) { + // get rid of last layer + return; + } + + assert(0==offsets[0]); + int top = (i>offsets[5]) ? 5: 0; + while (i>=offsets[++top]){}; + assert(top<10); + auto bottom = top-1; + if (bottom==3 || bottom==6) { + // do not have UP... (9 we got rid already) + return; + } + assert(i>=offsets[bottom]); + assert(i iphicut) + continue; + ++tot; + } + } + if (0==hist[top].nspills) assert(tot>=nmin); + // look in spill bin as well.... + } + + __global__ + void getDoubletsFromHisto(siPixelRecHitsHeterogeneousProduct::HitsOnGPU const * hhp, float phiCut) { + auto const & hh = *hhp; + doubletsFromHisto(hh.iphi_d,hh.hist_d,hh.hitsLayerStart_d,phiCut); + } + +} // namespace end + +#endif // RecoLocalTracker_SiPixelRecHits_plugins_gpuPixelDouplets_h diff --git a/Validation/RecoTrack/python/TrackValidation_cff.py b/Validation/RecoTrack/python/TrackValidation_cff.py index b6251473fcab4..dff571862bb46 100644 --- a/Validation/RecoTrack/python/TrackValidation_cff.py +++ b/Validation/RecoTrack/python/TrackValidation_cff.py @@ -709,9 +709,14 @@ def _uniqueFirstLayers(layerList): ### Pixel tracking only mode (placeholder for now) -tpClusterProducerPixelTrackingOnly = tpClusterProducer.clone( - pixelClusterSrc = "siPixelClustersPreSplitting" + +tpClusterProducerHeterogeneousPixelTrackingOnly = tpClusterProducerHeterogeneous.clone( + pixelClusterSrc = "siPixelClustersPreSplitting" ) +tpClusterProducerPixelTrackingOnly = tpClusterProducer.clone() +# Need to use the modifier to customize because the exact EDProducer type depends on the modifier +gpu.toModify(tpClusterProducerPixelTrackingOnly, src = "tpClusterProducerHeterogeneousPixelTrackingOnly") + quickTrackAssociatorByHitsPixelTrackingOnly = quickTrackAssociatorByHits.clone( cluster2TPSrc = "tpClusterProducerPixelTrackingOnly" ) @@ -739,12 +744,18 @@ def _uniqueFirstLayers(layerList): tracksValidationTruthPixelTrackingOnly.replace(quickTrackAssociatorByHits, quickTrackAssociatorByHitsPixelTrackingOnly) tracksValidationTruthPixelTrackingOnly.replace(trackingParticleRecoTrackAsssociation, trackingParticlePixelTrackAsssociation) tracksValidationTruthPixelTrackingOnly.replace(VertexAssociatorByPositionAndTracks, PixelVertexAssociatorByPositionAndTracks) + +_tracksValidationTruthPixelTrackingOnlyGPU = tracksValidationTruthPixelTrackingOnly.copy() +_tracksValidationTruthPixelTrackingOnlyGPU.insert(0, tpClusterProducerHeterogeneousPixelTrackingOnly) +gpu.toReplaceWith(tracksValidationTruthPixelTrackingOnly, _tracksValidationTruthPixelTrackingOnlyGPU) + tracksValidationPixelTrackingOnly = cms.Sequence( tracksValidationTruthPixelTrackingOnly + trackValidatorPixelTrackingOnly ) + ### Lite mode (only generalTracks and HP) trackValidatorLite = trackValidator.clone( label = ["generalTracks", "cutsRecoTracksHp"]