Convolution of two functions yielding Inf/NaN

[Laura_Munteanu]

Hello,

I’m trying to fit a waveform with the convolution of 2 functions. I’ve attempted to use TF1Convolution in the following way:

TF1* tf1_A = new TF1("tf1_A", "[2]*(x>[0])*TMath::Exp(-1*(x>[0])*(x-[0])*3/[1])*TMath::Power((x>[0])*(x-[0])/[1],3)*TMath::Sin((x>[0])*(x-[0])/[1])");
tf1_A->SetParameter(0, 9000.);
tf1_A->SetParameter(1, 600.);
tf1_A->SetParameter(2,1);

TF1* tf1_L = new TF1("tf1_L", "gaus(0)");
tf1_L->SetParameter(0, 1);
tf1_L->SetParameter(1, 9000);
tf1_L->SetParameter(2, sigma_in);

TF1Convolution* tf1conv_AL = new TF1Convolution("tf1_L", "tf1_A", 0.0, tmax, true);
tf1conv_AL->SetNofPointsFFT(10000);
 TF1   *f = new TF1("f",*tf1conv_AL, 0.01, tmax, tf1conv_AL->GetNpar());
   f->SetParameters(9000, 600, 1, 1, 9000, sigma_in);

TCanvas* c_tf = new TCanvas();
//tf1_A->Draw();
//tf1_L->Draw("same");
f->Draw();

with sigma_in=16 and tmax = 40000. The two separate functions are drawn correctly, but when I try to draw the convolution I get the following message:

Warning in <TCanvas::ResizePad>: Inf/NaN propagated to the pad. Check drawn objects.
Warning in <TCanvas::ResizePad>: c1 height changed from 0 to 10

Could you help me figure out what is wrong with the way I’m doing the convolution?

Thanks,
Laura

[bellenot]

@moneta can most probably help you

[moneta]

Hi,

I am not having a problem using code similar to yours. Can you tell me exactly what are the function parameters. so I can run exactly your code and reproduce your problem ?
I suspect is a numerical issue in the function calculation

Lorenzo

[Laura_Munteanu]

Hi Lorenzo,

The functions I’m using so far are:

• a Gaussian meant to describe the longitudianal diffusion of a cluster of electrons (described by tf1_L in my code) with a width of sigma_in = 16 (ns).
• the product of a sine function, an exponential decay and a 3rd degree polynomial to describe the electronics shaping effects. This is a function of (x-(par[0])/par[1] where par[0] is supposed to be the offset (~9000 ns in my case) and par[1] is the time constant of the decay and is of the order of 600 (ns). I also make sure that the function is only defined for x>par[0] by putting in that condition. par[2] is just a normalization constant so far, I put it equal to 1, so that it’s comparable to the gaussian.
  (I will further have to convolute the result of this convolution with a third function containing TMath::Erf).
  These functions should be convoluted in the range 0…40000 (ns).

Otherwise I’m just using the parameters in the code above.

Laura

[moneta]

Hi,
Can you please attach a running code showing this problem ? This would be helpful to understand the problem

Lorenzo

[Laura_Munteanu]

Hi,

Sure, here is the code. If I don’t do the convolution, the two functions get drawn as they should.
testROOTTF1Convolution.C (1.3 KB)

Laura

[moneta]

Thank you for posting the code. I can reproduce it now and I will investigate it

Lorenzo

[moneta]

Hi,

The problem you are inverting the function in the constructor. You should do:

	TF1Convolution* tf1conv_AL = new TF1Convolution("tf1_A", "tf1_L", 0.01, tmax, true);

where passing Ttf1_A as first function. Also I would use a smaller range for the convolution eventually since the functions are zero outside the peak

Lorenzo

[Laura_Munteanu]

Thank you very much, this works! However, shouldn’t convolution be commutative? What’s the driving factor in which function I should put first?

[moneta]

Hi,
Yes the convolution is commutative if you define it in the full -[-inf, +inf] range. Otherwise one needs to correct change the range when the order of the functions is reversed (e.g. if the gaussian is the first function)

Lorenzo