diff --git a/bin/desi_zfind.py b/bin/desi_zfind.py
index db78da6e2..29cf08a5f 100755
--- a/bin/desi_zfind.py
+++ b/bin/desi_zfind.py
@@ -11,7 +11,7 @@
 from desispec import io
 from desispec.interpolation import resample_flux
 from desispec.log import get_logger
-from desispec.zfind import RedMonsterZfind
+from desispec.zfind.redmonster import RedMonsterZfind
 
 import optparse
 
diff --git a/py/desispec/zfind/__init__.py b/py/desispec/zfind/__init__.py
new file mode 100644
index 000000000..680631d4f
--- /dev/null
+++ b/py/desispec/zfind/__init__.py
@@ -0,0 +1 @@
+from .zfind import ZfindBase
\ No newline at end of file
diff --git a/py/desispec/zfind.py b/py/desispec/zfind/redmonster.py
similarity index 70%
rename from py/desispec/zfind.py
rename to py/desispec/zfind/redmonster.py
index b963fef16..20536f612 100644
--- a/py/desispec/zfind.py
+++ b/py/desispec/zfind/redmonster.py
@@ -1,3 +1,5 @@
+from __future__ import division, absolute_import
+
 import os
 
 import numpy as np
@@ -5,67 +7,9 @@
 from redmonster.physics.zfitter import Zfitter
 from redmonster.physics.zpicker import Zpicker
 
+from desispec.zfind import ZfindBase
 from desispec.interpolation import resample_flux
 
-class ZfindBase(object):
-    def __init__(self, wave, flux, ivar, R=None, results=None):
-        """
-        Base class of classification / redshift finders.
-        
-        Args:
-            wave : 1D[nwave] wavelength grid [Angstroms]
-            flux : 2D[nspec, nwave] flux [erg/s/cm2/A]
-            ivar : 2D[nspec, nwave] inverse variance of flux
-            
-        Optional:
-            R : 1D[nspec] list of resolution objects
-            results : ndarray with keys such as z, zerr, zwarn (see below)
-                all results.dtype.names are added to this object
-            
-        Subclasses should perform classification and redshift fitting
-        upon initialization and set the following member variables:
-            nspec : number of spectra
-            nwave : number of wavelegths (may be resampled from input)
-            z     : 1D[nspec] best fit redshift
-            zerr  : 1D[nspec] redshift uncertainty estimate
-            zwarn : 1D[nspec] integer redshift warning bitmask (details TBD)
-            type  : 1D[nspec] classification [GALAXY, QSO, STAR, ...]
-            subtype : 1D[nspec] sub-classification 
-            wave  : 1D[nwave] wavelength grid used; may be resampled from input
-            flux  : 2D[nspec, nwave] flux used; may be resampled from input
-            ivar  : 2D[nspec, nwave] ivar of flux
-            model : 2D[nspec, nwave] best fit model
-            
-            chi2?
-            zbase?
-            
-        For the purposes of I/O, it is possible to create a ZfindBase
-        object that contains only the results, without the input
-        wave, flux, ivar, or output model.
-        """
-        #- Inputs
-        if flux is not None:
-            nspec, nwave = flux.shape
-            self.nspec = nspec
-            self.nwave = nwave
-            self.wave = wave
-            self.flux = flux
-            self.ivar = ivar
-            self.R = R
-        
-        #- Outputs to fill
-        if results is None:
-            self.model = np.zeros((nspec, nwave), dtype=flux.dtype)
-            self.z = np.zeros(nspec)
-            self.zerr = np.zeros(nspec)
-            self.zwarn = np.zeros(nspec, dtype=int)
-            self.type = np.zeros(nspec, dtype='S20')
-            self.subtype = np.zeros(nspec, dtype='S20')
-        else:
-            for key in results.dtype.names:
-                self.__setattr__(key.lower(), results[key])
-    
-
 class RedMonsterZfind(ZfindBase):
     def __init__(self, wave, flux, ivar, R=None, dloglam=1e-4):
         """
diff --git a/py/desispec/zfind/zfind.py b/py/desispec/zfind/zfind.py
new file mode 100644
index 000000000..1b7dd6f35
--- /dev/null
+++ b/py/desispec/zfind/zfind.py
@@ -0,0 +1,60 @@
+import numpy as np
+
+class ZfindBase(object):
+    def __init__(self, wave, flux, ivar, R=None, results=None):
+        """
+        Base class of classification / redshift finders.
+        
+        Args:
+            wave : 1D[nwave] wavelength grid [Angstroms]
+            flux : 2D[nspec, nwave] flux [erg/s/cm2/A]
+            ivar : 2D[nspec, nwave] inverse variance of flux
+            
+        Optional:
+            R : 1D[nspec] list of resolution objects
+            results : ndarray with keys such as z, zerr, zwarn (see below)
+                all results.dtype.names are added to this object
+            
+        Subclasses should perform classification and redshift fitting
+        upon initialization and set the following member variables:
+            nspec : number of spectra
+            nwave : number of wavelegths (may be resampled from input)
+            z     : 1D[nspec] best fit redshift
+            zerr  : 1D[nspec] redshift uncertainty estimate
+            zwarn : 1D[nspec] integer redshift warning bitmask (details TBD)
+            type  : 1D[nspec] classification [GALAXY, QSO, STAR, ...]
+            subtype : 1D[nspec] sub-classification 
+            wave  : 1D[nwave] wavelength grid used; may be resampled from input
+            flux  : 2D[nspec, nwave] flux used; may be resampled from input
+            ivar  : 2D[nspec, nwave] ivar of flux
+            model : 2D[nspec, nwave] best fit model
+            
+            chi2?
+            zbase?
+            
+        For the purposes of I/O, it is possible to create a ZfindBase
+        object that contains only the results, without the input
+        wave, flux, ivar, or output model.
+        """
+        #- Inputs
+        if flux is not None:
+            nspec, nwave = flux.shape
+            self.nspec = nspec
+            self.nwave = nwave
+            self.wave = wave
+            self.flux = flux
+            self.ivar = ivar
+            self.R = R
+        
+        #- Outputs to fill
+        if results is None:
+            self.model = np.zeros((nspec, nwave), dtype=flux.dtype)
+            self.z = np.zeros(nspec)
+            self.zerr = np.zeros(nspec)
+            self.zwarn = np.zeros(nspec, dtype=int)
+            self.type = np.zeros(nspec, dtype='S20')
+            self.subtype = np.zeros(nspec, dtype='S20')
+        else:
+            for key in results.dtype.names:
+                self.__setattr__(key.lower(), results[key])
+