Logo ROOT   6.14/05
Reference Guide
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
pict1_FitAwmi.C.png
Processing /mnt/build/workspace/root-makedoc-v614/rootspi/rdoc/src/v6-14-00-patches/tutorials/spectrum/FitAwmi.C...
created -9.76 9.97356 2
created -9.28 34.9074 7
created -8.8 49.8678 10
created -8.32 29.9207 6
created -7.84 29.9207 6
created -7.36 14.9603 3
created -6.88 49.8678 10
created -6.4 4.98678 1
created -5.92 34.9074 7
created -5.44 49.8678 10
created -4.96 19.9471 4
created -4.48 14.9603 3
created -4 39.8942 8
created -3.52 24.9339 5
created -3.04 39.8942 8
created -2.56 44.881 9
created -2.08 34.9074 7
created -1.6 44.881 9
created -1.12 49.8678 10
created -0.64 39.8942 8
created -0.16 39.8942 8
created 0.32 34.9074 7
created 0.8 34.9074 7
created 1.28 49.8678 10
created 1.76 19.9471 4
created 2.24 34.9074 7
created 2.72 49.8678 10
created 3.2 9.97356 2
created 3.68 4.98678 1
created 4.16 9.97356 2
created 4.64 14.9603 3
created 5.12 9.97356 2
created 5.6 39.8942 8
created 6.08 34.9074 7
created 6.56 49.8678 10
created 7.04 49.8678 10
created 7.52 34.9074 7
created 8 44.881 9
created 8.48 14.9603 3
created 8.96 9.97356 2
created 9.44 49.8678 10
the total number of created peaks = 41 with sigma = 0.08
the total number of found peaks = 41 with sigma = 0.0800011 (+-2.92618e-05)
fit chi^2 = 5.62566e-06
found -8.8 (+-0.000261155) 49.8678 (+-0.160445) 10.0001 (+-0.0010625)
found -6.88 (+-0.00025935) 49.8673 (+-0.160316) 10 (+-0.00106164)
found -5.44 (+-0.000260823) 49.8677 (+-0.160421) 10.0001 (+-0.00106234)
found -1.12 (+-0.000261668) 49.868 (+-0.160484) 10.0002 (+-0.00106276)
found 1.28 (+-0.000260823) 49.8677 (+-0.160421) 10.0001 (+-0.00106234)
found 2.72 (+-0.000260364) 49.8676 (+-0.160389) 10.0001 (+-0.00106212)
found 6.56 (+-0.00026165) 49.868 (+-0.160483) 10.0002 (+-0.00106275)
found 7.04 (+-0.00026165) 49.868 (+-0.160483) 10.0002 (+-0.00106275)
found 9.44 (+-0.000258193) 49.8677 (+-0.160249) 10.0001 (+-0.0010612)
found -2.56 (+-0.00027578) 44.8812 (+-0.152246) 9.00015 (+-0.0010082)
found -1.6 (+-0.000276029) 44.8813 (+-0.152264) 9.00017 (+-0.00100832)
found 8 (+-0.000274897) 44.8809 (+-0.152187) 9.0001 (+-0.00100781)
found -4 (+-0.00029144) 39.8941 (+-0.143475) 8.00008 (+-0.000950115)
found -3.04 (+-0.00029256) 39.8944 (+-0.143543) 8.00014 (+-0.000950565)
found -0.64 (+-0.000293209) 39.8946 (+-0.143583) 8.00018 (+-0.000950831)
found -0.16 (+-0.000292777) 39.8945 (+-0.143555) 8.00015 (+-0.00095065)
found 5.6 (+-0.000291487) 39.8941 (+-0.143479) 8.00009 (+-0.000950145)
found -9.28 (+-0.000312358) 34.9076 (+-0.134254) 7.00012 (+-0.000889051)
found -5.92 (+-0.000311874) 34.9075 (+-0.134231) 7.00011 (+-0.0008889)
found -2.08 (+-0.000313829) 34.9079 (+-0.13433) 7.00018 (+-0.000889559)
found 0.32 (+-0.000313332) 34.9077 (+-0.134303) 7.00015 (+-0.000889376)
found 0.800001 (+-0.000313642) 34.9078 (+-0.13432) 7.00017 (+-0.000889491)
found 2.24 (+-0.000312994) 34.9077 (+-0.134285) 7.00014 (+-0.000889263)
found 6.08 (+-0.000313817) 34.9079 (+-0.13433) 7.00018 (+-0.000889555)
found 7.52 (+-0.000313979) 34.908 (+-0.134339) 7.00019 (+-0.000889615)
found -8.32 (+-0.00033902) 29.9211 (+-0.124368) 6.00016 (+-0.000823589)
found -7.84 (+-0.000337406) 29.9207 (+-0.124293) 6.00009 (+-0.000823088)
found -3.52 (+-0.00037208) 24.9344 (+-0.11356) 5.00016 (+-0.000752016)
found -4.96 (+-0.000415623) 19.9475 (+-0.101562) 4.00013 (+-0.000672559)
found 1.76 (+-0.000417146) 19.9477 (+-0.10161) 4.00017 (+-0.000672877)
found 8.48 (+-0.000480215) 14.9607 (+-0.0879634) 3.00011 (+-0.000582509)
found -7.36 (+-0.000482821) 14.961 (+-0.0880258) 3.00016 (+-0.000582923)
found -4.48 (+-0.000481231) 14.9608 (+-0.0879859) 3.00012 (+-0.000582659)
found 4.64 (+-0.000476892) 14.9603 (+-0.0878815) 3.00004 (+-0.000581967)
found 3.19999 (+-0.000589486) 9.97399 (+-0.0718482) 2.00011 (+-0.000475792)
found 5.12 (+-0.000591164) 9.97399 (+-0.071871) 2.00011 (+-0.000475943)
found 8.96 (+-0.000592051) 9.9741 (+-0.0718871) 2.00013 (+-0.00047605)
found -9.76 (+-0.000588168) 9.97379 (+-0.0718179) 2.00007 (+-0.000475591)
found 4.16 (+-0.000585436) 9.97363 (+-0.071777) 2.00004 (+-0.00047532)
found -6.4 (+-0.000849942) 4.98758 (+-0.0509449) 1.00017 (+-0.000337366)
found 3.68 (+-0.000833855) 4.98692 (+-0.050801) 1.00004 (+-0.000336414)
#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;
}
Author

Definition in file FitAwmi.C.