Skip to content

Commit 6859d10

Browse files
ferdymercurydpiparo
authored andcommitted
[hist] Implement calling TSpectrum2 and TSpectrum3 Background from TH2 or TH3
Since it's implemented for the 1D version, and for the 2D versions you can call it via raw pointers but not directly via TH2 or TH3. Functions were marked before as to-do, so no one could be calling them directly since they were already doing nothing and raising error. Hence, function signature and name were adapted to match the 2D and 3D needed parameters. Since these are legacy classes, the code from the implemented 1D version was copy-pasted and adapted, but no code modernization was performed, such as avoiding new/delete etc. Fixes https://root-forum.cern.ch/t/tspectrum2-bac-error-function-not-yet-implemented/64427 I've been hit myself with this, too.
1 parent d2c87b2 commit 6859d10

File tree

10 files changed

+212
-70
lines changed

10 files changed

+212
-70
lines changed

hist/hist/inc/TH2.h

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -130,9 +130,9 @@ class TH2 : public TH1 {
130130
virtual void SetShowProjectionX(Int_t nbins=1); // *MENU*
131131
virtual void SetShowProjectionY(Int_t nbins=1); // *MENU*
132132
virtual void SetShowProjectionXY(Int_t nbinsY=1, Int_t nbinsX=1); // *MENU*
133-
TH1 *ShowBackground(Int_t niter=20, Option_t *option="same") override;
134-
Int_t ShowPeaks(Double_t sigma=2, Option_t *option="", Double_t threshold=0.05) override; // *MENU*
135-
void Smooth(Int_t ntimes=1, Option_t *option="") override; // *MENU*
133+
virtual TH1 *ShowBackground2D(Int_t nIterX = 20, Int_t nIterY = 20, Option_t *option = "same");
134+
Int_t ShowPeaks(Double_t sigma = 2, Option_t *option = "", Double_t threshold = 0.05) override; // *MENU*
135+
void Smooth(Int_t ntimes = 1, Option_t *option = "") override; // *MENU*
136136

137137
ClassDefOverride(TH2,5) //2-Dim histogram base class
138138
};

hist/hist/inc/TH3.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -144,6 +144,7 @@ class TH3 : public TH1, public TAtt3D {
144144
void SetBinContent(Int_t bin, Int_t, Double_t content) override { SetBinContent(bin, content); }
145145
void SetBinContent(Int_t binx, Int_t biny, Int_t binz, Double_t content) override { SetBinContent(GetBin(binx, biny, binz), content); }
146146
virtual void SetShowProjection(const char *option="xy",Int_t nbins=1); // *MENU*
147+
virtual TH1 *ShowBackground3D(Int_t nIterX = 20, Int_t nIterY = 20, Int_t nIterZ = 20, Option_t *option = "same");
147148

148149
protected:
149150

hist/hist/src/TH2.cxx

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2636,13 +2636,13 @@ void TH2::SetShowProjectionXY(Int_t nbinsY, Int_t nbinsX)
26362636
////////////////////////////////////////////////////////////////////////////////
26372637
/// This function calculates the background spectrum in this histogram.
26382638
/// The background is returned as a histogram.
2639-
/// to be implemented (may be)
26402639

2641-
TH1 *TH2::ShowBackground(Int_t niter, Option_t *option)
2640+
TH1 *TH2::ShowBackground2D(Int_t nIterX, Int_t nIterY, Option_t *option)
26422641
{
26432642

2644-
return (TH1 *)gROOT->ProcessLineFast(TString::Format("TSpectrum2::StaticBackground((TH1*)0x%zx,%d,\"%s\")",
2645-
(size_t)this, niter, option).Data());
2643+
return (TH1 *)gROOT->ProcessLineFast(
2644+
TString::Format("TSpectrum2::StaticBackground((TH1*)0x%zx,%d,%d,\"%s\")", (size_t)this, nIterX, nIterY, option)
2645+
.Data());
26462646
}
26472647

26482648

hist/hist/src/TH3.cxx

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4883,3 +4883,15 @@ TH3D operator/(TH3D const &h1, TH3D const &h2)
48834883
hnew.SetDirectory(nullptr);
48844884
return hnew;
48854885
}
4886+
4887+
////////////////////////////////////////////////////////////////////////////////
4888+
/// This function calculates the background spectrum in this histogram.
4889+
/// The background is returned as a histogram.
4890+
4891+
TH1 *TH3::ShowBackground3D(Int_t nIterX, Int_t nIterY, Int_t nIterZ, Option_t *option)
4892+
{
4893+
4894+
return (TH1 *)gROOT->ProcessLineFast(
4895+
TString::Format("TSpectrum3::StaticBackground((TH1*)0x%zx,%d,%d,%d,\"%s\")", (size_t)this, nIterX, nIterY, nIterZ, option)
4896+
.Data());
4897+
}

