TGraph crashes for big data

ROOT
,
1

Dear all,

I have some simulations in C++ where I heavily make use of ROOT libraries and produce a ROOT file for data analysis. My ROOT file has around 800,000 points, and I’m using RDataFrame to read everything, plot the spacial coordinates with TGraph2d, and other computed variables with TGraph. Funny enough, TGraph2D does its job without any problems for the 800,000 points, but TGraph always crashes. I have tried cutting my RVec’s with (much) less data, and then it starts plotting stuff, but it struggles.

Is that a known issue/limitation with TGraph? Are there known methods to fix it?

Many thanks in advance.

ROOT Version: 6.30
Platform: GNU/Linux Debian 12
Compiler: Cling

2

Hi Leonardo,

Sorry to hear you are experiencing this issue.
Just to understand: can we reproduce the issue with a RDF generating 800k points and plotting them in a TGraph?
Have you checked you checked you are not hitting a limitation of the system, e.g. memory?

Best,
D

1 Like
3

This works fine for me:

void graph800000() {
   auto g = new TGraph();
   for (int i=0; i<800000; i++) {
      g->AddPoint(i,i);
   }
   g->Draw();
}
1 Like
4

Dear Danilo,

Yes of course. I’m sorry that I haven’t immediately shared a reproducible example. Here’s my code to read the file:

#include <iostream>

#include "TKey.h"
#include "TMath.h"
#include "TFile.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TString.h"
#include "TCanvas.h"
#include "TGraph2D.h"
#include "ROOT/RVec.hxx"
#include "ROOT/RDataFrame.hxx"

// Constants
const Double_t e = 1.60217663e-19;  // Electron charge [C]
const Double_t c = 299792458.0;     // Speed of light [m.s-1]
const Double_t me = 9.1093837e-31;  // Electron mass [kg]
const Double_t mp = 1.67262192e-27; // Proton mass [kg] 
const Double_t Me = 5.97216787e+24; // Earth's mass [kg]
const Double_t Mj = 1.89812460e+27; // Jupiter's mass [kg]
const Double_t Re = 6378100.0;      // Earth Radius [m]
const Double_t Rj = 71492000.0;     // jupiter radius [m]

class Planet{

private:
    Int_t id;
    TString name;
    TString color;
    Double_t mass;
    Double_t radius;

public:
    Int_t getID() {return id;}
    TString getName() {return name;};
    TString getColor() {return color;}
    Double_t getMass() {return mass;}
    Double_t getRadius() {return radius;}

    void setID(Int_t ID) {id = ID;}
    void setName(TString Name) {name = Name;}
    void setColor(TString Color) {color = Color;}
    void setMass(Double_t Mass) {mass = Mass;}
    void setRadius(Double_t Radius) {radius = Radius;}
};

class Particle{

private:
    Int_t pid;
    TString name;
    Double_t mass;
    Double_t charge;
    bool relativistic;

public:
    Int_t getPID() {return pid;}
    TString getName() {return name;}
    Double_t getCharge() {return charge;}
    Double_t getMass() {return mass;}
    bool isRelativistic() {return relativistic;}

    void setPID(Int_t p) {pid = p;}
    void setName(TString n) {name = n;}
    void setCharge(Double_t q) {charge = q;}
    void setMass(Double_t m) {mass = m;}
    void setRelativistic(bool rel) {relativistic = rel;}
};

Planet getPlanet(Int_t id);
Particle getParticle(Int_t pid, Int_t rel);

ROOT::RVec<TString> getTrees(TString file);

Double_t norm(ROOT::RVec<Double_t> vector);

ROOT::RVec<Double_t> magneticMoment(Particle particle,
                                    ROOT::RVec<ROOT::RVec<Double_t>> velocity,
                                    ROOT::RVec<ROOT::RVec<Double_t>> field);

