TF1 and lambda change the values of a vector

Newbie
1

ROOT Version: 6.18/04
Platform: 18.04
Compiler: cc 7.5.0 (Ubuntu 7.5.0-3ubuntu1~18.04)

Hello co-rooters,

Here’s a weird one.
I’m defining a few TF1 objects using lambda functions.
Here’s my MWE.


#include "TF1.h"
#include "TH1D.h"
#include "TPad.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TGraph.h"
#include "TString.h"
#include "TGraphErrors.h"
#include "Math/Functor.h"
#include "Math/WrappedFunction.h"
#include "Math/IFunction.h"
#include "Math/Integrator.h"


#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>

using namespace std;

void cumulative_Gn0(){

    ROOT::Math::IntegratorOneDimOptions::SetDefaultIntegrator("AdaptiveSingular");
    
    //Let's do the thresholds
    double E_low  = 0.001;
    double E_high = 550.;
    double a      = 5.e-4;
    double b      = 1.70;
    
    TF1 *f_Porter_Thomas = new TF1("f_Porter_Thomas", "ROOT::Math::chisquared_pdf(x,1,0)",0,10);
    TF1 *f_threshold_gGn0 = new TF1("f_threshold_gGn0", "[0]*pow(x,[1])", E_low, E_high);
    f_threshold_gGn0->SetParameters(a, b);
    
    int    mobs                 = 26;
    double sum_Gn0              = 0.;
    double Gn0_average          = 3.63193;    
    double cmeas                = a/Gn0_average;    
    int    number_of_iterations = 2;
    double integral_tolerance   = 1.e-5;
    int    energy_points        = 1.e4;
    
    TF1 *f_alpha = new TF1("f_alpha", "[0]*pow(x,[1])" , E_low, E_high);
    f_alpha->SetParameters(cmeas, b);
    
    TF1 *f_p   = new TF1("f_p"  , [=](double *x, double *p){ return f_Porter_Thomas->Integral(0., f_alpha->Eval(x[0])); }, E_low, E_high, 0);
    TF1 *f_1mp = new TF1("f_1mp", [=](double *x, double *p) {return 1.-f_Porter_Thomas->Integral(0., f_alpha->Eval(x[0])); }, E_low, E_high, 0);
    TF1 *f_q   = new TF1("f_q"  , "2*sqrt( 0.5*f_alpha )*TMath::Exp( -0.5*f_alpha )/sqrt(TMath::Pi())", E_low, E_high);

    double c_vector[number_of_iterations],
           gn0_vector[number_of_iterations],
           M_vector[number_of_iterations],
           D_vector[number_of_iterations],
           SGn0_vector[number_of_iterations],
           S_vector[number_of_iterations];
    int    iteration_vector[number_of_iterations];

    TF1** f_alpha_vector = new TF1*[number_of_iterations];
    TF1** f_p_vector     = new TF1*[number_of_iterations];
    TF1** f_1mp_vector   = new TF1*[number_of_iterations];
    TF1** f_q_vector     = new TF1*[number_of_iterations];
    
    double c, gn0, M, D, SGn0, S;
    int    iteration;
    
    
    iteration     = 0;
    gn0_vector[0] = Gn0_average;
    c_vector[0]   = a/Gn0_average;
    
    f_alpha  ->SetParameters(c_vector[0], b);
    f_alpha_vector[0] = new TF1("falpha0", "f_alpha", E_low, E_high);
    f_p_vector[0] = new TF1("fp0"  ,[=](double *x, double *p){ return f_Porter_Thomas->Integral(0., f_alpha->Eval(x[0])); }, E_low, E_high, 0);
    f_1mp_vector[0] = new TF1("f_1mp", [=](double *x, double *p) {return 1.-f_Porter_Thomas->Integral(0., f_alpha->Eval(x[0])); }, E_low, E_high, 0);
    f_q_vector[0] = new TF1("fq0", "2*sqrt(0.5*f_alpha)*TMath::Exp(-0.5*f_alpha)/sqrt(TMath::Pi())", E_low, E_high);
    
    M_vector[0]    = mobs*( 1. + f_p_vector[0]->Integral(E_low, E_high, integral_tolerance)/((E_high-E_low)-f_p_vector[0]->Integral(E_low, E_high, integral_tolerance)) );
    D_vector[0]    = (E_high-E_low)/M_vector[0];
    SGn0_vector[0] = gn0_vector[0]*( f_p_vector[0]->Integral(E_low, E_high, integral_tolerance)-f_q_vector[0]->Integral(E_low, E_high, integral_tolerance) )/D_vector[0];
	S_vector[0]    = 0.;
	
	cout << "SGn0_vector[0] = " << SGn0_vector[0] << endl;
    //gn0_vector[0] = (sum_Gn0 + SGn0_vector[0]) / M_vector[0];
	cout << "SGn0_vector[0] = " << SGn0_vector[0] << endl;
	cout << "SGn0_vector[0] = " << gn0_vector[0]*( f_p_vector[0]->Integral(E_low, E_high, integral_tolerance)-f_q_vector[0]->Integral(E_low, E_high, integral_tolerance) )/D_vector[0] << endl;

}

The output of that codes is the following

SGn0_vector[0] = -122.11
SGn0_vector[0] = -122.11
SGn0_vector[0] = -122.11

If I uncomment the line gn0_vector[0] = (sum_Gn0 + SGn0_vector[0]) / M_vector[0]; my output changes to the following:


SGn0_vector[0] = -122.11
SGn0_vector[0] = -122.11
SGn0_vector[0] = 43.188

Not sure if it’s related to the way the TF1s are defined, but I can’t seem to be able to explain the behaviour.

Any ideas are more than welcome!

2

Hello,

A preliminary observation about your post: according to your scores on the forum, you have quite some experience with scientific coding and ROOT, perhaps your post are really useful to the community and having them in the newbie community (where they disappear after a while) is a pity! Of course, the decision is yours :slight_smile:

I cannot reproduce your numbers: on my computer I get

SGn0_vector[0] = -94.4302
SGn0_vector[0] = -94.4302
SGn0_vector[0] = -94.4302

Could you doublecheck? Then we can dive further. Thanks!

Cheers,
D