// @(#)root/hist:$Name: $:$Id: TSpectrum2.cxx,v 1.12 2006/02/03 21:55:39 pcanal Exp $
// Author: Miroslav Morhac 17/01/2006
/////////////////////////////////////////////////////////////////////////////
// THIS CLASS CONTAINS ADVANCED SPECTRA PROCESSING FUNCTIONS. //
// //
// ONE-DIMENSIONAL BACKGROUND ESTIMATION FUNCTIONS //
// TWO-DIMENSIONAL BACKGROUND ESTIMATION FUNCTIONS //
// ONE-DIMENSIONAL SMOOTHING FUNCTIONS //
// TWO-DIMENSIONAL SMOOTHING FUNCTIONS //
// ONE-DIMENSIONAL DECONVOLUTION FUNCTIONS //
// TWO-DIMENSIONAL DECONVOLUTION FUNCTIONS //
// ONE-DIMENSIONAL PEAK SEARCH FUNCTIONS //
// TWO-DIMENSIONAL PEAK SEARCH FUNCTIONS //
// //
// These functions were written by: //
// Miroslav Morhac //
// Institute of Physics //
// Slovak Academy of Sciences //
// Dubravska cesta 9, 842 28 BRATISLAVA //
// SLOVAKIA //
// //
// email:fyzimiro@savba.sk, fax:+421 7 54772479 //
// //
// The original code in C has been repackaged as a C++ class by R.Brun // //
// //
// The algorithms in this class have been published in the following //
// references: //
// [1] M.Morhac et al.: Background elimination methods for //
// multidimensional coincidence gamma-ray spectra. Nuclear //
// Instruments and Methods in Physics Research A 401 (1997) 113- //
// 132. //
// //
// [2] M.Morhac et al.: Efficient one- and two-dimensional Gold //
// deconvolution and its application to gamma-ray spectra //
// decomposition. Nuclear Instruments and Methods in Physics //
// Research A 401 (1997) 385-408. //
// //
// [3] M.Morhac et al.: Identification of peaks in multidimensional //
// coincidence gamma-ray spectra. Nuclear Instruments and Methods in //
// Research Physics A 443(2000), 108-125. //
// //
// These NIM papers are also available as Postscript files from: //
//
/*
ftp://root.cern.ch/root/SpectrumDec.ps.gz
ftp://root.cern.ch/root/SpectrumSrc.ps.gz
ftp://root.cern.ch/root/SpectrumBck.ps.gz
*/
//
/////////////////////////////////////////////////////////////////////////////
/////////////////////NEW FUNCTIONS January 2006
//
//______________________________________________________________________________
#include "TSpectrum2.h"
#include "TPolyMarker.h"
#include "TMath.h"
#define PEAK_WINDOW 1024
ClassImp(TSpectrum2)
//______________________________________________________________________________
TSpectrum2::TSpectrum2() :TNamed("Spectrum", "Miroslav Morhac peak finder")
{
Int_t n = 100;
fMaxPeaks = n;
fPosition = new Float_t[n];
fPositionX = new Float_t[n];
fPositionY = new Float_t[n];
fResolution = 1;
fHistogram = 0;
fNPeaks = 0;
}
//______________________________________________________________________________
TSpectrum2::TSpectrum2(Int_t maxpositions, Float_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder")
{
// maxpositions: maximum number of peaks
// resolution: determines resolution of the neighboring peaks
// default value is 1 correspond to 3 sigma distance
// between peaks. Higher values allow higher resolution
// (smaller distance between peaks.
// May be set later through SetResolution.
Int_t n = TMath::Max(maxpositions, 100);
fMaxPeaks = n;
fPosition = new Float_t[n];
fPositionX = new Float_t[n];
fPositionY = new Float_t[n];
fHistogram = 0;
fNPeaks = 0;
SetResolution(resolution);
}
//______________________________________________________________________________
TSpectrum2::~TSpectrum2()
{
delete [] fPosition;
delete [] fPositionX;
delete [] fPositionY;
delete fHistogram;
}
//______________________________________________________________________________
const char *TSpectrum2::Background(const TH1 * h, int number_of_iterations,
Option_t * option)
{
/////////////////////////////////////////////////////////////////////////////
// ONE-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION //
// This function calculates background spectrum from source in h. //
// The result is placed in the vector pointed by spectrum pointer. //
// //
// Function parameters: //
// spectrum: pointer to the vector of source spectrum //
// size: length of spectrum and working space vectors //
// number_of_iterations, for details we refer to manual //
// //
/////////////////////////////////////////////////////////////////////////////
Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn"
, h->GetName(), number_of_iterations, option);
return 0;
}
//______________________________________________________________________________
void TSpectrum2::Print(Option_t *) const
{
// Print the array of positions
printf("\nNumber of positions = %d\n",fNPeaks);
for (Int_t i=0;i<fNPeaks;i++) {
printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]);
}
}
//______________________________________________________________________________
Int_t TSpectrum2::Search(const TH1 * hin, Double_t sigma,
Option_t * option, Double_t threshold)
{
/////////////////////////////////////////////////////////////////////////////
// ONE-DIMENSIONAL PEAK SEARCH FUNCTION //
// This function searches for peaks in source spectrum in hin //
// The number of found peaks and their positions are written into //
// the members fNpeaks and fPositionX. //
// //
// Function parameters: //
// hin: pointer to the histogram of source spectrum //
// sigma: sigma of searched peaks, for details we refer to manual //
// Note that sigma is in number of bins //
// threshold: (default=0.05) peaks with amplitude less than //
// threshold*highest_peak are discarded. //
// //
// if option is not equal to "goff" (goff is the default), then //
// a polymarker object is created and added to the list of functions of //
// the histogram. The histogram is drawn with the specified option and //
// the polymarker object drawn on top of the histogram. //
// The polymarker coordinates correspond to the npeaks peaks found in //
// the histogram. //
// A pointer to the polymarker object can be retrieved later via: //
// TList *functions = hin->GetListOfFunctions(); //
// TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker") //
// //
/////////////////////////////////////////////////////////////////////////////
if (hin == 0)
return 0;
Int_t dimension = hin->GetDimension();
if (dimension != 2) {
Error("Search", "Must be a 2-d histogram");
return 0;
}
Int_t sizex = hin->GetXaxis()->GetNbins();
Int_t sizey = hin->GetYaxis()->GetNbins();
Int_t i, j, binx,biny, npeaks;
Float_t ** source = new float *[sizex];
Float_t ** dest = new float *[sizex];
for (i = 0; i < sizex; i++) {
source[i] = new float[sizey];
dest[i] = new float[sizey];
for (j = 0; j < sizey; j++) {
source[i][j] = (Float_t) hin->GetBinContent(i + 1, j + 1);
}
}
//npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, kTRUE, 3, kTRUE, 10);
//the smoothing option is used for 1-d but not for 2-d histograms
npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, kTRUE, 3, kFALSE, 3);
//The logic in the loop should be improved to use the fact
//that fPositionX,Y give a precise position inside a bin.
//The current algorithm takes the center of the bin only.
for (i = 0; i < npeaks; i++) {
binx = 1 + Int_t(fPositionX[i] + 0.5);
biny = 1 + Int_t(fPositionY[i] + 0.5);
fPositionX[i] = hin->GetXaxis()->GetBinCenter(binx);
fPositionY[i] = hin->GetYaxis()->GetBinCenter(biny);
}
for (i = 0; i < sizex; i++) {
delete [] source[i];
delete [] dest[i];
}
delete [] source;
delete [] dest;
if (strstr(option, "goff"))
return npeaks;
if (!npeaks) return 0;
TPolyMarker * pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker");
if (pm) {
hin->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new TPolyMarker(npeaks, fPositionX, fPositionY);
hin->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1.3);
((TH1*)hin)->Draw(option);
return npeaks;
}
//______________________________________________________________________________
void TSpectrum2::SetResolution(Float_t resolution)
{
// resolution: determines resolution of the neighboring peaks
// default value is 1 correspond to 3 sigma distance
// between peaks. Higher values allow higher resolution
// (smaller distance between peaks.
// May be set later through SetResolution.
if (resolution > 1)
fResolution = resolution;
else
fResolution = 1;
}
//_____________________________________________________________________________
//_____________________________________________________________________________
/////////////////////NEW FUNCTIONS JANUARY 2006
//______________________________________________________________________________
const char *TSpectrum2::Background(float **spectrum,
Int_t ssizex, Int_t ssizey,
Int_t numberIterationsX,
Int_t numberIterationsY,
Int_t direction,
Int_t filterType)
{
/////////////////////////////////////////////////////////////////////////////
// TWO-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION - RECTANGULAR RIDGES //
// This function calculates background spectrum from source spectrum. //
// The result is placed to the array pointed by spectrum pointer. //
// //
// Function parameters: //
// spectrum-pointer to the array of source spectrum //
// ssizex-x length of spectrum //
// ssizey-y length of spectrum //
// numberIterationsX-maximal x width of clipping window //
// numberIterationsY-maximal y width of clipping window //
// for details we refer to manual //
// direction- direction of change of clipping window //
// - possible values=kBackIncreasingWindow //
// kBackDecreasingWindow //
// filterType-determines the algorithm of the filtering //
// -possible values=kBackSuccessiveFiltering //
// kBackOneStepFiltering //
// //
// //
/////////////////////////////////////////////////////////////////////////////
//
//
Background
estimation
Goal:
Separation of useful information (peaks) from useless information (background)
•
method is based on Sensitive
Nonlinear Iterative Peak (SNIP) clipping algorithm [1]
•
there exist two algorithms for the
estimation of new value in the channel “
”
Algorithm
based on Successive Comparisons
It
is an extension of one-dimensional SNIP algorithm to another dimension. For
details we refer to [2].
Algorithm
based on One Step Filtering
New
value in the estimated channel is calculated as
.
where
p = 1, 2, …, number_of_iterations.
Function:
const
char*
TSpectrum2::Background
(float
**spectrum, int ssizex, int
ssizey, int numberIterationsX, int
numberIterationsY, int direction, int
filterType)
This
function calculates background spectrum from the source spectrum. The result
is placed in the matrix pointed by spectrum pointer. One can also switch the
direction of the change of the clipping window and to select one of the two
above given algorithms. On successful completion it returns 0. On error it
returns pointer to the string describing error.
Parameters:
spectrum-pointer
to the matrix of source spectrum
ssizex, ssizey-lengths
of the spectrum matrix
numberIterationsX, numberIterationsYmaximal
widths of clipping
window,
direction-
direction of change of clipping window
- possible
values=kBackIncreasingWindow
kBackDecreasingWindow
filterType-type
of the clipping algorithm,
-possible values=kBack SuccessiveFiltering
kBackOneStepFiltering
References:
[1]
C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the
quantitative analysis of PIXE spectra in geoscience applications. NIM, B34
(1988), 396-402.
[2]
M. Morháč, J. Kliman, V.
Matoušek, M. Veselský, I. Turzo.:
Background elimination methods for multidimensional gamma-ray spectra. NIM,
A401 (1997) 113-132.
int i, x, y, sampling, r1, r2;
float a, b, p1, p2, p3, p4, s1, s2, s3, s4;
if (ssizex <= 0 || ssizey <= 0)
return "Wrong parameters";
if (numberIterationsX < 1 || numberIterationsY < 1)
return "Width of Clipping Window Must Be Positive";
if (ssizex < 2 * numberIterationsX + 1
|| ssizey < 2 * numberIterationsY + 1)
return ("Too Large Clipping Window");
float **working_space = new float *[ssizex];
for (i = 0; i < ssizex; i++)
working_space[i] = new float[ssizey];
sampling =
(int) TMath::Max(numberIterationsX, numberIterationsY);
if (direction == kBackIncreasingWindow) {
if (filterType == kBackSuccessiveFiltering) {
for (i = 1; i <= sampling; i++) {
r1 = (int) TMath::Min(i, numberIterationsX), r2 =
(int) TMath::Min(i, numberIterationsY);
for (y = r2; y < ssizey - r2; y++) {
for (x = r1; x < ssizex - r1; x++) {
a = spectrum[x][y];
p1 = spectrum[x - r1][y - r2];
p2 = spectrum[x - r1][y + r2];
p3 = spectrum[x + r1][y - r2];
p4 = spectrum[x + r1][y + r2];
s1 = spectrum[x][y - r2];
s2 = spectrum[x - r1][y];
s3 = spectrum[x + r1][y];
s4 = spectrum[x][y + r2];
b = (p1 + p2) / 2.0;
if (b > s2)
s2 = b;
b = (p1 + p3) / 2.0;
if (b > s1)
s1 = b;
b = (p2 + p4) / 2.0;
if (b > s4)
s4 = b;
b = (p3 + p4) / 2.0;
if (b > s3)
s3 = b;
s1 = s1 - (p1 + p3) / 2.0;
s2 = s2 - (p1 + p2) / 2.0;
s3 = s3 - (p3 + p4) / 2.0;
s4 = s4 - (p2 + p4) / 2.0;
b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
p3 +
p4) / 4.0;
if (b < a && b > 0)
a = b;
working_space[x][y] = a;
}
}
for (y = r2; y < ssizey - r2; y++) {
for (x = r1; x < ssizex - r1; x++) {
spectrum[x][y] = working_space[x][y];
}
}
}
}
else if (filterType == kBackOneStepFiltering) {
for (i = 1; i <= sampling; i++) {
r1 = (int) TMath::Min(i, numberIterationsX), r2 =
(int) TMath::Min(i, numberIterationsY);
for (y = r2; y < ssizey - r2; y++) {
for (x = r1; x < ssizex - r1; x++) {
a = spectrum[x][y];
b = -(spectrum[x - r1][y - r2] +
spectrum[x - r1][y + r2] + spectrum[x + r1][y -
r2]
+ spectrum[x + r1][y + r2]) / 4 +
(spectrum[x][y - r2] + spectrum[x - r1][y] +
spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
if (b < a && b > 0)
a = b;
working_space[x][y] = a;
}
}
for (y = i; y < ssizey - i; y++) {
for (x = i; x < ssizex - i; x++) {
spectrum[x][y] = working_space[x][y];
}
}
}
}
}
else if (direction == kBackDecreasingWindow) {
if (filterType == kBackSuccessiveFiltering) {
for (i = sampling; i >= 1; i--) {
r1 = (int) TMath::Min(i, numberIterationsX), r2 =
(int) TMath::Min(i, numberIterationsY);
for (y = r2; y < ssizey - r2; y++) {
for (x = r1; x < ssizex - r1; x++) {
a = spectrum[x][y];
p1 = spectrum[x - r1][y - r2];
p2 = spectrum[x - r1][y + r2];
p3 = spectrum[x + r1][y - r2];
p4 = spectrum[x + r1][y + r2];
s1 = spectrum[x][y - r2];
s2 = spectrum[x - r1][y];
s3 = spectrum[x + r1][y];
s4 = spectrum[x][y + r2];
b = (p1 + p2) / 2.0;
if (b > s2)
s2 = b;
b = (p1 + p3) / 2.0;
if (b > s1)
s1 = b;
b = (p2 + p4) / 2.0;
if (b > s4)
s4 = b;
b = (p3 + p4) / 2.0;
if (b > s3)
s3 = b;
s1 = s1 - (p1 + p3) / 2.0;
s2 = s2 - (p1 + p2) / 2.0;
s3 = s3 - (p3 + p4) / 2.0;
s4 = s4 - (p2 + p4) / 2.0;
b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 +
p3 +
p4) / 4.0;
if (b < a && b > 0)
a = b;
working_space[x][y] = a;
}
}
for (y = r2; y < ssizey - r2; y++) {
for (x = r1; x < ssizex - r1; x++) {
spectrum[x][y] = working_space[x][y];
}
}
}
}
else if (filterType == kBackOneStepFiltering) {
for (i = sampling; i >= 1; i--) {
r1 = (int) TMath::Min(i, numberIterationsX), r2 =
(int) TMath::Min(i, numberIterationsY);
for (y = r2; y < ssizey - r2; y++) {
for (x = r1; x < ssizex - r1; x++) {
a = spectrum[x][y];
b = -(spectrum[x - r1][y - r2] +
spectrum[x - r1][y + r2] + spectrum[x + r1][y -
r2]
+ spectrum[x + r1][y + r2]) / 4 +
(spectrum[x][y - r2] + spectrum[x - r1][y] +
spectrum[x + r1][y] + spectrum[x][y + r2]) / 2;
if (b < a && b > 0)
a = b;
working_space[x][y] = a;
}
}
for (y = i; y < ssizey - i; y++) {
for (x = i; x < ssizex - i; x++) {
spectrum[x][y] = working_space[x][y];
}
}
}
}
}
for (i = 0; i < ssizex; i++)
delete[]working_space[i];
delete[]working_space;
return 0;
}
//
Example 1– script Back_gamma64.c
:
Fig.
1 Original two-dimensional gamma-gamma-ray spectrum
Fig.
2 Background estimated from data from Fig. 1 using decreasing clipping window with
widths 4, 4 and algorithm based on successive comparisons. The estimate
includes not only continuously changing background but also one-dimensional
ridges.
Fig.
3 Resulting peaks after subtraction of the estimated background (Fig. 2) from original
two-dimensional gamma-gamma-ray spectrum (Fig. 1).
Script:
// Example to illustrate the background estimator (class
TSpectrum).
// To execute this example, do
// root > .x Back_gamma64.C
#include <TSpectrum>
void Back_gamma64() {
Int_t i, j;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new float[nbinsy];
TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
back=(TH2F*) f->Get("back1;1");
TCanvas *Background = new
TCanvas("Background","Estimation of background with increasing
window",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = back->GetBinContent(i +
1,j + 1);
}
}
s->Background(source,nbinsx,nbinsy,4,4,kBackDecreasingWindow,kBackSuccessiveFiltering);
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
back->SetBinContent(i + 1,j + 1, source[i][j]);
}
back->Draw("SURF");
}
Example 2– script Back_gamma256.c
:
Fig.
4 Original two-dimensional gamma-gamma-ray spectrum 256x256 channels
Fig.
5 Peaks after subtraction of the estimated background (increasing clipping
window with widths 8, 8 and algorithm based on successive comparisons) from original
two-dimensional gamma-gamma-ray spectrum (Fig. 4).
Script:
// Example to illustrate the background estimator (class
TSpectrum).
// To execute this example, do
// root > .x Back_gamma256.C
#include <TSpectrum>
void Back_gamma256() {
Int_t i, j;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new float[nbinsy];
TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
back=(TH2F*) f->Get("back2;1");
TCanvas *Background = new
TCanvas("Background","Estimation of background with increasing
window",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = back->GetBinContent(i +
1,j + 1);
}
}
s->Background(source,nbinsx,nbinsy,8,8,kBackIncreasingWindow,kBackSuccessiveFiltering);
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
back->SetBinContent(i + 1,j + 1, source[i][j]);
}
back->Draw("SURF");
}
Example 3– script Back_synt256.c
:
Fig.
6 Original two-dimensional synthetic spectrum 256x256 channels
Fig.
7 Peaks after subtraction of the estimated background (increasing clipping
window with widths 8, 8 and algorithm based on successive comparisons) from original
two-dimensional gamma-gamma-ray spectrum (Fig. 6). One can observe artifacts
(ridges) around the peaks due to the employed algorithm.
Fig.
8 Peaks after subtraction of the estimated background (increasing clipping
window with widths 8, 8 and algorithm based on one step filtering) from original
two-dimensional gamma-gamma-ray spectrum (Fig. 6). The artifacts from the
above given Fig. 7 disappeared.
Script:
// Example to illustrate the background estimator (class
TSpectrum).
// To execute this example, do
// root > .x Back_synt256.C
#include <TSpectrum>
void Back_synt256() {
Int_t i, j;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new float[nbinsy];
TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
back=(TH2F*) f->Get("back3;1");
TCanvas *Background = new
TCanvas("Background","Estimation of background with increasing
window",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = back->GetBinContent(i +
1,j + 1);
}
}
s->Background(source,nbinsx,nbinsy,8,8,kBackIncreasingWindow,kBackSuccessiveFiltering);//kBackOneStepFiltering
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
back->SetBinContent(i + 1,j + 1, source[i][j]);
}
back->Draw("SURF");
}
//_____________________________________________________________________________
const char* TSpectrum2::SmoothMarkov(float **source, Int_t ssizex, Int_t ssizey, Int_t averWindow)
{
/////////////////////////////////////////////////////////////////////////////
// TWO-DIMENSIONAL MARKOV SPECTRUM SMOOTHING FUNCTION
//
// This function calculates smoothed spectrum from source spectrum
// based on Markov chain method.
// The result is placed in the array pointed by source pointer.
//
// Function parameters:
// source-pointer to the array of source spectrum
// ssizex-x length of source
// ssizey-y length of source
// averWindow-width of averaging smoothing window
//
/////////////////////////////////////////////////////////////////////////////
//
Smoothing
Goal: Suppression of
statistical fluctuations
•
the algorithm is based on discrete
Markov chain, which has very simple invariant distribution
being defined from the normalization condition 
n
is the length of the smoothed spectrum and
is the probability of the change of
the peak position from channel i to the channel i+1. 
is the normalization constant so
that 
and m is a width of smoothing window. We have extended this
algortihm to two dimensions.
Function:
const char*
TSpectrum2::SmoothMarkov(float
**fSpectrum, int ssizex, int
ssizey, int averWindow)
This
function calculates smoothed spectrum from the source spectrum based on Markov
chain method. The result is placed in the vector pointed by source pointer. On
successful completion it returns 0. On error it returns pointer to the string
describing error.
Parameters:
fSpectrum-pointer
to the matrix of source spectrum
ssizex, ssizey
-lengths of the spectrum matrix
averWindow-width
of averaging smoothing window
Reference:
[1]
Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
(1996), 451.
int xmin, xmax, ymin, ymax, i, j, l;
double a, b, maxch;
double nom, nip, nim, sp, sm, spx, spy, smx, smy, plocha = 0;
if(averWindow <= 0)
return "Averaging Window must be positive";
float **working_space = new float* [ssizex];
for(i = 0; i < ssizex; i++)
working_space[i] = new float[ssizey];
xmin = 0;
xmax = ssizex - 1;
ymin = 0;
ymax = ssizey - 1;
for(i = 0, maxch = 0; i < ssizex; i++){
for(j = 0; j < ssizey; j++){
working_space[i][j] = 0;
if(maxch < source[i][j])
maxch = source[i][j];
plocha += source[i][j];
}
}
if(maxch == 0)
return 0;
nom = 0;
working_space[xmin][ymin] = 1;
for(i = xmin; i < xmax; i++){
nip = source[i][ymin] / maxch;
nim = source[i + 1][ymin] / maxch;
sp = 0,sm = 0;
for(l = 1; l <= averWindow; l++){
if((i + l) > xmax)
a = source[xmax][ymin] / maxch;
else
a = source[i + l][ymin] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a = TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
sp = sp + b;
if(i - l + 1 < xmin)
a = source[xmin][ymin] / maxch;
else
a = source[i - l + 1][ymin] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a = TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
sm = sm + b;
}
a = sp / sm;
a = working_space[i + 1][ymin] = a * working_space[i][ymin];
nom = nom + a;
}
for(i = ymin; i < ymax; i++){
nip = source[xmin][i] / maxch;
nim = source[xmin][i + 1] / maxch;
sp = 0,sm = 0;
for(l = 1; l <= averWindow; l++){
if((i + l) > ymax)
a = source[xmin][ymax] / maxch;
else
a = source[xmin][i + l] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a = TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
sp = sp + b;
if(i - l + 1 < ymin)
a = source[xmin][ymin] / maxch;
else
a = source[xmin][i - l + 1] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a = TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
sm = sm + b;
}
a = sp / sm;
a = working_space[xmin][i + 1] = a * working_space[xmin][i];
nom = nom + a;
}
for(i = xmin; i < xmax; i++){
for(j = ymin; j < ymax; j++){
nip = source[i][j + 1] / maxch;
nim = source[i + 1][j + 1] / maxch;
spx = 0,smx = 0;
for(l = 1; l <= averWindow; l++){
if(i + l > xmax)
a = source[xmax][j] / maxch;
else
a = source[i + l][j] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a = TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
spx = spx + b;
if(i - l + 1 < xmin)
a = source[xmin][j] / maxch;
else
a = source[i - l + 1][j] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a = TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
smx = smx + b;
}
spy = 0,smy = 0;
nip = source[i + 1][j] / maxch;
nim = source[i + 1][j + 1] / maxch;
for(l = 1; l <= averWindow; l++){
if(j + l > ymax)
a = source[i][ymax]/maxch;
else
a = source[i][j + l] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a = TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
spy = spy + b;
if(j - l + 1 < ymin)
a = source[i][ymin] / maxch;
else
a = source[i][j - l + 1] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a = TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
smy = smy + b;
}
a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy);
working_space[i + 1][j + 1] = a;
nom = nom + a;
}
}
for(i = xmin; i <= xmax; i++){
for(j = ymin; j <= ymax; j++){
working_space[i][j] = working_space[i][j] / nom;
}
}
for(i = 0;i < ssizex; i++){
for(j = 0; j < ssizey; j++){
source[i][j] = plocha * working_space[i][j];
}
}
for (i = 0; i < ssizex; i++)
delete[]working_space[i];
delete[]working_space;
return 0;
}
//
Example 4 – script Smooth.c
:
Fig. 9 Original noisy
spectrum. Fig. 10 Smoothed spectrum m=3
Peaks can hardly be
observed. Peaks become apparent.
Fig. 11 Smoothed spectrum
m=5 Fig.12 Smoothed spectrum m=7
Script:
// Example to illustrate the Markov smoothing (class
TSpectrum).
// To execute this example, do
// root > .x Smooth.C
#include <TSpectrum>
void Smooth() {
Int_t i, j;
Double_t nbinsx = 256;
Double_t nbinsy = 256;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
TH2F *smooth = new
TH2F("smooth","Background
estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
smooth=(TH2F*) f->Get("smooth1;1");
TCanvas *Smoothing = new
TCanvas("Smoothing","Markov smoothing",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = smooth->GetBinContent(i +
1,j + 1);
}
}
s->SmoothMarkov(source,nbinsx,nbinsx,3);//5,7
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
smooth->SetBinContent(i + 1,j + 1,
source[i][j]);
}
smooth->Draw("SURF");
}
//______________________________________________________________________________________________________________________________
const char *TSpectrum2::Deconvolution(float **source, float **resp,
Int_t ssizex, Int_t ssizey,
Int_t numberIterations,
Int_t numberRepetitions,
Double_t boost)
{
/////////////////////////////////////////////////////////////////////////////
// TWO-DIMENSIONAL DECONVOLUTION FUNCTION
// This function calculates deconvolution from source spectrum
// according to response spectrum
// The result is placed in the matrix pointed by source pointer.
//
// Function parameters:
// source-pointer to the matrix of source spectrum
// resp-pointer to the matrix of response spectrum
// ssizex-x length of source and response spectra
// ssizey-y length of source and response spectra
// numberIterations, for details we refer to manual
// numberRepetitions, for details we refer to manual
// boost, boosting factor, for details we refer to manual
//
/////////////////////////////////////////////////////////////////////////////
//
Deconvolution
Goal:
Improvement of the resolution in spectra, decomposition of multiplets
Mathematical formulation of
the 2-dimensional convolution system is
where h(i,j) is the impulse
response function, x, y are input and output matrices, respectively, 
are the lengths
of x and h matrices
•
let us assume that we know the
response and the output matrices (spectra) of the above given system.
•
the deconvolution represents
solution of the overdetermined system of linear equations, i.e., the
calculation of the matrix x.
•
from numerical stability point of
view the operation of deconvolution is extremely critical (ill-posed problem)
as well as time consuming operation.
•
the Gold deconvolution algorithm
proves to work very well even for 2-dimensional systems. Generalization of the
algorithm for 2-dimensional systems was presented in [1], [2].
•
for Gold deconvolution algorithm
as well as for boosted deconvolution algorithm we refer also to TSpectrum
Function:
const
char* TSpectrum2::Deconvolution(float **source,
const float
**resp, int ssizex,
int ssizey, int numberIterations,
int numberRepetitions,
double boost)
This
function calculates deconvolution from source spectrum according to response
spectrum using Gold deconvolution algorithm. The result is placed in the matrix
pointed by source pointer. On successful completion it returns 0. On error it
returns pointer to the string describing error. If desired after every
numberIterations one can apply boosting operation (exponential function with
exponent given by boost coefficient) and repeat it numberRepetitions times.
Parameters:
source-pointer
to the matrix of source spectrum
resp-pointer
to the matrix of response spectrum
ssizex, ssizey-lengths
of the spectrum matrix
numberIterations-number of iterations
numberRepetitions-number of repetitions
for boosted deconvolution. It must be
greater or equal to one.
boost-boosting coefficient, applies only
if numberRepetitions is greater than one.
Recommended range <1,2>.
References:
[1]
M. Morháč, J. Kliman, V.
Matoušek, M. Veselský, I. Turzo.:
Efficient one- and two-dimensional Gold deconvolution and its application to
gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.
[2]
Morháč M., Matoušek V., Kliman J., Efficient algorithm of multidimensional
deconvolution and its application to nuclear data processing, Digital Signal
Processing 13 (2003) 144.
int i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max,
i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet;
double lda, ldb, ldc, area, maximum = 0;
if (ssizex <= 0 || ssizey <= 0)
return "Wrong parameters";
if (numberIterations <= 0)
return "Number of iterations must be positive";
if (numberRepetitions <= 0)
return "Number of repetitions must be positive";
double **working_space = new double *[ssizex];
for (i = 0; i < ssizex; i++)
working_space[i] = new double[5 * ssizey];
area = 0;
lhx = -1, lhy = -1;
for (i = 0; i < ssizex; i++) {
for (j = 0; j < ssizey; j++) {
lda = resp[i][j];
if (lda != 0) {
if ((i + 1) > lhx)
lhx = i + 1;
if ((j + 1) > lhy)
lhy = j + 1;
}
working_space[i][j] = lda;
area = area + lda;
if (lda > maximum) {
maximum = lda;
positx = i, posity = j;
}
}
}
if (lhx == -1 || lhy == -1)
return ("Zero response data");
//calculate ht*y and write into p
for (i2 = 0; i2 < ssizey; i2++) {
for (i1 = 0; i1 < ssizex; i1++) {
ldc = 0;
for (j2 = 0; j2 <= (lhy - 1); j2++) {
for (j1 = 0; j1 <= (lhx - 1); j1++) {
k2 = i2 + j2, k1 = i1 + j1;
if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) {
lda = working_space[j1][j2];
ldb = source[k1][k2];
ldc = ldc + lda * ldb;
}
}
}
working_space[i1][i2 + ssizey] = ldc;
}
}
//calculate matrix b=ht*h
i1min = -(lhx - 1), i1max = lhx - 1;
i2min = -(lhy - 1), i2max = lhy - 1;
for (i2 = i2min; i2 <= i2max; i2++) {
for (i1 = i1min; i1 <= i1max; i1++) {
ldc = 0;
j2min = -i2;
if (j2min < 0)
j2min = 0;
j2max = lhy - 1 - i2;
if (j2max > lhy - 1)
j2max = lhy - 1;
for (j2 = j2min; j2 <= j2max; j2++) {
j1min = -i1;
if (j1min < 0)
j1min = 0;
j1max = lhx - 1 - i1;
if (j1max > lhx - 1)
j1max = lhx - 1;
for (j1 = j1min; j1 <= j1max; j1++) {
lda = working_space[j1][j2];
if (i1 + j1 < ssizex && i2 + j2 < ssizey)
ldb = working_space[i1 + j1][i2 + j2];
else
ldb = 0;
ldc = ldc + lda * ldb;
}
}
working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc;
}
}
//initialization in x1 matrix
for (i2 = 0; i2 < ssizey; i2++) {
for (i1 = 0; i1 < ssizex; i1++) {
working_space[i1][i2 + 3 * ssizey] = 1;
working_space[i1][i2 + 4 * ssizey] = 0;
}
}
//START OF ITERATIONS
for (repet = 0; repet < numberRepetitions; repet++) {
if (repet != 0) {
for (i = 0; i < ssizex; i++) {
for (j = 0; j < ssizey; j++) {
working_space[i][j + 3 * ssizey] =
TMath::Power(working_space[i][j + 3 * ssizey], boost);
}
}
}
for (lindex = 0; lindex < numberIterations; lindex++) {
for (i2 = 0; i2 < ssizey; i2++) {
for (i1 = 0; i1 < ssizex; i1++) {
ldb = 0;
j2min = i2;
if (j2min > lhy - 1)
j2min = lhy - 1;
j2min = -j2min;
j2max = ssizey - i2 - 1;
if (j2max > lhy - 1)
j2max = lhy - 1;
j1min = i1;
if (j1min > lhx - 1)
j1min = lhx - 1;
j1min = -j1min;
j1max = ssizex - i1 - 1;
if (j1max > lhx - 1)
j1max = lhx - 1;
for (j2 = j2min; j2 <= j2max; j2++) {
for (j1 = j1min; j1 <= j1max; j1++) {
ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey];
lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey];
ldb = ldb + lda * ldc;
}
}
lda = working_space[i1][i2 + 3 * ssizey];
ldc = working_space[i1][i2 + 1 * ssizey];
if (ldc * lda != 0 && ldb != 0) {
lda = lda * ldc / ldb;
}
else
lda = 0;
working_space[i1][i2 + 4 * ssizey] = lda;
}
}
for (i2 = 0; i2 < ssizey; i2++) {
for (i1 = 0; i1 < ssizex; i1++)
working_space[i1][i2 + 3 * ssizey] =
working_space[i1][i2 + 4 * ssizey];
}
}
}
for (i = 0; i < ssizex; i++) {
for (j = 0; j < ssizey; j++)
source[(i + positx) % ssizex][(j + posity) % ssizey] =
area * working_space[i][j + 3 * ssizey];
}
for (i = 0; i < ssizex; i++)
delete[]working_space[i];
delete[]working_space;
return 0;
}
//
Example 5 – script Decon.c
:
•
response function (usually peak)
should be shifted to the beginning of the coordinate system (see Fig. 13)
Fig.
13 2-dimensional response spectrum
Fig.
14 2-dimensional gamma-gamma-ray input spectrum (before deconvolution)
Fig.
15 Spectrum from Fig. 14 after deconvolution (1000 iterations)
Script:
// Example to illustrate the Gold deconvolution (class
TSpectrum2).
// To execute this example, do
// root > .x Decon.C
#include <TSpectrum2>
void Decon() {
Int_t i, j;
Double_t nbinsx = 256;
Double_t nbinsy = 256;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
TH2F *decon = new TH2F("decon","Gold
deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
decon=(TH2F*) f->Get("decon1;1");
Float_t ** response = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
response[i]=new
float[nbinsy];
TH2F *resp = new TH2F("resp","Response
matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
resp=(TH2F*) f->Get("resp1;1");
TCanvas *Deconvol = new
TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = decon->GetBinContent(i +
1,j + 1);
}
}
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
response[i][j] = resp->GetBinContent(i +
1,j + 1);
}
}
s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
decon->SetBinContent(i + 1,j + 1, source[i][j]);
}
decon->Draw("SURF");
}
Example 6 – script
Decon2.c :
Fig.
16 Response spectrum
Fig.
17 Original synthetic input spectrum (before deconvolution). It is composed of
17 peaks. 5 peaks are overlapping in the outlined multiplet and two peaks in
doublet.
Fig.
18 Spectrum from Fig. 17 after deconvolution (1000 iterations). Resolution is
improved but the peaks in multiplet remained unresolved.
Script:
// Example to illustrate the Gold deconvolution (class
TSpectrum2).
// To execute this example, do
// root > .x Decon2.C
#include <TSpectrum2>
void Decon2() {
Int_t i, j;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
TH2F *decon = new TH2F("decon","Gold
deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new TFile("spectra2\\TSpectrum2.root");
decon=(TH2F*) f->Get("decon2;1");
Float_t ** response = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
response[i]=new
float[nbinsy];
TH2F *resp = new TH2F("resp","Response
matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
resp=(TH2F*) f->Get("resp2;1");
TCanvas *Deconvol = new
TCanvas("Deconvolution","Gold
deconvolution",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = decon->GetBinContent(i +
1,j + 1);
}
}
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
response[i][j] = resp->GetBinContent(i +
1,j + 1);
}
}
s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
decon->SetBinContent(i + 1,j + 1, source[i][j]);
}
decon->Draw("SURF");
}
Example 7 – script
Decon2HR.c :
Fig.
19 Spectrum from Fig. 17 after boosted deconvolution (50 iterations repeated 20
times, boosting coefficient was 1.2). All the peaks in multiplet as well as in
doublet are completely decomposed.
Script:
// Example to illustrate boosted Gold deconvolution (class
TSpectrum2).
// To execute this example, do
// root > .x Decon2HR.C
//#include <TSpectrum2>
void Decon2HR() {
Int_t i, j;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
TH2F *decon = new TH2F("decon","Boosted
Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new TFile("spectra2\\TSpectrum2.root");
decon=(TH2F*) f->Get("decon2;1");
Float_t ** response = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
response[i]=new
float[nbinsy];
TH2F *resp = new TH2F("resp","Response
matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
resp=(TH2F*) f->Get("resp2;1");
TCanvas *Deconvol = new
TCanvas("Deconvolution","Gold
deconvolution",10,10,1000,700);
TSpectrum *s = new TSpectrum();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = decon->GetBinContent(i +
1,j + 1);
}
}
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
response[i][j] = resp->GetBinContent(i +
1,j + 1);
}
}
s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1);
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++)
decon->SetBinContent(i + 1,j + 1, source[i][j]);
}
decon->Draw("SURF");
}
//____________________________________________________________________________
Int_t TSpectrum2::SearchHighRes(float **source,float **dest, Int_t ssizex, Int_t ssizey,
Double_t sigma, Double_t threshold,
Bool_t backgroundRemove,Int_t deconIterations,
Bool_t markov, Int_t averWindow)
{
/////////////////////////////////////////////////////////////////////////////
// TWO-DIMENSIONAL HIGH-RESOLUTION PEAK SEARCH FUNCTION //
// This function searches for peaks in source spectrum //
// It is based on deconvolution method. First the background is //
// removed (if desired), then Markov spectrum is calculated //
// (if desired), then the response function is generated //
// according to given sigma and deconvolution is carried out. //
// //
// Function parameters: //
// source-pointer to the matrix of source spectrum //
// dest-pointer to the matrix of resulting deconvolved spectrum //
// ssizex-x length of source spectrum //
// ssizey-y length of source spectrum //
// sigma-sigma of searched peaks, for details we refer to manual //
// threshold-threshold value in % for selected peaks, peaks with //
// amplitude less than threshold*highest_peak/100 //
// are ignored, see manual //
// backgroundRemove-logical variable, set if the removal of //
// background before deconvolution is desired //
// deconIterations-number of iterations in deconvolution operation //
// markov-logical variable, if it is true, first the source spectrum //
// is replaced by new spectrum calculated using Markov //
// chains method. //
// averWindow-averanging window of searched peaks, for details //
// we refer to manual (applies only for Markov method) //
// //
/////////////////////////////////////////////////////////////////////////////
//
Peaks searching
Goal:
to identify automatically the peaks in spectrum with the presence of the
continuous background, one-fold coincidences (ridges) and statistical
fluctuations - noise.
The common problems connected with correct peak
identification in two-dimensional coincidence spectra are
- non-sensitivity to noise, i.e., only statistically
relevant peaks should be identified
- non-sensitivity of the algorithm to continuous
background
- non-sensitivity to one-fold coincidences (coincidences
peak – background in both dimensions) and their crossings
- ability to identify peaks close to the edges of the
spectrum region. Usually peak finders fail to detect them
- resolution, decomposition of doublets and multiplets.
The algorithm should be able to recognize close positioned peaks.
- ability to identify peaks with different sigma
Function:
Int_t TSpectrum2::SearchHighRes (float **source,float **dest, int ssizex, int ssizey, float sigma, double threshold,
bool backgroundRemove,int deconIterations,
bool markov,
int averWindow)
This
function searches for peaks in source spectrum. It is based on deconvolution
method. First the background is removed (if desired), then Markov smoothed
spectrum is calculated (if desired), then the response function is generated
according to given sigma and deconvolution is carried out. The order of peaks
is arranged according to their heights in the spectrum after background
elimination. The highest peak is the first in the list. On success it returns
number of found peaks.
Parameters:
source-pointer to the matrix of source
spectrum
dest-resulting spectrum after deconvolution
ssizex, ssizey-lengths of the source and
destination spectra
sigma-sigma of searched peaks
threshold- threshold
value in % for selected peaks, peaks with amplitude less than
threshold*highest_peak/100 are ignored
backgroundRemove- background_remove-logical variable, true if the removal of
background before deconvolution is desired
deconIterations-number of iterations in
deconvolution operation
markov-logical variable, if it is true,
first the source spectrum is replaced by new spectrum calculated using Markov
chains method
averWindow-width of
averaging smoothing window
References:
[1]
M.A. Mariscotti: A method for identification of peaks in the presence of
background and its application to spectrum analysis. NIM 50 (1967), 309-320.
[2]
M. Morháč, J. Kliman, V.
Matoušek, M. Veselský, I. Turzo.:Identification
of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000)
108-125.
[3]
Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376
(1996), 451.
int number_of_iterations = (int)(4 * sigma + 0.5);
int k, lindex, priz;
double lda, ldb, ldc, area, maximum;
int xmin, xmax, l, peak_index = 0, ssizex_ext = ssizex + 4 * number_of_iterations, ssizey_ext = ssizey + 4 * number_of_iterations, shift = 2 * number_of_iterations;
int ymin, ymax, i, j;
double a, b, ax, ay, maxch, plocha = 0;
double nom, nip, nim, sp, sm, spx, spy, smx, smy;
double p1, p2, p3, p4, s1, s2, s3, s4;
int x, y;
int lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity;
if (sigma < 1) {
Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1");
return 0;
}
if(threshold<=0||threshold>=100){
Error("SearchHighRes", "Invalid threshold, must be positive and less than 100");
return 0;
}
j = (int) (5.0 * sigma + 0.5);
if (j >= PEAK_WINDOW / 2) {
Error("SearchHighRes", "Too large sigma");
return 0;
}
if (markov == true) {
if (averWindow <= 0) {
Error("SearchHighRes", "Averanging window must be positive");
return 0;
}
}
if(backgroundRemove == true){
if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){
Error("SearchHighRes", "Too large clipping window");
return 0;
}
}
i = (int)(4 * sigma + 0.5);
i = 4 * i;
double **working_space = new double *[ssizex + i];
for (j = 0; j < ssizex + i; j++) {
Double_t *wsk = working_space[j] = new double[16 * (ssizey + i)];
for (Int_t k=0;k<16 * (ssizey + i);k++) wsk[k] = 0;
}
for(j = 0; j < ssizey_ext; j++){
for(i = 0; i < ssizex_ext; i++){
if(i < shift){
if(j < shift)
working_space[i][j + ssizey_ext] = source[0][0];
else if(j >= ssizey + shift)
working_space[i][j + ssizey_ext] = source[0][ssizey - 1];
else
working_space[i][j + ssizey_ext] = source[0][j - shift];
}
else if(i >= ssizex + shift){
if(j < shift)
working_space[i][j + ssizey_ext] = source[ssizex - 1][0];
else if(j >= ssizey + shift)
working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1];
else
working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift];
}
else{
if(j < shift)
working_space[i][j + ssizey_ext] = source[i - shift][0];
else if(j >= ssizey + shift)
working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1];
else
working_space[i][j + ssizey_ext] = source[i - shift][j - shift];
}
}
}
if(backgroundRemove == true){
for(i = 1; i <= number_of_iterations; i++){
for(y = i; y < ssizey_ext - i; y++){
for(x = i; x < ssizex_ext - i; x++){
a = working_space[x][y + ssizey_ext];
p1 = working_space[x - i][y + ssizey_ext - i];
p2 = working_space[x - i][y + ssizey_ext + i];
p3 = working_space[x + i][y + ssizey_ext - i];
p4 = working_space[x + i][y + ssizey_ext + i];
s1 = working_space[x][y + ssizey_ext - i];
s2 = working_space[x - i][y + ssizey_ext];
s3 = working_space[x + i][y + ssizey_ext];
s4 = working_space[x][y + ssizey_ext + i];
b = (p1 + p2) / 2.0;
if(b > s2)
s2 = b;
b = (p1 + p3) / 2.0;
if(b > s1)
s1 = b;
b = (p2 + p4) / 2.0;
if(b > s4)
s4 = b;
b = (p3 + p4) / 2.0;
if(b > s3)
s3 = b;
s1 = s1 - (p1 + p3) / 2.0;
s2 = s2 - (p1 + p2) / 2.0;
s3 = s3 - (p3 + p4) / 2.0;
s4 = s4 - (p2 + p4) / 2.0;
b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0;
if(b < a)
a = b;
working_space[x][y] = a;
}
}
for(y = i;y < ssizey_ext - i; y++){
for(x = i; x < ssizex_ext - i; x++){
working_space[x][y + ssizey_ext] = working_space[x][y];
}
}
}
for(j = 0;j < ssizey_ext; j++){
for(i = 0; i < ssizex_ext; i++){
if(i < shift){
if(j < shift)
working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext];
else if(j >= ssizey + shift)
working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext];
else
working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext];
}
else if(i >= ssizex + shift){
if(j < shift)
working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext];
else if(j >= ssizey + shift)
working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext];
else
working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext];
}
else{
if(j < shift)
working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext];
else if(j >= ssizey + shift)
working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext];
else
working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext];
}
}
}
}
for(j = 0; j < ssizey_ext; j++){
for(i = 0; i < ssizex_ext; i++){
working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext];
}
}
if(markov == true){
for(i = 0;i < ssizex_ext; i++){
for(j = 0; j < ssizey_ext; j++)
working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j];
}
xmin = 0;
xmax = ssizex_ext - 1;
ymin = 0;
ymax = ssizey_ext - 1;
for(i = 0, maxch = 0; i < ssizex_ext; i++){
for(j = 0; j < ssizey_ext; j++){
working_space[i][j] = 0;
if(maxch < working_space[i][j + 2 * ssizey_ext])
maxch = working_space[i][j + 2 * ssizey_ext];
plocha += working_space[i][j + 2 * ssizey_ext];
}
}
if(maxch == 0)
return 0;
nom=0;
working_space[xmin][ymin] = 1;
for(i = xmin; i < xmax; i++){
nip = working_space[i][ymin + 2 * ssizey_ext] / maxch;
nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch;
sp = 0,sm = 0;
for(l = 1;l <= averWindow; l++){
if((i + l) > xmax)
a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch;
else
a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a=TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
sp = sp + b;
if(i - l + 1 < xmin)
a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
else
a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a=TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
sm = sm + b;
}
a = sp / sm;
a = working_space[i + 1][ymin] = a * working_space[i][ymin];
nom = nom + a;
}
for(i = ymin; i < ymax; i++){
nip = working_space[xmin][i + 2 * ssizey_ext] / maxch;
nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch;
sp = 0,sm = 0;
for(l = 1; l <= averWindow; l++){
if((i + l) > ymax)
a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch;
else
a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch;
b = a - nip;
if(a + nip <= 0)
a=1;
else
a=TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
sp = sp + b;
if(i - l + 1 < ymin)
a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch;
else
a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch;
b = a - nim;
if(a + nim <= 0)
a = 1;
else
a=TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
sm = sm + b;
}
a = sp / sm;
a = working_space[xmin][i + 1] = a * working_space[xmin][i];
nom = nom + a;
}
for(i = xmin; i < xmax; i++){
for(j = ymin; j < ymax; j++){
nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch;
nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
spx = 0,smx = 0;
for(l = 1; l <= averWindow; l++){
if(i + l > xmax)
a = working_space[xmax][j + 2 * ssizey_ext] / maxch;
else
a = working_space[i + l][j + 2 * ssizey_ext] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a=TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
spx = spx + b;
if(i - l + 1 < xmin)
a = working_space[xmin][j + 2 * ssizey_ext] / maxch;
else
a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch;
b = a - nim;
if(a + nim <= 0)
a=1;
else
a=TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
smx = smx + b;
}
spy = 0,smy = 0;
nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch;
nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch;
for(l = 1; l <= averWindow; l++){
if(j + l > ymax)
a = working_space[i][ymax + 2 * ssizey_ext] / maxch;
else
a = working_space[i][j + l + 2 * ssizey_ext] / maxch;
b = a - nip;
if(a + nip <= 0)
a = 1;
else
a=TMath::Sqrt(a + nip);
b = b / a;
b = TMath::Exp(b);
spy = spy + b;
if(j - l + 1 < ymin)
a = working_space[i][ymin + 2 * ssizey_ext] / maxch;
else
a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch;
b=a-nim;
if(a + nim <= 0)
a = 1;
else
a=TMath::Sqrt(a + nim);
b = b / a;
b = TMath::Exp(b);
smy = smy + b;
}
a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy);
working_space[i + 1][j + 1] = a;
nom = nom + a;
}
}
for(i = xmin; i <= xmax; i++){
for(j = ymin; j <= ymax; j++){
working_space[i][j] = working_space[i][j] / nom;
}
}
for(i = 0; i < ssizex_ext; i++){
for(j = 0; j < ssizey_ext; j++){
working_space[i][j + ssizey_ext] = working_space[i][j] * plocha;
working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j];
}
}
}
//deconvolution starts
area = 0;
lhx = -1,lhy = -1;
positx = 0,posity = 0;
maximum = 0;
//generate response matrix
for(i = 0; i < ssizex_ext; i++){
for(j = 0; j < ssizey_ext; j++){
lda = (double)i - 3 * sigma;
ldb = (double)j - 3 * sigma;
lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma);
k=(int)(1000 * TMath::Exp(-lda));
lda = k;
if(lda != 0){
if((i + 1) > lhx)
lhx = i + 1;
if((j + 1) > lhy)
lhy = j + 1;
}
working_space[i][j] = lda;
area = area + lda;
if(lda > maximum){
maximum = lda;
positx = i,posity = j;
}
}
}
//read source matrix
for(i = 0;i < ssizex_ext; i++){
for(j = 0;j < ssizey_ext; j++){
working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]);
}
}
//calculate matrix b=ht*h
i = lhx - 1;
if(i > ssizex_ext)
i = ssizex_ext;
j = lhy - 1;
if(j>ssizey_ext)
j = ssizey_ext;
i1min = -i,i1max = i;
i2min = -j,i2max = j;
for(i2 = i2min; i2 <= i2max; i2++){
for(i1 = i1min; i1 <= i1max; i1++){
ldc = 0;
j2min = -i2;
if(j2min < 0)
j2min = 0;
j2max = lhy - 1 - i2;
if(j2max > lhy - 1)
j2max = lhy - 1;
for(j2 = j2min; j2 <= j2max; j2++){
j1min = -i1;
if(j1min < 0)
j1min = 0;
j1max = lhx - 1 - i1;
if(j1max > lhx - 1)
j1max = lhx - 1;
for(j1 = j1min; j1 <= j1max; j1++){
lda = working_space[j1][j2];
ldb = working_space[i1 + j1][i2 + j2];
ldc = ldc + lda * ldb;
}
}
k = (i1 + ssizex_ext) / ssizex_ext;
working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc;
}
}
//calculate at*y and write into p
i = lhx - 1;
if(i > ssizex_ext)
i = ssizex_ext;
j = lhy - 1;
if(j > ssizey_ext)
j = ssizey_ext;
i2min = -j,i2max = ssizey_ext + j - 1;
i1min = -i,i1max = ssizex_ext + i - 1;
for(i2 = i2min; i2 <= i2max; i2++){
for(i1=i1min;i1<=i1max;i1++){
ldc=0;
for(j2 = 0; j2 <= (lhy - 1); j2++){
for(j1 = 0; j1 <= (lhx - 1); j1++){
k2 = i2 + j2,k1 = i1 + j1;
if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){
lda = working_space[j1][j2];
ldb = working_space[k1][k2 + 14 * ssizey_ext];
ldc = ldc + lda * ldb;
}
}
}
k = (i1 + ssizex_ext) / ssizex_ext;
working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc;
}
}
//move matrix p
for(i2 = 0; i2 < ssizey_ext; i2++){
for(i1 = 0; i1 < ssizex_ext; i1++){
k = (i1 + ssizex_ext) / ssizex_ext;
ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext];
working_space[i1][i2 + 14 * ssizey_ext] = ldb;
}
}
//initialization in x1 matrix
for(i2 = 0; i2 < ssizey_ext; i2++){
for(i1 = 0; i1 < ssizex_ext; i1++){
working_space[i1][i2 + ssizey_ext] = 1;
working_space[i1][i2 + 2 * ssizey_ext] = 0;
}
}
//START OF ITERATIONS
for(lindex = 0; lindex < deconIterations; lindex++){
for(i2 = 0; i2 < ssizey_ext; i2++){
for(i1 = 0; i1 < ssizex_ext; i1++){
lda = working_space[i1][i2 + ssizey_ext];
ldc = working_space[i1][i2 + 14 * ssizey_ext];
if(lda > 0.000001 && ldc > 0.000001){
ldb=0;
j2min=i2;
if(j2min > lhy - 1)
j2min = lhy - 1;
j2min = -j2min;
j2max = ssizey_ext - i2 - 1;
if(j2max > lhy - 1)
j2max = lhy - 1;
j1min = i1;
if(j1min > lhx - 1)
j1min = lhx - 1;
j1min = -j1min;
j1max = ssizex_ext - i1 - 1;
if(j1max > lhx - 1)
j1max = lhx - 1;
for(j2 = j2min; j2 <= j2max; j2++){
for(j1 = j1min; j1 <= j1max; j1++){
k = (j1 + ssizex_ext) / ssizex_ext;
ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext];
lda = working_space[i1 + j1][i2 + j2 + ssizey_ext];
ldb = ldb + lda * ldc;
}
}
lda = working_space[i1][i2 + ssizey_ext];
ldc = working_space[i1][i2 + 14 * ssizey_ext];
if(ldc * lda != 0 && ldb != 0){
lda =lda * ldc / ldb;
}
else
lda=0;
working_space[i1][i2 + 2 * ssizey_ext] = lda;
}
}
}
for(i2 = 0; i2 < ssizey_ext; i2++){
for(i1 = 0; i1 < ssizex_ext; i1++)
working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext];
}
}
//looking for maximum
maximum=0;
for(i = 0; i < ssizex_ext; i++){
for(j = 0; j < ssizey_ext; j++){
working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext];
if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext])
maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext];
}
}
//searching for peaks in deconvolved spectrum
for(i = 1; i < ssizex_ext - 1; i++){
for(j = 1; j < ssizey_ext - 1; j++){
if(working_space[i][j] > working_space[i - 1][j] && working_space[i][j] > working_space[i - 1][j - 1] && working_space[i][j] > working_space[i][j - 1] && working_space[i][j] > working_space[i + 1][j - 1] && working_space[i][j] > working_space[i + 1][j] && working_space[i][j] > working_space[i + 1][j + 1] && working_space[i][j] > working_space[i][j + 1] && working_space[i][j] > working_space[i - 1][j + 1]){
if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){
if(working_space[i][j] > threshold * maximum / 100.0){
if(peak_index < fMaxPeaks){
for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){
a += (double)(k - shift) * working_space[k][j];
b += working_space[k][j];
}
ax=a/b;
if(ax < 0)
ax = 0;
if(ax >= ssizex)
ax = ssizex - 1;
for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){
a += (double)(k - shift) * working_space[i][k];
b += working_space[i][k];
}
ay=a/b;
if(ay < 0)
ay = 0;
if(ay >= ssizey)
ay = ssizey - 1;
if(peak_index == 0){
fPositionX[0] = ax;
fPositionY[0] = ay;
peak_index = 1;
}
else{
for(k = 0, priz = 0; k < peak_index && priz == 0; k++){
if(working_space[shift+(int)(ax+0.5)][15 * ssizey_ext + shift + (int)(ay+0.5)] > working_space[shift+(int)(fPositionX[k]+0.5)][15 * ssizey_ext + shift + (int)(fPositionY[k]+0.5)])
priz = 1;
}
if(priz == 0){
if(k < fMaxPeaks){
fPositionX[k] = ax;
fPositionY[k] = ay;
}
}
else{
for(l = peak_index; l >= k; l--){
if(l < fMaxPeaks){
fPositionX[l] = fPositionX[l - 1];
fPositionY[l] = fPositionY[l - 1];
}
}
fPositionX[k - 1] = ax;
fPositionY[k - 1] = ay;
}
if(peak_index < fMaxPeaks)
peak_index += 1;
}
}
}
}
}
}
}
//writing back deconvolved spectrum
for(i = 0; i < ssizex; i++){
for(j = 0; j < ssizey; j++){
dest[i][j] = working_space[i + shift][j + shift];
}
}
i = (int)(4 * sigma + 0.5);
i = 4 * i;
for (j = 0; j < ssizex + i; j++)
delete[]working_space[j];
delete[]working_space;
fNPeaks = peak_index;
return fNPeaks;
}
//
Examples of peak searching
method
SearchHighRes function provides users with the possibility
to vary the input parameters and with the access to the output deconvolved data
in the destination spectrum. Based on the output data one can tune the
parameters.
Example 8 – script Src.c:
Fig.
20 Two-dimensional spectrum with found peaks denoted by markers (
,
threshold=5%, 3 iterations steps in the deconvolution)
Fig.
21 Spectrum from Fig. 20 after background elimination and deconvolution
Script:
// Example to illustrate high resolution peak searching
function (class TSpectrum).
// To execute this example, do
// root > .x Src.C
#include <TSpectrum2>
void Src() {
Int_t i, j, nfound;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
Float_t ** dest = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
dest[i]=new
float[nbinsy];
TH2F *search = new TH2F("search","High
resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new TFile("spectra2\\TSpectrum2.root");
search=(TH2F*) f->Get("search4;1");
TCanvas *Searching = new
TCanvas("Searching","High resolution peak
searching",10,10,1000,700);
TSpectrum2 *s = new TSpectrum2();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = search->GetBinContent(i +
1,j + 1);
}
}
nfound = s->SearchHighRes(source, dest, nbinsx,
nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);
printf("Found %d candidate peaks\n",nfound);
for(i=0;i<nfound;i++)
printf("posx= %d, posy= %d, value=
%d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5),
(int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]);
}
Example 9 – script Src2.c:
Fig.
22 Two-dimensional noisy spectrum with found peaks denoted by markers (,
threshold=10%, 10 iterations steps in the deconvolution). One can observe that
the algorithm is insensitive to the crossings of one-dimensional ridges. It
identifies only two-cooincidence peaks.
Fig.
23 Spectrum from Fig. 22 after background elimination and deconvolution
Script:
// Example to illustrate high resolution peak searching
function (class TSpectrum).
// To execute this example, do
// root > .x Src2.C
#include <TSpectrum2>
void Src2() {
Int_t i, j, nfound;
Double_t nbinsx = 256;
Double_t nbinsy = 256;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
Float_t ** dest = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
dest[i]=new
float[nbinsy];
TH2F *search = new TH2F("search","High
resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
search=(TH2F*) f->Get("back3;1");
TCanvas *Searching = new
TCanvas("Searching","High resolution peak
searching",10,10,1000,700);
TSpectrum2 *s = new TSpectrum2();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = search->GetBinContent(i +
1,j + 1);
}
}
nfound = s->SearchHighRes(source, dest, nbinsx,
nbinsy, 2, 10, kTRUE, 10, kFALSE, 3);
printf("Found %d candidate peaks\n",nfound);
for(i=0;i<nfound;i++)
printf("posx= %d, posy= %d, value=
%d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5),
(int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]);
}
Example 10 – script Src3.c:
Fig.
24 Two-dimensional spectrum with 15 found peaks denoted by markers. Some peaks
are positioned close to each other. It is necessary to increase number of
iterations in the deconvolution. In next 3 Figs. we shall study the influence
of this parameter.
Fig.
25 Spectrum from Fig. 24 after deconvolution (# of iterations = 3). Number of
identified peaks = 13.
Fig.
26 Spectrum from Fig. 24 after deconvolution (# of iterations = 10). Number of
identified peaks = 13.
Fig.
27 Spectrum from Fig. 24 after deconvolution (# of iterations = 100). Number of
identified peaks = 15. Now the algorithm is able to decompose two doublets in
the spectrum.
Script:
// Example to illustrate high resolution peak searching
function (class TSpectrum).
// To execute this example, do
// root > .x Src3.C
#include <TSpectrum2>
void Src3() {
Int_t i, j, nfound;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
Float_t ** dest = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
dest[i]=new
float[nbinsy];
TH2F *search = new TH2F("search","High
resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
search=(TH2F*) f->Get("search1;1");
TCanvas *Searching = new
TCanvas("Searching","High resolution peak
searching",10,10,1000,700);
TSpectrum2 *s = new TSpectrum2();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = search->GetBinContent(i +
1,j + 1);
}
}
nfound = s->SearchHighRes(source, dest, nbinsx,
nbinsy, 2, 2, kFALSE, 3, kFALSE, 1);//3, 10, 100
printf("Found %d candidate peaks\n",nfound);
for(i=0;i<nfound;i++)
printf("posx= %d, posy= %d, value=
%d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5),
(int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]);
}
Example 11 – script Src4.c:
Fig.
28 Two-dimensional spectrum with peaks with different sigma denoted by markers (
,
threshold=5%, 10 iterations steps in the deconvolution, Markov smoothing with
window=3)
Fig.
29 Spectrum from Fig. 28 after smoothing and deconvolution.
Script:
// Example to illustrate high resolution peak searching
function (class TSpectrum).
// To execute this example, do
// root > .x Src4.C
#include <TSpectrum2>
void Src4() {
Int_t i, j, nfound;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
Float_t ** dest = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
dest[i]=new
float[nbinsy];
TH2F *search = new TH2F("search","High
resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
search=(TH2F*) f->Get("search2;1");
TCanvas *Searching = new
TCanvas("Searching","High resolution peak
searching",10,10,1000,700);
TSpectrum2 *s = new TSpectrum2();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = search->GetBinContent(i +
1,j + 1);
}
}
nfound = s->SearchHighRes(source, dest, nbinsx,
nbinsy, 3, 5, kFALSE, 10, kTRUE, 3);
printf("Found %d candidate peaks\n",nfound);
for(i=0;i<nfound;i++)
printf("posx= %d, posy= %d, value=
%d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5),
(int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]);
}
Example 12 – script Src5.c:
Fig.
30 Two-dimensional spectrum with peaks positioned close to the edges denoted by
markers (, threshold=5%, 10 iterations
steps in the deconvolution)
Fig.
31 Spectrum from Fig. 30 after deconvolution.
Script:
// Example to illustrate high resolution peak searching
function (class TSpectrum).
// To execute this example, do
// root > .x Src5.C
#include <TSpectrum2>
void Src5() {
Int_t i, j, nfound;
Double_t nbinsx = 64;
Double_t nbinsy = 64;
Double_t xmin = 0;
Double_t xmax = (Double_t)nbinsx;
Double_t ymin = 0;
Double_t ymax = (Double_t)nbinsy;
Float_t ** source = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
source[i]=new
float[nbinsy];
Float_t ** dest = new float *[nbinsx];
for (i=0;i<nbinsx;i++)
dest[i]=new
float[nbinsy];
TH2F *search = new TH2F("search","High
resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
TFile *f = new
TFile("spectra2\\TSpectrum2.root");
search=(TH2F*) f->Get("search3;1");
TCanvas *Searching = new
TCanvas("Searching","High resolution peak
searching",10,10,1000,700);
TSpectrum2 *s = new TSpectrum2();
for (i = 0; i < nbinsx; i++){
for (j = 0; j < nbinsy; j++){
source[i][j] = search->GetBinContent(i +
1,j + 1);
}
}
nfound = s->SearchHighRes(source, dest, nbinsx,
nbinsy, 2, 5, kFALSE, 10, kFALSE, 1);
printf("Found %d candidate peaks\n",nfound);
for(i=0;i<nfound;i++)
printf("posx= %d, posy= %d, value=
%d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5),
(int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]);
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.
