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)
face = bm.select_history.active
o = face.calc_center_median()

axis_src = face.normal
axis_src2 = face.calc_tangent_edge()
axis_dst = Vector((0, 0, 1))
axis_dst2 = Vector((0, 1, 0))

vec2 = axis_src * obj.matrix_world.inverted()
matrix_rotate = axis_dst.rotation_difference(vec2).to_matrix().to_4x4()

vec1 = axis_src2 * obj.matrix_world.inverted()
axis_dst2 = axis_dst2*matrix_rotate.inverted()
mat_tmp = axis_dst2.rotation_difference(vec1).to_matrix().to_4x4()
matrix_rotate = mat_tmp*matrix_rotate
matrix_translation = Matrix.Translation(mw * o) #

obj2 = context.scene.objects.get("Cube.001")
obj2.matrix_world = matrix_translation * matrix_rotate.to_4x4()