Bone Matrix Problem

I have this simplified blend file were I need to match the bone ‘shoulder’ to this other bone called ‘shoulder_extra_snap’:

So, both bones have constraints that rotate them. In order to make the rotations match, I’m calculating a rotational difference with matrix.to_euler. and then I’m applying that rotation to ‘shoulder’.

This is the script you will find inside of the blend file:

import bpy
from mathutils import Vector

Original_rot =[‘Armature’].pose.bones[‘shoulder’].matrix_basis.copy()[‘Armature’].pose.bones[‘shoulder’].rotation_euler = Vector((0,0,0))

a =[‘Armature’].pose.bones[‘shoulder’].matrix.copy()
b =[‘Armature’].pose.bones[‘shoulder_extra_snap’].matrix.copy()[‘Armature’].pose.bones[‘shoulder’].matrix = b

difEuler = Vector(b.to_euler()[:]) - Vector(a.to_euler()[:]) + Vector(Original_rot.to_euler()[:])[‘Armature’].pose.bones[‘shoulder’].rotation_euler = difEuler

This calculation works as long as I don’t rotate the ‘master’ bone, which is the parent of the other bones. Therefore I need to fix that.

I’m pretty sure the solution must be to multiply some matrix by some other matrix so that I can get the local rotation difference but also take the transforms of the hierarchy into account, but I don’t know anything about it! (Although I did read some answers regarding similar topics around here, but still I don’t get it and it’s not my area of expertise, haha)

Any help will be appreciated! thanks!