EdbPeak2 Class Reference

peak analyser for EdbH2 More...

#include <EdbCell2.h>

Inheritance diagram for EdbPeak2:

Collaboration diagram for EdbPeak2:

Public Member Functions
void	Delete ()
	EdbPeak2 ()
	EdbPeak2 (const EdbH2 &h)
float	EstimatePeakMeanPosition (int iv[2], int ir[2], float &x, float &y)
float	EstimatePeakVolume (int ipeak)
float	EstimatePeakVolumeSafe (int ipeak)
float	FindGlobalPeak (float &x, float &y, float ratio=0.1)
int	FindPeak (float &x, float &y)
int	FindPeak (float v[2])
int	FindPeak (int iv[2])
float	FindPeak9 (float &x, float &y)
void	Init (const EdbH2 &h, int npeaks=10)
void	InitPeaks (int npeaks)
float	Mean (int i=0) const
float	Mean3 (int i=0) const
float	Peak (int i=0) const
void	Print ()
float	ProbPeak ()
float	ProbPeak (float &x, float &y)
float	ProbPeak (int iv[2], int ir[2])
int	ProbPeaks (int npeak)
int	ProbPeaks (int npeak, int ir[2])
float	Smooth (Option_t *option="k5a")
int	WipePeak (int iv[2], int ir[2])
float	Xmean ()
float	Ymean ()
	~EdbPeak2 ()
Public Member Functions inherited from EdbH2
void	AddBin (int jcell, int n)
int	Bin (float x, float y) const
int	Bin (int iv[2]) const
int	Bin (int ix, int iy) const
int	Bin (int j) const
void	CleanCells ()
void	Copy (const EdbH2 &h)
void	Delete ()
int	DiscardHighCells (int nmax)
TH2F *	DrawH2 (const char *name="plot2d", const char *title="EdbH2plot2D")
TH1I *	DrawSpectrum (const char *name="plot1d", const char *title="EdbH2 DrawSpectrun")
	EdbH2 ()
	EdbH2 (const EdbH2 &h)
	EdbH2 (int nx, float minx, float maxx, int ny, float miny, float maxy)
int	Fill (float x, float y)
int	Fill (float x, float y, int n)
int	InitH2 (const EdbH2 &h)
int	InitH2 (int n[2], float min[2], float max[2])
int	InitH2 (int nx, float minx, float maxx, int ny, float miny, float maxy)
Long_t	Integral ()
Long_t	Integral (int iv[2], int ir[2])
int	IX (float x) const
int	IX (int jcell) const
int	IY (float y) const
int	IY (int jcell) const
int	Jcell (float v[2]) const
int	Jcell (float x, float y) const
int	Jcell (int ix, int iy) const
int	MaxBin ()
Float_t	Mean ()
int	Ncell () const
int	NX () const
int	NY () const
void	PrintStat ()
EdbH1 *	ProjectionX ()
EdbH1 *	ProjectionY ()
void	Set0 ()
void	SetBin (int ix, int iy, int n)
void	SetBin (int j, int n)
float	X (int ix) const
float	Xbin () const
float	Xj (int j) const
float	Xmax () const
float	XmaxA (float level=0)
float	Xmin () const
float	XminA (float level=0)
float	Y (int iy) const
float	Ybin () const
float	Yj (int j) const
float	Ymax () const
float	YmaxA (float level=0)
float	Ymin () const
float	YminA (float level=0)
	~EdbH2 ()

Public Attributes
TArrayF	eMean
TArrayF	eMean3
Float_t	eNorm
	the norm-factor in case of the smoothing applied More...
Int_t	eNpeaks
	number of found peaks More...
TArrayF	ePeak

Additional Inherited Members
Protected Attributes inherited from EdbH2
Float_t	eBin [2]
	bin size More...
Float_t	eMax [2]
	max More...
Float_t	eMin [2]
	min More...
Int_t	eN [2]
	divisions More...
Int_t *	eNC
	[eNcell] number of objects/cell More...
Int_t	eNcell
	eNx*eNy More...

Detailed Description

peak analyser for EdbH2

Constructor & Destructor Documentation

EdbPeak2() [1/2]

EdbPeak2::EdbPeak2 ( )

inline

{ InitPeaks(10); }

EdbPeak2() [2/2]

EdbPeak2::EdbPeak2 ( const EdbH2 &  h )

