Ytdyklly

I am working on a simulation with a high number of particles, so my code must be the fastest as possible. I am new to C++ and ROOT (I am learning at the same time I work on my code), so I don't know any techniques for optimization.

My code solves a system of differential equations using GSL libraries. I also need ROOT to obtain the initial values for my particles.

At first, the computational times weren't too big, but I needed to introduce a function that checks the position of the particle at any time, and it returns a value after a for loop that is being computed each time I call the function. This function is called a lot of times, since my steps are very small, and this function appear in all the equations of the system.

This is my code:

#include "Rtypes.h"

R__LOAD_LIBRARY(libgsl)

R__LOAD_LIBRARY(gsl) 



#include <iostream>

#include <math.h>

#include <iomanip>

#include <time.h> 

#include <stdio.h>

#include <vector>

#include <fstream>



#include <TTree.h>

#include <TFile.h>





#include <gsl/gsl_odeiv2.h>



using namespace std;



class DataHolder

{

    DataHolder()

    {

    ifstream inFile;



    inFile.open("data.txt");



    double v1, v2;



    while(inFile >> v1 >> v2){

        a.push_back(v1); 

        b.push_back(v2 - 30.4); 

    }

    }



    public:



        static DataHolder& getInstance()

        {

            static DataHolder d;

            return d;

        }

        vector<double> a, b;

        vector<double> alpha() {return a;};

        vector<double> field() {return b;};



};













double* E(double x, double y, double z)

{

    static double Efield[3] = {0.}; 





    static const double e_const = 1.6021766208e-19;

    static const double m_const = 1.883531594e-28;

    static const double c_const = 299792458;

    static const double a_const = 11659208.0 * 1e-10;

    static const double Gamma_magic = sqrt(1. + (1. / a_const));

    static const double beta_magic = sqrt(1. - (1. / (Gamma_magic*Gamma_magic)));

    static const double B_const = 1.4513;



    static const double R_const = Gamma_magic * m_const * beta_magic*c_const / (e_const*B_const);



    static double x_radial; 

    x_radial = sqrt(x * x + z * z) - R_const;



    static double theta;

    theta = atan2(z, x)*180/M_PI;



        if (atan2(z, x) < 0) theta = 360 + theta;



    static double r_local;

    r_local = hypot(x_radial, y);





    static double theta_local;

    theta_local = atan2(y, x_radial);



    static double normal = {-28.6*1.e3, 0, -0.89*1.e3, 0., 0.564*1.e3, 0., 0.135*1.e3, 0., 0.7894*1.e3, 0., 0.0031*1.e3, 0., -0.0523*1.e3}; 



    static int array_length = sizeof(normal)/sizeof(normal[0]);



    static const double f_q = 4.*(13+26)/360;

    static const double n_index = 0.108; 

    static const double r0_2_ = 0.045*0.045; 



    static const double R = - n_index / (2*R_const / (beta_magic * c_const) / B_const / r0_2_ * f_q)*1/normal[0]; 



    double E_r_local = 0.;

    double E_theta = 0;



    double Er, Ey;



        if ((theta > 16 && theta < 30) || (theta > 33 && theta < 58) || (theta > 106 && theta < 119) || (theta > 123 && theta < 149) || (theta > 196 && theta < 210) || (theta > 214 && theta < 239) || (theta > 286 && theta < 300) || (theta > 304 && theta < 329))

      {



        for(int n = 2; n < array_length + 2; n++){

          E_r_local = E_r_local - n * pow(r_local, n-1)/pow(0.045, n)*normal[n-2]*R*cos(n*theta_local);

          E_theta = E_theta + n * pow(r_local, n-1)/pow(0.045, n)*normal[n-2]*R*sin(n*theta_local);

        }



        Er = E_r_local*cos(theta_local) - E_theta*sin(theta_local);

        Ey = E_r_local*sin(theta_local) + E_theta*cos(theta_local);



        Efield[0] = Er * x / sqrt(x * x + z * z);

        Efield[1] = Ey;

        Efield[2] = Er * z / sqrt(x * x + z * z);

      }



    else {

        Efield[0] = 0;

        Efield[1] = 0;

        Efield[2] = 0;

    }



    return Efield;

}







double* B(double x, double y, double z)

