Monte-Carlo integration does not work properly after implementing pdfs

[minits]

TL;DR Summary

Hello everyone,

I calculated the total cross section of the u b > d t process by calculating the Matrix element and then did the phase space integration via the vegas Monte-Carlo algorithm. Everything worked but then I wanted to implement the pdfs of the u- and b-quark in order to obtain the hadronic cross section of this process. Now the calculated standard deviation is of the same order than the result of the integral and it depends on the given seed of my random number generator.

Does anyone have experience with such strange behavior in Monte-Carlo methods? I think it is a conceptional problem and I am just missing a key point in how to set up the integration instead of a error in my code itself. I use data files from LHAPDF and also checked that my variables give the desired outputs. Here is the relevant part of my code:

#include <iostream>
#include <cmath>
#include "vegas.h"
#include "ranlxd.h"
#include "eventn.cpp"
#include "LHAPDF/LHAPDF.h"

#define DIMENSION 4

extern "C" void ranlxd_(double x[],int & n){
  ranlxd(x,n);
}
extern "C" void rlxd_init_(int level,int seed){
  rlxd_init(level,seed);
}

double integrand(double x[4], double* wgt){
double p1[4],p2[4],pf[2][4],fwgt,jac; //declaration four-vectors and jacobian
double XA=x[0]; //assigning random values for x value of quark pdfs
double XB=x[1];

std::vector pdfa,pdfb; pdfa=LHAPDF::xfx(XA,173);
pdfb=LHAPDF::xfx(XB,173);
double upart=pdfa.at(8)/XA; //assigning the value of the pdf for the corresponding quark type and divide by the x since xfx returns x*f) double bpart=pdfb.at(11)/XB;

 //constants
const double sqrts = 13000;
const double s=sqrts*sqrts;
const double mf[2]={0,173};
const double g_w4 = pow(0.65323293034757990,4);
const double m_W2 = pow(80.419002445756163,2);
const double pi(3.14159265358979323846264338327950288419716939937510582097494459230);
const double prefac = 1/(4*pow(pi,2));

double z[2]; //assigning random values needed for the phase-space generator
z[0]=x[2];
z[1]=x[3];
eventn(2,sqrts,z,mf,pf,fwgt); //phase-space generator overwrites the values of the pf array (four-vectors in final state)
jac=fwgt/(2*s); // Jacobian of the transformation

//four-vectors in initial state
  p1[0]=sqrts/2;
  p1[1]=0;
  p1[2]=0;
  p1[3]=p1[0];
  p2[0]=sqrts/2;
  p2[1]=0;
  p2[2]=0;
  p2[3]=-p2[0];
  double k = *wgt;

//squared matrix element of u b > d t
double Msquared=g_w4*(s/2)*(pf[0][0]*pf[1][0]-pf[0][1]*pf[1][1]-pf[0][2]*pf[1][2]-pf[0][3]*pf[1][3])/ (pow(-2*(p1[0]*pf[0][0]-p1[1]*pf[0][1]-p1[2]*pf[0][2]-p1[3]*pf[0][3])-m_W2,2));

//integrand with jacobian and pdfs
double ds=upart*bpart*jac*Msquared*prefac; return ds;
 }

int main(int argc, char **argv){
LHAPDF::initPDFSet("NNPDF23_nlo_as_0119",LHAPDF::LHGRID,0);
double acc = 1.e-20;
double sigmas[1000];
double errors[100];
double uncertainties[1000];
double bins[1000];
int dim=DIMENSION;
int itmx=10; // iterations
int ncall=1000; // sampling points per iteration
int nprn = 1; // printlevel (0/1)
int doof = 0;
rlxd_init(1,42743); //seed for the random number generator
/* integrate over hyper cube [0,1]^DIMENSIONS */
vegas_(integrand,&acc,&dim,&ncall,&itmx,&nprn,&doof);
cout << "Integral = " << result_.s1<< endl;

Does anyone see a mistake I made when implementing the pdfs? Any help is highly appreciated!

 

[Future Bruno]

It is possible that you are not correctly implementing the PDFs in your Monte-Carlo integration. In order for your Monte-Carlo integration to work correctly, you must correctly account for the variance of the PDFs in your integrand. This can be done by incorporating the PDFs into the weight assigned to each sample point. Additionally, you must make sure that the PDFs are correctly normalized. If they are not correctly normalized, then your results will be incorrect.