FitAwmi

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

Author:
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Sunday, October 02, 2022 at 09:38 AM.

Definition of a helper function:

In [1]:
%%cpp -d

#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;
}

Arguments are defined.

In [2]:
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
created -9.88 19.9471 2
created -9.64 49.8678 5
created -9.4 39.8942 4
created -9.16 79.7885 8
created -8.92 99.7356 10
created -8.68 59.8413 6
created -8.44 19.9471 2
created -8.2 39.8942 4
created -7.96 79.7885 8
created -7.72 79.7885 8
created -7.48 59.8413 6
created -7.24 59.8413 6
created -7 29.9207 3
created -6.76 49.8678 5
created -6.52 29.9207 3
created -6.28 29.9207 3
created -6.04 69.8149 7
created -5.8 29.9207 3
created -5.56 69.8149 7
created -5.32 29.9207 3
created -5.08 69.8149 7
created -4.84 89.762 9
created -4.6 99.7356 10
created -4.36 29.9207 3
created -4.12 69.8149 7
created -3.88 19.9471 2
created -3.64 59.8413 6
created -3.4 19.9471 2
created -3.16 29.9207 3
created -2.92 9.97356 1
created -2.68 99.7356 10
created -2.44 19.9471 2
created -2.2 59.8413 6
created -1.96 29.9207 3
created -1.72 99.7356 10
created -1.48 29.9207 3
created -1.24 49.8678 5
created -1 89.762 9
created -0.76 9.97356 1
created -0.52 19.9471 2
created -0.28 69.8149 7
created -0.04 49.8678 5
created 0.2 29.9207 3
created 0.44 39.8942 4
created 0.68 19.9471 2
created 0.92 9.97356 1
created 1.16 69.8149 7
created 1.4 79.7885 8
created 1.64 79.7885 8
created 1.88 89.762 9
created 2.12 99.7356 10
created 2.36 19.9471 2
created 2.6 69.8149 7
created 2.84 19.9471 2
created 3.08 59.8413 6
created 3.32 9.97356 1
created 3.56 49.8678 5
created 3.8 19.9471 2
created 4.04 29.9207 3
created 4.28 99.7356 10
created 4.52 89.762 9
created 4.76 69.8149 7
created 5 79.7885 8
created 5.24 59.8413 6
created 5.48 29.9207 3
created 5.72 69.8149 7
created 5.96 99.7356 10
created 6.2 59.8413 6
created 6.44 79.7885 8
created 6.68 49.8678 5
created 6.92 49.8678 5
created 7.16 29.9207 3
created 7.4 39.8942 4
created 7.64 9.97356 1
created 7.88 29.9207 3
created 8.12 49.8678 5
created 8.36 69.8149 7
created 8.6 99.7356 10
created 8.84 69.8149 7
created 9.08 29.9207 3
created 9.32 39.8942 4
created 9.56 79.7885 8
created 9.8 79.7885 8
the total number of created peaks = 83 with sigma = 0.04

searching for candidate peaks positions

In [3]:
nfound = s->SearchHighRes(source, dest, nbins, 2., 2., kFALSE, 10000, kFALSE, 0);

filling in the initial estimates of the input parameters

In [4]:
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);

In [5]:
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

In [6]:
sigma /= TMath::Sqrt2(); sigmaErr /= TMath::Sqrt2();

convert "bin numbers" into "x-axis values"