{

    static double Bfield[3] = {0.}; 



    static const double B_const = 1.4513;







    auto& d = DataHolder::getInstance();

    vector<double> pos;

    vector<double> B_radial_data;

    pos=d.alpha();

    B_radial_data=d.field();



    static double beta = 359.9; 

    static int entry;

    static double B_radial, alpha;



    alpha = atan2(z, x)*180/M_PI;



    if (atan2(z, x) < 0) alpha = alpha + 360;



    if (alpha < beta) entry=0;



    while (alpha > pos.at(entry)) entry++;



    B_radial = 1.e-6*B_radial_data.at(entry-1)*B_const + (alpha-pos.at(entry-1)) * (B_radial_data.at(entry) - B_radial_data.at(entry-1))/(pos.at(entry) - pos.at(entry-1))*1.e-6*B_const; 

    beta = alpha;





    Bfield[0] = B_radial * x / sqrt(x * x + z * z);

    Bfield[1] = B_const;

    Bfield[2] = B_radial * z / sqrt(x * x + z * z);



    return Bfield;

}







double dot(double vect_A, double vect_B)

{

    double product = 0;



    for (int i = 0; i < 3; i++)

        product = product + vect_A[i] * vect_B[i];

    return product;

}



double* cross(double vect_A, double vect_B)

{

    static double cross_P[3] = {};



    cross_P[0] = vect_A[1] * vect_B[2] - vect_A[2] * vect_B[1];

    cross_P[1] = vect_A[2] * vect_B[0] - vect_A[0] * vect_B[2];

    cross_P[2] = vect_A[0] * vect_B[1] - vect_A[1] * vect_B[0];



    return cross_P;

}







struct const_type {

     double e, m, c, g, Gamma; 

};





int

func(double t, const double var, double dvar, void *params)

{

    (void)(t); /* avoid unused parameter warning */

    const_type *my_params_pointer = (const_type *)params; 

     double e = my_params_pointer->e;

     double m = my_params_pointer->m;

     double c = my_params_pointer->c;

     double g = my_params_pointer->g;



    double beta[3] = { var[0], var[1], var[2] };

    double Spin[3] = { var[3], var[4], var[5] };



    static double Gamma_Global;

    Gamma_Global = 1. / sqrt(1. - dot(beta, beta));





    dvar[0] = e / (m*Gamma_Global*c)*(E(var[6], var[7], var[8])[0] + c * cross(beta, B(var[6], var[7], var[8]))[0] - beta[0] * dot(beta, E(var[6], var[7], var[8]))); 

    dvar[1] = e / (m*Gamma_Global*c)*(E(var[6], var[7], var[8])[1] + c * cross(beta, B(var[6], var[7], var[8]))[1] - beta[1] * dot(beta, E(var[6], var[7], var[8]))); 

    dvar[2] = e / (m*Gamma_Global*c)*(E(var[6], var[7], var[8])[2] + c * cross(beta, B(var[6], var[7], var[8]))[2] - beta[2] * dot(beta, E(var[6], var[7], var[8]))); 

    dvar[3] = e / m * ((g / 2.0 - 1.0 + 1.0 / Gamma_Global)*cross(Spin, B(var[6], var[7], var[8]))[0] - (g / 2.0 - 1.0)*Gamma_Global / (Gamma_Global + 1.0)*cross(Spin, beta)[0] * dot(beta, B(var[6], var[7], var[8])) - 1 / c * (g / 2.0 - Gamma_Global / (Gamma_Global + 1.0))*(beta[0] * dot(Spin, E(var[6], var[7], var[8])) - E(var[6], var[7], var[8])[0] * dot(Spin, beta))); 

    dvar[4] = e / m * ((g / 2.0 - 1.0 + 1.0 / Gamma_Global)*cross(Spin, B(var[6], var[7], var[8]))[1] - (g / 2.0 - 1.0)*Gamma_Global / (Gamma_Global + 1.0)*cross(Spin, beta)[1] * dot(beta, B(var[6], var[7], var[8])) - 1 / c * (g / 2.0 - Gamma_Global / (Gamma_Global + 1.0))*(beta[1] * dot(Spin, E(var[6], var[7], var[8])) - E(var[6], var[7], var[8])[1] * dot(Spin, beta))); 

    dvar[5] = e / m * ((g / 2.0 - 1.0 + 1.0 / Gamma_Global)*cross(Spin, B(var[6], var[7], var[8]))[2] - (g / 2.0 - 1.0)*Gamma_Global / (Gamma_Global + 1.0)*cross(Spin, beta)[2] * dot(beta, B(var[6], var[7], var[8])) - 1 / c * (g / 2.0 - Gamma_Global / (Gamma_Global + 1.0))*(beta[2] * dot(Spin, E(var[6], var[7], var[8])) - E(var[6], var[7], var[8])[2] * dot(Spin, beta))); 

    dvar[6] = var[0] * c; 

    dvar[7] = var[1] * c; 

    dvar[8] = var[2] * c; 



    return GSL_SUCCESS;



}



