Align to vector using Vector.rotation_difference(v)
Altered the script from this answer to find the quaternion rotation difference between the local z axis of the edit mesh and the normal of its selected (active) face.
Run script with target object in edit mode, one face selected and active. An object in scene named "Cube.001" will be located at face centre and oriented by face normal.
import bpy
from mathutils import Matrix, Vector
import bmesh
context = bpy.context
obj = context.edit_object
mw = obj.matrix_world.copy()
bm = bmesh.from_edit_mesh(obj.data)
# for this example jmake a face active
face = bm.select_history.active
o = face.calc_center_median()
# calculate the axis dif in local coords
axis_src = face.normal
# local z-axis
axis_dst = Vector((0, 0, 1))
matrix_rotate = mw.to_3x3()
matrix_rotate = matrix_rotate * axis_src.rotation_difference(axis_dst).to_matrix()
matrix_translation = Matrix.Translation(mw * o) #
obj2 = context.scene.objects.get("Cube.001")
obj2.matrix_world = matrix_translation * matrix_rotate.to_4x4()
Result of script. Moved cube in edit mode, showing local coords.