In [7]:
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);
the total number of found peaks = 83 with sigma = 0.040002 (+-1.67053e-05)
fit chi^2 = 1.36343e-05
found -8.92 (+-0.000193828) 99.7342 (+-0.476136) 10.0004 (+-0.00165417)
found -4.6 (+-0.000193524) 99.7337 (+-0.476049) 10.0003 (+-0.00165387)
found 2.12 (+-0.000193334) 99.7334 (+-0.475998) 10.0003 (+-0.00165369)
found 5.96 (+-0.000193737) 99.7339 (+-0.476107) 10.0003 (+-0.00165407)
found 8.6 (+-0.000193836) 99.7342 (+-0.476138) 10.0004 (+-0.00165418)
found -2.68 (+-0.000192228) 99.7312 (+-0.475679) 10.0001 (+-0.00165258)
found -1.72 (+-0.000192865) 99.732 (+-0.475851) 10.0002 (+-0.00165318)
found 4.28 (+-0.000193524) 99.7336 (+-0.476049) 10.0003 (+-0.00165387)
found -4.84 (+-0.000204755) 89.7619 (+-0.451827) 9.00045 (+-0.00156972)
found 1.88 (+-0.000204856) 89.7622 (+-0.451855) 9.00048 (+-0.00156982)
found 4.52 (+-0.000204755) 89.7619 (+-0.451827) 9.00045 (+-0.00156972)
found -1 (+-0.000203213) 89.759 (+-0.451419) 9.00015 (+-0.0015683)
found -7.72 (+-0.000217054) 79.7881 (+-0.425955) 8.00037 (+-0.00147984)
found 1.4 (+-0.000217176) 79.7884 (+-0.425986) 8.0004 (+-0.00147994)
found 1.64 (+-0.000217393) 79.7889 (+-0.426041) 8.00045 (+-0.00148013)
found 5 (+-0.000216942) 79.7878 (+-0.425927) 8.00034 (+-0.00147974)
found 6.44 (+-0.000216685) 79.7873 (+-0.425863) 8.00029 (+-0.00147952)
found 9.8 (+-0.000215179) 79.7886 (+-0.425574) 8.00042 (+-0.00147851)
found -9.16 (+-0.000216966) 79.7881 (+-0.425938) 8.00037 (+-0.00147977)
found -7.96 (+-0.000216765) 79.7876 (+-0.425885) 8.00032 (+-0.00147959)
found 9.56 (+-0.000216765) 79.7876 (+-0.425885) 8.00032 (+-0.00147959)
found 4.76 (+-0.000232663) 69.8159 (+-0.398584) 7.00045 (+-0.00138474)
found 8.84 (+-0.000231978) 69.8148 (+-0.398438) 7.00035 (+-0.00138424)
found -6.04 (+-0.000230965) 69.8129 (+-0.398215) 7.00016 (+-0.00138346)
found -5.56 (+-0.000230965) 69.8129 (+-0.398215) 7.00016 (+-0.00138346)
found -5.08 (+-0.000231868) 69.8145 (+-0.398413) 7.00032 (+-0.00138415)
found -4.12 (+-0.000230712) 69.8127 (+-0.398165) 7.00013 (+-0.00138329)
found -0.279999 (+-0.000231083) 69.8132 (+-0.398243) 7.00018 (+-0.00138356)
found 1.16 (+-0.000231132) 69.8137 (+-0.398267) 7.00024 (+-0.00138364)
found 2.6 (+-0.00023046) 69.8124 (+-0.398114) 7.00011 (+-0.00138311)
found 5.72 (+-0.000231978) 69.8148 (+-0.398438) 7.00035 (+-0.00138424)
found 8.36 (+-0.000232355) 69.8153 (+-0.398517) 7.0004 (+-0.00138451)
found -8.68 (+-0.000250572) 59.8415 (+-0.368889) 6.00032 (+-0.00128158)
found 6.2 (+-0.000251779) 59.8431 (+-0.36911) 6.00048 (+-0.00128235)
found -7.48 (+-0.000251208) 59.842 (+-0.368999) 6.00037 (+-0.00128196)
found -7.24 (+-0.000250298) 59.8407 (+-0.368828) 6.00024 (+-0.00128137)
found 5.24 (+-0.000250602) 59.8412 (+-0.368887) 6.00029 (+-0.00128157)
found -3.64 (+-0.000249124) 59.8393 (+-0.368616) 6.00011 (+-0.00128063)
found -2.2 (+-0.000249411) 59.8396 (+-0.368666) 6.00014 (+-0.0012808)
found 3.08 (+-0.000248715) 59.8391 (+-0.368548) 6.00008 (+-0.0012804)
found -6.76 (+-0.000273844) 49.8668 (+-0.336637) 5.00016 (+-0.00116953)
found -1.24 (+-0.000275071) 49.8684 (+-0.336835) 5.00032 (+-0.00117022)
found -0.0400015 (+-0.000274739) 49.8679 (+-0.33678) 5.00027 (+-0.00117003)
found 6.68 (+-0.000275419) 49.8687 (+-0.336887) 5.00035 (+-0.0011704)
found 6.92 (+-0.000274346) 49.8674 (+-0.336716) 5.00021 (+-0.0011698)
found 8.12 (+-0.000274739) 49.8679 (+-0.33678) 5.00027 (+-0.00117003)
found -9.64 (+-0.000273777) 49.8668 (+-0.336629) 5.00016 (+-0.0011695)
found 3.56 (+-0.000272697) 49.866 (+-0.336471) 5.00008 (+-0.00116895)
found -9.4 (+-0.000308564) 39.8957 (+-0.301404) 4.00035 (+-0.00104713)
found 0.439999 (+-0.000306222) 39.8935 (+-0.301105) 4.00013 (+-0.00104609)
found 7.4 (+-0.000305645) 39.8933 (+-0.301039) 4.00011 (+-0.00104586)
found 9.32 (+-0.00030794) 39.8951 (+-0.301325) 4.00029 (+-0.00104685)
found -8.2 (+-0.000307527) 39.8949 (+-0.301276) 4.00027 (+-0.00104668)
found -4.36 (+-0.000358477) 29.9237 (+-0.261244) 3.00045 (+-0.000907602)
found -1.48 (+-0.000357837) 29.9231 (+-0.26118) 3.0004 (+-0.00090738)
found -7 (+-0.000356777) 29.9221 (+-0.261071) 3.00029 (+-0.000907)
found -6.52 (+-0.000355635) 29.9213 (+-0.260958) 3.00021 (+-0.00090661)
found -5.8 (+-0.000357713) 29.9229 (+-0.261165) 3.00037 (+-0.000907328)
found -5.32 (+-0.000357713) 29.9229 (+-0.261165) 3.00037 (+-0.000907328)
found -1.96 (+-0.000358173) 29.9234 (+-0.261213) 3.00043 (+-0.000907496)
found 0.199999 (+-0.000356069) 29.9215 (+-0.261) 3.00024 (+-0.000906756)
found 5.48 (+-0.000357411) 29.9226 (+-0.261134) 3.00035 (+-0.000907222)
found 7.16 (+-0.000356069) 29.9215 (+-0.261) 3.00024 (+-0.000906756)
found 9.08 (+-0.0003567) 29.9221 (+-0.261064) 3.0003 (+-0.000906978)
found -6.28 (+-0.000356263) 29.9218 (+-0.261022) 3.00027 (+-0.000906832)
found -3.16 (+-0.000353044) 29.92 (+-0.260716) 3.00008 (+-0.000905767)
found 4.04001 (+-0.000356481) 29.9224 (+-0.261051) 3.00032 (+-0.000906934)
found 7.88 (+-0.000354362) 29.9208 (+-0.260844) 3.00016 (+-0.000906213)
found -2.44 (+-0.00044089) 19.9504 (+-0.213434) 2.00043 (+-0.000741504)
found 2.36 (+-0.00044132) 19.9506 (+-0.213464) 2.00045 (+-0.000741607)
found -8.44 (+-0.000438383) 19.9488 (+-0.21326) 2.00027 (+-0.000740898)
found -3.88 (+-0.000439831) 19.9496 (+-0.213359) 2.00035 (+-0.000741242)
found -3.4 (+-0.000437748) 19.9485 (+-0.213218) 2.00024 (+-0.000740752)
found 2.84 (+-0.000439831) 19.9496 (+-0.213359) 2.00035 (+-0.000741242)
found 3.8 (+-0.000437276) 19.9482 (+-0.213185) 2.00021 (+-0.000740639)
found 0.679997 (+-0.000434909) 19.9474 (+-0.213034) 2.00014 (+-0.000740115)
found -9.88 (+-0.000435317) 19.9475 (+-0.213045) 2.00014 (+-0.000740152)
found -0.519994 (+-0.000436323) 19.9483 (+-0.213133) 2.00022 (+-0.000740458)
found -0.760013 (+-0.000624797) 9.97601 (+-0.150976) 1.0003 (+-0.000524516)
found 3.32 (+-0.000626459) 9.97598 (+-0.151025) 1.00029 (+-0.000524686)
found 7.64 (+-0.000622559) 9.97492 (+-0.150887) 1.00019 (+-0.000524204)
found -2.91999 (+-0.000626808) 9.97653 (+-0.151047) 1.00035 (+-0.00052476)
found 0.92001 (+-0.000623581) 9.97547 (+-0.150929) 1.00024 (+-0.000524351)

cleanup

In [8]:
delete pfit;
delete [] Amp;
delete [] FixAmp;
delete [] FixPos;
delete s;
delete [] dest;
delete [] source;
return;

Draw all canvases

In [9]:
gROOT->GetListOfCanvases()->Draw()