int main(int argc, char *argv)

{





    const  double e_const = 1.6021766208e-19;

    const  double m_const = 1.883531594e-28; 

    const  double m_GeV = m_const / (1.e9 * 1.782661907*1e-36); 

    const  double c_const = 299792458;

    const  double a_const = 11659208.0 * 1e-10;

    const  double g_const = 2.+2.*a_const; 

    const  double Gamma_magic = sqrt(1. + (1. / a_const)); 

    const  double beta_magic = sqrt(1. - (1. / (Gamma_magic*Gamma_magic)));

    const  double B_const = 1.4513;

    const  double R_const = Gamma_magic * m_const * beta_magic*c_const / (e_const*B_const);

    const double pi = M_PI;









    struct const_type my_const = { e_const , m_const , c_const, g_const};



    const  double S_0 = 1.0;



    size_t dimension = 9; 



    double eps_abs = 1.e-12;    /* absolute error requested */

    double eps_rel = 1.e-12;    /* relative error requested */



    /* define the type of routine for making steps: */

    const gsl_odeiv2_step_type *type_ptr = gsl_odeiv2_step_rk4;



    /*

    allocate/initialize the stepper, the control function, and the

    evolution function.

    */

    gsl_odeiv2_step *step_ptr = gsl_odeiv2_step_alloc(type_ptr, dimension);

    gsl_odeiv2_control *control_ptr = gsl_odeiv2_control_y_new(eps_abs, eps_rel);

    gsl_odeiv2_evolve *evolve_ptr = gsl_odeiv2_evolve_alloc(dimension);



    gsl_odeiv2_system sys = { func, NULL, dimension, &my_const };



    double var[9]; 



    double t, t_next; 



    double tmin = 0.;           /* starting t value */

    double tmax = 1.e-4;            /* final t value */

    double delta_t = 1.e-9; 



    double h = 1.e-12; 



    TFile *fileinput = new TFile("fileinput.root", "READ"); 

    TTree *b_tree = (TTree*)(fileinput->Get("file_tree"));



    int NPart = b_tree->GetEntries(); 

    cout << " Number of particles = " << NPart;

    char name[20], file[200];

    /*

    vector<double> *x_0 = 0;

    vector<double> *y_0 = 0;

    vector<double> *Px_0 = 0;

    vector<double> *Py_0 = 0;

    vector<double> *Pz_0 = 0;

    vector<double> *Sx_0 = 0;

    vector<double> *Sy_0 = 0;

    vector<double> *Sz_0 = 0;

    */



    double x_0, y_0, Px_0, Py_0, Pz_0, Sx_0, Sy_0, Sz_0;

    b_tree->SetBranchAddress("x_radial", &x_0);

    b_tree->SetBranchAddress("y", &y_0);

    b_tree->SetBranchAddress("P_x", &Px_0);

    b_tree->SetBranchAddress("P_y", &Py_0);

    b_tree->SetBranchAddress("P_z", &Pz_0);

    b_tree->SetBranchAddress("S_x", &Sx_0);

    b_tree->SetBranchAddress("S_y", &Sy_0);

    b_tree->SetBranchAddress("S_z", &Sz_0);



    int node_number = atoi(argv[1]);



    int nJobs=16*16-1; 



    int Quotient = NPart/nJobs;

    int Reminder = NPart%nJobs;



    sprintf(file, "result_%i.root", node_number+1);



    TFile *fileout = new TFile(file, "RECREATE"); 



    for (int k = node_number * Quotient; k < NPart - (2*nJobs - node_number - 2)/nJobs * (nJobs - node_number - 1)*Quotient - (2*nJobs - node_number - 1)/nJobs*Reminder; k++)

    {

        b_tree->GetEntry(k);

        sprintf(name, "Part%i",k+1);





        double Gamma_0 = sqrt(1 + (Px_0*Px_0 + Py_0*Py_0 + Pz_0*Pz_0)*1.e-6/(m_GeV * m_GeV));



        var[0] = Pz_0/(Gamma_0 * m_GeV)*1.e-3*sin(89*pi/180) + Px_0/(Gamma_0 * m_GeV)*1.e-3*cos(89*pi/180); 

        var[1] = Py_0/(Gamma_0 * m_GeV)*1.e-3;

        var[2] = -Px_0/(Gamma_0 * m_GeV)*1.e-3*sin(89*pi/180) + Pz_0/(Gamma_0 * m_GeV)*1.e-3*cos(89*pi/180); 

        var[3] = Sz_0*sin(89*pi/180) + Sx_0*cos(89*pi/180);

        var[4] = Sy_0;

        var[5] = sqrt(S_0*S_0 - var[3]*var[3] - var[4]*var[4]); 

        var[6] = R_const + x_0*1.e-3; 

        var[7] = y_0*1.e-3;

        var[8] = 0.;



        double beta_module = sqrt(var[0]*var[0]+ var[1] * var[1]+ var[2] * var[2]);

        double S_module = sqrt(var[3] * var[3] + var[4] * var[4] + var[5] * var[5]);

        double x_e = sqrt(var[6]*var[6]+var[8]*var[8])-R_const;

        double x_radial = sqrt(var[6] * var[6] + var[8] * var[8]) - R_const;

        double Projection = 1. / (beta_module * S_module) * (var[0] * var[3] + var[1] * var[4] + var[2] * var[5]);

        double Vertical_focusing = E(var[6], var[7], var[8])[1];



        double S_module0 = S_module; 





        TTree *tree = new TTree(name, name); 

        tree->Branch("t", &t, "t/D");

        tree->Branch("beta_y", &var[1], "beta_y/D");

        tree->Branch("x", &var[6], "x/D");

        tree->Branch("y", &var[7], "y/D");

        tree->Branch("z", &var[8], "z/D");

        tree->Branch("x_e", &x_e, "x_e/D"); 

        tree->Branch("Projection", &Projection, "Projection/D"); 

    tree->Branch("Vertical_focusing", &Vertical_focusing, "Vertical_focusing/D"); 



        t = tmin; 



        tree->GetEntry(0);

        tree->Fill();



        printf("%.5e %.5e  %.5e  %.5e  %.5e  %.5e  %.5e  %.5e  %.5e %.5en", t, var[0], var[1], var[2], var[3], var[4], var[5], var[6], var[7], var[8] );   /* initial values */



        int i = 1;      



        /* step to tmax from tmin */

        for (t_next = tmin + delta_t; t_next <= tmax; t_next += delta_t)

        {



            while (t < t_next)  

            {

                gsl_odeiv2_evolve_apply(evolve_ptr, control_ptr, step_ptr,

                    &sys, &t, t_next, &h, var); 

            }



            beta_module = sqrt(var[0]*var[0]+ var[1] * var[1]+ var[2] * var[2]);







            S_module = sqrt(var[3] * var[3] + var[4] * var[4] + var[5] * var[5]);



            var[3] = var[3] * S_module0 / S_module;

            var[4] = var[4] * S_module0 / S_module;

            var[5] = var[5] * S_module0 / S_module;



            S_module = sqrt(var[3] * var[3] + var[4] * var[4] + var[5] * var[5]); 



            Projection = 1 / (S_module*beta_module) * (var[0] * var[3] + var[1] * var[4] + var[2] * var[5]);

        Vertical_focusing = E(var[6], var[7], var[8])[1];



        x_radial = sqrt(var[6] * var[6] + var[8] * var[8]) - R_const;

        x_e = (t*1.e9*x_e + x_radial) * 1 / ((t + delta_t)*1.e9); /







        if(sqrt(x_radial * x_radial + var[7] * var[7]) > 0.045) 

          {cout << "Radius = " << sqrt(x_radial * x_radial + var[7] * var[7]) << " Particle = " << k << " time = " << t <<  endl;

        tree->Fill(); 

        break; 

          }









            tree->GetEntry(i);

            tree->Fill();

            i++;

        }

        tree->Write();





        }



    fileout->Close();



    /* all done; free up the gsl_odeiv2 stuff */

    gsl_odeiv2_evolve_free(evolve_ptr);

    gsl_odeiv2_control_free(control_ptr);

    gsl_odeiv2_step_free(step_ptr);





    return 0;

}

I am pretty sure that the high computation time is due to the functions E and B, but some experts will find mistakes at other places.

Let me know if something is unclear or you need more information.