inline

: EdbH2( h ) { InitPeaks(10); }

~EdbPeak2()

EdbPeak2::~EdbPeak2 ( )

inline

{}

Member Function Documentation

Delete()

void EdbPeak2::Delete

(

)

{
  ((EdbH2*)this)->Delete();
  eNpeaks=0;
}

EstimatePeakMeanPosition()

float EdbPeak2::EstimatePeakMeanPosition	(	int	iv[2],
		int	ir[2],
		float &	x,
		float &	y
	)

{
  double x0=0, y0=0, volume=0;
  for(int ix=iv[0]-ir[0]; ix<=iv[0]+ir[0]; ix++) 
    for(int iy=iv[1]-ir[1]; iy<=iv[1]+ir[1]; iy++) 
      { 
    x0     += Bin(ix,iy)*X(ix);
    y0     += Bin(ix,iy)*Y(iy);
    volume += Bin(ix,iy);
      }
  if(volume>0) { x =  x0/volume;  y = y0/volume; }
  return volume;
}

EstimatePeakVolume()

float EdbPeak2::EstimatePeakVolume

(

int

ipeak

)

{
  if(ipeak>(eNpeaks-1)||ipeak<0)  return 0;
  return ePeak[ipeak] + 8*eMean3[ipeak] - 9*eMean[ipeak];
}

EstimatePeakVolumeSafe()

float EdbPeak2::EstimatePeakVolumeSafe

(

int

ipeak

)

add to the background one standard deviation for the safety

{
  if(ipeak>(eNpeaks-1)||ipeak<0)  return 0;
  return ePeak[ipeak] + 8*eMean3[ipeak] - 9*(eMean[ipeak] + Sqrt(eMean[ipeak]));
}

FindGlobalPeak()

float EdbPeak2::FindGlobalPeak	(	float &	x,
		float &	y,
		float	ratio = 0.1
	)

to find the center of wide multibin peak select ratio*Nbin of highest bins and calculate CoG

{
  int n = Ncell();
  TArrayI ind(n);
  Sort(n, eNC , ind.GetArray(), 0);  // sort encreasing order
  Int_t nwipe = Min( (Int_t)((1.-ratio)*n), n-1 );
  double mean=0;
  for(int i=0; i<nwipe; i++) {
    mean += eNC[ind[i]];
    eNC[ind[i]] = 0;
  }
  mean /= nwipe;
  x = Xmean();
  y = Ymean();
  float peakvolume = Integral() - (n-nwipe)*mean;
  Log(3,"EdbPeak2::FindGlobalPeak","of volume %8.2f at %10.2f %10.2f using %d bins", 
      peakvolume,x,y, n-nwipe );
  return peakvolume;
}

FindPeak() [1/3]

int EdbPeak2::FindPeak	(	float &	x,
		float &	y
	)

{
  int iv[2];
  int peak=FindPeak(iv);
  x = X(iv[0]);
  y = Y(iv[1]);
  return peak;
}

FindPeak() [2/3]

int EdbPeak2::FindPeak

(

float

v[2]

)

{
  int iv[2];
  int peak=FindPeak(iv);
  v[0] = X(iv[0]);
  v[1] = Y(iv[1]);
  return peak;
}

FindPeak() [3/3]

int EdbPeak2::FindPeak

(

int

iv[2]

)

{
  int peak=0;
  for(int ix=0; ix<eN[0]; ix++)
    for(int iy=0; iy<eN[1]; iy++) {
      int j = Jcell(ix,iy);
      if(eNC[j]>peak) {
        peak =eNC[j];
        iv[0] = ix;
        iv[1] = iy;
      }
    }
  Log(4,"EdbPeak2::FindPeak","of %d at %d %d", peak, iv[0],iv[1] );
  return peak;
}

FindPeak9()

float EdbPeak2::FindPeak9	(	float &	x,
		float &	y
	)

particular case of 3x3 peak neiboring, return peak volume and estimate mean position

{
  int iv[2], ir[2]={1,1};
  FindPeak(iv);
  return EstimatePeakMeanPosition(iv, ir, x, y);
}

Init()

void EdbPeak2::Init	(	const EdbH2 &	h,
		int	npeaks = 10
	)

{
  Delete();
  ((EdbH2*)this)->Copy(h);
  InitPeaks(npeaks);
}

