Set matrix rotation value (not add to it)
Similarly to the answer to Convert mouse movement into rotation value (not in bge), customized to make the all selected objects rotate about the active object.

Need to subtract the pivot point of the rotation, which gives the "pivot arm" vector for the rotation.  In the previous answer, where just the active object was being rotated it had been hard coded to (0, 0, 0), since that is the result of subtracting an objects global location from itself.

import bpy, math
from bpy_extras.view3d_utils import (
        location_3d_to_region_2d,
        region_2d_to_vector_3d)

from mathutils import Matrix, Vector

class CustomRotate(bpy.types.Operator):
    bl_idname = "transform.custom_rotate"
    bl_label = "Custom Rotate"
    pivot_point = Vector((0, 0))
    v2d = Vector((0, 0))
    axis = Vector((1, 0, 0))
    mw = Matrix()

    def modal(self, context, event):
        ob = context.object
        if event.type in ['LEFTMOUSE','ENTER']:
            print('finished')
            return{'FINISHED'}

        if event.type in ['ESC','RIGHTMOUSE']: 
            ob.matrix_world = self.mw
            return {'CANCELLED'}

        if event.value == 'PRESS' and event.type == self.axis:
            # two times make it local
            axis = Vector([int(c == event.type) for c in 'XYZ'])
            self.axis = ob.matrix_world * axis  - ob.matrix_world * Vector()
            return {'RUNNING_MODAL'}                
        elif event.value == 'PRESS' and event.type in "XYZ":
            # to rotate about global axes
            self.axis = event.type
            # local
            return {'RUNNING_MODAL'}

        v2d = Vector((event.mouse_region_x, event.mouse_region_y)) - self.pivot_point
        if not v2d.length:
            return {'PASS_THROUGH'}
        #Get rotation angle and update object rotation
        a = self.v2d.angle_signed(v2d)
        rot_mat = Matrix.Rotation(a, 4, self.axis)
        t = ob.matrix_world.translation.copy()
        #Print angle in degrees to help with testing...
        #print('alpha angle:',math.degrees(a))
        for o in context.selected_objects:            
            mw = o.matrix_world
            mw.translation -=  t            
            mw = rot_mat * mw 
            mw.translation  += t
            o.matrix_world = mw

        self.v2d = v2d        

        return {'RUNNING_MODAL'}

    def invoke(self, context, event):        
        region = context.region
        rv3d = context.region_data

        ob = context.object
        self.pivot_point = location_3d_to_region_2d(region, rv3d, ob.matrix_world.translation)
        self.v2d = Vector((event.mouse_region_x, event.mouse_region_y)) - self.pivot_point
        z_axis = Vector((0, 0, -1))

        self.axis = rv3d.view_matrix.transposed().to_3x3() * z_axis
        self.mw = ob.matrix_world.copy()
        context.window_manager.modal_handler_add(self)

        return {'RUNNING_MODAL'}

def register():
    bpy.utils.register_class(CustomRotate)             
def unregister():
    bpy.utils.unregister_class(CustomRotate) 

if __name__ == "__main__":
    register()