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’:

https://pasteall.org/blend/9f468cbe96f145eb9d77553f3ebf1b89

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 = bpy.data.objects[‘Armature’].pose.bones[‘shoulder’].matrix_basis.copy()

bpy.data.objects[‘Armature’].pose.bones[‘shoulder’].rotation_euler = Vector((0,0,0))
bpy.context.evaluated_depsgraph_get()

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

bpy.data.objects[‘Armature’].pose.bones[‘shoulder’].matrix = b

difEuler = Vector(b.to_euler()[:]) - Vector(a.to_euler()[:]) + Vector(Original_rot.to_euler()[:])

bpy.data.objects[‘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!