Index out of bounds on TMatrixDSym using RDataFrame in PyRoot

[bthornbe]

Hello, I am getting an error when trying to fill a TMatrixDSym in a .Define() in RDataFrames in PyRoot. The error is just the following two lines over and over:

Error in <operator=>: matrices not compatible
Error in <TMatrixTRow_const(const TMatrixTSym &,Int_t)>: row index out of bounds

However, the program finishes running after the errors print.
Here is a reproducer for the error:

import ROOT

fullname = "root:path_to_a_root_file"

tname = "my_tree's_name"

mystring= \
'''
TArrayD array(9);   
for (int i=0; i<nTracks; i++) { 
    for (int j=0; j<3; j++) {
        for (int k=0; k<3; k++) {
            array[j+3*k] += Momenta.at(i).at(j)*Momenta.at(i).at(k)/(sqrt(Tracks[i].Mag2())*Denominator);
        }
    }
} 
TMatrixDSym s(3, array.GetArray());
return s;
'''
dFrame= ROOT.ROOT.RDataFrame(tname, fullname).Define("nTracks", "Tracks.size()") \
        .Filter("nTracks > 0") \
        .Define("CutHT", "double cutht=0; for (int i=0; i<Jets.size(); i++) if (Jets[i].Pt()>30 and abs(Jets[i].eta())<2.4) cutht+=Jets[i].Pt(); return cutht") \
        .Filter("CutHT>500") \
        .Define("Momenta", "vector<vector<double>> p; for (int i=0; i<nTracks; i++) {p.emplace_back(); p[i].push_back(Tracks[i].x()); p[i].push_back(Tracks[i].y()); p[i].push_back(Tracks[i].z());} return p;") \
        .Define("Denominator", "double denom=0; for (int i=0; i<nTracks; i++) denom += sqrt(Tracks[i].Mag2()); return denom;") \
        .Define("SphericityTensor", mystring) \
        .Define("Val", "return (SphericityTensor[0][0])") \


model = ROOT.RDF.TH1DModel("Val", "Val", 50, -1000., 1000.)
hist = dFrame.Histo1D(model, "Val")

can = ROOT.TCanvas("canName", "canTitle")
file = ROOT.TFile('reproducerHists', 'RECREATE')
hist.Write()

I defined mystring separately so that it’s syntax is more clear because I think that is probably where the problem is. I can’t find anything that would cause an index to be out of bounds, though. Does anyone have any ideas?

Thanks in advance for sharing your knowledge!

---

Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

---

[eguiraud]

Hi @bthornbe ,
ah this is tricky, here’s what I think is happening.
When you do return s, RDF tries to assign s to a variable that stores the result of your computation for that entry. Initially, that variable will be a default-constructed TMatrixDSym.
Unfortunately, however, it is not possible to assign a TMatrixDSym with a certain shape to a TMatrixDSym with another shape!

So basically this happens:

~ root -l
root [0] TMatrixDSym rdf_var; // this is the internal RDF variable, default-constructed
root [1] TMatrixDSym s(3, 4); // this is your `s`
root [2] rdf_var = s // does not work!
Error in <operator=>: matrices not compatible

So you can’t just return s. Maybe @moneta can suggest what to use instead.

Cheers,
Enrico

[bthornbe]

What is the default constructor of TMatrixDSym? I am struggling to find clear documentation on it.

[eguiraud]

It’s what you get if you call TMatrixDSym m; – it constructs an empty matrix (zero rows and zero columns). The problem is that assigning another non-empty TMatrixDSym to an empty matrix is an error, so I am not sure TMatrixDSym can be used as the result of a Define.

[bthornbe]

eguiraud was right about the cause of the problem, and I’m going to post my workaround in case anyone else finds this thread while searching for a solution to a similar problem. The reason I wanted SphericityTensor in the form of a TMatrixDSym was that I wanted to use TMatrixDSymEigen to find its eigenvalues. For that reason, the reproducer below has different output in order to show how the workaround works.

Instead of returning a TMatrixDSym in the definition of SphericityTensor, I am returning a 2-d c++ vector, and then I make a TMatrixDSym out of that vector within the definition of EigenVals. TVectorD appears to have the same type of issue with the default constructor, so I also have a loop in that definition that takes the output of TMatrixDSymEigen.GetEigenValues() from a TVectorD and puts it in a c++ vector.

import ROOT

fullname = "root:path_to_a_root_file"

tname = "my_tree's_name"

tensorString = \
    '''
    vector<vector<double>> s{{0,0,0},{0,0,0},{0,0,0}};    
    for (int i=0; i<nTracks; i++) 
    { 
        for (int j=0; j<3; j++) 
        {
            for (int k=0; k<3; k++) 
            {
                s.at(j).at(k) += Momenta.at(i).at(j)*Momenta.at(i).at(k)/(sqrt(Tracks[i].Mag2())*Denominator);
            }
        }
    } 
    return s;
    '''
eigenString = \
'''
double a[9];
for (int i = 0; i < 3; i++)
{
    for (int j = 0; j < 3; j++)
    {
        a[i+3*j] = SphericityTensor.at(i).at(j);
    }
}
TMatrixDSym s(3, a);
TMatrixDSymEigen eigen(s); 
TVectorD eigenVals = eigen.GetEigenValues();
vector <double> eigenVals2;
for (int i = 0; i < 3; i++)
{
        eigenVals2.push_back(eigenVals[i]);
}
return eigenVals2;
'''

dFrame= ROOT.ROOT.RDataFrame(tname, fullname).Define("nTracks", "Tracks.size()") \
        .Filter("nTracks > 0") \
        .Define("CutHT", "double cutht=0; for (int i=0; i<Jets.size(); i++) if (Jets[i].Pt()>30 and abs(Jets[i].eta())<2.4) cutht+=Jets[i].Pt(); return cutht") \
        .Filter("CutHT>500") \
        .Define("Momenta", "vector<vector<double>> p; for (int i=0; i<nTracks; i++) {p.emplace_back(); p[i].push_back(Tracks[i].x()); p[i].push_back(Tracks[i].y()); p[i].push_back(Tracks[i].z());} return p;") \
        .Define("Denominator", "double denom=0; for (int i=0; i<nTracks; i++) denom += sqrt(Tracks[i].Mag2()); return denom;") \
        .Define("SphericityTensor", tensorString) \
        .Define("EigenVals", eigenString) \
        .Define("NewVal", "return (EigenVals[0])")


model = ROOT.RDF.TH1DModel("NewVal", "NewVal", 50, -1000., 1000.)
hist = dFrame.Histo1D(model, "NewVal")

can = ROOT.TCanvas("canName", "canTitle")
file = ROOT.TFile('reproducerHists', 'RECREATE')
hist.Write()

[moneta]

Hi,
Given the fact you can early build a TMatrix or a TVector from a simple C array or std::vector, I would use them to return your data values, and then build the TMatric,TVector from those

Cheers

Lorenzo