From 8c2709492fc8f2a5a36b737d01e9d4e4dfa9dea1 Mon Sep 17 00:00:00 2001 From: cmsbuild Date: Wed, 12 Sep 2018 09:47:46 +0200 Subject: [PATCH] Add new concurrent vertex algo (cms-patatrack#158) Add a new concurrent vertex algorithm based on DBSCAN (that reuses parts of the pixel cluster algorithm). It still needs to be tuned to reach at least the performance of current DivisiveClustering algorithm used at HLT. --- .../python/RecoPixelVertexing_cff.py | 7 +- .../python/PixelVertexes_cfi.py | 3 + .../PixelVertexFinding/test/BuildFile.xml | 34 ++- .../test/gpuVertexFinder_t.cu | 207 ++++++++++++++++++ 4 files changed, 238 insertions(+), 13 deletions(-) create mode 100644 RecoPixelVertexing/PixelVertexFinding/test/gpuVertexFinder_t.cu diff --git a/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py b/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py index 34ee6fadb04de..5f541dd19a412 100644 --- a/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py +++ b/RecoPixelVertexing/Configuration/python/RecoPixelVertexing_cff.py @@ -4,7 +4,6 @@ # # for STARTUP ONLY use try and use Offline 3D PV from pixelTracks, with adaptive vertex # -#from RecoPixelVertexing.PixelVertexFinding.PixelVertexes_cff import * -from RecoVertex.PrimaryVertexProducer.OfflinePixel3DPrimaryVertices_cfi import * -recopixelvertexingTask = cms.Task(pixelTracksTask,pixelVertices) -recopixelvertexing = cms.Sequence(recopixelvertexingTask) +from RecoPixelVertexing.PixelVertexFinding.PixelVertexes_cff import * +# from RecoVertex.PrimaryVertexProducer.OfflinePixel3DPrimaryVertices_cfi import * +recopixelvertexing = cms.Sequence(pixelTracksSequence*pixelVertices) diff --git a/RecoPixelVertexing/PixelVertexFinding/python/PixelVertexes_cfi.py b/RecoPixelVertexing/PixelVertexFinding/python/PixelVertexes_cfi.py index 77a9f367b9d9b..ea9e4b1e4e037 100644 --- a/RecoPixelVertexing/PixelVertexFinding/python/PixelVertexes_cfi.py +++ b/RecoPixelVertexing/PixelVertexFinding/python/PixelVertexes_cfi.py @@ -20,3 +20,6 @@ ) +from Configuration.ProcessModifiers.gpu_cff import gpu +from RecoPixelVertexing.PixelVertexFinding.pixelVertexHeterogeneousProducer_cfi import pixelVertexHeterogeneousProducer as _pixelVertexHeterogeneousProducer +gpu.toReplaceWith(pixelVertices, _pixelVertexHeterogeneousProducer) diff --git a/RecoPixelVertexing/PixelVertexFinding/test/BuildFile.xml b/RecoPixelVertexing/PixelVertexFinding/test/BuildFile.xml index 0f4f4dee63832..ad1f03999fbea 100644 --- a/RecoPixelVertexing/PixelVertexFinding/test/BuildFile.xml +++ b/RecoPixelVertexing/PixelVertexFinding/test/BuildFile.xml @@ -1,9 +1,25 @@ - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/RecoPixelVertexing/PixelVertexFinding/test/gpuVertexFinder_t.cu b/RecoPixelVertexing/PixelVertexFinding/test/gpuVertexFinder_t.cu new file mode 100644 index 0000000000000..f47c4362503ae --- /dev/null +++ b/RecoPixelVertexing/PixelVertexFinding/test/gpuVertexFinder_t.cu @@ -0,0 +1,207 @@ +#include +#include +#include +#include + +#include "RecoPixelVertexing/PixelVertexFinding/src/gpuClusterTracks.h" +using namespace gpuVertexFinder; +#include + + +struct Event { + std::vector zvert; + std::vector itrack; + std::vector ztrack; + std::vector eztrack; + std::vector ivert; +}; + +struct ClusterGenerator { + + explicit ClusterGenerator(float nvert, float ntrack) : + rgen(-13.,13), errgen(0.005,0.025), clusGen(nvert), trackGen(ntrack), gauss(0.,1.) + {} + + void operator()(Event & ev) { + + int nclus = clusGen(reng); + ev.zvert.resize(nclus); + ev.itrack.resize(nclus); + for (auto & z : ev.zvert) { + z = 3.5f*gauss(reng); + } + + ev.ztrack.clear(); + ev.eztrack.clear(); + ev.ivert.clear(); + for (int iv=0; iv rgen; + std::uniform_real_distribution errgen; + std::poisson_distribution clusGen; + std::poisson_distribution trackGen; + std::normal_distribution gauss; + + +}; + + +#include + +int main() { + + if (cuda::device::count() == 0) { + std::cerr << "No CUDA devices on this system" << "\n"; + exit(EXIT_FAILURE); + } + + auto current_device = cuda::device::current::get(); + + auto zt_d = cuda::memory::device::make_unique(current_device, 64000); + auto ezt2_d = cuda::memory::device::make_unique(current_device, 64000); + auto zv_d = cuda::memory::device::make_unique(current_device, 256); + auto wv_d = cuda::memory::device::make_unique(current_device, 256); + auto chi2_d = cuda::memory::device::make_unique(current_device, 256); + + auto izt_d = cuda::memory::device::make_unique(current_device, 64000); + auto nn_d = cuda::memory::device::make_unique(current_device, 64000); + auto iv_d = cuda::memory::device::make_unique(current_device, 64000); + + auto nv_d = cuda::memory::device::make_unique(current_device, 1); + + auto onGPU_d = cuda::memory::device::make_unique(current_device, 1); + + OnGPU onGPU; + + onGPU.zt = zt_d.get(); + onGPU.ezt2 = ezt2_d.get(); + onGPU.zv = zv_d.get(); + onGPU.wv = wv_d.get(); + onGPU.chi2 = chi2_d.get(); + onGPU.nv = nv_d.get(); + onGPU.izt = izt_d.get(); + onGPU.nn = nn_d.get(); + onGPU.iv = iv_d.get(); + + + cuda::memory::copy(onGPU_d.get(), &onGPU, sizeof(OnGPU)); + + + Event ev; + + for (int nav=30;nav<80;nav+=20){ + + ClusterGenerator gen(nav,10); + + for (int i=8; i<20; ++i) { + + auto kk=i/4; // M param + + gen(ev); + + std::cout << ev.zvert.size() << ' ' << ev.ztrack.size() << std::endl; + + cuda::memory::copy(onGPU.zt,ev.ztrack.data(),sizeof(float)*ev.ztrack.size()); + cuda::memory::copy(onGPU.ezt2,ev.eztrack.data(),sizeof(float)*ev.eztrack.size()); + + float eps = 0.1f; + + std::cout << "M eps " << kk << ' ' << eps << std::endl; + + if ( (i%4) == 0 ) + cuda::launch(clusterTracks, + { 1, 1024 }, + ev.ztrack.size(), onGPU_d.get(),kk,eps, + 0.02f,12.0f + ); + + if ( (i%4) == 1 ) + cuda::launch(clusterTracks, + { 1, 1024 }, + ev.ztrack.size(), onGPU_d.get(),kk,eps, + 0.02f,9.0f + ); + + if ( (i%4) == 2 ) + cuda::launch(clusterTracks, + { 1, 1024 }, + ev.ztrack.size(), onGPU_d.get(),kk,eps, + 0.01f,9.0f + ); + + if ( (i%4) == 3 ) + cuda::launch(clusterTracks, + { 1, 1024 }, + ev.ztrack.size(), onGPU_d.get(),kk,0.7f*eps, + 0.01f,9.0f + ); + + + + uint32_t nv; + cuda::memory::copy(&nv, onGPU.nv, sizeof(uint32_t)); + float zv[nv]; + float wv[nv]; + float chi2[nv]; + int32_t nn[nv]; + cuda::memory::copy(&zv, onGPU.zv, nv*sizeof(float)); + cuda::memory::copy(&wv, onGPU.wv, nv*sizeof(float)); + cuda::memory::copy(&chi2, onGPU.chi2, nv*sizeof(float)); + cuda::memory::copy(&nn, onGPU.nn, nv*sizeof(int32_t)); + for (auto j=0U; j0) chi2[j]/=float(nn[j]); + + { + auto mx = std::minmax_element(wv,wv+nv); + std::cout << "min max error " << 1./std::sqrt(*mx.first) << ' ' << 1./std::sqrt(*mx.second) << std::endl; + } + { + auto mx = std::minmax_element(chi2,chi2+nv); + std::cout << "min max chi2 " << *mx.first << ' ' << *mx.second << std::endl; + } + + + float dd[nv]; + uint32_t ii=0; + for (auto zr : zv) { + auto md=500.0f; + for (auto zs : ev.ztrack) { + auto d = std::abs(zr-zs); + md = std::min(d,md); + } + dd[ii++] = md; + } + assert(ii==nv); + if (i==6) { + for (auto d:dd) std::cout << d << ' '; + std::cout << std::endl; + } + auto mx = std::minmax_element(dd,dd+nv); + float rms=0; + for (auto d:dd) rms+=d*d; rms = std::sqrt(rms)/(nv-1); + std::cout << "min max rms " << *mx.first << ' ' << *mx.second << ' ' << rms << std::endl; + + } // loop on events + } // lopp on ave vert + + return 0; +}