TF1 as sum of many exponential functions (with single parameter)

[Piet]

Dear experts.

I would like to make a TF1 function that is the sum of many exponentials, which have all the same variable lambda that I would like to vary in a fit. This is the function I would like to implement:

with lambda > 0 and with H the Heavyside function. So the terms of the sum should be 0 for x<i and exp(-lambda x) for x>=i. For illustration:

Now, after some googling, I based myself on an example of the root forum back in 2008. I attach a sample program at the bottom of this topic that you can run by root -l SumOfNExp.C+. However I am experiencing some problems / instabilities and I also have some questions.

1. Even though an exp(-10x) drops quicker to zero than for instance exp(-1x), I obtain NAN in the calculation of the evaluation of the function f. I tried understand which term lead to the infinity and I found out it was quite random and depends on the value of lambda. For instance for the calculation of f(0) I printed out the the evaluation at 0 for the various components:

 f68(0) = 0 cumulative sum = 0
 f69(0) = 0 cumulative sum = 0
 f70(0) = 0 cumulative sum = 0
 f71(0) = nan cumulative sum = nan
 f72(0) = nan cumulative sum = nan
 f73(0) = nan cumulative sum = nan

2. Playing with lambda, I can find values for which I do not have this problem, but I would like to use lambda as a fit variable, so this is looking for some trouble i.m.h.o.
3. Will I be able to use this function to fit, and by defining the individual terms in a loop, will they all share the same variable [3] or will these parameters be all different? In the latter case, how can I lock them?

std::vector<TF1 *> func_vec;
for(int i=0; i<nbx; ++i) {

    std::string fncname = "f";
    fncname = fncname+std::to_string(i);

    TF1 * f = new TF1(fncname.c_str(),"(x<[0])*[1] + (x>=[0] && x<[0]+10)*[2]*exp(-1*(x-[0])*[3]) + (x>=[0]+10)*[1]",0,nbx-1);
    f->FixParameter(0,i);
    f->FixParameter(1,0);
    f->FixParameter(2,lumivec[i]);
    f->SetParameter(3,10); // 10 non funziona                                                                                                                                                                                                
    func_vec.push_back(f);
}

4. Is this the most efficient way to solve my problem, or are there other tricks to obtain the function I need that are more efficient and computationally effective?

Thanks a lot
Kind regards
Piet Verwilligen

SumOfNExp.C (5.3 KB)

---

ROOT Version: 6.32.12
Platform: macosx64
Compiler: apple clang version 17.0.0

---

[devajith]

Hi @Piet,Thank you for your question.

Perhaps @jonas could help you with this.

Cheers,
Dev

[Piet]

Thanks, I would really appreciate some help on this point. Maybe my email was a bit too long, so I summarize it here in a single line:

My TF1 is a finite sum of converging functions but I obtain a NAN value when I evaluate the function, which I believe should not happen.

Thanks
Piet

[couet]

I took the first function producing a Nan according to the ouput your code generates. That’s number 71. I did the following:

root [1] auto f = new TF1("f","(x<72) * 0 + (x>=72 && x<82) * 100 * exp( - ( x - 72) *  10) + (x>=82) * 0",0.,71.)
(TF1 *) 0x7faaec8151b0
root [2] f->Eval(0)
(double) nan
root [3] 

So this single function evaluated at 0 (like you do in your code) produces a Nan.

[Wile_E_Coyote]

TF1 *f = new TF1(fncname.c_str(), "(x<[0])*[1] + (x>=[0] && x<[0]+10.)*[2]*exp((x>=[0] && x<[0]+10.)*([0]-x)*[3]) + (x>=[0]+10.)*[1]", 0., nbx - 1.);

[Piet]

@Wile_E_Coyote

TF1 f = new TF1(fncname.c_str(), "(x<[0])[1] + (x>=[0] && x<[0]+10.)[2]exp((x>=[0] && x<[0]+10.)([0]-x)[3]) + (x>=[0]+10.)*[1]", 0., nbx - 1.);

Dear @Wile_E_Coyote so what you propose is some additional protection such that only in the range [0]<=x<[0]+10 the exponent is evaluated? Do I understand that even if the criterion (x>=[0] && x<[0]+10.) is not met, the exponent is still evaluated resulting in 0 * nan = nan? Implementing this fix it worked for me, thanks a lot.
Kind regards
Piet

[Piet]

I have a remaining side-question though:

Will I be able to use this function to fit, and by defining the individual terms in a loop, will they all share the same variable [3] or will these parameters be all different? In the latter case, how can I lock them?

Thanks
Piet

[Wile_E_Coyote]

You define the “fsum” with 0 parameters (you set npar = 0 in its constructor).
Try: fsum->Print();

You then use “fFuncList[i]->EvalPar(x,p)” so, ALL functions will be evaluated with exactly the same set of parameters (“p”).

As the “fsum” has no parameters, I assume “p” will always be 0 (so ALL functions will use their own preset parameters).