diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt
index 47a74b5..0ce797c 100644
--- a/sample/CMakeLists.txt
+++ b/sample/CMakeLists.txt
@@ -2,6 +2,7 @@ SET(tolink ${GSL_LIBRARIES} CosmoTool ${CosmoTool_LIBS})
 include_directories(${CMAKE_SOURCE_DIR}/src)
 include_directories(${FFTW3_INCLUDE_DIRS}
 ${EIGEN3_INCLUDE_DIRS} ${NETCDF_INCLUDE_PATH} ${GSL_INCLUDE_PATH})
+include_directories(${CMAKE_SOURCE_DIR}/sample)
 
 IF(SHARP_INCLUDE_PATH)
   include_directories(BEFORE ${SHARP_INCLUDE_PATH})
diff --git a/sample/simple3DFilter.cpp b/sample/simple3DFilter.cpp
index 1932e97..c333f2e 100644
--- a/sample/simple3DFilter.cpp
+++ b/sample/simple3DFilter.cpp
@@ -1,3 +1,4 @@
+#include "omptl/algorithm"
 #include <cassert>
 #include "yorick.hpp"
 #include "sphSmooth.hpp"
@@ -106,7 +107,7 @@ int main(int argc, char **argv)
 
   array3_type out_den_1(boost::extents[Nres][Nres][Nres]);
   array4_type out_v3d_1(boost::extents[Nres][Nres][Nres][3]);
-  array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
+//  array3_type out_rad_1(boost::extents[Nres][Nres][Nres]);
  
   cout << "Weighing..." << endl;
 
@@ -174,7 +175,7 @@ int main(int argc, char **argv)
                     out_v3d_1[rx][ry][rz][1] = 0;
                     out_v3d_1[rx][ry][rz][2] = 0;
                     out_den_1[rx][ry][rz] = 0;
-                    out_rad_1[rx][ry][rz] = -1;
+                    //out_rad_1[rx][ry][rz] = -1;
                     continue;
                   }
 
@@ -184,9 +185,7 @@ int main(int argc, char **argv)
                 else
                   smooth1.fetchNeighbours(c);
 
-                float val;
-
-                out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen();
+                //out_rad_1[rx][ry][rz] = smooth1.getSmoothingLen();
                 out_v3d_1[rx][ry][rz][0] = smooth1.computeSmoothedValue(c, getVelocity<0>);
                 out_v3d_1[rx][ry][rz][1] = smooth1.computeSmoothedValue(c, getVelocity<1>);
                 out_v3d_1[rx][ry][rz][2] = smooth1.computeSmoothedValue(c, getVelocity<2>);
@@ -196,7 +195,7 @@ int main(int argc, char **argv)
       }
   }
       
-  hdf5_write_array(out_f, "radii", out_rad_1);
+  //hdf5_write_array(out_f, "radii", out_rad_1);
   hdf5_write_array(out_f, "velocity", out_v3d_1);
   hdf5_write_array(out_f, "density", out_den_1);
 
diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp
index bf1eaa0..c21da8e 100644
--- a/src/mykdtree.hpp
+++ b/src/mykdtree.hpp
@@ -40,6 +40,10 @@ knowledge of the CeCILL license and that you accept its terms.
 #include "config.hpp"
 #include "bqueue.hpp"
 
+#ifndef __KD_TREE_ACTIVE_CELLS
+#define __KD_TREE_ACTIVE_CELLS 1
+#endif
+
 namespace CosmoTool {
 
   template<int N, typename CType = ComputePrecision> 
@@ -52,7 +56,9 @@ namespace CosmoTool {
   template<int N, typename ValType, typename CType = ComputePrecision>
   struct KDCell
   {
+#if __KD_TREE_ACTIVE_CELLS == 1
     bool active;
+#endif
     ValType val;
     typename KDDef<N,CType>::KDCoordinates coord;
   };
diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc
index e276cad..2424a06 100644
--- a/src/mykdtree.tcc
+++ b/src/mykdtree.tcc
@@ -184,7 +184,9 @@ namespace CosmoTool {
     int axis = level % N;
     CoordType d2 = 0;
 
+#if __KD_TREE_ACTIVE_CELLS == 1
     if (node->value->active) 
+#endif
       {
 	for (int j = 0; j < 3; j++)
 	  {
@@ -229,6 +231,7 @@ namespace CosmoTool {
     uint32_t swapId = Ncells-1;
     uint32_t i = 0;
 
+#if __KD_TREE_ACTIVE_CELLS == 1
     while (!cells[swapId]->active && swapId > 0)
       swapId--;
 
@@ -244,6 +247,7 @@ namespace CosmoTool {
 	  }
 	i++;
       }
+#endif
     return swapId+1;
   }
 
@@ -251,7 +255,7 @@ namespace CosmoTool {
   void  KD_default_cell_splitter<N,ValType,CType>::operator()(KDCell<N,ValType,CType> **cells, uint32_t Ncells, uint32_t& split_index, int axis, typename KDDef<N,CType>::KDCoordinates minBound, typename KDDef<N,CType>::KDCoordinates maxBound)
   {
     CellCompare<N,ValType,CType> compare(axis);
-    std::sort(cells, cells+Ncells, compare);
+    omptl::sort(cells,cells+Ncells,compare); // std::sort(cells, cells+Ncells, compare);
     split_index = Ncells/2;
   }
 
@@ -266,10 +270,13 @@ namespace CosmoTool {
     if (Ncells == 0)
       return 0;
 
+    Node *node;
     int axis = depth % N;
-    Node *node = &nodes[lastNode++];
     uint32_t mid;
     coords tmpBound;
+
+#pragma omp critical
+    node = &nodes[lastNode++];
     
     // Isolate the environment
     splitter(cell0, Ncells, mid, axis, minBound, maxBound);
@@ -282,12 +289,20 @@ namespace CosmoTool {
     tmpBound[axis] = node->value->coord[axis];
 
     depth++;
-    node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
+#pragma omp task private(tmpBound)
+    {
+      node->children[0] = buildTree(cell0, mid, depth, minBound, tmpBound);
+    }
     
     memcpy(tmpBound, minBound, sizeof(coords));
     tmpBound[axis] = node->value->coord[axis];
-    node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth, 
+#pragma omp task private(tmpBound)
+    {
+      node->children[1] = buildTree(cell0+mid+1, Ncells-mid-1, depth, 
 				  tmpBound, maxBound);
+    }
+
+#pragma omp taskwait
 
 #ifdef __KD_TREE_NUMNODES
     node->numNodes = (node->children[0] != 0) ? node->children[0]->numNodes : 0;
@@ -304,7 +319,9 @@ namespace CosmoTool {
     uint32_t numActive = 0;
     for (uint32_t i = 0; i < lastNode; i++)
       {
+#if __KD_TREE_ACTIVE_CELLS == 1
 	if (nodes[i].value->active)
+#endif
 	  numActive++;
       }
     return numActive;