InitPeaks()

void EdbPeak2::InitPeaks

(

int

npeaks

)

{
  eNpeaks=0;
  ePeak.Set(npeaks);  ePeak.Reset(0);
  eMean3.Set(npeaks); eMean3.Reset(0);
  eMean.Set(npeaks);  eMean.Reset(0);
  eNorm = 1.;  // default norm factor;
}

Mean()

float EdbPeak2::Mean ( int  i = 0 ) const

inline

{ if(i>=0&&i<eNpeaks) return eMean[i];  else return -1; }

Mean3()

float EdbPeak2::Mean3 ( int  i = 0 ) const

inline

{ if(i>=0&&i<eNpeaks) return eMean3[i]; else return -1; }

Peak()

float EdbPeak2::Peak ( int  i = 0 ) const

inline

{ if(i>=0&&i<eNpeaks) return ePeak[i];  else return -1; }

Print()

void EdbPeak2::Print

(

)

{
  printf("%d peaks:\n", eNpeaks);
  for(int i=0; i<eNpeaks; i++) 
      printf("%d :  %f %f %f \n", i, ePeak[i], eMean3[i], eMean[i]);
}

ProbPeak() [1/3]

float EdbPeak2::ProbPeak

(

)

{
  int   iv[2];
  FindPeak(iv);
  int   ir[2]={1,1};              // 3x3 neigbouring
  return ProbPeak(iv,ir);
}

ProbPeak() [2/3]

float EdbPeak2::ProbPeak	(	float &	x,
		float &	y
	)

{
  int   iv[2];
  FindPeak(iv);
  x = X(iv[0]);
  y = Y(iv[1]);
  int   ir[2]={1,1};              // 3x3 neigbouring
  return ProbPeak(iv,ir);
}

ProbPeak() [3/3]

float EdbPeak2::ProbPeak	(	int	iv[2],
		int	ir[2]
	)

{
  float prob=0;
  int   npeak      = Bin(iv);
  int   nbinPeak   = (2*ir[0]+1)*(2*ir[1]+1);
  float meanNeib   = 1.*(Integral(iv,ir) - npeak)/(nbinPeak-1);
  float meanNoPeak = 1.*(Integral() - Integral(iv,ir))/(Ncell()-nbinPeak);
  Log(3,"ProbPeak","found at (%3d %3d): %4d  %6.3f  %6.3f", iv[0],iv[1],npeak, meanNeib,meanNoPeak);
  ePeak[eNpeaks]  = npeak;
  eMean3[eNpeaks] = meanNeib;
  eMean[eNpeaks]  = meanNoPeak;
  prob = npeak - meanNoPeak;                       // ToDo: this is not normalised value
  eNpeaks++;
  return prob;
}

ProbPeaks() [1/2]

int EdbPeak2::ProbPeaks

(

int

npeak

)

{
  int ir[2]={1,1};
  return ProbPeaks(npeaks,ir);
}

ProbPeaks() [2/2]

int EdbPeak2::ProbPeaks	(	int	npeak,
		int	ir[2]
	)

{
  int   iv[2], ic=0;
  for(int i=0; i<npeaks; i++) {
    FindPeak(iv);
    ProbPeak(iv,ir);
    WipePeak(iv,ir);
    ic++;
  }
  return ic;
}

Smooth()

float EdbPeak2::Smooth

(

Option_t *

option = "k5a"

)

Smooth bin contents of this 2-d histogram using kernel algorithms
similar to the ones used in the raster graphics community.
Bin contents in the active range are replaced by their smooth values.
This code is copied from root_5.22 mainly to not depend on the newest root versions
Adopted to fedra by VT
to speedup processing do all calculations as integer, so the histogram do not normalised at output
return value - norm-factor of the kernel

