Finding the peak positions

[Sandy]

Dear @couet,

Is it possible to get the X, Y positions of the peak of a fit (it may be Gaussian or may not be)?

Does this example work for finding peak from fit too?
https://root.cern.ch/root/html534/tutorials/spectrum/peaks.C.html

[Wile_E_Coyote]

Try to play with this:

Trying to find peak areas using example FitAwmi.C

@axel The standard version of the ${ROOTSYS}/tutorials/spectrum/peaks.C script fits “peaks’ heights”, while the version attached below can fit “peaks’ areas” instead (and you can run this version with ROOT 5 or ROOT 6). peaks.C (4.3 KB)

and that:

Trying to find peak areas using example FitAwmi.C

The attached “SpectrumFit.cxx” macro fits 1D spectra using either the AWMI method (an algorithm without matrix inversion) or the Stiefel-Hestens method (a conjugate gradient algorithm) from the TSpectrumFit class. The TSpectrum class is used to find peaks. To try this macro, in a ROOT (5 or 6) prompt, do … … for the “AWMI” method … .x SpectrumFit.cxx ... or ... .x SpectrumFit.cxx++ SpectrumFit(); another random set of peaks (integer areas) SpectrumFit(0, kFALSE); another random set of peaks (…

[Sandy]

I meant that I have a fit already. Which is not a Gaussian. Can I get the X, Y centers of the its peak.
Like in 1D we have GetMean() for Gaussian which gives X mean.

Likewise for 2D fit (not Gaussian), is there any which give the value of X mean and Y mean?

[Wile_E_Coyote]

I’m sure your own fit function has some “parameters”. So you need to use the fitted parameters’ values to calculate what you want (it may well be that one of the “parameters” is simply the “mean”).

TF1, TF2 and TF3 functions do not have any “GetMean” method. You are mixing here “functions” with “histograms” and / or “graphs”.

On the other hand, TF1, TF2 and TF3 functions do provide various “Mean...”, “Moment...” and “CentralMoment...” methods, which you may find convenient.

[Sandy]

The problem with MeanX or MeanY is it does not take function with the form

fit= new TF2 (“fit”,func,0.,1.,0.,1.)
where “func” is a pre-defined function.

[Wile_E_Coyote]

Sorry, no idea where you got “MeanX” and “MeanY” from (i.e. it seems to me that you again try to apply some methods to an inappropriate object).

[Sandy]

Suppose I have a TF2. How do I get Mean, RMS both in X and Y direction?

[Wile_E_Coyote]

If by “RMS” you mean the “standard deviation” then you could use the “sqrt(second central moment)” to calculate it.

If by “mean”, you mean the “expected value” then you could use the “first (raw) moment” to calculate it.

[Sandy]

Does that mean

RMS = SQRT( myfunction->TF2::CentralMoment2());

and expected value = myfunction->TF2::Moment2();

[Wile_E_Coyote]

Try to play, e.g. with this function:

{
  TF2 *f = new TF2("f", "[0] * TMath::Gaus(x, [1], [2]) * TMath::Gaus(y, [3], [4])", -10., 10., -20., 20.);
  // f->Print();
  f->SetParameters(10., 2., 3., 4., 5.);
  f->Draw("colz");
}

[Sandy]

Are you asking to use TF2:moment2() or TF2:CentralMoment2() with your bi Gaussian?

[Wile_E_Coyote]

Yes, do try to apply various class methods to some easy well defined function, with precisely known parameters’ settings. This way you will learn how to use them properly (and what they really return).

[Sandy]

The following is giving error (Error: Can’t call TF2::CentralMoment2(-10.,10.,-20.,20.) in current scope):

  TF2 *f = new TF2("f", "[0] * TMath::Gaus(x, [1], [2]) * TMath::Gaus(y, [3], [4])", -10., 10., -20., 20.);
  f->SetParameters(10., 2., 3., 4., 5.);
  f->Draw("Surf");

  Double_t RMS = SQRT(f->CentralMoment2(-10., 10., -20., 20.));
  Double_t MeanY = f->Moment2(-10., 10., -20., 20.);

[Wile_E_Coyote]

See: TF2::CentralMoment2

See: TF2::Moment2

I don’t know where you take the “SQRT” from but, I would use TMath::Sqrt there, thought it’s just my private preference.

[Sandy]

Sorry.
Now when I do this:

  Double_t RMS1 = TMath::Sqrt(f->CentralMoment2(1,-10., 10., 1, -20., 20., 0.000001));
  Double_t Mean1 = f->Moment2(1,-10., 10., 1, -20., 20., 0.000001);
  Double_t RMS2 = TMath::Sqrt(f->CentralMoment2(2,-10., 10., 2, -20., 20., 0.000001));
  Double_t Mean2 = f->Moment2(2,-10., 10., 2, -20., 20., 0.000001);

I get the following:
RMS1 =-nan
Mean1 = 7.84092
RMS2 =14.7077
Mean2 = 512.375

I don’t know among those which is for X direction and which is for Y direction?
And what are the errors on the calculated values?