I have found assigning to obj.rotation_quaternion does nothing, whereas assigning to obj.rotation_euler does, but I'm not sure if this will help.

Firstly, I'm not clear you are getting the normal of the selected face as opposed to the first face.

Quaternions are tricky. I have found assigning to obj.rotation_quaternion does nothing, whereas assigning to obj.rotation_euler does, but I'm not sure I don't think I set rotation mode as you do at the end.

I don't think there is any value in setting the rotation of the first object to the selected face, this apart from not working (perhaps because of the above) doesn't help align the cone.

And I've got in the habit of using a bmesh, but you might be able to do it directly from the data. But once you have the rotation Quaternion, you have to multiply the object's existing rotation by the inverse of the Quaternion. Well, it worked for me.

# This example assumes we have a mesh object in edit-mode

import bpy
import bmesh
import mathutils

# Get the active mesh (assumes starting in edit mode)
obj = bpy.context.edit_object
me = obj.data
coneradius=1     # Im not sure how you plan to scale the cone.
Up=mathutils.Vector((0,0,1))

# Get a BMesh representation
bm = bmesh.from_edit_mesh(me)
selected_faces=[face for face in bm.faces if face.select]
if len(selected_faces)>0:
    norm=selected_faces[0].normal
    my_location=selected_faces[0].calc_center_median()
    # convert this to a rotation quaternion    
    Qrot=norm.rotation_difference(Up)

    #Add the cone in object mode
    bpy.ops.object.mode_set(mode='OBJECT')
    bpy.ops.mesh.primitive_cone_add(location=my_location)
    cone=bpy.context.object
    cone.rotation_mode='QUATERNION'
    cone.rotation_quaternion*=Qrot.inverted()
    
    # Im not sure how you plan to scale the cone and translate it.