void dataAnalysis(){

    TString fileName = "lorentz_Sph10.root";

    ROOT::RVec<TString> trees = getTrees(fileName);
    ROOT::RDataFrame rdf(trees[0], "lorentz_Sph10.root");

    auto planetID = rdf.Take<ROOT::RVec<Int_t>>("planet").GetValue()[0][0];
    auto particleID = rdf.Take<ROOT::RVec<Int_t>>("particle").GetValue()[0][0];
    auto relativistic = rdf.Take<ROOT::RVec<Int_t>>("relativistic").GetValue()[0][0];

    Planet planet = getPlanet(planetID);
    Particle particle = getParticle(particleID, relativistic);

    auto t = rdf.Take<ROOT::RVec<Double_t>>("t").GetValue()[0];
    
    auto x = rdf.Take<ROOT::RVec<Double_t>>("x").GetValue()[0];
    auto y = rdf.Take<ROOT::RVec<Double_t>>("y").GetValue()[0];
    auto z = rdf.Take<ROOT::RVec<Double_t>>("z").GetValue()[0];

    auto vx = rdf.Take<ROOT::RVec<Double_t>>("vx").GetValue()[0];
    auto vy = rdf.Take<ROOT::RVec<Double_t>>("vy").GetValue()[0];
    auto vz = rdf.Take<ROOT::RVec<Double_t>>("vz").GetValue()[0];

    auto Bx = rdf.Take<ROOT::RVec<Double_t>>("Bx").GetValue()[0];
    auto By = rdf.Take<ROOT::RVec<Double_t>>("By").GetValue()[0];
    auto Bz = rdf.Take<ROOT::RVec<Double_t>>("Bz").GetValue()[0];

    ROOT::RVec<ROOT::RVec<Double_t>> velocity{vx, vy, vz};
    ROOT::RVec<ROOT::RVec<Double_t>> bField{Bx, By, Bz};

    ROOT::RVec<Double_t> moment = magneticMoment(particle, velocity, bField);

    gStyle->SetCanvasPreferGL(kTRUE);

    auto c0 = new TCanvas("c0", "Particle Trajectories", 1200, 800);
    auto g0 = new TGraph2D(t.size());
    for(Int_t i=0; i<t.size(); i++){
        g0->SetPoint(i,x[i]/planet.getRadius(),
                       y[i]/planet.getRadius(),
                       z[i]/planet.getRadius());
    }
    g0->SetLineColor(kBlue);
    g0->SetLineWidth(2);
    g0->SetTitle("Adiabatic Motion of Particles - " + planet.getName());
    if(planetID == 3){
        g0->GetXaxis()->SetTitle("X [Re]");
        g0->GetYaxis()->SetTitle("Y [Re]");
        g0->GetZaxis()->SetTitle("Z [Re]");
    }
    else if(planetID == 5){
        g0->GetXaxis()->SetTitle("X [Rj]");
        g0->GetYaxis()->SetTitle("Y [Rj]");
        g0->GetZaxis()->SetTitle("Z [Rj]");
    }    
    g0->Draw("Line");
    c0->Modified();
    c0->Update();

    auto c1 = new TCanvas("c1", "Magnetic Moment", 1200, 800);
    c1->SetGrid(1,1);
    auto g1 = new TGraph();
    for(Int_t i=0; i<t.size(); i++) g1->SetPoint(i,t[i],moment[i]);
    g1->SetLineColor(kBlue);
    g1->SetLineWidth(2);
    g1->GetXaxis()->SetTitle("Time [seconds]");
    g1->GetYaxis()->SetTitle("Magnetic Moment [seconds]");
    g1->SetTitle("Magnetic Moment [J.T-1]");
    g1->Draw("AC");
    c1->Modified();
    c1->Update();
}

Planet getPlanet(Int_t id){
    Planet planet;
    if(id == 3){
        planet.setID(id);
        planet.setName("Earth");
        planet.setMass(Me);
        planet.setRadius(Re);
    }
    else if(id == 5){
        planet.setID(id);
        planet.setName("Jupiter");
        planet.setMass(Mj);
        planet.setRadius(Rj);
    }
    return planet;
}

