diff --git a/azint.cpp b/azint.cpp
index 7d9b339..96fae8a 100644
--- a/azint.cpp
+++ b/azint.cpp
@@ -80,121 +80,6 @@ void rotation_matrix(float rot[3][3], Poni poni)
     matrix_multiplication(rot, tmp, rot1);
 }
 
-void
-generateQXQYMatrix(const Poni& poni,
-		   const py::array_t<float>& pixel_corners,
-		   int n_splitting, 
-		   std::vector<RListMatrix>& segments,
-		   const int8_t* mask,
-		   int nX, const float* qx_bins,
-		   int nY, const float* qy_bins)
-{
-  float rot[3][3];
-  rotation_matrix(rot, poni);
-  // h, w, corner index [A, B, C, D], coordinates [z, y, x]
-  // A D
-  // B C
-  auto pc = pixel_corners.unchecked<4>();
-  auto shape = pixel_corners.shape();
-
-#pragma omp parallel for schedule(static)
-  
-  for (int i=0; i<shape[0]; i++)
-    {
-      int rank = omp_get_thread_num();
-
-      for (int j=0; j<shape[1]; j++)
-	{
-	  int pixel_index = i*shape[1] + j;
-	  if (mask[pixel_index]) continue;
-
-	  // coordinate of A1/A2/A3
-	  float A1 = pc(i, j, 0, 1);
-	  float A2 = pc(i, j, 0, 2);
-	  float A3 = pc(i, j, 0, 0);
-	  
-	  float BA1 = pc(i, j, 1, 1) - A1;
-	  float BA2 = pc(i, j, 1, 2) - A2;
-	  float BA3 = pc(i, j, 1, 0) - A3;
-          
-	  float DA1 = pc(i, j, 3, 1) - A1;
-	  float DA2 = pc(i, j, 3, 2) - A2;
-	  float DA3 = pc(i, j, 3, 0) - A3;
-          
-	  for (int k=0; k<n_splitting; k++)
-	    {
-	      float delta1 = (k + 0.5) / n_splitting;
-	      for (int l=0; l<n_splitting; l++)
-		{
-		  float delta2 = (l + 0.5) / n_splitting;
-		  float p[] =
-		    {
-		      A1 + delta1 * BA1 + delta2 * DA1 - poni.poni1,
-		      A2 + delta1 * BA2 + delta2 * DA2 - poni.poni2,
-		      A3 + delta1 * BA3 + delta2 * DA3 + poni.dist
-                    };
-		  float pos[3];
-		  dot(pos, rot, p);
-
-		  float r = sqrtf(pos[0]*pos[0] + pos[1]*pos[1]);
-		  float tth = atan2f(r, pos[2]);
-		  float phi = atan2f(-pos[0], -pos[1]);
-		  float q=4.0e-10 * M_PI / poni.wavelength * sinf(0.5*tth);
-		  float qxValue=q * sinf(phi);
-		  float qyValue=q * cosf(phi);
-
-		  /*
-		    float phiDeg=180.0*phi/M_PI;
-		  if (((flag & 1) && phiDeg>90.0) ||
-		      ((flag & 2) && phiDeg<45.0 && phiDeg>0.0) ||
-		      ((flag & 4) && phiDeg<0.0))
-		    {
-		      std::cout<<"Qxx = "<<phiDeg<<"["<<q<<"]"
-			       <<" -> "<<qxValue<<":"<<qyValue
-			       <<" == "<<sinf(phi)<<" "<<cosf(phi)<<" xx "
-			       <<q*cosf(phi)-qyValue<<" "<<q*sinf(phi)-qxValue<<" --- "
-			       <<sqrtf(qxValue*qxValue+qyValue*qyValue)-q
-			       <<std::endl;
-		      if (phiDeg>90.0)
-			flag&=6;
-		      else if (phiDeg<0.0)
-			flag&=3;
-		      else 
-			flag&=5;
-		    }
-		  */
-		  auto  qxIter= std::lower_bound
-		    (qx_bins,qx_bins+nX+1,qxValue);
-		  auto  qyIter= std::lower_bound
-		    (qy_bins,qy_bins+nY+1,qyValue);
-		  
-		  int qXindex = std::distance(qx_bins, qxIter) - 1;
-		  int qYindex = std::distance(qy_bins, qyIter) - 1;
-		  if (qXindex < 0 || qXindex >= nX ||
-		      qYindex < 0 || qYindex >= nY)
-		    continue;
-		  
-		  int binIndex=qXindex+nX*qYindex;
-		  auto& row = segments[rank].rows[binIndex];
-
-		  // sum duplicate entries
-		  if (!row.empty() && row.back().col==pixel_index)
-		    {
-		      row.back().value+=
-			1.0f / (n_splitting * n_splitting);
-                    }
-		  else
-		    {
-		      row.emplace_back(pixel_index,
-				       1.0f/(n_splitting*n_splitting));
-		    }
-		}
-	    }
-	}
-    }
-  return;
-}
-
 void generate_matrix(const Poni& poni,
                      const py::array_t<float>& pixel_corners,
                      int n_splitting, 
@@ -212,12 +97,10 @@ void generate_matrix(const Poni& poni,
     // B C
     auto pc = pixel_corners.unchecked<4>();
     auto shape = pixel_corners.shape();
-
-#pragma omp parallel for schedule(static)
     
+    #pragma omp parallel for schedule(static)
     for (int i=0; i<shape[0]; i++) {
-      int rank = omp_get_thread_num();
-
+        int rank = omp_get_thread_num();
         for (int j=0; j<shape[1]; j++) {
             int pixel_index = i*shape[1] + j;
             if (mask[pixel_index]) {
@@ -250,7 +133,7 @@ void generate_matrix(const Poni& poni,
                     
                     float r = sqrtf(pos[0]*pos[0] + pos[1]*pos[1]);
                     float tth = atan2f(r, pos[2]);
-		    
+                    
                     float radial_coord = 0.0f; 
                     switch(output_unit) {
                         case Unit::q:
@@ -302,114 +185,68 @@ void generate_matrix(const Poni& poni,
     }
 }
 
-Sparse::Sparse
-(
- py::object py_poni,
- py::array_t<float> pixel_corners,
- int n_splitting, 
- py::array_t<int8_t> mask,
- py::array_t<float, py::array::c_style | py::array::forcecast> qxBins,
- py::array_t<float, py::array::c_style | py::array::forcecast> qyBins
- )
-{    
-
-  Poni poni;
-  poni.dist = py_poni.attr("dist").cast<float>();
-  poni.poni1 = py_poni.attr("poni1").cast<float>();
-  poni.poni2 = py_poni.attr("poni2").cast<float>();
-  poni.rot1 = py_poni.attr("rot1").cast<float>();
-  poni.rot2 = py_poni.attr("rot2").cast<float>();
-  poni.rot3 = py_poni.attr("rot3").cast<float>();
-  poni.wavelength = py_poni.attr("wavelength").cast<float>();
-
-  int max_threads = omp_get_max_threads();
-
-
-  std::vector<RListMatrix> segments;
-  int nX = qxBins.size() - 1;
-  int nY = qyBins.size() - 1;
-  int nRows= nX*nY;
-  std::cout<<"size == "<<max_threads<<" "<<nRows<<std::endl;
-  segments.resize(max_threads, nRows);
-  std::cout<<"MAx == "<<max_threads<<std::endl;
-  generateQXQYMatrix(poni,
-		     pixel_corners, 
-		     n_splitting, 
-		     segments, 
-		     mask.data(),
-		     nX,qxBins.data(),
-		     nY,qyBins.data());
-    
-  tocsr(segments,nRows,col_idx,row_ptr,values);
-}
-
 Sparse::Sparse(py::object py_poni,
                py::array_t<float> pixel_corners,
                int n_splitting, 
                py::array_t<int8_t> mask,
                const std::string& unit,
-               py::array_t<float, py::array::c_style | py::array::forcecast>
-	       radial_bins,
-               std::optional<py::array_t<float, py::array::c_style |
-	       py::array::forcecast> > phi_bins)
+               py::array_t<float, py::array::c_style | py::array::forcecast> radial_bins,
+               std::optional<py::array_t<float, py::array::c_style | py::array::forcecast> > phi_bins)
 {
-  Poni poni;
-  poni.dist = py_poni.attr("dist").cast<float>();
-  poni.poni1 = py_poni.attr("poni1").cast<float>();
-  poni.poni2 = py_poni.attr("poni2").cast<float>();
-  poni.rot1 = py_poni.attr("rot1").cast<float>();
-  poni.rot2 = py_poni.attr("rot2").cast<float>();
-  poni.rot3 = py_poni.attr("rot3").cast<float>();
-  poni.wavelength = py_poni.attr("wavelength").cast<float>();
-  
-  Unit output_unit;
-  if (unit == "q") {
-    output_unit = Unit::q;
-  }
-  else {
-    output_unit = Unit::tth;
-  }
-
-  int max_threads =  omp_get_max_threads();
-
-  std::vector<RListMatrix> segments;
-  int nradial_bins = radial_bins.size() - 1;
-  
-  int nrows, nphi_bins;
-  float* phi_data;
-  // 1D integration
-  if (!phi_bins.has_value()) {
-    nrows = nradial_bins;
-    nphi_bins = 0;
-    phi_data = nullptr;
-  }
-  // 2D integration
-  else {
-    nphi_bins = phi_bins.value().size() - 1;
-    nrows = nphi_bins * nradial_bins;
-    phi_data = phi_bins.value().mutable_data();
-  }
-
-  segments.resize(max_threads, nrows);
-  std::cout<<"Tread SIZE / Rows == "<<max_threads<<" "<<nrows<<std::endl;  
-  generate_matrix(poni, 
-		  pixel_corners, 
-		  n_splitting, 
-		  segments, 
-		  mask.data(),
-		  output_unit,
-		  nradial_bins, radial_bins.data(),
-		  nphi_bins, phi_data);
-  
-  tocsr(segments, nrows, col_idx, row_ptr, values);
+    Poni poni;
+    poni.dist = py_poni.attr("dist").cast<float>();
+    poni.poni1 = py_poni.attr("poni1").cast<float>();
+    poni.poni2 = py_poni.attr("poni2").cast<float>();
+    poni.rot1 = py_poni.attr("rot1").cast<float>();
+    poni.rot2 = py_poni.attr("rot2").cast<float>();
+    poni.rot3 = py_poni.attr("rot3").cast<float>();
+    poni.wavelength = py_poni.attr("wavelength").cast<float>();
+    
+    Unit output_unit;
+    if (unit == "q") {
+        output_unit = Unit::q;
+    }
+    else {
+        output_unit = Unit::tth;
+    }
+    
+    int max_threads = omp_get_max_threads();
+    std::vector<RListMatrix> segments;
+    int nradial_bins = radial_bins.size() - 1;
+    
+    int nrows, nphi_bins;
+    float* phi_data;
+    // 1D integration
+    if (!phi_bins.has_value()) {
+        nrows = nradial_bins;
+        nphi_bins = 0;
+        phi_data = nullptr;
+    }
+    // 2D integration
+    else {
+        nphi_bins = phi_bins.value().size() - 1;
+        nrows = nphi_bins * nradial_bins;
+        phi_data = phi_bins.value().mutable_data();
+    }
+    
+    segments.resize(max_threads, nrows);
+    generate_matrix(poni, 
+                    pixel_corners, 
+                    n_splitting, 
+                    segments, 
+                    mask.data(),
+                    output_unit,
+                    nradial_bins, radial_bins.data(),
+                    nphi_bins, phi_data);
+    
+    tocsr(segments, nrows, col_idx, row_ptr, values);
 }
 
 Sparse::Sparse(std::vector<int>&& c,
                std::vector<int>&& r,
                std::vector<float>&& v,
                std::vector<float>&& vc,
-               std::vector<float>&& vc2) :
-  col_idx(c), row_ptr(r), values(v), values_corrected(vc), values_corrected2(vc2)
+               std::vector<float>&& vc2) : col_idx(c), row_ptr(r), values(v), values_corrected(vc), values_corrected2(vc2)
 {
 }
 
@@ -511,7 +348,6 @@ py::array_t<float> Sparse::spmv_corrected2(py::array x)
 PYBIND11_MODULE(_azint, m) {
     py::class_<Sparse>(m, "Sparse")
         .def(py::init<py::object, py::array_t<float>, int, py::array_t<int8_t>, std::string, py::array_t<float>, std::optional<py::array_t<float> > >())
-        .def(py::init<py::object, py::array_t<float>, int, py::array_t<int8_t>, py::array_t<float>, py::array_t<float> >())
         .def("set_correction", &Sparse::set_correction)
         .def("spmv", &Sparse::spmv)
         .def("spmv_corrected", &Sparse::spmv_corrected)
diff --git a/azint.hpp b/azint.hpp
index f4de815..c2991bb 100644
--- a/azint.hpp
+++ b/azint.hpp
@@ -7,8 +7,8 @@ namespace py = pybind11;
 
 enum class Unit
 {
-  q,
-  tth
+    q,
+    tth
 };
 
 
@@ -28,72 +28,38 @@ struct RListMatrix
 
 struct Poni
 {
-  float dist;
-  float poni1;
-  float poni2;
-  float rot1;
-  float rot2;
-  float rot3;
-  float wavelength;
+    float dist;
+    float poni1;
+    float poni2;
+    float rot1;
+    float rot2;
+    float rot3;
+    float wavelength;
 };
 
 class Sparse
 {
 public:
-
-  Sparse(py::object py_poni,
-	 py::array_t<float> pixel_corners,
-	 int n_splitting, 
-	 py::array_t<int8_t> mask,
-	 const std::string& unit,
-	 py::array_t<float, py::array::c_style | py::array::forcecast> radial_bins,
-	 std::optional<py::array_t
-	 <float, py::array::c_style | py::array::forcecast> > phi_bins);
-
-  Sparse(py::object py_poni,
-	 py::array_t<float> pixel_corners,
-	 int n_splitting, 
-	 py::array_t<int8_t> mask,
-	 py::array_t<float, py::array::c_style | py::array::forcecast>,
-	 py::array_t<float, py::array::c_style | py::array::forcecast>);
-
-  Sparse(std::vector<int>&& c,
-	 std::vector<int>&& r,
-	 std::vector<float>&& v,
-	 std::vector<float>&& vc,
-	 std::vector<float>&& vc2);
-  void set_correction(py::array_t<float> corrections);
-  py::array_t<float> spmv(py::array x);
-  py::array_t<float> spmv_corrected(py::array x);
-  py::array_t<float> spmv_corrected2(py::array x);
-  // sparse csr matrix
-  std::vector<int> col_idx;
-  std::vector<int> row_ptr;
-  std::vector<float> values;
-  std::vector<float> values_corrected;
-  std::vector<float> values_corrected2;
-};
-
-class SparseQX
-{
-public:
-
-  SparseQX(py::object py_poni,
-	 py::array_t<float> pixel_corners,
-	 int n_splitting, 
-	 py::array_t<int8_t> mask,
-	 py::array_t<float, py::array::c_style | py::array::forcecast>,
-	 py::array_t<float, py::array::c_style | py::array::forcecast>);
-
-
-  void set_correction(py::array_t<float> corrections);
-  py::array_t<float> spmv(py::array x);
-  py::array_t<float> spmv_corrected(py::array x);
-  py::array_t<float> spmv_corrected2(py::array x);
-  // sparse csr matrix
-  std::vector<int> col_idx;
-  std::vector<int> row_ptr;
-  std::vector<float> values;
-  std::vector<float> values_corrected;
-  std::vector<float> values_corrected2;
+    Sparse(py::object py_poni,
+           py::array_t<float> pixel_corners,
+           int n_splitting, 
+           py::array_t<int8_t> mask,
+           const std::string& unit,
+           py::array_t<float, py::array::c_style | py::array::forcecast> radial_bins,
+           std::optional<py::array_t<float, py::array::c_style | py::array::forcecast> > phi_bins);
+    Sparse(std::vector<int>&& c,
+           std::vector<int>&& r,
+           std::vector<float>&& v,
+           std::vector<float>&& vc,
+           std::vector<float>&& vc2);
+    void set_correction(py::array_t<float> corrections);
+    py::array_t<float> spmv(py::array x);
+    py::array_t<float> spmv_corrected(py::array x);
+    py::array_t<float> spmv_corrected2(py::array x);
+    // sparse csr matrix
+    std::vector<int> col_idx;
+    std::vector<int> row_ptr;
+    std::vector<float> values;
+    std::vector<float> values_corrected;
+    std::vector<float> values_corrected2;
 };
diff --git a/azint/azint.py b/azint/azint.py
index 64253d5..ef7c7a9 100644
--- a/azint/azint.py
+++ b/azint/azint.py
@@ -82,19 +82,18 @@ def transform(poni, p1, p2, p3):
     phi = np.arctan2(pos[0], pos[1])
     return tth, phi
 
+
 def setup_radial_bins(poni, radial_bins, unit, tth):
     if unit == 'q':
         # calculate auto range min/max radial bins
         if not isinstance(radial_bins, Iterable):
             # q = 4pi/lambda sin( 2theta / 2 ) in A-1
-            print('theta = ',tth)
             q = 4.0e-10 * np.pi / poni.wavelength * np.sin(0.5*tth)
             radial_bins = np.linspace(np.amin(q), np.amax(q), radial_bins+1)
         else:
             radial_bins = np.asarray(radial_bins)
                 
     elif unit == '2th':
-        print('YY: ',tth)
         # radial bins in 2th are in deg
         if not isinstance(radial_bins, Iterable):
             radial_bins = np.rad2deg(np.linspace(np.amin(tth), np.amax(tth), radial_bins+1))
@@ -104,35 +103,18 @@ def setup_radial_bins(poni, radial_bins, unit, tth):
     return radial_bins
     
     
-def setup_azimuth_bins(azimuth_bins,unit=None):
-
+def setup_azimuth_bins(azimuth_bins):
     if azimuth_bins is None:
         return None
-    print('AZA == ',azimuth_bins)
+    
     if not isinstance(azimuth_bins, Iterable):
         azimuth_bins = np.linspace(0, 360, azimuth_bins+1)
     else:
         azimuth_bins = np.asarray(azimuth_bins)
     return azimuth_bins
 
-def calcQXQYbins(poni,qxqy,tth,phi):
-
-    if qxqy is None or qxqy<1:
-        return None,None
-
-    print('QX/QY bins == ',qxqy)
-
-    ## This is 3(?) values:
-    q = 4.0e-10 * np.pi / poni.wavelength * np.sin(0.5*tth)
-    qmin=np.amin(q)
-    qmax=np.amax(q)
-    
-    qx=np.linspace(-qmax,qmax,qxqy+1)
-    qy=np.linspace(-qmax,qmax,qxqy+1)
-    return qx,qy
     
 def setup_corrections(poni, solid_angle, polarization_factor, p1, p2, tth, phi):
-    ## Build all the corrections assuming
     corrections = np.ones(p1.size, dtype=np.float32)
     if solid_angle:
         solid_angle = poni.dist / np.sqrt(poni.dist**2 + p1*p1 + p2*p2)
@@ -159,8 +141,7 @@ def __init__(self,
                  mask: np.ndarray = None, 
                  solid_angle: bool = True,
                  polarization_factor: Optional[float] = None,
-                 error_model: Optional[str] = None,
-                 qxqy: Optional[Union[int,Sequence]] = None):
+                 error_model: Optional[str] = None):
         """
         Args:
             poni: Name of Poni file or instance of Poni
@@ -202,10 +183,6 @@ def __init__(self,
         self.radial_axis = 0.5*(radial_bins[1:] + radial_bins[:-1])
         azimuth_bins = setup_azimuth_bins(azimuth_bins)
         self.azimuth_axis = 0.5*(azimuth_bins[1:] + azimuth_bins[:-1]) if azimuth_bins is not None else None
-
-        qxBins,qyBins = calcQXQYbins(poni,qxqy,tth,phi)
-        self.qxBins=qxBins
-        self.qyBins=qyBins
         
         shape = pixel_centers.shape[:2]
         self.input_size = np.prod(shape)
@@ -219,30 +196,15 @@ def __init__(self,
         
         self.sparse_matrix = Sparse(poni, poni.det.pixel_corners, n_splitting, 
                                     mask, unit, radial_bins, azimuth_bins)
-            
         self.norm = self.sparse_matrix.spmv(np.ones(shape[0]*shape[1], dtype=np.float32))
         corrections = setup_corrections(poni, solid_angle, polarization_factor, p1, p2, tth, phi)
         self.sparse_matrix.set_correction(corrections)
-
-        self.qxqy_shape = None
-        self.qxqy = None
-        if (qxqy):
-            self.qxqy=qxqy
-            self.qxqy_shape = [len(qxBins) - 1, len(qyBins) - 1]
-            self.sparseQXQY = Sparse(poni,
-                                     poni.det.pixel_corners,
-                                     n_splitting,
-                                     mask,
-                                     qxBins,
-                                     qyBins)
-            self.sparseQXQY.set_correction(corrections)
         
         
     def integrate(self, 
                   img: np.ndarray, 
                   mask: Optional[np.ndarray] = None,
-                  normalized: Optional[bool] = True) -> tuple[np.ndarray, np.ndarray,
-                                                              Optional[np.ndarray],Optional[np.ndarray]]:
+                  normalized: Optional[bool] = True) -> tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]:
         """
         Calculates the azimuthal integration of the input image
         
@@ -258,7 +220,7 @@ def integrate(self,
             the norm if normalized is False
         """
         img = np.ascontiguousarray(img)
-
+        
         if img.size != self.input_size:
             raise RuntimeError('Size of image is wrong!\nExpected %d\nActual size %d' %(self.input_size, img.size))
         if mask is None:
@@ -267,32 +229,20 @@ def integrate(self,
             inverted_mask = 1 - mask
             img = img*inverted_mask
             norm = self.sparse_matrix.spmv(inverted_mask.reshape(-1))
-            
+                
         signal = self.sparse_matrix.spmv_corrected(img).reshape(self.output_shape)
+        norm = norm.reshape(self.output_shape)
         
-        outQXQY=None
-        if (self.qxqy):
-            normQXQY = self.sparseQXQY.spmv(inverted_mask.reshape(-1))
-            normQXQY=normQXQY.reshape(self.qxqy_shape)
-            qxqySignal=self.sparseQXQY.spmv_corrected(img). \
-              reshape(self.qxqy_shape)
-            outQXQY=np.divide(qxqySignal, normQXQY,
-                              out=np.zeros_like(qxqySignal),
-                              where=normQXQY!=0.0)
-##            print('Sum == ',np.sum(outQXQY))
-            
-        norm = norm.reshape(self.output_shape)        
         errors = None
         if self.error_model:
             # poisson error model
             errors = np.sqrt(self.sparse_matrix.spmv_corrected2(img)).reshape(self.output_shape)
             if normalized:
                 errors = np.divide(errors, norm, out=np.zeros_like(errors), where=norm!=0.0)
-
-
+        
         if normalized:
             result = np.divide(signal, norm, out=np.zeros_like(signal), where=norm!=0.0)
-            return result, errors, outQXQY
+            return result, errors
         else:
-            return signal, errors, norm, outQXQY
+            return signal, errors, norm