Double sum in Rcpp

Question

I have to calculate this double sum in Rcpp:

. Here, "t" and "m" are variables. "beta", "c" and "p" are the parameters. For which I wrote the following code. The first code is the functor which I am going to integrate. The second code is used for integrating. The third code is for calculating the double sum. I want to know if I am doing it correctly. Any help would be much appreciated:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>


using namespace Numer;
using namespace Rcpp;
using namespace std;


// this is the functor we are going to integrate
class hnr: public MFunc
{
  private:

   double beta;
   double p;
   double c;

   public:
      hnr(double beta_, double p_, double c_):  beta(beta_),  p(p_),  c(c_) {};

  double operator()(Constvec& x)
   {    

    return    R::dexp(x[0],1/beta,0)) * (1/(pow(x[1]+c,p))) ;


   }
 };


// [[Rcpp::export]]
double integrate_test3(  double beta, double p, double c,  Eigen::VectorXd lower, Eigen::VectorXd upper)
{
  // this is the function to integrate the above function.
  // the first observation of "lower" and "upper" vector
  // is for m; the second ones are for t
  hnr f(  beta,  p, c);
   double err_est;
   int err_code;
   int maxeval = 5000;
   double res = integrate ( f, lower, upper,err_est,err_code,maxeval);
   return res;

}

 // This is the final double sum
// [[Rcpp::export]]
 double partB2(std::vector<double> t, std::vector<double> m,         std::vector<double> theta){
   // extracting the parameter estimates


   double beta=theta[0];
   double c=theta[1];
   double p=theta[2];

  int n= ts.size();
  std::vector<double> temp2(n);

  for(int i=0; i <n; i++){
      std::vector<double> newts();
      newts.insert(newts.end(),t.begin(),t.begin()+i);

      std::vector<double> newms();
      newms.insert(newms.end(),m.begin(),m.begin()+i );

      int l= newts.size();
      double temp2in=0;
      for(int j=0; j <l; j++){

         Eigen::VectorXd lower(2);
         lower << 0, (t[i]-newts[j]);
         Eigen::VectorXd upper(2);
        upper << 9, (t[i+1]-newts[j]);
        temp2in += integrate_test3(  beta, p,  c,  lower,   upper);
      }

      temp2[i]= temp2in;
  }
  double partC = std::accumulate(temp2.begin(), temp2.end(), 0.0);
  return partC;
}

Am I doing it correctly?