wrong implementation of gdml matrices?

[ManuelHu]

I think I found a quite serious mis-implementation of GDML matrices. Luckily, this does not affect all use cases (that is why we never have seen this in our optical use case)

see this short notebook example, in which I would expect a different result:

---

from pyg4ometry import gdml, geant4
import numpy as np

reg = geant4.Registry()

mv = gdml.MatrixFromVectors([1, 2, 3], [100, 200, 300], "testmatrix", reg)

mv.values

[Expression : matrix_expr_idx0_val_testmatrix = expr_matrix_expr_idx0_val_testmatrix : 1*eV,
 Expression : matrix_expr_idx1_val_testmatrix = expr_matrix_expr_idx1_val_testmatrix : 100,
 Expression : matrix_expr_idx2_val_testmatrix = expr_matrix_expr_idx2_val_testmatrix : 2*eV,
 Expression : matrix_expr_idx3_val_testmatrix = expr_matrix_expr_idx3_val_testmatrix : 200,
 Expression : matrix_expr_idx4_val_testmatrix = expr_matrix_expr_idx4_val_testmatrix : 3*eV,
 Expression : matrix_expr_idx5_val_testmatrix = expr_matrix_expr_idx5_val_testmatrix : 300]

--> this has the right order

mv.eval()

array([[9.99999566e-07, 1.00000000e+02, 1.99999913e-06],
       [2.00000000e+02, 2.99999870e-06, 3.00000000e+02]])

--> this does not have the right structure at all? It is just the values array arbitrarily folded in half?

mv.eval().reshape(-1).reshape(3, 2)

array([[9.99999566e-07, 1.00000000e+02],
       [1.99999913e-06, 2.00000000e+02],
       [2.99999870e-06, 3.00000000e+02]])

--> this would be the right result, with the two columns in GDML organized in columns of the evaled matrix?

this is not just the transposition of the "wrong" matrix above, but a complete re-ordering of all values. (if eval() returns a matrix like this, or the transposed value, could still be discussed; especially for 1D "matrices")

[stewartboogert]

@ManuelHu looks bad! For the accelerator work we do I feel that this is unlikely well tested. What about creating the matrix from gdml.Matrix constructor opposed to the helper function gdml.MatrixFromVectors.

[ManuelHu]

I think that gdml.MatrixFromVectors is producing the right matrices. I had always verified that the material property tables had been read correctly.

The problem only appears when trying to read the matrix from python. reading and writing from Geant4 is luckily fine, I suppose. I would just guess that nobody did try pyg4ometry for this task until now...

---

I will still verify that this is really what Geant4 expects, but I am very sure that Geant4 that this interleaved format of the column vectors:

mv2 = gdml.Matrix("testmatrix2", 2, ["1*eV", 100, "2*eV", 200, "3*eV", 300], reg)

which gives exactly the same result:

np.all(mv.eval() == mv2.eval())

np.True_

---

https://github.com/Geant4/geant4/blob/df176550b3e44cdd70406ff2b581896426b5bcbe/source/persistency/gdml/src/G4GDMLEvaluator.cc#L122-L131

see also this reading code, which exactly reads such an interleaved value string