Logo ROOT  

Reference Guide
 
Loading...
Searching...
No Matches
FitAwmi.C File Reference
Tutorials » Spectrum tutorials

Detailed Description

View in nbviewer Open in SWAN
This macro fits the source spectrum using the AWMI algorithm from the "TSpectrumFit" class ("TSpectrum" class is used to find peaks).

To try this macro, in a ROOT (5 or 6) prompt, do:

root > .x FitAwmi.C

or:

root > .x FitAwmi.C++
root > FitAwmi(); // re-run with another random set of peaks
created -9.76 39.8942 8
created -9.28 44.881 9
created -8.8 29.9207 6
created -8.32 39.8942 8
created -7.84 19.9471 4
created -7.36 44.881 9
created -6.88 24.9339 5
created -6.4 49.8678 10
created -5.92 9.97356 2
created -5.44 9.97356 2
created -4.96 9.97356 2
created -4.48 39.8942 8
created -4 9.97356 2
created -3.52 14.9603 3
created -3.04 19.9471 4
created -2.56 44.881 9
created -2.08 34.9074 7
created -1.6 4.98678 1
created -1.12 14.9603 3
created -0.64 34.9074 7
created -0.16 24.9339 5
created 0.32 24.9339 5
created 0.8 44.881 9
created 1.28 44.881 9
created 1.76 34.9074 7
created 2.24 39.8942 8
created 2.72 49.8678 10
created 3.2 19.9471 4
created 3.68 39.8942 8
created 4.16 19.9471 4
created 4.64 9.97356 2
created 5.12 34.9074 7
created 5.6 49.8678 10
created 6.08 44.881 9
created 6.56 44.881 9
created 7.04 24.9339 5
created 7.52 19.9471 4
created 8 9.97356 2
created 8.48 14.9603 3
created 8.96 9.97356 2
created 9.44 24.9339 5
the total number of created peaks = 41 with sigma = 0.08
the total number of found peaks = 41 with sigma = 0.0800011 (+-2.19075e-05)
fit chi^2 = 2.83639e-06
found -6.4 (+-0.000184667) 49.8675 (+-0.11387) 10.0001 (+-0.000754069)
found 2.72 (+-0.000185291) 49.8677 (+-0.113916) 10.0001 (+-0.00075437)
found 5.6 (+-0.00018571) 49.8679 (+-0.113947) 10.0002 (+-0.000754577)
found -9.28 (+-0.000195712) 44.8811 (+-0.108097) 9.00014 (+-0.000715835)
found -7.36 (+-0.000195126) 44.8809 (+-0.108057) 9.00009 (+-0.000715571)
found -2.56 (+-0.000195354) 44.881 (+-0.108072) 9.00011 (+-0.000715675)
found 0.800001 (+-0.000195682) 44.8811 (+-0.108095) 9.00014 (+-0.000715823)
found 1.28 (+-0.000195913) 44.8812 (+-0.10811) 9.00016 (+-0.000715927)
found 6.08 (+-0.00019619) 44.8814 (+-0.10813) 9.00019 (+-0.000716056)
found 6.56 (+-0.000195682) 44.8811 (+-0.108095) 9.00014 (+-0.000715823)
found -9.76 (+-0.000207601) 39.8942 (+-0.101908) 8.00009 (+-0.000674851)
found -8.32 (+-0.000207252) 39.8942 (+-0.101894) 8.0001 (+-0.000674763)
found -4.48 (+-0.000206177) 39.8939 (+-0.101832) 8.00004 (+-0.000674348)
found 2.24 (+-0.000208086) 39.8946 (+-0.101946) 8.00017 (+-0.000675103)
found 3.68 (+-0.000206967) 39.8941 (+-0.101877) 8.00008 (+-0.000674649)
found -2.08 (+-0.000221345) 34.9075 (+-0.0953061) 7.0001 (+-0.000631134)
found -0.64 (+-0.000221418) 34.9074 (+-0.0953059) 7.00008 (+-0.000631133)
found 1.76 (+-0.000222723) 34.9078 (+-0.0953763) 7.00017 (+-0.000631599)
found 5.12 (+-0.000221794) 34.9076 (+-0.0953283) 7.00012 (+-0.000631281)
found -8.8 (+-0.000240901) 29.9211 (+-0.0883174) 6.00017 (+-0.000584854)
found -6.88 (+-0.000264648) 24.9345 (+-0.0806532) 5.00019 (+-0.0005341)
found -0.16 (+-0.000263484) 24.9342 (+-0.080606) 5.00012 (+-0.000533788)
found 7.04 (+-0.000263569) 24.9342 (+-0.0806101) 5.00013 (+-0.000533815)
found 0.320001 (+-0.000263801) 24.9343 (+-0.080619) 5.00014 (+-0.000533874)
found 9.44 (+-0.000259722) 24.9339 (+-0.0804756) 5.00007 (+-0.000532924)
found -7.84 (+-0.000296229) 19.9477 (+-0.0721499) 4.00017 (+-0.00047779)
found 3.2 (+-0.000296401) 19.9478 (+-0.0721558) 4.00018 (+-0.000477828)
found 4.16 (+-0.000294356) 19.9474 (+-0.0720908) 4.0001 (+-0.000477398)
found -3.04 (+-0.000294948) 19.9475 (+-0.0721092) 4.00012 (+-0.00047752)
found 7.52 (+-0.000293699) 19.9472 (+-0.0720691) 4.00007 (+-0.000477255)
found -3.52 (+-0.000339558) 14.9604 (+-0.0624233) 3.00006 (+-0.000413378)
found -1.12 (+-0.000339762) 14.9605 (+-0.0624308) 3.00008 (+-0.000413428)
found 8.48 (+-0.000338622) 14.9603 (+-0.0624013) 3.00004 (+-0.000413232)
found -5.92001 (+-0.000419635) 9.97404 (+-0.0510324) 2.00012 (+-0.000337946)
found -4 (+-0.000419763) 9.97399 (+-0.0510328) 2.00011 (+-0.000337949)
found 4.64 (+-0.000420009) 9.97399 (+-0.0510363) 2.00011 (+-0.000337972)
found 8 (+-0.000418074) 9.97378 (+-0.0510035) 2.00007 (+-0.000337754)
found 8.96 (+-0.000418576) 9.97383 (+-0.051012) 2.00008 (+-0.000337811)
found -5.44 (+-0.000416) 9.97363 (+-0.0509701) 2.00004 (+-0.000337533)
found -4.96 (+-0.000419009) 9.97393 (+-0.051021) 2.0001 (+-0.00033787)
found -1.60001 (+-0.000598224) 4.98723 (+-0.0361259) 1.0001 (+-0.000239232)
#include "TROOT.h"
#include "TMath.h"
#include "TRandom.h"
#include "TH1.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TSpectrum.h"
#include "TSpectrumFit.h"
#include "TPolyMarker.h"
#include "TList.h"
#include <iostream>
TH1F *FitAwmi_Create_Spectrum(void) {
Int_t nbins = 1000;
Double_t xmin = -10., xmax = 10.;
delete gROOT->FindObject("h"); // prevent "memory leak"
TH1F *h = new TH1F("h", "simulated spectrum", nbins, xmin, xmax);
h->SetStats(kFALSE);
TF1 f("f", "TMath::Gaus(x, [0], [1], 1)", xmin, xmax);
// f.SetParNames("mean", "sigma");
gRandom->SetSeed(0); // make it really random
// create well separated peaks with exactly known means and areas
// note: TSpectrumFit assumes that all peaks have the same sigma
Double_t sigma = (xmax - xmin) / Double_t(nbins) * Int_t(gRandom->Uniform(2., 6.));
Int_t npeaks = 0;
while (xmax > (xmin + 6. * sigma)) {
npeaks++;
xmin += 3. * sigma; // "mean"
f.SetParameters(xmin, sigma);
Double_t area = 1. * Int_t(gRandom->Uniform(1., 11.));
h->Add(&f, area, ""); // "" ... or ... "I"
std::cout << "created "
<< xmin << " "
<< (area / sigma / TMath::Sqrt(TMath::TwoPi())) << " "
<< area << std::endl;
xmin += 3. * sigma;
}
std::cout << "the total number of created peaks = " << npeaks
<< " with sigma = " << sigma << std::endl;
return h;
}
void FitAwmi(void) {
TH1F *h = FitAwmi_Create_Spectrum();
TCanvas *cFit = ((TCanvas *)(gROOT->GetListOfCanvases()->FindObject("cFit")));
if (!cFit) cFit = new TCanvas("cFit", "cFit", 10, 10, 1000, 700);
else cFit->Clear();
h->Draw("L");
Int_t i, nfound, bin;
Int_t nbins = h->GetNbinsX();
Double_t *source = new Double_t[nbins];
Double_t *dest = new Double_t[nbins];
for (i = 0; i < nbins; i++) source[i] = h->GetBinContent(i + 1);
TSpectrum *s = new TSpectrum(); // note: default maxpositions = 100
// searching for candidate peaks positions
nfound = s->SearchHighRes(source, dest, nbins, 2., 2., kFALSE, 10000, kFALSE, 0);
// filling in the initial estimates of the input parameters
Bool_t *FixPos = new Bool_t[nfound];
Bool_t *FixAmp = new Bool_t[nfound];
for(i = 0; i < nfound; i++) FixAmp[i] = FixPos[i] = kFALSE;
Double_t *Pos, *Amp = new Double_t[nfound]; // ROOT 6
Pos = s->GetPositionX(); // 0 ... (nbins - 1)
for (i = 0; i < nfound; i++) {
bin = 1 + Int_t(Pos[i] + 0.5); // the "nearest" bin
Amp[i] = h->GetBinContent(bin);
}
TSpectrumFit *pfit = new TSpectrumFit(nfound);
pfit->SetFitParameters(0, (nbins - 1), 1000, 0.1, pfit->kFitOptimChiCounts,
pfit->kFitAlphaHalving, pfit->kFitPower2,
pfit->kFitTaylorOrderFirst);
pfit->SetPeakParameters(2., kFALSE, Pos, FixPos, Amp, FixAmp);
// pfit->SetBackgroundParameters(source[0], kFALSE, 0., kFALSE, 0., kFALSE);
pfit->FitAwmi(source);
Double_t *Positions = pfit->GetPositions();
Double_t *PositionsErrors = pfit->GetPositionsErrors();
Double_t *Amplitudes = pfit->GetAmplitudes();
Double_t *AmplitudesErrors = pfit->GetAmplitudesErrors();
Double_t *Areas = pfit->GetAreas();
Double_t *AreasErrors = pfit->GetAreasErrors();
delete gROOT->FindObject("d"); // prevent "memory leak"
TH1F *d = new TH1F(*h); d->SetNameTitle("d", ""); d->Reset("M");
for (i = 0; i < nbins; i++) d->SetBinContent(i + 1, source[i]);
Double_t x1 = d->GetBinCenter(1), dx = d->GetBinWidth(1);
Double_t sigma, sigmaErr;
pfit->GetSigma(sigma, sigmaErr);
// current TSpectrumFit needs a sqrt(2) correction factor for sigma
sigma /= TMath::Sqrt2(); sigmaErr /= TMath::Sqrt2();
// convert "bin numbers" into "x-axis values"
sigma *= dx; sigmaErr *= dx;
std::cout << "the total number of found peaks = " << nfound
<< " with sigma = " << sigma << " (+-" << sigmaErr << ")"
<< std::endl;
std::cout << "fit chi^2 = " << pfit->GetChi() << std::endl;
for (i = 0; i < nfound; i++) {
bin = 1 + Int_t(Positions[i] + 0.5); // the "nearest" bin
Pos[i] = d->GetBinCenter(bin);
Amp[i] = d->GetBinContent(bin);
// convert "bin numbers" into "x-axis values"
Positions[i] = x1 + Positions[i] * dx;
PositionsErrors[i] *= dx;
Areas[i] *= dx;
AreasErrors[i] *= dx;
std::cout << "found "
<< Positions[i] << " (+-" << PositionsErrors[i] << ") "
<< Amplitudes[i] << " (+-" << AmplitudesErrors[i] << ") "
<< Areas[i] << " (+-" << AreasErrors[i] << ")"
<< std::endl;
}
d->SetLineColor(kRed); d->SetLineWidth(1);
d->Draw("SAME L");
TPolyMarker *pm = ((TPolyMarker*)(h->GetListOfFunctions()->FindObject("TPolyMarker")));
if (pm) {
h->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new TPolyMarker(nfound, Pos, Amp);
h->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1);
// cleanup
delete pfit;
delete [] Amp;
delete [] FixAmp;
delete [] FixPos;
delete s;
delete [] dest;
delete [] source;
return;
}
d
#define d(i)
Definition RSha256.hxx:102
f
#define f(i)
Definition RSha256.hxx:104
h
#define h(i)
Definition RSha256.hxx:106
Bool_t
bool Bool_t
Definition RtypesCore.h:63
Int_t
int Int_t
Definition RtypesCore.h:45
kFALSE
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
Double_t
double Double_t
Definition RtypesCore.h:59
kRed
@ kRed
Definition Rtypes.h:66
TRangeDynCast
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Definition TCollection.h:358
dest
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t dest
Definition TGWin32VirtualXProxy.cxx:164
x1
Option_t Option_t TPoint TPoint const char x1
Definition TGWin32VirtualXProxy.cxx:70
xmin
float xmin
Definition THbookFile.cxx:95
xmax
float xmax
Definition THbookFile.cxx:95
gROOT
#define gROOT
Definition TROOT.h:414
gRandom
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
TCanvas
The Canvas class.
Definition TCanvas.h:23
TF1
1-Dim function class
Definition TF1.h:233
TH1F
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:621
TPolyMarker
A PolyMarker is defined by an array on N points in a 2-D space.
Definition TPolyMarker.h:31
TRandom::SetSeed
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:615
TRandom::Uniform
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition TRandom.cxx:682
TSpectrumFit
Advanced 1-dimensional spectra fitting functions.
Definition TSpectrumFit.h:18
TSpectrum
Advanced Spectra Processing.
Definition TSpectrum.h:18
TSpectrum::SearchHighRes
Int_t SearchHighRes(Double_t *source, Double_t *destVector, Int_t ssize, Double_t sigma, Double_t threshold, bool backgroundRemove, Int_t deconIterations, bool markov, Int_t averWindow)
One-dimensional high-resolution peak search function.
Definition TSpectrum.cxx:2129
TSpectrum::GetPositionX
Double_t * GetPositionX() const
Definition TSpectrum.h:58
sigma
const Double_t sigma
Definition h1analysisProxy.h:11
TMath::Sqrt2
constexpr Double_t Sqrt2()
Definition TMath.h:86
TMath::Sqrt
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
TMath::TwoPi
constexpr Double_t TwoPi()
Definition TMath.h:44
Author

Definition in file FitAwmi.C.