Smear a 2d histogram

ROOT
1

Hi! I would like to smear a 2D histogram (i.e. smear both axis in the same way), using a gaussian function.
In other words, the two axis are energy signals from a certain detector and I would like to reduce the energy resolution in both axis by convoluting with a guassian of a given sigma, in order to “simulate” what will be the result with detectors with a worse energy resolution

Is there any easy way to do it ?

Andrea

2

Hi,

We don’t have a function to perform this smearing. But it is not difficult to do. You could smear the event using as x the bin centre or a random position within the bin. If there is no correlation between the 2 coordinates you can do each axis independently.

Best Regards

Lorenzo

3

This could do the job. It is brute force and a bit slow, but anyway.

/**
 * \brief A gaussian filter on a 2D histogram. 
 * \param h The histogram to be smeared. It is modified, so copy first the original one if you do not want to lose it!
 * \param sX the sigma of the x-Gaussian
 * \param sY the sigma of the y-Gaussian
 * \note x and y are treated independently
 */
void Gaus2DFilter(TH2F* const h,const Double_t sX,const Double_t sY)
{
	if(!h) return;
	
	TH2F* h2 = (TH2F*)h->Clone();
	h->Reset();
	double content;
	double x, y, w, z;
	const UInt_t nX = h->GetNbinsX();
	const UInt_t nY = h->GetNbinsY();
	UInt_t kmin, kmax, lmin, lmax;
	
	for(UInt_t i=1;i<=nX;++i)
	{
		x = h->GetXaxis()->GetBinCenter(i);
		for(UInt_t j=1;j<=nY;++j)
		{
			y = h->GetYaxis()->GetBinCenter(j);
			
			
			content = 0;
			
			//For bin i,j, we look in the old histogram at all bins in a nearby region and calculate the probability of being smeared in the current one
			
			kmin = h2->GetXaxis()->FindBin(x-3*sX);//For being faster, we reduce the range to 3-sigma. Use kmin = 1 if you want the whole range
			kmax = h2->GetXaxis()->FindBin(x+3*sX);//For being faster, we reduce the range to 3-sigma. Use kmax = nX if you want the whole range
			
			for(UInt_t k=kmin;k<=kmax;++k)
			{
				w = h2->GetXaxis()->GetBinCenter(k);
				lmin = h2->GetYaxis()->FindBin(y-3*sY);//For being faster, we reduce the range to 3-sigma. Use lmin = 1 if you want the whole range
				lmax = h2->GetYaxis()->FindBin(y+3*sY);//For being faster, we reduce the range to 3-sigma. Use lmax = nY if you want the whole range
				
				for(UInt_t l=lmin;l<=lmax;++l)
				{
					z = h2->GetYaxis()->GetBinCenter(l);
					//We approximate the integral of the gaussian function between binLowEdge and binUpEdge by the integrand times the bin width. This bin width and the normalization is divided after the end of the loop. The Error Function could be used instead.
					content+=h2->GetBinContent(k,l)*TMath::Gaus(x,w,sX,false)*TMath::Gaus(y,z,sY,false);
				}
			}
			content=content*h->GetXaxis()->GetBinWidth(i)*h->GetYaxis()->GetBinWidth(j)/2./TMath::Pi()/sX/sY;
			h->SetBinContent(i,j,content);
		}
	}	
}