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;
+}