diff --git a/cscripts/csphagen.py b/cscripts/csphagen.py
index 282fd5a..2b7ee40 100755
--- a/cscripts/csphagen.py
+++ b/cscripts/csphagen.py
@@ -363,6 +363,116 @@ class csphagen(ctools.csobservation):
         # Return
         return
 
+    def _compute_region_separation(self, pnt_dir):
+        """ Compute the separation angle for reflected off regions in radians.
+
+        Returns
+        -------
+        angle : `float`
+            Separation angle of two off regions (rad)
+        """
+        # Initialise the result
+        separation = -1
+
+        # Compute offset of reflected regions to pointing position
+        offset = pnt_dir.dist_deg(self._src_dir)
+
+        if self._srcshape == 'CIRCLE':
+
+            # Compute apparent diameter of the circle
+            separation = 2.0 * self._rad / offset
+
+        elif self._srcshape == 'RECT':
+
+            # Compute the opening angle towards combinations of rectangle corners.
+            # This method overestimates the real need of space between the
+            # rectangles, so the method may be optimised to gain more off regions!
+            # Anyway, it is assured that the off regions will never overlap.
+
+            # Get the corners
+            cs = [self._src_reg[0].get_corner(icorner) for icorner in range(4)]
+
+            # Auxillary, local method for opening angle computation
+            def __compute_posang(pnt_dir,a,b):
+                """ Compute the difference in position angle to the pointing in
+                    degrees.
+                """
+                posang_a = pnt_dir.posang_deg(a) %360
+                posang_b = pnt_dir.posang_deg(b) %360
+                posang = abs(posang_a - posang_b)
+                return posang
+
+            # Compute the 6 opening angles
+            combinations = [[0,1], [0,2], [0,3], [1,2], [1,3], [2,3]]
+            angles = [__compute_posang(pnt_dir, cs[i], cs[j]) for i,j in combinations]
+
+            # The desired separation is the maximum opening angle
+            separation = max(angles) *gammalib.deg2rad
+
+
+            ## Different approach:
+            ## Find the separation angle for which two rectangles to not overlap
+            ## Issue: Scanning not optimal.
+            # posang    = pnt_dir.posang_deg(self._src_dir)
+            # dphi_incr = 0.1
+            # dphi      = dphi_incr
+            # print("Src ctr direction: ra={} dec={}".format(self._src_dir.ra_deg(), self._src_dir.dec_deg()))
+            # print(self._src_reg)
+            # while dphi<360:
+
+            #     ctr_dir = pnt_dir.clone()
+            #     ctr_dir.rotate_deg(posang + dphi, offset)
+            #     print("Query ctr direction: ra={} dec={}".format(ctr_dir.ra_deg(), ctr_dir.dec_deg()))
+
+            #     # Create test rectangle
+            #     reg = gammalib.GSkyRegionRect(ctr_dir
+            #                                  ,self._reg_width
+            #                                  ,self._reg_height
+            #                                  ,self._reg_posang + dphi
+            #                                  )
+            #     print(reg)
+
+            #     # Check test rectangle corner containment in reference rectangle
+            #     is_corner_contained = self._src_reg.overlaps(reg)
+            #     # is_corner_contained = False
+            #     # for icorner in range(4):
+            #     #     corner_dir = reg.get_corner(icorner)
+            #     #     is_corner_contained |= self._src_reg.contains(corner_dir)
+
+            #     # Check if any corner was contained
+            #     if not is_corner_contained:
+            #         diameter = dphi
+            #         break
+
+            #     # Increment dphi
+            #     dphi += dphi_incr
+
+
+            ## Very different approach: Compute the width along the circular
+            ## region placement path.
+            ## Issue: Long regions will overlap in the inner part.
+            # def l(w,h,phi):
+            #     """ Super complicated method to compute the height along dec-axis
+            #         of a rotated rectangle.
+            #     """
+            #     phi = np.deg2rad(phi) % (2*np.pi)
+            #     phi_c = np.arctan2(w,h)
+            #     if phi>=(0.5*np.pi) and (phi<np.pi):
+            #         phi = np.pi-phi
+            #     elif (phi>=np.pi) and (phi<(1.5*np.pi)):
+            #         phi -=np.pi
+            #     elif (phi>=(1.5*np.pi)) and (phi<(2*np.pi)):
+            #         phi = 2*np.pi - phi
+
+            #     if (0<=phi) and (phi<phi_c):
+            #         result = h/2/np.cos(phi)
+            #     else:
+            #         result = w/2/np.sin(phi)
+            #     return 2*result
+
+        # Return
+        return separation
+
     def _reflected_regions(self, obs):
         """
         Calculate list of reflected regions for a single observation (pointing)
@@ -402,7 +512,7 @@ class csphagen(ctools.csobservation):
         # ... otherwise
         else:
             posang = pnt_dir.posang_deg(self._src_dir)
-            if self._srcshape == 'CIRCLE':
+            if (self._srcshape == 'CIRCLE') or (self._srcshape == 'RECT'):
 
                 # Determine number of background regions to skip
                 N_skip  = self['bkgregskip'].integer()
@@ -412,7 +522,7 @@ class csphagen(ctools.csobservation):
                 # Compute the angular separation of reflected regions wrt
                 # camera center. The factor 1.05 ensures background regions
                 # do not overlap due to numerical precision issues
-                alpha = 1.05 * 2.0 * self._rad / offset
+                alpha = 1.05 * self._compute_region_separation(pnt_dir)
 
                 # Compute number of reflected regions by dividing the angular
                 # separation by 2 pi.
@@ -437,7 +547,15 @@ class csphagen(ctools.csobservation):
                     while dphi <= dphi_max:
                         ctr_dir = pnt_dir.clone()
                         ctr_dir.rotate_deg(posang + dphi, offset)
-                        region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
+                        if self._srcshape == 'CIRCLE':
+                            region = gammalib.GSkyRegionCircle(ctr_dir, self._rad)
+                        elif self._srcshape == 'RECT':
+                            # Adjust the posang of the rectangle correspondingly
+                            region = gammalib.GSkyRegionRect(ctr_dir
+                                                            ,self._reg_width
+                                                            ,self._reg_height
+                                                            ,self._reg_posang + dphi
+                                                            )
                         if self._has_exclusion:
                             if self._excl_reg.overlaps(region):