Particle getParticle(Int_t pid, Int_t rel){
    Particle particle;
    if(pid == +11){
        particle.setPID(+11);
        particle.setName("electron");
        particle.setMass(me);
        particle.setCharge(-e);
    }
    else if(pid == -11){
        particle.setPID(-11);
        particle.setName("positron");
        particle.setMass(me);
        particle.setCharge(+e);
    }
    else if(pid == +2212){
        particle.setPID(+2212);
        particle.setName("proton");
        particle.setMass(mp);
        particle.setCharge(+e);
    }
    if(rel == 0) particle.setRelativistic(false);
    else particle.setRelativistic(true);

    return particle;
}

ROOT::RVec<TString> getTrees(TString fileName){
    ROOT::RVec<TString> trees;
    std::unique_ptr<TFile> file(TFile::Open(fileName));
    for (auto&& keyAsObj : *file->GetListOfKeys()){
        auto key = (TKey*) keyAsObj;
        trees.push_back(key->GetName());
    }
    return trees;
}

Double_t norm(ROOT::RVec<Double_t> vector){
    return TMath::Sqrt(TMath::Power(vector[0],2) +
                       TMath::Power(vector[1],2) +
                       TMath::Power(vector[2],2));
}

ROOT::RVec<Double_t> magneticMoment(Particle particle,
                                    ROOT::RVec<ROOT::RVec<Double_t>> velocity,
                                    ROOT::RVec<ROOT::RVec<Double_t>> field)
{
    ROOT::RVec<Double_t> vx = velocity[0];
    ROOT::RVec<Double_t> vy = velocity[1];
    ROOT::RVec<Double_t> vz = velocity[2];

    ROOT::RVec<Double_t> Bx = field[0];
    ROOT::RVec<Double_t> By = field[1];
    ROOT::RVec<Double_t> Bz = field[2];

    ROOT::RVec<Double_t> moment(vx.size());

    for(Int_t i=0; i<vx.size(); i++){

        ROOT::RVec<Double_t> v{vx[i], vy[i], vz[i]};
        ROOT::RVec<Double_t> B{Bx[i], By[i], Bz[i]};

        Double_t vMag = norm(v);
        Double_t bMag = norm(B);

        Double_t theta = TMath::ACos(ROOT::VecOps::Dot(v,B)/(vMag*bMag));

        Double_t vPerpendicular = vMag*TMath::Sin(theta);
        Double_t vParallel = vMag*TMath::Cos(theta);

        if(particle.isRelativistic()){
            moment[i] = particle.getMass()*TMath::Power(vPerpendicular,2)/(2*bMag);
        }
        else{
            Double_t gamma = TMath::Sqrt(1/(1-TMath::Power(vMag,2)/TMath::Power(c,2)));
            moment[i] = gamma*particle.getMass()*TMath::Power(vPerpendicular,2)/(2*bMag);
        }
    }
    return moment;
}

In regard to the ROOT file, I see that I cannot attach it here due to maximum size limitation. I therefore added the file in my Drive: https://drive.google.com/file/d/1PQVAx7aY0Kkw5iLHwWArPIPef5zUAf8F/view?usp=drive_link

Being the file in the same folder as the script, it should just run as

root -l dataAnalysis.C

Finally, in regard to reaching a limitation of the system, I’d guess (definitely not sure) that this is not the case, since I can do the same plot in python’s matplotlib without any issues.

Many thanks for the support.

5

In my case, on Mac, it does not crash but I see the memory usage going up to 745.3 MB and then root.exe is marked as “Not responding” by the activity monitor… It looks like more a infinite loop.

1 Like
6

@couet thanks.

I think by “crash” I mean the same as you’re experiencing. The graph window opens but no data is plotted, and even my terminal freezes when I try to cancel the job (I need to force close it).

When you mention an “infinite loop”, would that be an issue from my own code possibly?

And thanks for the feedback so far.

7

Problem found,

you need:

g1->Draw("AL");

The option “C” does a smoothing on the data points. Given the large number of points the algorithm does not succeed an never finishes or crashes. In your case smoothing, on a such large number of points, makes no real sense. Option “L” should be use, it simply draw a line.

1 Like
8

Oh my…super thanks!

After applying the correction, it opened the graph really quick now, with no issues.

That’s what I get for not reading some options in detail when using a function :sweat_smile:

Many thanks @couet.

1 Like
9

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.