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