minimum-phase phase via Hilbert transform returned values
I computed and compared the minimum phase HRIR and the original one.
This is my final code:
def min_phase_coversion(HRIR):
    '''

    :param HRIR: the desired HRIR impulse response to convert into minimum phase
    :return: the minimum phase version of the original HRIR
    '''

    HRIR_fft = fft(HRIR,44100)

    #computing magnitude, tested with sinusoid, it works.
    xf = fftfreq(len(HRIR_fft), 1/44100)

    fft_magnitude = np.abs(HRIR_fft/len(HRIR_fft))

    phase=np.angle(HRIR_fft)

    minimum_phase = np.imag(-hilbert(np.log(np.abs(HRIR_fft))))

    plt.figure()
    plt.plot(xf, fft_magnitude) #computing and plotting magnitude
    plt.plot(xf, phase, label='original phase')
    plt.plot(xf,minimum_phase, label='minimum phase')
    plt.legend()
    plt.plot()
    plt.grid()
    plt.show()

    min_phase_HRIR = irfft(np.abs(HRIR_fft)*np.exp(1j*minimum_phase), 44100) # |H(w)|*e^(jPhi(w))

    HRIR_original = irfft(HRIR_fft, 44100)


    plt.figure()
    plt.plot(min_phase_HRIR[0:512], label=&quot;Min_phase&quot;, linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.plot(HRIR_original[0:512], label=&quot;original phase&quot;, linewidth=0.5, marker='o', markersize=1) #ricostruisce bene
    plt.legend()
    plt.show()


    plt.figure()
    min_phase_fft = fft(min_phase_HRIR, 44100)
    min_phase_mangnitude =  np.abs(min_phase_fft/len(min_phase_fft))
    plt.plot(xf,fft_magnitude)
    plt.plot(xf,min_phase_mangnitude)
    plt.show()

    return min_phase_HRIR

The results I get seem reasonable looking at the plots:

and the magnitude is the same as I expected (the original one and the min phase one are perfectly overlapping):

However, I am not totally sure of the phase behaviour (as mentioned in my question) and its returned values.