Help with TF2 integral

I’m not sure if It’s the right place to ask this question.
I’m trying to write a TF2 function for integration purposes. It is something like this:

Integral = \int_0^umx du u f(u) \int_fmin(u)^fmax(u) dE crossection(E)

If the energy integral had some fix limits, I could just create the TF2 and integrate on those limits. The problem here is that for each u I’ll have different boundaries for the energy function.

I could write it in python like this:

	def integrand(E,u):
		fu = fCrossInterp(u)
		integrand = crossSection(element, m_A, E) * u * fu

		return integrand

	# Calculate the intersection uInt of eMin and eMax given a specific rIndex
	uInt = EminEmaxIntersection(element, m_X)

	uLow = 0
	uHigh = min(uInt, V_gal) 
	eLow = lambda u: eMin(u, m_X)
	eHigh = lambda u: eMax(element, m_X,  u)
	integral = integrate.dblquad(integrand, uLow, uHigh, eLow, eHigh)[0]
	return integral

Sadly, I can’t use python since the code here will be used to generate some events for a geant4 simulation. In this case, is there anything similar to solve this for a TF2 function or anyother thing in cling?
Worst case scenario, I’ll turn it in a numerical integration.

Anyway, thanks in advance.
Bests.

Hi,

I don’t see the problem, you can define a similar function in C++ and use it directly with the Integrator class or via TF1. You would need to compute as an Integral of an Integral. I think Python, as seeing in scipy.integrate.dblquad — SciPy v1.6.3 Reference Guide does the same internally.
See the reference doc, ROOT: ROOT::Math::IntegratorOneDim Class Reference.
If you need I can provide you a C++ example for it

Cheers

Lorenzo

Hi Lorenzo,

I’ll try to read the reference. I would like an example since I’m not used to c++ class syntax.

Thank you!
bests,
Ricardo

So, I’ve tried some functions. For instance, first I tried to compute

Int_0^2 dx Int_0^2 dy x*y

using this code:

double func2(const double *x){

	return x[0]*x[1];

}

void integral(){

ROOT::Math::Functor f2(& func2, 2);
double a[2] = {0,2};
double b[2] = {1,2};

ROOT::Math::IntegratorMultiDim ig1(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig1.SetFunction(f2);
double val = ig1.Integral(a,b);

cout<<val<<endl;
}

It worked. Them, I tried to change the y integration limits to

Blockquote int_fmin(x)^fmax(x)

like in the code, to simulate the lambda function in the python code:

**Python code**

eLow = lambda u: eMin(u, m_X)
	eHigh = lambda u: eMax(element, m_X,  u)
	integral = integrate.dblquad(integrand, uLow, uHigh, eLow, eHigh)[0]
**cling**

double fmin(x){
	
	return x+1;
	
}

double fmax(x){

	return 2*x+2;

}



void integral(){

ROOT::Math::Functor f2(& func2, 2);
double a[2] = {0,fmin(x[0])};
double b[2] = {1,fmax(x[0])};

ROOT::Math::IntegratorMultiDim ig1(ROOT::Math::IntegrationMultiDim::kADAPTIVE);
ig1.SetFunction(f2);
double val = ig1.Integral(a,b);

cout<<val<<endl;
}

this didn’t worked. If you could provide any help, I would appreciate.

Thank you,
Bests
Ricardo

Ah, I’ve figured it

I can do something like this right:

double func(double* x, double* p) {
   return x[0];
}

double func2(double* x, double* p) {
    double b = x[0];
    double xmin = x[1]+1;
    double xmax = 2*x[1] + 2;
    TF1 f("func", func, xmin, xmax);;
    return x[1]*f.Integral(xmin, xmax);
}

and then, integrate the TF1 func2 in the range I’m intrested right? I’ll test it too

Thank you
Ricardo

Hi,
Yes, this could work fine

Cheers

Lorenzo

It worked,

thank you for the patience

Bests,
Ricardo