Hi Antonio,
You have a problem when filling your histogram. The bin numbering convention
is the following:
bin 0 : bin with underflows
bin 1 : first normal bin
bin N : last normal bin if you have created the histogram with N bins
bin N+1 : bin with overflows
In your code, modify your two lines filling the histogram with:
TH1F *hpx = new TH1F("hpx","Peaks",npts,-0.5,npts-0.5);
for (Int_t j=0 ; j<npts ; j++) hpx->SetBinContent(j+1,y[j]) ;
There was effectively an offset of one bin in the returned peak positions.
I have fixed this problem in the development version in CVS.
Note that the TSpectrum class has essentially functions expecting arrays
as arguments. You can use these functions directly. The function
TSpectrum::Search
is just a more convenient function when you have an histogram. In addition
it can display the position of the peaks.
Rene Brun
Antonio Kanaan wrote:
>
> Hi All,
>
> I am trying to use TSpectrum to produce a list of peaks in spectra. I
> produced a data file containing two peaks. This is my code's output:
>
> root [0]
> Processing search.cxx...
> Peaks found 2
> Peak at = 6.5
> Peak at = 19.5
> (int)0
> root [1]
>
> Almost fine weren't for the peak positions. They are actually at 7.5
> and 20.5. I must be missing something but don't know what.
>
> I copy here the code that does the peak finding using TSpectrum. It
> reads the original x,y table (which is attached as a comment to the
> end of the code), converts it to a histogram (I must confess I am
> still struggling to understand the reason why everything in root uses
> histograms instead of arrays of values, but that's for another SOS
> message) then uses TSpectrum::Search to find the peaks. Histogram,
> x-y table and peaks are shown on the canvas. The red triangles show
> the found peaks which are off the peak by 1 unit.
>
> Any suggestions will be appreciated,
>
> Antonio
>
> ***************************************************************
> Antônio Kanaan
> Departamento de Física - Universidade Federal de Santa Catarina
> CP 476 -- CEP 88040-900 -- Florianópolis -- SC -- Brasil
> http://www.astro.ufsc.br/~kanaan e-mail: kanaan@astro.ufsc.br
> Phone: 55,48,3319069 FAX : 55,48,3319946
> ***************************************************************
>
> #include "TCanvas.h"
> #include "TGraph.h"
>
> #include <iostream>
> #include <fstream>
> using namespace std;
>
> const Int_t MAXNPT = 140; /* number of points */
> Float_t x[MAXNPT],y[MAXNPT],z ; /* data arrays */
>
> search()
> {
> Int_t n,npts;
>
> TCanvas *c1 = new TCanvas("c1","the fit canvas",800,600);
>
> // fills up the x and y reading file "
> ifstream infile("./peak.1") ;
> n = 0 ;
> while ( (infile >> x[n] >> y[n++]) )
> npts = n;
>
> TH1F *hpx = new TH1F("hpx","Peaks",npts,0,npts);
> for (Int_t j=0 ; j<npts ; j++) hpx->SetBinContent(j,y[j]) ;
> hpx->Draw() ;
>
> TSpectrum *sp1 = new TSpectrum();
> sp1->Search(hpx,1.0,"") ;
> Int_t npeaks = sp1->GetNPeaks() ;
> cout << "Peaks found " << npeaks << endl ;
>
> Float_t *peaks ;
> peaks = sp1->GetPositionX() ;
> for (j=0 ; j<npeaks ; j++)
> cout << "Peak at = " << peaks[j] << endl ;
>
> // creates a TGraph object to display the computed curves
> TGraph *gr = new TGraph (npts, x, y);
> gr->Draw("C*") ;
>
> return(0) ;
>
> }
>
> /*
>
> File peak.1 to be read by the above code:
>
> 0 10
> 1 10
> 2 10
> 3 10
> 4 10
> 5 20
> 6 30
> 7 50
> 8 100
> 9 50
> 10 30
> 11 20
> 12 10
> 13 10
> 14 10
> 15 10
> 16 10
> 17 10
> 18 20
> 19 30
> 20 50
> 21 100
> 22 50
> 23 30
> 24 20
> 25 10
> 26 10
> 27 10
> 28 10
> 29 10
> 30 10
>
> */
This archive was generated by hypermail 2b29 : Tue Jan 01 2002 - 17:50:47 MET