diff --git a/azint.cpp b/azint.cpp
index 96fae8a..7d9b339 100644
--- a/azint.cpp
+++ b/azint.cpp
@@ -80,6 +80,121 @@ 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, 
@@ -97,10 +212,12 @@ 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]) {
@@ -133,7 +250,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:
@@ -185,68 +302,114 @@ 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);
-    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);
+  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);
 }
 
 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)
 {
 }
 
@@ -348,6 +511,7 @@ 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 c2991bb..f4de815 100644
--- a/azint.hpp
+++ b/azint.hpp
@@ -7,8 +7,8 @@ namespace py = pybind11;
 
 enum class Unit
 {
-    q,
-    tth
+  q,
+  tth
 };
 
 
@@ -28,38 +28,72 @@ 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(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;
+
+  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;
 };
diff --git a/azint/azint.py b/azint/azint.py
index ef7c7a9..64253d5 100644
--- a/azint/azint.py
+++ b/azint/azint.py
@@ -82,18 +82,19 @@ 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))
@@ -103,18 +104,35 @@ def setup_radial_bins(poni, radial_bins, unit, tth):
     return radial_bins
     
     
-def setup_azimuth_bins(azimuth_bins):
+def setup_azimuth_bins(azimuth_bins,unit=None):
+
     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)
@@ -141,7 +159,8 @@ def __init__(self,
                  mask: np.ndarray = None, 
                  solid_angle: bool = True,
                  polarization_factor: Optional[float] = None,
-                 error_model: Optional[str] = None):
+                 error_model: Optional[str] = None,
+                 qxqy: Optional[Union[int,Sequence]] = None):
         """
         Args:
             poni: Name of Poni file or instance of Poni
@@ -183,6 +202,10 @@ 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)
@@ -196,15 +219,30 @@ 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]]:
+                  normalized: Optional[bool] = True) -> tuple[np.ndarray, np.ndarray,
+                                                              Optional[np.ndarray],Optional[np.ndarray]]:
         """
         Calculates the azimuthal integration of the input image
         
@@ -220,7 +258,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:
@@ -229,20 +267,32 @@ 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
+            return result, errors, outQXQY
         else:
-            return signal, errors, norm
+            return signal, errors, norm, outQXQY