- #1

minits

- 13

- 0

- 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:

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

C:

```
#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;
```

Last edited by a moderator: