Fitting with gaus+landau

ROOT
1

Hi everyone! I would like to fit a histogram with a sum of Gaus and Landau, using the following code:

[code]#include
#include “TFile.h”
#include “TTree.h”
#include “TRandom.h”
#include “TMath.h”
#include “TF1.h”
#include “TH1D.h”

void Fitting(const char *a,double b,double c)
{
TFile f(a);
TCanvas c1 =(TCanvas)f.Get(“c1”);
TF1 *l = new TF1(“l”,“gaus(0)+landau(3)”,-1.,1.);
TH1D h1 = (TH1D)c1->FindObject(“h”);
c1->Draw();
h1->Fit(“l”,"","",b,c);
}
[/code]

(a is the filename, and b and c are the fitting limits). If I write just “gaus”, it works perfectly fine, but “gaus(0)+landau(3)” absolutely refuses to work (same goes for “[0]*TMath::Gaus(x,[1],[2])”, and other variants).

Thanks in advance!

2

If you want to simply ADD both functions, which is not the usual case, your code should work if you set correctly the borders of the fit, specify initial parameters, … (Play first with the Fit Panel).

If you want to CONVOLUTE both functions, which is the usual case for non-symmetric gaussian functions, then continue reading:

For simple applications, try with:
http://root.cern.ch/root/html/tutorials/fit/langaus.C.html

For more robust fit routines, use Roofit (you have to install fftw3 and reinstall/extend ROOT installation):
http://root.cern.ch/drupal/content/roofit

See this example code modified/adapted from the on-line tutorials to read from TH1 objects.
http://root.cern.ch/root/html/tutorials/roofit/index.html

#include "RooRealVar.h"
// #include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooLandau.h"
#include "RooFFTConvPdf.h"
// #include "RooGenericPdf.h"
#include "RooPlot.h"
#include "RooCurve.h"
#include "RooDataHist.h"

//assuming your histogram is the variable "px".
Double_t mean=px->GetXaxis()->GetBinCenter(maxBin);
			Double_t sigma=(rightX-leftX)/2.35;

			// S e t u p   c o m p o n e n t   p d f s 
			// ---------------------------------------

			// Construct observable
			RooRealVar t("t","t",px->GetXaxis()->GetBinLowEdge(1),px->GetXaxis()->GetBinUpEdge(nbins));

			// Define fit range
			t.setRange((TString)"ROI_"+i2s(i),minX,maxX);
	
			// Construct gauss(t,mg,sg)
			RooRealVar mg("mg","mg",0) ;
			RooRealVar sg("sg","sg",sigma,0.1*sigma,5.*sigma) ;
			RooGaussian gauss("gauss","gauss",t,mg,sg) ;

			// Construct landau(t,ml,sl) ;
			RooRealVar ml("ml","mean landau",mean,mean-sigma,mean+sigma) ;
			RooRealVar sl("sl","sigma landau",0.04,0.,0.2) ;
			RooLandau landau("lx","lx",t,ml,sl) ;

			// C o n s t r u c t   c o n v o l u t i o n   p d f 
			// ---------------------------------------

			// Set #bins to be used for FFT sampling
			t.setBins(5000,"cache") ; 

			// Construct landau (x) gauss
			RooFFTConvPdf lxg("lxg","landau (X) gauss",t,landau,gauss) ;

			// S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f 
			// ----------------------------------------------------------------------

			RooDataHist* data = new RooDataHist("dh","dh",t,Import(*px)) ;

			// Fit gxlx to data
			lxg.fitTo(*data,Range((TString)"ROI_1));
			
			// Plot data, landau pdf, landau (X) gauss pdf
			RooPlot* frame = t.frame(Title((TString)"FitProjection")) ;
			data->plotOn(frame) ;
			RooPlot* fitLine=lxg.plotOn(frame) ;

				lxg.paramOn(frame);
				data->statOn(frame);
			//~ landau.plotOn(frame,LineStyle(kDashed)) ;

			TF1* flxg = lxg.asTF (RooArgList(t)) ;//NOT NORMALIZED
				
			Double_t modX=flxg->GetMaximumX();
			
			RooCurve* fitCurve = (RooCurve*)fitLine->findObject("lxg_Norm[t]_Range[fit_nll_lxg_dh]_NormRange[fit_nll_lxg_dh]");
			Double_t maxpico=fitCurve?fitCurve->getYAxisMax():0; 
			Double_t maxpicoU=flxg->Eval(modX);//unnormalized
			Double_t extrems[2];
			
			extrems[0]=flxg->GetX(maxpicoU/2.0,minX,modX);
			extrems[1]=flxg->GetX(maxpicoU/2.0,modX,maxX);
			
			// Draw frame on canvas
			frame->Draw();
		
			TLine* l;
			l = new TLine(extrems[0],maxpico/2.0,extrems[1],maxpico/2.0);
			l->SetLineWidth(1);
			l->SetLineColor(1);
			l->Draw("same");

			l = new TLine(extrems[0],px->GetMinimum(),extrems[1],px->GetMinimum());
			l->SetLineWidth(1);
			l->SetLineColor(1);
			l->Draw("same");
			
			
			l = new TLine(modX,px->GetMinimum(),modX,maxpico);
			l->SetLineWidth(0.5);
			l->SetLineColor(1);
			l->Draw("same");
			l=NULL;