2D Lorentzian Fit

[Tawfik_Gaballah]

Hello, I’m trying to make a 2D lorentzian fit of a TH2D histogram, but the fit doesn’t seem to be right. Any idea how I can improve it?
I have attached the code that I’m using here

#include <iostream>
#include <iomanip>
#include <TF2.h>
#include <vector>

double Lorentzian2Dr(double *x, double *par) {
  double r = sqrt(pow((x[0] - par[2]), 2.) + pow((x[1] - par[3]), 2.)); 
	
  return (0.5*par[0]*par[1]/TMath::Pi()) /
    TMath::Max( 1.e-10, r*r
		+ .25*par[1]*par[1]);


	
}



void Graph2D(){




  TH2D *gxyz = new TH2D("","",30,0,1500,30,0,1500);
  TFile* myfile = new TFile("output_CO_2D_2_test.root", "READ");
  TFile* proj = new TFile("2DGraph.root", "RECREATE");
  TH3D* hxyz=(TH3D*)myfile->Get("LaBr3_Tdiff_3D_0_8");
  hxyz->Write();		 
  Int_t nx = hxyz->GetNbinsX();
  Int_t ny = hxyz->GetNbinsY();
  Int_t nz = hxyz->GetNbinsZ();
  double Constant=1000;
  double MeanX= hxyz->GetMean(2);
  double SigmaX= hxyz->GetStdDev(2);
  double MeanY= hxyz->GetMean(3);
  double SigmaY= hxyz->GetStdDev(3);



  Int_t nn = 0;
  Double_t x, y, z, nw, nevts;
  Double_t avgz = 0.0;
  for (Int_t j = 1; j <= nx; j++) { 
    x = hxyz->GetXaxis()->GetBinCenter(j);   
    for (Int_t i = 1; i <= ny; i++) {
      y = hxyz->GetYaxis()->GetBinCenter(i);
      double nevts = 0;
      double avgz = 0.0;
      for (Int_t k = 1; k <= nz; k++) {
	z = hxyz->GetZaxis()->GetBinCenter(k);
	nw = hxyz->GetBinContent(j,i,k);
	avgz += (z*nw);
	nevts += nw;
      }
      if (avgz != 0) {
	gxyz->Fill(x,y);
	gxyz->SetBinContent(j,i,avgz/nevts);
      }
    } 
  }  			 
  TCanvas *c = new TCanvas("c","c",1000,600);

  double maxResponse= gxyz->GetMaximumBin();
  double Max= gxyz->GetBinContent(maxResponse);
  double xMax= gxyz->GetXaxis()->GetBinCenter(maxResponse);

  double yMax= gxyz->GetYaxis()->GetBinCenter(maxResponse);
  TF2 *fLorentz = 0;  
  double par[] = {Max, -2, xMax, yMax}; 
  fLorentz = new TF2("fLorentz",Lorentzian2Dr,
		     0, 1500, 
		     0, 1500,4
		     );
  fLorentz->SetParameters(par); 
  // fLorentz->FixParameter(2, xMax); 
  // fLorentz->FixParameter(3, yMax);
				
  int fitResult = gxyz->Fit("fLorentz"); 				
  fLorentz->GetParameters(par); 

  gxyz->Fit("fLorentz");
		       	  
  gxyz->Draw("lego2");
  proj->Close();
  myfile->Close();







}

[Eddy_Offermann]

Things I notice:

1. at y=0 and 1400 you have many entries with 0, exclude them from the fit

2. par[0] controls the amplitude of the Lorentzian but remember that the
   maximum of your function is par[0]*2/pi/par[1], so par[0]=max might
   not be a good starting value.

3. if you supply the data, others can check better

[Tawfik_Gaballah]

Thanks for the reply

1. I set the bin error=0 for y=0, y=1400 to avoid including them in the fit. Using this condition:

gxyz[l][m]->Fill(x,y);
if(y==0 || y==1400){
  gxyz[l][m]->SetBinError(j,i,0);}

