So you compute the matrix that best converts white balanced raw data to reference xyz/Lab values resulting from the target being lit by illuminant E?

Yes. I specifically converted the Babel color spectral data to XYZ using NO illuminant spectrum. That results in the cleanest XYZ-E reference for my purposes. I can make those XYZ colors available if you'd need them.

Then perform a chromatic adaptation to D65 for a/sRGB?

Yes, although I convert to ICC D50 using straight multiplication for adaptation. Then let colorsync figure out conversion to whatever the user selects as outputspace. I may change this later to incorporate colorclipping control. The trick then is to adapt the outputmatrix to illuminant E by simply dividing out its original illuminant.

I prefer to use simple division/multiplication for illuminant conversions because during several steps of the process it is already done that way. i.e. whitebalancing the raw data is done by multiplication, Lab to XYZ is generally defined with multiplication. Using bradford transforms

has a disputable theoretical advantage if you allow 2 degree observer logic to assess a 10 degree observer problem for extreme conversions. That advantage is however not of the precision that would radically change the desired results, especially in smaller conversions.

Since the matrix varies with (the inverse of) color temperature, how accurate does it come out compared to trying to estimate it for the correct illuminant as I did? Maybe I should try it out.

Some thoughts:

the algorithm trying to find the matrix uses reference XYZ values. If those values are given with D50 or D65 illuminant, then the algorithm should end up finding D50 or D65 illuminant, since you already whitebalanced the data.

If you have the spectral data of the patches, and the spectral data of the original scene illuminant, then you could theoretically adapt the patches using spectral conversion, though I doubt that the entire process recognises that kind of precision, and the objective of colormatching is to generalise the matrixes so they can easily be adapted under normal PCS circumstances.

That is: we're trying to find matrixes that are optimised for use in XYZ and straight illuminant conversions.

Note also that CCT is not an equivalence. That is an unfortunate mistake perhaps introduced by RAW converters. CCT is a one-way conversion, and very sensitive, it may be insightful to try and figure out the size of the differences in XYZ vs the corresponding size of the differences in CCT (and then try to translate that to a perceptually uniform error space). But that would be (euphemistically) left as an exercise for the reader...