hist/spectrum/inc/TSpectrum.h

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -52,11 +52,11 @@ static Int_t fgIterations; ///< Maximum number of decon iterations (de
5252
TSpectrum();
5353
TSpectrum(Int_t maxpositions, Double_t resolution=1); // resolution is *NOT USED*
5454
~TSpectrum() override;
55-
virtual TH1 *Background(const TH1 *hist,Int_t niter=20, Option_t *option="");
55+
virtual TH1 *Background(const TH1 *hist, Int_t nIter = 20, Option_t *option = "");
5656
TH1 *GetHistogram() const {return fHistogram;}
5757
Int_t GetNPeaks() const {return fNPeaks;}
58-
Double_t *GetPositionX() const {return fPositionX;}
59-
Double_t *GetPositionY() const {return fPositionY;}
58+
Double_t *GetPositionX() const { return fPositionX; }
59+
Double_t *GetPositionY() const { return fPositionY; }
6060
void Print(Option_t *option="") const override;
6161
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05);
6262
static void SetAverageWindow(Int_t w=3); //set average window

hist/spectrum/inc/TSpectrum2.h

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,8 @@ class TH1;
1717

1818
class TSpectrum2 : public TNamed {
1919
protected:
20-
Int_t fMaxPeaks; ///< Maximum number of peaks to be found
21-
Int_t fNPeaks; ///< number of peaks found
20+
Int_t fMaxPeaks; ///< Maximum number of peaks to be found
21+
Int_t fNPeaks; ///< number of peaks found
2222
Double_t *fPosition; ///< [fNPeaks] array of current peak positions
2323
Double_t *fPositionX; ///< [fNPeaks] X position of peaks
2424
Double_t *fPositionY; ///< [fNPeaks] Y position of peaks
@@ -38,11 +38,11 @@ static Int_t fgIterations; ///< Maximum number of decon iterations (def
3838
TSpectrum2();
3939
TSpectrum2(Int_t maxpositions, Double_t resolution=1); // resolution is *NOT USED*
4040
~TSpectrum2() override;
41-
virtual TH1 *Background(const TH1 *hist, Int_t niter=20, Option_t *option="");
41+
virtual TH1 *Background(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Option_t *option = "");
4242
TH1 *GetHistogram() const {return fHistogram;}
4343
Int_t GetNPeaks() const {return fNPeaks;}
44-
Double_t *GetPositionX() const {return fPositionX;}
45-
Double_t *GetPositionY() const {return fPositionY;}
44+
Double_t *GetPositionX() const {return fPositionX;}
45+
Double_t *GetPositionY() const {return fPositionY;}
4646
void Print(Option_t *option="") const override;
4747
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05);
4848
static void SetAverageWindow(Int_t w=3); //set average window
@@ -55,8 +55,8 @@ static Int_t fgIterations; ///< Maximum number of decon iterations (def
5555
const char *Deconvolution(Double_t **source, Double_t **resp, Int_t ssizex, Int_t ssizey,Int_t numberIterations, Int_t numberRepetitions, Double_t boost);
5656
Int_t SearchHighRes(Double_t **source,Double_t **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);
5757

58-
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05);
59-
static TH1 *StaticBackground(const TH1 *hist,Int_t niter=20, Option_t *option="");
58+
static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05);
59+
static TH1 *StaticBackground(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Option_t *option = "");
6060

6161
ClassDefOverride(TSpectrum2,1) //Peak Finder, background estimator, Deconvolution for 2-D histograms
6262
};

hist/spectrum/inc/TSpectrum3.h

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -37,20 +37,21 @@ class TSpectrum3 : public TNamed {
3737
TSpectrum3();
3838
TSpectrum3(Int_t maxpositions, Double_t resolution=1); // resolution is *NOT USED*
3939
~TSpectrum3() override;
40-
virtual const char *Background(const TH1 *hist, Int_t niter, Option_t *option="goff");
40+
virtual TH1 *Background(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Int_t nIterZ = 20, Option_t *option = "goff");
4141
const char *Background(Double_t ***spectrum, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t numberIterationsX,Int_t numberIterationsY, Int_t numberIterationsZ, Int_t direction,Int_t filterType);
4242
const char *Deconvolution(Double_t ***source, const Double_t ***resp, Int_t ssizex, Int_t ssizey, Int_t ssizez,Int_t numberIterations, Int_t numberRepetitions, Double_t boost);
4343
TH1 *GetHistogram() const {return fHistogram;}
4444
Int_t GetNPeaks() const {return fNPeaks;}
45-
Double_t *GetPositionX() const {return fPositionX;}
46-
Double_t *GetPositionY() const {return fPositionY;}
47-
Double_t *GetPositionZ() const {return fPositionZ;}
48-
void Print(Option_t *option="") const override;
45+
Double_t *GetPositionX() const {return fPositionX;}
46+
Double_t *GetPositionY() const {return fPositionY;}
47+
Double_t *GetPositionZ() const {return fPositionZ;}
48+
void Print(Option_t *option="") const override;
4949
virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05);
5050
Int_t SearchFast(const Double_t ***source, Double_t ***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez, Double_t sigma, Double_t threshold, Bool_t markov, Int_t averWindow);
5151
Int_t SearchHighRes(const Double_t ***source,Double_t ***dest, Int_t ssizex, Int_t ssizey, Int_t ssizez, Double_t sigma, Double_t threshold, Bool_t backgroundRemove,Int_t deconIterations, Bool_t markov, Int_t averWindow);
5252
void SetResolution(Double_t resolution=1); // *NOT USED*
5353
const char *SmoothMarkov(Double_t ***source, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t averWindow);
54+
static TH1 *StaticBackground(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Int_t nIterZ = 20, Option_t *option = "");
5455