2. I added this relation to calculate par [0]

double  A= Max*TMath::Pi()/-2*2;
  double par[] = {A, -2, xMax, yMax}; 

_Here is how the fit looks now:

3. Here is the root file for reference:
   output_CO_2D_test.root (153.2 KB)

[Eddy_Offermann]

The input file you added contains a branch called “LaBr3_Tdiff_3D_0_9”
instead of “LaBr3_Tdiff_3D_0_8”. So my results will probably differ.

The main problem in your code is that you do a “Fill” and a “SetBinContent”.
As long as you are filling the same cell, the last statement overwrites the
first one, but in many cases that is not happening and therefore you are
creating lots of bins with content 1.

I believe the following code will give you a good fit:

#include <iostream>
#include <iomanip>
#include <TF2.h>
#include <vector>

double Lorentzian2Dr(double *x, double *par) {
  double r = sqrt(pow((x[0] - par[2]), 2.) + pow((x[1] - par[3]), 2.));

  return (0.5*par[0]*par[1]/TMath::Pi()) /
    TMath::Max( 1.e-10, r*r
                + .25*par[1]*par[1]);
}

void Graph2D(){

  TH2D *gxyz = new TH2D("","",50,0,1400,50,0,1400);
  TFile* myfile = new TFile("output_CO_2D_test.root", "READ");
  TFile* proj = new TFile("2DGraph.root", "RECREATE");
  TH3D* hxyz=(TH3D*)myfile->Get("LaBr3_Tdiff_3D_0_9");
  hxyz->Write();

  Int_t nx = hxyz->GetNbinsX();
  Int_t ny = hxyz->GetNbinsY();
  Int_t nz = hxyz->GetNbinsZ();

  double MeanX= hxyz->GetMean(2);
  double SigmaX= hxyz->GetStdDev(2);
  double MeanY= hxyz->GetMean(3);
  double SigmaY= hxyz->GetStdDev(3);

  Int_t nn = 0;
  for (Int_t j = 1; j <= nx; j++) {
    Double_t x = hxyz->GetXaxis()->GetBinCenter(j);
    for (Int_t i = 1; i <= ny; i++) {
      Double_t y = hxyz->GetYaxis()->GetBinCenter(i);
      Double_t nevts = 0;
      Double_t avgz = 0.0;
      for (Int_t k = 1; k <= nz; k++) {
        Double_t z = hxyz->GetZaxis()->GetBinCenter(k);
        Double_t nw = hxyz->GetBinContent(j,i,k);
        avgz += (z*nw);
        nevts += nw;
      }
      if (avgz != 0) {
        gxyz->SetBinContent(j,i,avgz/nevts);
      } else {
        gxyz->SetBinContent(j,i,0.);
      }
    }
  }

  TCanvas *c = new TCanvas("c","c",1000,600);

  double maxResponse= gxyz->GetMaximumBin();
  double Max= gxyz->GetBinContent(maxResponse);
  double xMax= gxyz->GetXaxis()->GetBinCenter(maxResponse);

  double yMax= gxyz->GetYaxis()->GetBinCenter(maxResponse);
  TF2 *fLorentz = 0;
  double par[] = {-Max, 50, xMax, yMax};
  fLorentz = new TF2("fLorentz",Lorentzian2Dr,
                     0, 1500,
                     0, 1500,4
                     );
  fLorentz->SetParameters(par);
  // fLorentz->FixParameter(2, xMax); 
  // fLorentz->FixParameter(3, yMax);

  int fitResult = gxyz->Fit("fLorentz");
  fLorentz->GetParameters(par);

  gxyz->Fit("fLorentz");

  gxyz->Draw("lego2");
  proj->Close();
  myfile->Close();
}

Another consideration is the error in the bin, currently they are all equal .

[Tawfik_Gaballah]

I used this code and the fits look much better now, thank you so much.