I am trying to find out how to use either TSpectrum or TSpectrumFit to find the peak areas and uncertainties for the peaks that are found using the TSpectrum::Search() function or equivalent.
I have been trying to use the example script FitAwmi.C (from root.cern.ch/root/htmldoc/TSpectrumFit.html) as a starting point but I am encountering problems that I can’t seem to work around.
I have an edited version of FitAwmi.C with my edits marked and attached a root file (pointed to by the edited version of FitAwmi.C) which illustrates the issues I am experiencing.
When I run the attached script I get the error message [quote]Error in TSpectrumFit::SetPeakParameters: Invalid peak position, must be in the range fXmin, fXmax[/quote]
I believe I am getting the message because my spectrum starts at a bin numbered 420 rather than 0. How do I fix this.
The marked peak centroid is around bin 420 rather than 480 (I expect that this is related to problem 1)
The output of the commands
pfit->GetAreas();
pfit->GetAreasErrors();
cout << CalcArea[0] << " - " << CalcAreaError[0] << endl;
is given as [quote]3.42942e+07 - 3.42944e+07[/quote]
These numbers appear to be wrong as they greatly exceed the number of counts in the histogram.
Separately - what commands do I use to get the peak areas and peak errors for peaks identified by TSpectrum or TSpectrumFit if GetArea() or GetAreasErrors() are not correct?
I have gone over those examples (peaks.C and peaks2.C) very carefully over the last week - they do not appear to give any indication of the areas of the peaks found or of the errors associated with those peak areas.
The scripts run without difficulties - they generate no errors - they fit the curves beautifully - they just do not appear to get me any closer to the answers that I am looking for.
I am sorry if the answer is just considered obvious - I do not understand how these tools can be used to give me peak areas and associated errors.
The obvious commands (such as GetAreas()) do not produce the output I am looking for.
Rene advised me to seek answers to my questions on the areas of the peaks found by TSpectrum::Search and the errors associated with those peak areas in the tutorial scripts peaks.C and peaks2.C.
I had gone over those examples (peaks.C and peaks2.C) very carefully before asking my original question and again since Rene’s advice - they do not appear to give any indication of the areas of the peaks found or of the errors associated with those peak areas.
The scripts run without difficulties - they generate no errors - they fit the curves beautifully - they just do not appear to get me any closer to the answers that I am looking for.
I am sorry if the answer is just considered obvious - I do not understand how these tools can be used to give me peak areas and associated errors.
I’m facing the same problem as you had, did you find the answer if yes then please help me to find the error in position and calculate area and error in the area using TSpectrum.
@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)
Ok. I need to check also how it will look like in the reference guide.
I may need to comment both lines as you suggested in order to have a clean presentation in the guide.
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 (integer heights)
… for the “Stiefel-Hestens” method …
.x SpectrumFit.cxx(1) ... or ... .x SpectrumFit.cxx++(1)
SpectrumFit(1); another random set of peaks (integer areas)
SpectrumFit(1, kFALSE); another random set of peaks (integer heights)
I feel that the SpectrumFit.cxx macro is worth using for my data analysis. I would like to know, whether there is a possibility of removing "Compton Background’’ before fitting spectra using the two methods that you have used. I know this is possible using the TSpectrum class.
There are various “${ROOTSYS}/tutorials/spectrum” which deal with some “background”.
Just replace the built-in histograms / spectra with your own histogram / spectrum and then play with them, e.g. changing background estimation methods / parameters.
Once you are satisfied with the result, you can add the chosen procedure to your macro.
I have modified the SpectrumFit.cxx macro to search for peaks in the spectrum of 152Eu source. With the modified macro( SpectrumFit_v2.cxx) I am able to search and fit the peaks in “AB01” histogram in the (hist.root) file. However, the quality of the fit is extremely poor. This is apparent from the errors in areas and peak-heights.
Will you be able to check the code and tell my whether I am doing anything wrong?
The TSpectrumFit, used in your macro, assumes that all peaks have (almost) the same sigma.
That’s not really your case … the fit procedure says that it found peaks with sigma = 1.88 +-0.66.
Let’s take three example peaks (“manually fitted”) … at around x= 245, sigma = 1.6, at around x = 2815, sigma = 2.3 and at around x = 5230, sigma = 3.8 so it really changes a lot.
This may introduce significant uncertainties in fitted peaks’ amplitudes and areas.
Well, the returned errors of peaks’ areas (0.8 - 2.6 %) seem “justified” to me (taking into account “sigma problems”). @moneta However, the returned errors of peaks’ amplitudes (bigger than the returned amplitudes themselves by a factor 1.1 - 3.3) look really ridiculous to me and I have no explanation for them.
You may try to play with the other macro (it will try to fit “independent” gaussians to all found peaks):