{
 
   if(NX()<1||NY()<1)  return 1.;
 
   Int_t k5a[5][5] =  { { 0, 0, 1, 0, 0 }, 
            { 0, 2, 2, 2, 0 }, 
            { 1, 2, 5, 2, 1 }, 
            { 0, 2, 2, 2, 0 }, 
            { 0, 0, 1, 0, 0 } };   // norm = 25
   Int_t k5b[5][5] =  { { 0, 1, 2, 1, 0 },
            { 1, 2, 4, 2, 1 },
            { 2, 4, 8, 4, 2 },
            { 1, 2, 4, 2, 1 },
            { 0, 1, 2, 1, 0 } };
   Int_t k3a[3][3] =  { { 0, 1, 0 },
            { 1, 2, 1 },
            { 0, 1, 0 } };      // norm = 6
   Int_t k3b[3][3] =  { { 1, 2, 1 },        //gauss smoothing
            { 2, 4, 2 },
            { 1, 2, 1 } };      // norm = 16
 
   Int_t ksize_x, ksize_y;
   TString opt = option;
   opt.ToLower();
   Int_t *kernel = &k5a[0][0];
   if (opt.Contains("k5a")) {
     kernel = &k5a[0][0];
     ksize_x = ksize_y = 5;
   }
   else if (opt.Contains("k5b")) {
     kernel = &k5b[0][0];
     ksize_x = ksize_y = 5;
   }
   else if (opt.Contains("k3a")) {
      kernel = &k3a[0][0];
      ksize_x = ksize_y = 3;
   } 
   else if (opt.Contains("k3b")) {
      kernel = &k3b[0][0];
      ksize_x = ksize_y = 3;
   } 
   else return 1.;
 
   // Kernel tail sizes (kernel sizes must be odd for this to work!) 
   Int_t x_push = (ksize_x-1)/2;
   Int_t y_push = (ksize_y-1)/2; 
 
   Int_t norm=0;
   for(int i=0; i<ksize_x; i++)
     for(int j=0; j<ksize_y; j++) norm += kernel[i*ksize_y +j];
 
   int nx = NX(), ny= NY();
 
   EdbH2 buf( *((EdbH2*)this) );
 
   Long_t content;
   int bin=0, k=0;
 
   for (int i=0; i<nx; i++){             // input hist
     for (int j=0; j<=ny; j++) { 
       
       content = 0;
       for (int n=0; n<ksize_x; n++) {          // kernel
     for (int m=0; m<ksize_y; m++) { 
       Int_t xb = i+(n-x_push);
       Int_t yb = j+(m-y_push); 
       bin = buf.Bin(xb,yb);                      // out of dimensions - return 0;
       if (!bin) continue;
       k = kernel[n*ksize_y +m];
       if (!k) continue;
       content += k*bin;
     }
       }
       //SetBin( i, j, Nint(1.*content/norm) );
       SetBin( i, j, content );
     } 
   } 
 
   eNorm = (float)norm; 
   return eNorm;
}

WipePeak()

int EdbPeak2::WipePeak	(	int	iv[2],
		int	ir[2]
	)

{
  int mean = (int)Mean();
  int nbin=0;
  for(int ix=iv[0]-ir[0]; ix<=iv[0]+ir[0]; ix++) 
    for(int iy=iv[1]-ir[1]; iy<=iv[1]+ir[1]; iy++) 
      if(Jcell(ix,iy)>-1)
    {eNC[Jcell(ix,iy)] = mean; nbin++;}
  return nbin;
}

Xmean()

float EdbPeak2::Xmean

(

)

{
  Double_t mean=0, sum=0;
  for(int ix=0; ix<eN[0]; ix++) {
    float x = X(ix);
    for(int iy=0; iy<eN[1]; iy++) {
      int j = Jcell(ix,iy);
      mean += x*eNC[j];
      sum  += eNC[j];
    }
  }
  mean /= sum;
  return mean;
}

Ymean()

float EdbPeak2::Ymean

(

)

{
  Double_t mean=0, sum=0;
  for(int iy=0; iy<eN[1]; iy++) {
    float y = Y(iy);
    for(int ix=0; ix<eN[0]; ix++) {
      int j = Jcell(ix,iy);
      mean += y*eNC[j];
      sum  += eNC[j];
    }
  }
  mean /= sum;
  return mean;
}

Member Data Documentation

eMean

TArrayF EdbPeak2::eMean

eMean3

TArrayF EdbPeak2::eMean3

eNorm

Float_t EdbPeak2::eNorm

the norm-factor in case of the smoothing applied

eNpeaks

Int_t EdbPeak2::eNpeaks

number of found peaks

ePeak

TArrayF EdbPeak2::ePeak

---

The documentation for this class was generated from the following files:

• /home/antonio/fedra_doxygen/src/libEmath/EdbCell2.h
• /home/antonio/fedra_doxygen/src/libEmath/EdbCell2.cpp