5556
ClassDefOverride(TSpectrum3,1) //Peak Finder, Background estimator, Markov smoothing and Deconvolution for 3-D histograms
5657
};

hist/spectrum/src/TSpectrum.cxx

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -111,8 +111,8 @@ void TSpectrum::SetDeconIterations(Int_t n)
111111
/// #### Parameters:
112112
///
113113
/// - h: input 1-d histogram
114-
/// - numberIterations, (default value = 20)
115-
/// Increasing numberIterations make the result smoother and lower.
114+
/// - nIter, (default value = 20)
115+
/// Increasing number of iterations makes the result smoother and lower.
116116
/// - option: may contain one of the following options:
117117
///
118118
/// - to set the direction parameter
@@ -141,13 +141,13 @@ void TSpectrum::SetDeconIterations(Int_t n)
141141
/// as the input histogram h, but only bins from `binmin` to `binmax` will be filled
142142
/// with the estimated background.
143143

144-
TH1 *TSpectrum::Background(const TH1 * h, Int_t numberIterations,
144+
TH1 *TSpectrum::Background(const TH1 * h, Int_t nIter,
145145
Option_t * option)
146146
{
147147
if (h == nullptr) return nullptr;
148148
Int_t dimension = h->GetDimension();
149-
if (dimension > 1) {
150-
Error("Search", "Only implemented for 1-d histograms");
149+
if (dimension != 1) {
150+
Error("Background", "Only implemented for 1-d histograms");
151151
return nullptr;
152152
}
153153
TString opt = option;
@@ -180,7 +180,7 @@ TH1 *TSpectrum::Background(const TH1 * h, Int_t numberIterations,
180180
for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first);
181181

182182
//find background (source is input and in output contains the background
183-
Background(source,size,numberIterations, direction, filterOrder,smoothing,
183+
Background(source,size,nIter, direction, filterOrder,smoothing,
184184
smoothWindow,compton);
185185

186186
//create output histogram containing background

hist/spectrum/src/TSpectrum2.cxx

Lines changed: 74 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -44,8 +44,9 @@
4444
#include "TSpectrum2.h"
4545
#include "TPolyMarker.h"
4646
#include "TList.h"
47-
#include "TH1.h"
47+
#include "TH2.h"
4848
#include "TMath.h"
49+
#include "TVirtualPad.h"
4950
#define PEAK_WINDOW 1024
5051

5152
Int_t TSpectrum2::fgIterations = 3;
@@ -122,28 +123,15 @@ void TSpectrum2::SetDeconIterations(Int_t n)
122123
///
123124
/// Function parameters:
124125
/// - h: input 2-d histogram
125-
/// - numberIterations, (default value = 20)
126-
/// Increasing numberIterations make the result smoother and lower.
126+
/// - nIterX, nIterY, (default value = 20), iterations for X and Y
127+
/// Increasing number of iterations make the result smoother and lower.
127128
/// - option: may contain one of the following options
128-
/// - to set the direction parameter
129-
/// "BackIncreasingWindow". By default the direction is BackDecreasingWindow
130-
/// - filterOrder-order of clipping filter. Possible values:
131-
/// - "BackOrder2" (default)
132-
/// - "BackOrder4"
133-
/// - "BackOrder6"
134-
/// - "BackOrder8"
135-
/// - "nosmoothing"- if selected, the background is not smoothed
136-
/// By default the background is smoothed.
137-
/// - smoothWindow-width of smoothing window. Possible values:
138-
/// - "BackSmoothing3" (default)
139-
/// - "BackSmoothing5"
140-
/// - "BackSmoothing7"
141-
/// - "BackSmoothing9"
142-
/// - "BackSmoothing11"
143-
/// - "BackSmoothing13"
144-
/// - "BackSmoothing15"
145-
/// - "Compton" if selected the estimation of Compton edge
146-
/// will be included.
129+
/// - direction of change of clipping window
130+
/// - possible values=kBackIncreasingWindow
131+
/// kBackDecreasingWindow
132+
/// - filterType-determines the algorithm of the filtering
133+
/// - possible values=kBackSuccessiveFiltering
134+
/// kBackOneStepFiltering
147135
/// - "same" : if this option is specified, the resulting background
148136
/// histogram is superimposed on the picture in the current pad.
149137
///
@@ -153,12 +141,67 @@ void TSpectrum2::SetDeconIterations(Int_t n)
153141
/// as the input histogram h, but only bins from binmin to binmax will be filled
154142
/// with the estimated background.
155143

156-
TH1 *TSpectrum2::Background(const TH1 * h, Int_t number_of_iterations,
157-
Option_t * option)
144+
TH1 *TSpectrum2::Background(const TH1 *h, Int_t nIterX, Int_t nIterY, Option_t *option)
158145
{
159-
Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn"
160-
, h->GetName(), number_of_iterations, option);
161-
return nullptr;
146+
if (h == nullptr)
147+
return nullptr;
148+
Int_t dimension = h->GetDimension();
149+
if (dimension != 2) {
150+
Error("Background", "Only implemented for 2-d histograms");
151+
return nullptr;
152+
}
153+
TString opt = option;
154+
opt.ToLower();
155+
156+
// set options
157+
Int_t direction = kBackDecreasingWindow;
158+
if (opt.Contains("backincreasingwindow"))
159+
direction = kBackIncreasingWindow;
160+
Int_t filterType = kBackSuccessiveFiltering;
161+
if (opt.Contains("backonestepfiltering"))
162+
filterType = kBackOneStepFiltering;
163+
Int_t firstX = h->GetXaxis()->GetFirst();
164+
Int_t lastX = h->GetXaxis()->GetLast();
165+
Int_t sizeX = lastX - firstX + 1;
166+
Int_t firstY = h->GetYaxis()->GetFirst();
167+
Int_t lastY = h->GetYaxis()->GetLast();
168+
Int_t sizeY = lastY - firstY + 1;
169+
Int_t i, j;
170+
Double_t **source = new Double_t *[sizeX];
171+
for (i = 0; i < sizeX; i++) {
172+
source[i] = new Double_t[sizeY];
173+
for (j = 0; j < sizeY; j++)
174+
source[i][j] = h->GetBinContent(i + firstX, j + firstY);
175+
}
176+
177+
// find background (source is input and in output contains the background
178+
Background(source, sizeX, sizeY, nIterX, nIterY, direction, filterType);
179+
180+
// create output histogram containing background
181+
// only bins in the range of the input histogram are filled
182+
Int_t nch = strlen(h->GetName());
183+
char *hbname = new char[nch + 20];
184+
snprintf(hbname, nch + 20, "%s_background", h->GetName());
185+
TH2 *hb = (TH2 *)h->Clone(hbname);
186+
hb->Reset();
187+
hb->GetListOfFunctions()->Delete();
188+
for (i = 0; i < sizeX; i++)
189+
for (j = 0; j < sizeY; j++)
190+
hb->SetBinContent(i + firstX, j + firstY, source[i][j]);
191+
hb->SetEntries(sizeX * sizeY);
192+
193+
// if option "same is specified, draw the result in the pad
194+
if (opt.Contains("same")) {
195+
if (gPad)
196+
delete gPad->GetPrimitive(hbname);
197+
hb->Draw("same");
198+
}
199+
for (i = 0; i < sizeX; i++) {
200+
delete[] source[i];
201+
}
202+
delete[] source;
203+
delete[] hbname;
204+
return hb;
162205
}
163206

164207
////////////////////////////////////////////////////////////////////////////////
@@ -1706,7 +1749,7 @@ Int_t TSpectrum2::SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex
17061749
}
17071750

17081751
////////////////////////////////////////////////////////////////////////////////
1709-
/// static function (called by TH1), interface to TSpectrum2::Search
1752+
/// static function (called by TH2), interface to TSpectrum2::Search
17101753

17111754
Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold)
17121755
{
@@ -1715,10 +1758,10 @@ Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option
17151758
}
17161759

17171760
////////////////////////////////////////////////////////////////////////////////
1718-
/// static function (called by TH1), interface to TSpectrum2::Background
1761+
/// static function (called by TH2), interface to TSpectrum2::Background
17191762

1720-
TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option)
1763+
TH1 *TSpectrum2::StaticBackground(const TH1 *hist, Int_t nIterX, Int_t nIterY, Option_t *option)
17211764
{
17221765
TSpectrum2 s;
1723-
return s.Background(hist,niter,option);
1766+
return s.Background(hist, nIterX, nIterY, option);
17241767
}

0 commit comments

Comments
 (0)