Sum of TF1 in a loop

Dear,

In a loop of events, I want to calculate a Gaussian for each event, using the corresponding mean and width of the event, but also summing all the gaussians after each event. This is what I tried (naively), which doesn’t work:

Gaus_i = ROOT.TF1(“Gaus_i”,“gaus”,100, 200)
Gaus_tot = ROOT.TF1(“Gaus_tot”,“gaus”,100, 200)
i = 0

for iev in chain:
i+= 1
Gaus_i.Clear()
Gaus_i.SetParameters(i, iev.mean, iev.width)
Gaus_tot = ROOT.TF1(“Gaus_tot”,“Gaus_tot + Gaus_i”,100, 200)

Is there any solution to this problem? Btw is the use of TF1 here the best way, or should I use something else?

Many thanks!

Hi,

It is not clear to me exactly what you wnat to do it. It looks like you want to build a parametrize function as sum of many gaussian, each one for each single event. Is this really what you need ?
If you just want just to evaluate a Gaussian function on one event using a TF1 is not the optimal way to do it. You should use for example the TMath::Gaus function and do something like

gaus_tot = 0.
for iev in chain: 
   i+= 1
   gaus_tot += ROOT.TMath.Gaus(x, iev.mean, iev.width, true)  ## use true if you want a normalized gaussian

but it is not clear to me on what you are evaluating on, what is x ?

Lorenzo

Hi Lorenzo,

My apologies for the vague description, but I want to construct a PDF that is the sum of gaussians with as mean the value of the parameter I’m considering and the resolution as the standard deviation. Then I want to compare this resulting distribution to a certain histogram.

is the x here the amplitude of the Gaussian? Would it make sense if I equate it to the weight of the event?

Many thanks for you help!

Hi,
Apologies, but still I don’t fully understand what you want to constrict. Do you want to build a PDF as sum of gaussian kernels, something like this

For this we have already a class in ROOT doing that
https://root.cern.ch/doc/master/classTKDE.html

Lorenzo

Hi Lorenzo,

The piece of code you wrote is exactly what I want to calculate, I want to know that specifi distribution which is the sum of the gaussians, however I don’t know how to plot the resulting distribution, TMath returns just a number right? Because I want to compare it to a histogram

Hi,
You might want to create a function from that and plot. Maybe something like this

def my_function(x, par) : 
  gaus_tot = 0.
  for iev in chain: 
      i+= 1
      gaus_tot += ROOT.TMath.Gaus(x[0], iev.mean, iev.width, true) 
      return gaus_tot

f1 = ROOT.TF1('f1',my_function, 100, 200 , 0)
f1.Draw()

Hi Lorenzo,

This solution looks promising (and efficient), do you need to pass a variable here? i.e. “par” in the definition
I didn’t manage to run it

In the meantime I figured out a way by basically converting TF1’s to TH1’s then add the histograms:

for iev in chain1:
Gaus_i.Clear()
Gaus_i = ROOT.TF1(“Gaus_i”,“gaus”,100, 200)
Gaus_i.SetParameters(1., iev.mean, iev.sigma)
Hist_i = Gaus_i.GetHistogram()
Hist_tot.Add(Hist_i, iev.weight)

And this works, however on the edges of the histogram, peaks are appearing, which is due to accumulating events that doesn’t fall out the range. Any idea how to fix that?

Hi,

Sorry my code above was just an example. Attached here is a fully running code

For your solution, I would need a script to run, please provide me in case

Lorenzo

rootForum_39276.py (576 Bytes)