Logo ROOT  
Reference Guide
FitAwmi.C File Reference

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.88 89.762 9
created -9.64 69.8149 7
created -9.4 59.8413 6
created -9.16 9.97356 1
created -8.92 59.8413 6
created -8.68 9.97356 1
created -8.44 9.97356 1
created -8.2 69.8149 7
created -7.96 89.762 9
created -7.72 49.8678 5
created -7.48 29.9207 3
created -7.24 99.7356 10
created -7 49.8678 5
created -6.76 59.8413 6
created -6.52 19.9471 2
created -6.28 69.8149 7
created -6.04 39.8942 4
created -5.8 19.9471 2
created -5.56 59.8413 6
created -5.32 49.8678 5
created -5.08 89.762 9
created -4.84 69.8149 7
created -4.6 39.8942 4
created -4.36 59.8413 6
created -4.12 79.7885 8
created -3.88 89.762 9
created -3.64 69.8149 7
created -3.4 69.8149 7
created -3.16 29.9207 3
created -2.92 59.8413 6
created -2.68 19.9471 2
created -2.44 99.7356 10
created -2.2 29.9207 3
created -1.96 19.9471 2
created -1.72 49.8678 5
created -1.48 99.7356 10
created -1.24 49.8678 5
created -1 49.8678 5
created -0.76 99.7356 10
created -0.52 79.7885 8
created -0.28 89.762 9
created -0.04 59.8413 6
created 0.2 29.9207 3
created 0.44 29.9207 3
created 0.68 29.9207 3
created 0.92 79.7885 8
created 1.16 89.762 9
created 1.4 69.8149 7
created 1.64 59.8413 6
created 1.88 39.8942 4
created 2.12 79.7885 8
created 2.36 89.762 9
created 2.6 19.9471 2
created 2.84 49.8678 5
created 3.08 69.8149 7
created 3.32 69.8149 7
created 3.56 99.7356 10
created 3.8 59.8413 6
created 4.04 29.9207 3
created 4.28 9.97356 1
created 4.52 79.7885 8
created 4.76 99.7356 10
created 5 79.7885 8
created 5.24 49.8678 5
created 5.48 59.8413 6
created 5.72 89.762 9
created 5.96 69.8149 7
created 6.2 39.8942 4
created 6.44 79.7885 8
created 6.68 79.7885 8
created 6.92 79.7885 8
created 7.16 29.9207 3
created 7.4 89.762 9
created 7.64 59.8413 6
created 7.88 9.97356 1
created 8.12 69.8149 7
created 8.36 89.762 9
created 8.6 29.9207 3
created 8.84 19.9471 2
created 9.08 49.8678 5
created 9.32 39.8942 4
created 9.56 99.7356 10
created 9.8 59.8413 6
the total number of created peaks = 83 with sigma = 0.04
the total number of found peaks = 83 with sigma = 0.040002 (+-1.39171e-05)
fit chi^2 = 1.04592e-05
found 3.56 (+-0.000169685) 99.7339 (+-0.417001) 10.0003 (+-0.00144873)
found 4.76 (+-0.000169933) 99.7347 (+-0.417078) 10.0004 (+-0.00144899)
found -7.24 (+-0.000169162) 99.7326 (+-0.416847) 10.0002 (+-0.00144819)
found -2.44 (+-0.000168758) 99.7318 (+-0.416732) 10.0001 (+-0.00144779)
found -1.48 (+-0.000169404) 99.7331 (+-0.416916) 10.0003 (+-0.00144843)
found -0.759999 (+-0.000169668) 99.7339 (+-0.416997) 10.0003 (+-0.00144871)
found 9.56 (+-0.00016939) 99.7331 (+-0.416913) 10.0003 (+-0.00144842)
found -7.96 (+-0.000178886) 89.7606 (+-0.395609) 9.00032 (+-0.00137441)
found -3.88 (+-0.000179177) 89.7614 (+-0.39569) 9.0004 (+-0.00137469)
found -0.28 (+-0.000179081) 89.7611 (+-0.395663) 9.00037 (+-0.0013746)
found 1.16 (+-0.000179177) 89.7614 (+-0.39569) 9.0004 (+-0.00137469)
found 2.36 (+-0.000178528) 89.7601 (+-0.395519) 9.00026 (+-0.0013741)
found 5.72 (+-0.000178993) 89.7609 (+-0.395638) 9.00034 (+-0.00137451)
found 8.36 (+-0.000178621) 89.7601 (+-0.395539) 9.00026 (+-0.00137417)
found -9.88 (+-0.000178777) 89.7594 (+-0.395545) 9.00019 (+-0.00137419)
found -5.08 (+-0.000178886) 89.7606 (+-0.395609) 9.00032 (+-0.00137441)
found 7.4 (+-0.000178526) 89.7598 (+-0.395513) 9.00024 (+-0.00137407)
found -4.12 (+-0.000190199) 79.7884 (+-0.373099) 8.0004 (+-0.00129621)
found -0.52 (+-0.000190582) 79.7894 (+-0.373196) 8.0005 (+-0.00129654)
found 5 (+-0.000190166) 79.7884 (+-0.373092) 8.0004 (+-0.00129618)
found 6.68 (+-0.000190313) 79.7886 (+-0.373127) 8.00042 (+-0.0012963)
found 6.92 (+-0.000189696) 79.7873 (+-0.372977) 8.00029 (+-0.00129578)
found 0.920001 (+-0.000189787) 79.7876 (+-0.373001) 8.00032 (+-0.00129586)
found 2.12 (+-0.000189946) 79.7878 (+-0.373038) 8.00034 (+-0.00129599)
found 4.52 (+-0.000189387) 79.7873 (+-0.372918) 8.00029 (+-0.00129558)
found 6.44 (+-0.000189855) 79.7876 (+-0.373014) 8.00032 (+-0.00129591)
found -9.64 (+-0.000203547) 69.8153 (+-0.34905) 7.0004 (+-0.00121266)
found -4.84 (+-0.000203261) 69.8148 (+-0.348989) 7.00035 (+-0.00121244)
found -3.64 (+-0.000203668) 69.8156 (+-0.349077) 7.00043 (+-0.00121275)
found -3.4 (+-0.00020287) 69.814 (+-0.348904) 7.00027 (+-0.00121215)
found 1.4 (+-0.000203547) 69.8153 (+-0.34905) 7.0004 (+-0.00121266)
found 3.32 (+-0.000203765) 69.8159 (+-0.349099) 7.00045 (+-0.00121283)
found 5.96 (+-0.000203261) 69.8148 (+-0.348989) 7.00035 (+-0.00121244)
found -8.2 (+-0.00020254) 69.814 (+-0.348848) 7.00026 (+-0.00121195)
found -6.28 (+-0.000202246) 69.8129 (+-0.348772) 7.00016 (+-0.00121169)
found 3.08 (+-0.000203199) 69.8145 (+-0.348973) 7.00032 (+-0.00121239)
found 8.12 (+-0.00020254) 69.814 (+-0.348848) 7.00027 (+-0.00121195)
found -0.0400019 (+-0.000219609) 59.8415 (+-0.323115) 6.00032 (+-0.00112255)
found 3.8 (+-0.00021972) 59.8417 (+-0.323137) 6.00035 (+-0.00112263)
found 7.64 (+-0.00021899) 59.841 (+-0.323012) 6.00027 (+-0.0011222)
found 9.8 (+-0.000218106) 59.8425 (+-0.322901) 6.00042 (+-0.00112181)
found -9.4 (+-0.000218747) 59.8404 (+-0.322964) 6.00021 (+-0.00112203)
found -6.76 (+-0.00021882) 59.8402 (+-0.322967) 6.00019 (+-0.00112204)
found -4.36 (+-0.000219694) 59.8415 (+-0.323128) 6.00032 (+-0.0011226)
found -2.92 (+-0.000218448) 59.8396 (+-0.322898) 6.00013 (+-0.0011218)
found 1.64 (+-0.000219567) 59.8412 (+-0.323103) 6.00029 (+-0.00112251)
found 5.48 (+-0.000219987) 59.842 (+-0.323184) 6.00037 (+-0.00112279)
found -8.92 (+-0.000217482) 59.8388 (+-0.322736) 6.00005 (+-0.00112124)
found -5.56 (+-0.00021882) 59.8402 (+-0.322967) 6.00019 (+-0.00112204)
found -7.72 (+-0.000240922) 49.8685 (+-0.295019) 5.00032 (+-0.00102494)
found -7 (+-0.000241683) 49.8695 (+-0.295139) 5.00043 (+-0.00102536)
found -1.24 (+-0.000241499) 49.8692 (+-0.29511) 5.0004 (+-0.00102526)
found -5.32 (+-0.000241551) 49.8692 (+-0.295117) 5.0004 (+-0.00102528)
found -0.999998 (+-0.000241499) 49.8693 (+-0.29511) 5.0004 (+-0.00102526)
found 5.24 (+-0.00024141) 49.869 (+-0.295093) 5.00037 (+-0.0010252)
found -1.72 (+-0.000240755) 49.8684 (+-0.294997) 5.00032 (+-0.00102487)
found 2.84 (+-0.000240335) 49.8677 (+-0.294927) 5.00024 (+-0.00102462)
found 9.08 (+-0.00023979) 49.8668 (+-0.294839) 5.00016 (+-0.00102432)
found -6.04 (+-0.000269166) 39.8946 (+-0.263849) 4.00024 (+-0.000916654)
found -4.6 (+-0.000270296) 39.8957 (+-0.263991) 4.00035 (+-0.000917146)
found 1.88 (+-0.000270482) 39.8959 (+-0.264016) 4.00037 (+-0.000917233)
found 6.2 (+-0.000270686) 39.8962 (+-0.264043) 4.0004 (+-0.000917326)
found 9.32 (+-0.00027059) 39.8962 (+-0.264032) 4.0004 (+-0.000917289)
found -2.2 (+-0.000312226) 29.9224 (+-0.228643) 3.00032 (+-0.000794344)
found 7.16 (+-0.00031401) 29.9237 (+-0.228815) 3.00045 (+-0.00079494)
found 8.6 (+-0.000312023) 29.9221 (+-0.228622) 3.0003 (+-0.000794269)
found -7.48 (+-0.000313413) 29.9232 (+-0.228756) 3.0004 (+-0.000794734)
found -3.16 (+-0.000313041) 29.9226 (+-0.228716) 3.00035 (+-0.000794596)
found 0.199998 (+-0.000311774) 29.9216 (+-0.228591) 3.00024 (+-0.000794162)
found 4.04 (+-0.000310656) 29.921 (+-0.228491) 3.00019 (+-0.000793813)
found 0.44 (+-0.00031078) 29.9208 (+-0.228493) 3.00016 (+-0.00079382)
found 0.680003 (+-0.000312274) 29.9221 (+-0.228642) 3.00029 (+-0.00079434)
found 2.6 (+-0.000385452) 19.9498 (+-0.186889) 2.00037 (+-0.00064928)
found -6.52 (+-0.000385228) 19.9496 (+-0.186871) 2.00035 (+-0.000649221)
found -2.68 (+-0.000386156) 19.9504 (+-0.186938) 2.00043 (+-0.000649451)
found -5.8 (+-0.000383961) 19.9488 (+-0.186785) 2.00027 (+-0.00064892)
found -1.96 (+-0.00038299) 19.9483 (+-0.186719) 2.00022 (+-0.000648693)
found 8.84 (+-0.000382991) 19.9482 (+-0.186719) 2.00021 (+-0.000648693)
found -9.16 (+-0.000549414) 9.97625 (+-0.132303) 1.00032 (+-0.000459641)
found -8.68001 (+-0.000543769) 9.97494 (+-0.132111) 1.00019 (+-0.000458976)
found 7.88 (+-0.000550052) 9.97652 (+-0.132327) 1.00035 (+-0.000459724)
found 4.28001 (+-0.000548015) 9.976 (+-0.132257) 1.0003 (+-0.00045948)
found -8.43999 (+-0.000544388) 9.97522 (+-0.132135) 1.00022 (+-0.000459058)
#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->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->SetMarkerSize(1);
// cleanup
delete pfit;
delete [] Amp;
delete [] FixAmp;
delete [] FixPos;
delete s;
delete [] dest;
delete [] source;
return;
}
#define d(i)
Definition: RSha256.hxx:102
#define f(i)
Definition: RSha256.hxx:104
#define h(i)
Definition: RSha256.hxx:106
static const double x1[5]
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
@ kRed
Definition: Rtypes.h:64
float xmin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
#define gROOT
Definition: TROOT.h:415
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
The Canvas class.
Definition: TCanvas.h:31
void Clear(Option_t *option="")
Remove all primitives from the canvas.
Definition: TCanvas.cxx:715
1-Dim function class
Definition: TF1.h:211
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
A PolyMarker is defined by an array on N points in a 2-D space.
Definition: TPolyMarker.h:31
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:635
Advanced 1-dimensional spectra fitting functions.
Definition: TSpectrumFit.h:18
Double_t GetChi() const
Definition: TSpectrumFit.h:121
void SetPeakParameters(Double_t sigma, Bool_t fixSigma, const Double_t *positionInit, const Bool_t *fixPosition, const Double_t *ampInit, const Bool_t *fixAmp)
This function sets the following fitting parameters of peaks:
Double_t * GetAmplitudesErrors() const
Definition: TSpectrumFit.h:117
void FitAwmi(Double_t *source)
This function fits the source spectrum.
Double_t * GetAreasErrors() const
Definition: TSpectrumFit.h:119
@ kFitTaylorOrderFirst
Definition: TSpectrumFit.h:82
void GetSigma(Double_t &sigma, Double_t &sigmaErr)
This function gets the sigma parameter and its error.
Double_t * GetAreas() const
Definition: TSpectrumFit.h:118
Double_t * GetAmplitudes() const
Definition: TSpectrumFit.h:116
void SetFitParameters(Int_t xmin, Int_t xmax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
This function sets the following fitting parameters:
Double_t * GetPositionsErrors() const
Definition: TSpectrumFit.h:123
Double_t * GetPositions() const
Definition: TSpectrumFit.h:122
Advanced Spectra Processing.
Definition: TSpectrum.h:18
const Double_t sigma
static constexpr double s
constexpr Double_t Sqrt2()
Definition: TMath.h:89
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
constexpr Double_t TwoPi()
Definition: TMath.h:45
#define dest(otri, vertexptr)
Definition: triangle.c:1040
Author

Definition in file FitAwmi.C.