diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp
index e87b376..9d83d5b 100644
--- a/src/mykdtree.hpp
+++ b/src/mykdtree.hpp
@@ -7,20 +7,19 @@
 
 namespace CosmoTool {
 
-  template<int N> 
+  template<int N, typename CType = ComputePrecision> 
   struct KDDef
   {
-    typedef float CoordType;
+    typedef CType CoordType;
     typedef float KDCoordinates[N];
-    static const int NumCubes = 1 << (N);
   };
   
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType = ComputePrecision>
   struct KDCell
   {
     bool active;
     ValType val;
-    typename KDDef<N>::KDCoordinates coord;
+    typename KDDef<N,CType>::KDCoordinates coord;
   };
 
   class NotEnoughCells: public Exception
@@ -30,114 +29,113 @@ namespace CosmoTool {
     ~NotEnoughCells() throw () {}
   };
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType = ComputePrecision>
   struct KDTreeNode
   {
-    KDCell<N,ValType> *value;
-    KDTreeNode<N,ValType> *children[2];
-    typename KDDef<N>::KDCoordinates minBound, maxBound;
+    KDCell<N,ValType,CType> *value;
+    KDTreeNode<N,ValType,CType> *children[2];
+    typename KDDef<N,CType>::KDCoordinates minBound, maxBound;
   };
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType = ComputePrecision>
   class RecursionInfoCells
   {
   public:
 
-    typename KDDef<N>::KDCoordinates x;
-    typename KDDef<N>::CoordType r, r2;
-    KDCell<N, ValType> **cells;
-    typename KDDef<N>::CoordType *distances;
+    typename KDDef<N,CType>::KDCoordinates x;
+    typename KDDef<N,CType>::CoordType r, r2;
+    KDCell<N, ValType,CType> **cells;
+    typename KDDef<N,CType>::CoordType *distances;
     uint32_t currentRank;
     uint32_t numCells;
   };
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType = ComputePrecision>
   class RecursionMultipleInfo
   {
   public:
-    const typename KDDef<N>::KDCoordinates& x;
-    BoundedQueue< KDCell<N,ValType> *, typename KDDef<N>::CoordType> queue;
+    const typename KDDef<N,CType>::KDCoordinates& x;
+    BoundedQueue< KDCell<N,ValType,CType> *, typename KDDef<N,CType>::CoordType> queue;
     int traversed;
 
-    RecursionMultipleInfo(const typename KDDef<N>::KDCoordinates& rx,
-			  KDCell<N,ValType> **cells,
+    RecursionMultipleInfo(const typename KDDef<N,CType>::KDCoordinates& rx,
+			  KDCell<N,ValType,CType> **cells,
 			  uint32_t numCells)
       : x(rx), queue(cells, numCells, INFINITY),traversed(0)
     {      
     }
   };
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType = ComputePrecision>
   class KDTree
   {
   public:
-    static const int NumCubes = KDDef<N>::NumCubes;
-
-  public:
-    typedef typename KDDef<N>::CoordType CoordType;
+    typedef typename KDDef<N,CType>::CoordType CoordType;
     typedef typename KDDef<N>::KDCoordinates coords;
+    typedef KDCell<N,ValType,CType> Cell;
+    typedef KDTreeNode<N,ValType,CType> Node;
 
-    KDTree(KDCell<N,ValType> *cells, uint32_t Ncells);
+    KDTree(Cell *cells, uint32_t Ncells);
     ~KDTree();
 
     uint32_t getIntersection(const coords& x, CoordType r, 
-			     KDCell<N, ValType> **cells,
+			     Cell **cells,
 			     uint32_t numCells)
       throw (NotEnoughCells);
     uint32_t getIntersection(const coords& x, CoordType r, 
-			     KDCell<N, ValType> **cells,
+			     Cell **cells,
 			     CoordType *distances,
 			     uint32_t numCells)
       throw (NotEnoughCells);
 
-    KDCell<N, ValType> *getNearestNeighbour(const coords& x);
+    Cell *getNearestNeighbour(const coords& x);
 
     void getNearestNeighbours(const coords& x, uint32_t N,
-			      KDCell<N, ValType> **cells);
+			      Cell **cells);
     void getNearestNeighbours(const coords& x, uint32_t N,
-			      KDCell<N, ValType> **cells,
+			      Cell **cells,
 			      CoordType *distances);
 
-    KDTreeNode<N,ValType> *getRoot() { return root; }
+    Node *getRoot() { return root; }
 
     void optimize();
 
-    KDTreeNode<N,ValType> *getAllNodes() { return nodes; }
+    Node *getAllNodes() { return nodes; }
     uint32_t getNumNodes() const { return lastNode; }
 
     uint32_t countActives() const;
 
   protected:
-    KDTreeNode<N, ValType> *nodes;
+    Node *nodes;
     uint32_t numNodes;
     uint32_t lastNode;
 
-    KDTreeNode<N, ValType> *root;
-    KDCell<N, ValType> **sortingHelper;
+    Node *root;
+    Cell **sortingHelper;
 
-    KDTreeNode<N, ValType> *buildTree(KDCell<N,ValType> **cell0,
-				       uint32_t N,
-				       uint32_t depth,
-				       coords minBound,
-				       coords maxBound);
+    Node *buildTree(Cell **cell0,
+		    uint32_t N,
+		    uint32_t depth,
+		    coords minBound,
+		    coords maxBound);
     
-    void recursiveIntersectionCells(RecursionInfoCells<N,ValType>& info,
-				    KDTreeNode<N,ValType> *node,
+    void recursiveIntersectionCells(RecursionInfoCells<N,ValType, CType>& info,
+				    Node *node,
 				    int level)
       throw (NotEnoughCells);
 
-    CoordType computeDistance(KDCell<N,ValType> *cell, const coords& x);
-    void recursiveNearest(KDTreeNode<N, ValType> *node,
+    CoordType computeDistance(Cell *cell, const coords& x);
+    void recursiveNearest(Node *node,
 			  int level,
 			  const coords& x,
 			  CoordType& R2,
-			  KDCell<N,ValType>*& cell);
-    void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType>& info, KDTreeNode<N,ValType> *node,
+			  Cell*& cell);
+    void recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
 				  int level);    
   };
 
-  template<int N, class T>
-  uint32_t gatherActiveCells(KDCell<N,T> **cells, uint32_t numCells);
+  template<int N, typename T, typename CType>
+  uint32_t gatherActiveCells(KDCell<N,T,CType> **cells, uint32_t numCells);
 
 };
 
diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc
index 66ec29e..183e48e 100644
--- a/src/mykdtree.tcc
+++ b/src/mykdtree.tcc
@@ -5,7 +5,7 @@
 
 namespace CosmoTool {
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType>
   class CellCompare
   {
   public:
@@ -14,7 +14,7 @@ namespace CosmoTool {
       rank = k;
     }
 
-    bool operator()(const KDCell<N,ValType> *a, const KDCell<N,ValType> *b) const
+    bool operator()(const KDCell<N,ValType,CType> *a, const KDCell<N,ValType,CType> *b) const
     {
       return (a->coord[rank] < b->coord[rank]);
     }
@@ -22,27 +22,27 @@ namespace CosmoTool {
     int rank;
   };
 
-  template<int N, typename ValType>
-  KDTree<N,ValType>::~KDTree()
+  template<int N, typename ValType, typename CType>
+  KDTree<N,ValType,CType>::~KDTree()
   {
   }
 
-  template<int N, typename ValType>
-  KDTree<N,ValType>::KDTree(KDCell<N,ValType> *cells, uint32_t Ncells)
+  template<int N, typename ValType, typename CType>
+  KDTree<N,ValType,CType>::KDTree(Cell *cells, uint32_t Ncells)
   {
     
     numNodes = Ncells;
-    nodes = new KDTreeNode<N,ValType>[numNodes];
+    nodes = new Node[numNodes];
 
-    sortingHelper = new KDCell<N,ValType> *[Ncells];
+    sortingHelper = new Cell *[Ncells];
     for (uint32_t i = 0; i < Ncells; i++)
 	sortingHelper[i] = &cells[i];    
 
     optimize();
   }
 
-  template<int N, typename ValType>
-  void KDTree<N,ValType>::optimize()
+  template<int N, typename ValType, typename CType>
+  void KDTree<N,ValType,CType>::optimize()
   {
     coords absoluteMin, absoluteMax;
 
@@ -62,13 +62,13 @@ namespace CosmoTool {
     std::cout << "   done." << std::endl;
   }
 
-  template<int N, typename ValType>
-  uint32_t KDTree<N,ValType>::getIntersection(const coords& x, CoordType r, 
-					       KDCell<N, ValType> **cells,
-					       uint32_t numCells)
+  template<int N, typename ValType, typename CType>
+  uint32_t KDTree<N,ValType,CType>::getIntersection(const coords& x, CoordType r, 
+						    KDTree<N,ValType,CType>::Cell **cells,
+						    uint32_t numCells)
     throw (NotEnoughCells)
   {
-    RecursionInfoCells<N,ValType> info;
+    RecursionInfoCells<N,ValType,CType> info;
 
     memcpy(info.x, x, sizeof(x));
     info.r = r;
@@ -82,11 +82,11 @@ namespace CosmoTool {
     return info.currentRank;
   }
 
-  template<int N, typename ValType>
-  uint32_t KDTree<N,ValType>::getIntersection(const coords& x, CoordType r, 
-					      KDCell<N, ValType> **cells,
-					      CoordType *distances,
-					      uint32_t numCells)
+  template<int N, typename ValType, typename CType>
+  uint32_t KDTree<N,ValType,CType>::getIntersection(const coords& x, CoordType r, 
+						    Cell **cells,
+						    CoordType *distances,
+						    uint32_t numCells)
     throw (NotEnoughCells)
   {
     RecursionInfoCells<N,ValType> info;
@@ -103,10 +103,10 @@ namespace CosmoTool {
     return info.currentRank;
   }
 
-  template<int N, typename ValType>
-  void KDTree<N,ValType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType>& info, 
-						      KDTreeNode<N,ValType> *node,
-						      int level)
+  template<int N, typename ValType, typename CType>
+  void KDTree<N,ValType,CType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info, 
+							   Node *node,
+							   int level)
     throw (NotEnoughCells)
   {
     int axis = level % N;
@@ -147,9 +147,9 @@ namespace CosmoTool {
       }
   }
 
-  template<int N, typename ValType>
-  uint32_t gatherActiveCells(KDCell<N,ValType> **cells,
-			 uint32_t Ncells)
+  template<int N, typename ValType, typename CType>
+  uint32_t gatherActiveCells(KDCell<N,ValType,CType> **cells,
+			     uint32_t Ncells)
   {
     uint32_t swapId = Ncells-1;
     uint32_t i = 0;
@@ -172,24 +172,24 @@ namespace CosmoTool {
     return swapId+1;
   }
 
-  template<int N, typename ValType>
-  KDTreeNode<N, ValType> *KDTree<N,ValType>::buildTree(KDCell<N,ValType> **cell0,
-							 uint32_t Ncells,
-							 uint32_t depth,
-							 coords minBound,
-							 coords maxBound)
+  template<int N, typename ValType, typename CType>
+  KDTreeNode<N,ValType,CType> *KDTree<N,ValType,CType>::buildTree(Cell **cell0,
+								    uint32_t Ncells,
+								    uint32_t depth,
+								    coords minBound,
+								    coords maxBound)
   {
     if (Ncells == 0)
       return 0;
 
     int axis = depth % N;
-    KDTreeNode<N,ValType> *node = &nodes[lastNode++];
+    Node *node = &nodes[lastNode++];
     uint32_t mid = Ncells/2;
     coords tmpBound;
     
     // Isolate the environment
     {
-      CellCompare<N,ValType> compare(axis);
+      CellCompare<N,ValType,CType> compare(axis);
       std::sort(cell0, cell0+Ncells, compare);
     }
 
@@ -211,8 +211,8 @@ namespace CosmoTool {
     return node;
   }
 
-  template<int N, typename ValType>
-  uint32_t KDTree<N,ValType>::countActives() const
+  template<int N, typename ValType, typename CType>
+  uint32_t KDTree<N,ValType,CType>::countActives() const
   {
     uint32_t numActive = 0;
     for (uint32_t i = 0; i < lastNode; i++)
@@ -223,9 +223,9 @@ namespace CosmoTool {
     return numActive;
   }
 
-  template<int N, typename ValType>
-  typename KDDef<N>::CoordType
-  KDTree<N,ValType>::computeDistance(KDCell<N, ValType> *cell, const coords& x)
+  template<int N, typename ValType, typename CType>
+  typename KDDef<N,CType>::CoordType
+  KDTree<N,ValType,CType>::computeDistance(Cell *cell, const coords& x)
   {
     CoordType d2 = 0;
 
@@ -237,18 +237,18 @@ namespace CosmoTool {
     return d2;
   }
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType>
   void
-  KDTree<N,ValType>::recursiveNearest(
-				       KDTreeNode<N,ValType> *node,
+  KDTree<N,ValType,CType>::recursiveNearest(
+				       Node *node,
 				       int level,
 				       const coords& x,
 				       CoordType& R2,
-				       KDCell<N,ValType> *& best)
+				       Cell *& best)
   {
     CoordType d2 = 0;
     int axis = level % N;
-    KDTreeNode<N,ValType> *other, *go;
+    Node *other, *go;
 
     if (x[axis] < node->value->coord[axis])
       {
@@ -311,26 +311,26 @@ namespace CosmoTool {
       }
   }
 
- template<int N, typename ValType>
- KDCell<N, ValType> *
- KDTree<N,ValType>::getNearestNeighbour(const coords& x)
+ template<int N, typename ValType, typename CType>
+ KDCell<N,ValType,CType> *
+ KDTree<N,ValType,CType>::getNearestNeighbour(const coords& x)
  {
    CoordType R2 = INFINITY;
-   KDCell<N,ValType> *best = 0;
+   Cell *best = 0;
 
    recursiveNearest(root, 0, x, R2, best);   
    
    return best;
  }
 
-  template<int N, typename ValType>
+  template<int N, typename ValType, typename CType>
   void
-  KDTree<N,ValType>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType>& info, KDTreeNode<N,ValType> *node,
+  KDTree<N,ValType,CType>::recursiveMultipleNearest(RecursionMultipleInfo<N,ValType,CType>& info, Node *node,
 					       int level)
   {
     CoordType d2 = 0;
     int axis = level % N;
-    KDTreeNode<N,ValType> *other, *go;
+    Node *other, *go;
 
     if (info.x[axis] < node->value->coord[axis])
       {
@@ -380,9 +380,9 @@ namespace CosmoTool {
       }    
   }
 
- template<int N, typename ValType>
- void KDTree<N,ValType>::getNearestNeighbours(const coords& x, uint32_t N2,
-					       KDCell<N, ValType> **cells)
+ template<int N, typename ValType, typename CType>
+ void KDTree<N,ValType,CType>::getNearestNeighbours(const coords& x, uint32_t N2,
+						    Cell **cells)
  {
    RecursionMultipleInfo<N,ValType> info(x, cells, N2);
    
@@ -391,13 +391,13 @@ namespace CosmoTool {
 
    recursiveMultipleNearest(info, root, 0);
 
-   std::cout << "Traversed = " << info.traversed << std::endl;
+   //   std::cout << "Traversed = " << info.traversed << std::endl;
  }
 
- template<int N, typename ValType>
- void KDTree<N,ValType>::getNearestNeighbours(const coords& x, uint32_t N2,
-					      KDCell<N, ValType> **cells,
-					      CoordType *distances)
+ template<int N, typename ValType, typename CType>
+ void KDTree<N,ValType,CType>::getNearestNeighbours(const coords& x, uint32_t N2,
+						    Cell **cells,
+						    CoordType *distances)
  {
    RecursionMultipleInfo<N,ValType> info(x, cells, N2);
    
@@ -405,9 +405,9 @@ namespace CosmoTool {
      cells[i] = 0;
 
    recursiveMultipleNearest(info, root, 0);
-   memcpy(distances, info.getPriorities(), sizeof(CoordType)*N2);
+   memcpy(distances, info.queue.getPriorities(), sizeof(CoordType)*N2);
 
-   std::cout << "Traversed = " << info.traversed << std::endl;
+   //   std::cout << "Traversed = " << info.traversed << std::endl;
  }
 
 };