diff --git a/hist/hist/inc/TH2.h b/hist/hist/inc/TH2.h index 612c4c31b6bfc..e7d8101aecc58 100644 --- a/hist/hist/inc/TH2.h +++ b/hist/hist/inc/TH2.h @@ -130,9 +130,9 @@ class TH2 : public TH1 { virtual void SetShowProjectionX(Int_t nbins=1); // *MENU* virtual void SetShowProjectionY(Int_t nbins=1); // *MENU* virtual void SetShowProjectionXY(Int_t nbinsY=1, Int_t nbinsX=1); // *MENU* - TH1 *ShowBackground(Int_t niter=20, Option_t *option="same") override; - Int_t ShowPeaks(Double_t sigma=2, Option_t *option="", Double_t threshold=0.05) override; // *MENU* - void Smooth(Int_t ntimes=1, Option_t *option="") override; // *MENU* + virtual TH1 *ShowBackground2D(Int_t nIterX = 20, Int_t nIterY = 20, Option_t *option = "same"); + Int_t ShowPeaks(Double_t sigma = 2, Option_t *option = "", Double_t threshold = 0.05) override; // *MENU* + void Smooth(Int_t ntimes = 1, Option_t *option = "") override; // *MENU* ClassDefOverride(TH2,5) //2-Dim histogram base class }; diff --git a/hist/hist/inc/TH3.h b/hist/hist/inc/TH3.h index d05262975566d..92b6b67400574 100644 --- a/hist/hist/inc/TH3.h +++ b/hist/hist/inc/TH3.h @@ -144,6 +144,7 @@ class TH3 : public TH1, public TAtt3D { void SetBinContent(Int_t bin, Int_t, Double_t content) override { SetBinContent(bin, content); } void SetBinContent(Int_t binx, Int_t biny, Int_t binz, Double_t content) override { SetBinContent(GetBin(binx, biny, binz), content); } virtual void SetShowProjection(const char *option="xy",Int_t nbins=1); // *MENU* + virtual TH1 *ShowBackground3D(Int_t nIterX = 20, Int_t nIterY = 20, Int_t nIterZ = 20, Option_t *option = "same"); protected: diff --git a/hist/hist/src/TH2.cxx b/hist/hist/src/TH2.cxx index 5c8af92e19239..5792ac8f27906 100644 --- a/hist/hist/src/TH2.cxx +++ b/hist/hist/src/TH2.cxx @@ -2636,13 +2636,13 @@ void TH2::SetShowProjectionXY(Int_t nbinsY, Int_t nbinsX) //////////////////////////////////////////////////////////////////////////////// /// This function calculates the background spectrum in this histogram. /// The background is returned as a histogram. -/// to be implemented (may be) -TH1 *TH2::ShowBackground(Int_t niter, Option_t *option) +TH1 *TH2::ShowBackground2D(Int_t nIterX, Int_t nIterY, Option_t *option) { - return (TH1 *)gROOT->ProcessLineFast(TString::Format("TSpectrum2::StaticBackground((TH1*)0x%zx,%d,\"%s\")", - (size_t)this, niter, option).Data()); + return (TH1 *)gROOT->ProcessLineFast( + TString::Format("TSpectrum2::StaticBackground((TH1*)0x%zx,%d,%d,\"%s\")", (size_t)this, nIterX, nIterY, option) + .Data()); } diff --git a/hist/hist/src/TH3.cxx b/hist/hist/src/TH3.cxx index 6e611ab3efe50..6fd8244756202 100644 --- a/hist/hist/src/TH3.cxx +++ b/hist/hist/src/TH3.cxx @@ -4883,3 +4883,15 @@ TH3D operator/(TH3D const &h1, TH3D const &h2) hnew.SetDirectory(nullptr); return hnew; } + +//////////////////////////////////////////////////////////////////////////////// +/// This function calculates the background spectrum in this histogram. +/// The background is returned as a histogram. + +TH1 *TH3::ShowBackground3D(Int_t nIterX, Int_t nIterY, Int_t nIterZ, Option_t *option) +{ + + return (TH1 *)gROOT->ProcessLineFast( + TString::Format("TSpectrum3::StaticBackground((TH1*)0x%zx,%d,%d,%d,\"%s\")", (size_t)this, nIterX, nIterY, nIterZ, option) + .Data()); +} diff --git a/hist/spectrum/inc/TSpectrum.h b/hist/spectrum/inc/TSpectrum.h index 127015d7b2ee8..5d2d3af37641e 100644 --- a/hist/spectrum/inc/TSpectrum.h +++ b/hist/spectrum/inc/TSpectrum.h @@ -52,11 +52,11 @@ static Int_t fgIterations; ///< Maximum number of decon iterations (de TSpectrum(); TSpectrum(Int_t maxpositions, Double_t resolution=1); // resolution is *NOT USED* ~TSpectrum() override; - virtual TH1 *Background(const TH1 *hist,Int_t niter=20, Option_t *option=""); + virtual TH1 *Background(const TH1 *hist, Int_t nIter = 20, Option_t *option = ""); TH1 *GetHistogram() const {return fHistogram;} Int_t GetNPeaks() const {return fNPeaks;} - Double_t *GetPositionX() const {return fPositionX;} - Double_t *GetPositionY() const {return fPositionY;} + Double_t *GetPositionX() const { return fPositionX; } + Double_t *GetPositionY() const { return fPositionY; } void Print(Option_t *option="") const override; virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05); static void SetAverageWindow(Int_t w=3); //set average window diff --git a/hist/spectrum/inc/TSpectrum2.h b/hist/spectrum/inc/TSpectrum2.h index 9d35caade0997..01d25802e9706 100644 --- a/hist/spectrum/inc/TSpectrum2.h +++ b/hist/spectrum/inc/TSpectrum2.h @@ -17,8 +17,8 @@ class TH1; class TSpectrum2 : public TNamed { protected: - Int_t fMaxPeaks; ///< Maximum number of peaks to be found - Int_t fNPeaks; ///< number of peaks found + Int_t fMaxPeaks; ///< Maximum number of peaks to be found + Int_t fNPeaks; ///< number of peaks found Double_t *fPosition; ///< [fNPeaks] array of current peak positions Double_t *fPositionX; ///< [fNPeaks] X position of peaks Double_t *fPositionY; ///< [fNPeaks] Y position of peaks @@ -38,11 +38,11 @@ static Int_t fgIterations; ///< Maximum number of decon iterations (def TSpectrum2(); TSpectrum2(Int_t maxpositions, Double_t resolution=1); // resolution is *NOT USED* ~TSpectrum2() override; - virtual TH1 *Background(const TH1 *hist, Int_t niter=20, Option_t *option=""); + virtual TH1 *Background(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Option_t *option = ""); TH1 *GetHistogram() const {return fHistogram;} Int_t GetNPeaks() const {return fNPeaks;} - Double_t *GetPositionX() const {return fPositionX;} - Double_t *GetPositionY() const {return fPositionY;} + Double_t *GetPositionX() const {return fPositionX;} + Double_t *GetPositionY() const {return fPositionY;} void Print(Option_t *option="") const override; virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="", Double_t threshold=0.05); static void SetAverageWindow(Int_t w=3); //set average window @@ -55,8 +55,8 @@ static Int_t fgIterations; ///< Maximum number of decon iterations (def const char *Deconvolution(Double_t **source, Double_t **resp, Int_t ssizex, Int_t ssizey,Int_t numberIterations, Int_t numberRepetitions, Double_t boost); 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); - static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05); - static TH1 *StaticBackground(const TH1 *hist,Int_t niter=20, Option_t *option=""); + static Int_t StaticSearch(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05); + static TH1 *StaticBackground(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Option_t *option = ""); ClassDefOverride(TSpectrum2,1) //Peak Finder, background estimator, Deconvolution for 2-D histograms }; diff --git a/hist/spectrum/inc/TSpectrum3.h b/hist/spectrum/inc/TSpectrum3.h index 61f8edc05c26f..88ac1050e31c8 100644 --- a/hist/spectrum/inc/TSpectrum3.h +++ b/hist/spectrum/inc/TSpectrum3.h @@ -37,20 +37,21 @@ class TSpectrum3 : public TNamed { TSpectrum3(); TSpectrum3(Int_t maxpositions, Double_t resolution=1); // resolution is *NOT USED* ~TSpectrum3() override; - virtual const char *Background(const TH1 *hist, Int_t niter, Option_t *option="goff"); + virtual TH1 *Background(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Int_t nIterZ = 20, Option_t *option = "goff"); 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); 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); TH1 *GetHistogram() const {return fHistogram;} Int_t GetNPeaks() const {return fNPeaks;} - Double_t *GetPositionX() const {return fPositionX;} - Double_t *GetPositionY() const {return fPositionY;} - Double_t *GetPositionZ() const {return fPositionZ;} - void Print(Option_t *option="") const override; + Double_t *GetPositionX() const {return fPositionX;} + Double_t *GetPositionY() const {return fPositionY;} + Double_t *GetPositionZ() const {return fPositionZ;} + void Print(Option_t *option="") const override; virtual Int_t Search(const TH1 *hist, Double_t sigma=2, Option_t *option="goff", Double_t threshold=0.05); 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); 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); void SetResolution(Double_t resolution=1); // *NOT USED* const char *SmoothMarkov(Double_t ***source, Int_t ssizex, Int_t ssizey, Int_t ssizez, Int_t averWindow); + static TH1 *StaticBackground(const TH1 *hist, Int_t nIterX = 20, Int_t nIterY = 20, Int_t nIterZ = 20, Option_t *option = ""); ClassDefOverride(TSpectrum3,1) //Peak Finder, Background estimator, Markov smoothing and Deconvolution for 3-D histograms }; diff --git a/hist/spectrum/src/TSpectrum.cxx b/hist/spectrum/src/TSpectrum.cxx index 8051b653b0aca..bb5455d5d9aab 100644 --- a/hist/spectrum/src/TSpectrum.cxx +++ b/hist/spectrum/src/TSpectrum.cxx @@ -111,8 +111,8 @@ void TSpectrum::SetDeconIterations(Int_t n) /// #### Parameters: /// /// - h: input 1-d histogram -/// - numberIterations, (default value = 20) -/// Increasing numberIterations make the result smoother and lower. +/// - nIter, (default value = 20) +/// Increasing number of iterations makes the result smoother and lower. /// - option: may contain one of the following options: /// /// - to set the direction parameter @@ -141,13 +141,13 @@ void TSpectrum::SetDeconIterations(Int_t n) /// as the input histogram h, but only bins from `binmin` to `binmax` will be filled /// with the estimated background. -TH1 *TSpectrum::Background(const TH1 * h, Int_t numberIterations, +TH1 *TSpectrum::Background(const TH1 * h, Int_t nIter, Option_t * option) { if (h == nullptr) return nullptr; Int_t dimension = h->GetDimension(); - if (dimension > 1) { - Error("Search", "Only implemented for 1-d histograms"); + if (dimension != 1) { + Error("Background", "Only implemented for 1-d histograms"); return nullptr; } TString opt = option; @@ -180,7 +180,7 @@ TH1 *TSpectrum::Background(const TH1 * h, Int_t numberIterations, for (i = 0; i < size; i++) source[i] = h->GetBinContent(i + first); //find background (source is input and in output contains the background - Background(source,size,numberIterations, direction, filterOrder,smoothing, + Background(source,size,nIter, direction, filterOrder,smoothing, smoothWindow,compton); //create output histogram containing background diff --git a/hist/spectrum/src/TSpectrum2.cxx b/hist/spectrum/src/TSpectrum2.cxx index 2cf6f6e5b5acf..b0a9e588955ca 100644 --- a/hist/spectrum/src/TSpectrum2.cxx +++ b/hist/spectrum/src/TSpectrum2.cxx @@ -44,8 +44,9 @@ #include "TSpectrum2.h" #include "TPolyMarker.h" #include "TList.h" -#include "TH1.h" +#include "TH2.h" #include "TMath.h" +#include "TVirtualPad.h" #define PEAK_WINDOW 1024 Int_t TSpectrum2::fgIterations = 3; @@ -122,28 +123,15 @@ void TSpectrum2::SetDeconIterations(Int_t n) /// /// Function parameters: /// - h: input 2-d histogram -/// - numberIterations, (default value = 20) -/// Increasing numberIterations make the result smoother and lower. +/// - nIterX, nIterY, (default value = 20), iterations for X and Y +/// Increasing number of iterations make the result smoother and lower. /// - option: may contain one of the following options -/// - to set the direction parameter -/// "BackIncreasingWindow". By default the direction is BackDecreasingWindow -/// - filterOrder-order of clipping filter. Possible values: -/// - "BackOrder2" (default) -/// - "BackOrder4" -/// - "BackOrder6" -/// - "BackOrder8" -/// - "nosmoothing"- if selected, the background is not smoothed -/// By default the background is smoothed. -/// - smoothWindow-width of smoothing window. Possible values: -/// - "BackSmoothing3" (default) -/// - "BackSmoothing5" -/// - "BackSmoothing7" -/// - "BackSmoothing9" -/// - "BackSmoothing11" -/// - "BackSmoothing13" -/// - "BackSmoothing15" -/// - "Compton" if selected the estimation of Compton edge -/// will be included. +/// - direction of change of clipping window +/// - possible values=kBackIncreasingWindow +/// kBackDecreasingWindow +/// - filterType-determines the algorithm of the filtering +/// - possible values=kBackSuccessiveFiltering +/// kBackOneStepFiltering /// - "same" : if this option is specified, the resulting background /// histogram is superimposed on the picture in the current pad. /// @@ -153,12 +141,67 @@ void TSpectrum2::SetDeconIterations(Int_t n) /// as the input histogram h, but only bins from binmin to binmax will be filled /// with the estimated background. -TH1 *TSpectrum2::Background(const TH1 * h, Int_t number_of_iterations, - Option_t * option) +TH1 *TSpectrum2::Background(const TH1 *h, Int_t nIterX, Int_t nIterY, Option_t *option) { - Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn" - , h->GetName(), number_of_iterations, option); - return nullptr; + if (h == nullptr) + return nullptr; + Int_t dimension = h->GetDimension(); + if (dimension != 2) { + Error("Background", "Only implemented for 2-d histograms"); + return nullptr; + } + TString opt = option; + opt.ToLower(); + + // set options + Int_t direction = kBackDecreasingWindow; + if (opt.Contains("backincreasingwindow")) + direction = kBackIncreasingWindow; + Int_t filterType = kBackSuccessiveFiltering; + if (opt.Contains("backonestepfiltering")) + filterType = kBackOneStepFiltering; + Int_t firstX = h->GetXaxis()->GetFirst(); + Int_t lastX = h->GetXaxis()->GetLast(); + Int_t sizeX = lastX - firstX + 1; + Int_t firstY = h->GetYaxis()->GetFirst(); + Int_t lastY = h->GetYaxis()->GetLast(); + Int_t sizeY = lastY - firstY + 1; + Int_t i, j; + Double_t **source = new Double_t *[sizeX]; + for (i = 0; i < sizeX; i++) { + source[i] = new Double_t[sizeY]; + for (j = 0; j < sizeY; j++) + source[i][j] = h->GetBinContent(i + firstX, j + firstY); + } + + // find background (source is input and in output contains the background + Background(source, sizeX, sizeY, nIterX, nIterY, direction, filterType); + + // create output histogram containing background + // only bins in the range of the input histogram are filled + Int_t nch = strlen(h->GetName()); + char *hbname = new char[nch + 20]; + snprintf(hbname, nch + 20, "%s_background", h->GetName()); + TH2 *hb = (TH2 *)h->Clone(hbname); + hb->Reset(); + hb->GetListOfFunctions()->Delete(); + for (i = 0; i < sizeX; i++) + for (j = 0; j < sizeY; j++) + hb->SetBinContent(i + firstX, j + firstY, source[i][j]); + hb->SetEntries(sizeX * sizeY); + + // if option "same is specified, draw the result in the pad + if (opt.Contains("same")) { + if (gPad) + delete gPad->GetPrimitive(hbname); + hb->Draw("same"); + } + for (i = 0; i < sizeX; i++) { + delete[] source[i]; + } + delete[] source; + delete[] hbname; + return hb; } //////////////////////////////////////////////////////////////////////////////// @@ -1706,7 +1749,7 @@ Int_t TSpectrum2::SearchHighRes(Double_t **source, Double_t **dest, Int_t ssizex } //////////////////////////////////////////////////////////////////////////////// -/// static function (called by TH1), interface to TSpectrum2::Search +/// static function (called by TH2), interface to TSpectrum2::Search Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold) { @@ -1715,10 +1758,10 @@ Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option } //////////////////////////////////////////////////////////////////////////////// -/// static function (called by TH1), interface to TSpectrum2::Background +/// static function (called by TH2), interface to TSpectrum2::Background -TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option) +TH1 *TSpectrum2::StaticBackground(const TH1 *hist, Int_t nIterX, Int_t nIterY, Option_t *option) { TSpectrum2 s; - return s.Background(hist,niter,option); + return s.Background(hist, nIterX, nIterY, option); } diff --git a/hist/spectrum/src/TSpectrum3.cxx b/hist/spectrum/src/TSpectrum3.cxx index 3ee4baf60899f..0d4c75e8f6b84 100644 --- a/hist/spectrum/src/TSpectrum3.cxx +++ b/hist/spectrum/src/TSpectrum3.cxx @@ -48,8 +48,9 @@ */ #include "TSpectrum3.h" -#include "TH1.h" +#include "TH3.h" #include "TMath.h" +#include "TVirtualPad.h" #define PEAK_WINDOW 1024 @@ -108,17 +109,92 @@ TSpectrum3::~TSpectrum3() /// This function calculates background spectrum from source in h. /// The result is placed in the vector pointed by spectrum pointer. /// -/// Function parameters: -/// - spectrum: pointer to the vector of source spectrum -/// - size: length of spectrum and working space vectors -/// - number_of_iterations, for details we refer to manual - -const char *TSpectrum3::Background(const TH1 * h, Int_t number_of_iterations, - Option_t * option) +/// Function parameters: +/// - h: input 3-d histogram +/// - nIterX, nIterY, nIterZ, iterations for X and Y and Z axes +/// Increasing number of iterations make the result smoother and lower. +/// - option: may contain one of the following options +/// - direction of change of clipping window +/// - possible values=kBackIncreasingWindow +/// kBackDecreasingWindow +/// - filterType-determines the algorithm of the filtering +/// - possible values=kBackSuccessiveFiltering +/// kBackOneStepFiltering +/// - "same" : if this option is specified, the resulting background +/// histogram is superimposed on the picture in the current pad. +/// Default is "goff" ie no graphics output + +TH1 *TSpectrum3::Background(const TH1 *h, Int_t nIterX, Int_t nIterY, Int_t nIterZ, Option_t *option) { - Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn" - , h->GetName(), number_of_iterations, option); - return nullptr; + if (h == nullptr) + return nullptr; + Int_t dimension = h->GetDimension(); + if (dimension != 3) { + Error("Background", "Only implemented for 3-d histograms"); + return nullptr; + } + TString opt = option; + opt.ToLower(); + + // set options + Int_t direction = kBackDecreasingWindow; + if (opt.Contains("backincreasingwindow")) + direction = kBackIncreasingWindow; + Int_t filterType = kBackSuccessiveFiltering; + if (opt.Contains("backonestepfiltering")) + filterType = kBackOneStepFiltering; + Int_t firstX = h->GetXaxis()->GetFirst(); + Int_t lastX = h->GetXaxis()->GetLast(); + Int_t sizeX = lastX - firstX + 1; + Int_t firstY = h->GetYaxis()->GetFirst(); + Int_t lastY = h->GetYaxis()->GetLast(); + Int_t sizeY = lastY - firstY + 1; + Int_t firstZ = h->GetZaxis()->GetFirst(); + Int_t lastZ = h->GetZaxis()->GetLast(); + Int_t sizeZ = lastZ - firstZ + 1; + Int_t i, j, k; + Double_t ***source = new Double_t **[sizeX]; + for (i = 0; i < sizeX; i++) { + source[i] = new Double_t *[sizeY]; + for (j = 0; j < sizeY; j++) { + source[i][j] = new Double_t[sizeZ]; + for (k = 0; k < sizeZ; k++) + source[i][j][k] = h->GetBinContent(i + firstX, j + firstY, k + firstZ); + } + } + + // find background (source is input and in output contains the background + Background(source, sizeX, sizeY, sizeZ, nIterX, nIterY, nIterZ, direction, filterType); + + // create output histogram containing background + // only bins in the range of the input histogram are filled + Int_t nch = strlen(h->GetName()); + char *hbname = new char[nch + 20]; + snprintf(hbname, nch + 20, "%s_background", h->GetName()); + TH3 *hb = (TH3 *)h->Clone(hbname); + hb->Reset(); + hb->GetListOfFunctions()->Delete(); + for (i = 0; i < sizeX; i++) + for (j = 0; j < sizeY; j++) + for (k = 0; k < sizeZ; k++) + hb->SetBinContent(i + firstX, j + firstY, k + firstZ, source[i][j][k]); + hb->SetEntries(sizeX * sizeY * sizeZ); + + // if option "same is specified, draw the result in the pad + if (opt.Contains("same")) { + if (gPad) + delete gPad->GetPrimitive(hbname); + hb->Draw("same"); + } + for (i = 0; i < sizeX; i++) { + for (j = 0; j < sizeY; j++) { + delete[] source[i][j]; + } + delete[] source[i]; + } + delete[] source; + delete[] hbname; + return hb; } //////////////////////////////////////////////////////////////////////////////// @@ -4334,3 +4410,12 @@ Int_t TSpectrum3::SearchFast(const Double_t***source, Double_t***dest, Int_t ssi fNPeaks = peak_index; return fNPeaks; } + +//////////////////////////////////////////////////////////////////////////////// +/// static function (called by TH3), interface to TSpectrum3::Background + +TH1 *TSpectrum3::StaticBackground(const TH1 *hist, Int_t nIterX, Int_t nIterY, Int_t nIterZ, Option_t *option) +{ + TSpectrum3 s; + return s.Background(hist, nIterX, nIterY, nIterZ, option); +}