Fitting a Landau Convoluted Gaussian + Gaussian with RooFit

Roofit and RooStats
1

I’ve been trying to fit a landau convoluted gaussian plus another gaussian to a histogram using RooFit. The best fit I can achieve is attached. I’ve tried adjusting all the parameters and cannot get a fully accurate fit. It seems there are “tipping points” for some of these parameters, where if the parameter is changed past the point, the fit will change dramatically. Any suggestions?

Here is my code:

[code]#ifndef CINT
#include “RooGlobalFunc.h”
#endif
#include “RooRealVar.h”
#include “RooDataSet.h”
#include “RooGaussian.h”
#include “RooLandau.h”
#include “RooFFTConvPdf.h”
#include “RooPlot.h”
#include “TCanvas.h”
#include “TAxis.h”
#include "TH1.h"
using namespace RooFit;

void gl()
{

  RooRealVar t("t", "t", -30, 100);


  RooRealVar ml("ml", "mean landau", 15, -30, 100);
  RooRealVar sl("sl", "sigma landau", 0.8, 0.1, 10);
  RooLandau landau("lx", "lx", t, ml, sl);


  RooRealVar mg("mg","mg", 10, -30, 100);
  RooRealVar sg("sg","sg", 2.3, 0.1, 10);
  RooGaussian gauss("gauss", "gauss", t, mg, sg);


  RooRealVar mg2("mg2","mg2",5, -30, 100);
  RooRealVar sg2("sg2", "sg2",5, 0.1, 10);
  RooGaussian gauss2("gauss2", "gauss2", t, mg, sg);


  t.setBins(10000, "cache");


  TFile* f = TFile::Open("HGC_514_Reco_CM.root");
  f->cd("hgcaltbrechitsplotter_highgain_correlation_cm");
  TH1F* h =(TH1F*)f->FindObjectAny("Ski_1_Channel_26");
  RooDataHist* data = new RooDataHist("data", "data", t, Import(*h));


  RooFFTConvPdf lxg("lxg", "landau (X) gauss", t, landau, gauss);
  RooRealVar coeff("coeff", "coeff", .5, 0.01, 1);
  RooAddPdf lxgpg("lxgpg", "Landau Convoluted Gaussian + Gaussian", lxg, gauss2, coeff);

  lxgpg.fitTo(*data);


  RooPlot* frame = t.frame(Title("landau (x) gauss convolution"));


  data->plotOn(frame);
  lxgpg.plotOn(frame, LineColor(kBlue));


  TCanvas* c1 = new TCanvas("gausslandau", "gausslandau", 1000, 1000);
  gPad->SetLeftMargin(0.15);
  frame->GetYaxis()->SetTitleOffset(1.4);
  frame->Draw();


  c1->Print("Plots/Ski_1_Channel_26.pdf");

}
[/code]

– Ryan

2

Hi,

It is possible that a small change of the initial parameter can cause the fit to fail because you are starting in a region too far way from the minimum. This is not surprising and the effect depends on the model.
If you have the impression something is wrong, please attach the ROOT file so we can reproduce the problem.

Lorenzo

3

I ended up getting it to fit by using a Landau convoluted Crystal Ball instead of the Landau convoluted Gaussian.

However, when I fit the function, I’d like to plot the final function as well as both the convoluted function and the added gaussian. I’m able to plot all of them, but it seems the scale is incorrect for the landau convoluted crystal ball and the added gaussian. I’ve attached the resulting plot. Blue is the final fit, dashed red is the gaussian, and dashed orange is the landau convoluted crystal ball.

Is there a way to “shrink” these functions so that they seem to fit the curve?

Here is my code in case something needs to be changed:

[code]#ifndef CINT
#include “RooGlobalFunc.h”
#endif
#include “RooRealVar.h”
#include “RooDataSet.h”
#include “RooGaussian.h”
#include “RooLandau.h”
#include “RooCBShape.h”
#include “RooFFTConvPdf.h”
#include “RooPlot.h”
#include “TCanvas.h”
#include “TAxis.h”
#include “TH1.h”
#include
using namespace RooFit;

void gausscrystalfit()
{
//setting filenames to open
std::string FileName(“ryantest8308”);
std::string DirectoryName(“hgcaltbrechitsplotter_ryan”);
std::string HistName(“Ski_2_Channel_5_ADC1”);

//declaring variables and functions
RooRealVar t("t", "t", -20, 80);
t.setBins(10000, "cache");

RooRealVar ml("ml", "mean landau", 10, -100, 100);
RooRealVar sl("sl", "sigma landau", 2, 0.1, 10);
RooLandau landau("lx", "lx", t, ml, sl);

RooRealVar mc("mc","mc", 5, -100, 100);
RooRealVar sc("sc","sc", 2, 0.1, 10);
RooRealVar alpha("alpha", "alpha", 1, 0, 100);
RooRealVar nn("nn", "nn", 1, 0, 100);
RooCBShape crystal("crystal", "crystal", t, mc, sc, alpha, nn);

RooRealVar mg("mg","mg",0, -1000, 100);
RooRealVar sg("sg", "sg", 0.85, 0.1, 5);
RooGaussian gauss("gauss", "gauss", t, mg, sg);



//opening files and transferring histo contents to RooDataSet
std::stringstream fname;
fname << "Files/" << FileName << ".root";
TFile* f = TFile::Open(fname.str().c_str());
f->cd(DirectoryName.c_str());
TH1F* h =(TH1F*)f->FindObjectAny(HistName.c_str());
RooDataHist* data = new RooDataHist("data", "data", t, Import(*h));



//convoluting and adding functions
RooFFTConvPdf lxc("lxc", "Landau (X) CBall", t, landau, crystal);
RooRealVar coeff("coeff", "coeff", .5, 0.01, 1);
RooAddPdf lxcpg("lxcpg", "Landau Convoluted Crystal Ball + Gaussian", lxc, gauss, coeff);



//plotting and fitting functions to the data
RooPlot* frame = t.frame(Title("ADC Count"));
data->plotOn(frame);
lxcpg.fitTo(*data);
gauss.plotOn(frame, LineStyle(kDashed), LineColor(kRed));
lxc.plotOn(frame, LineStyle(kDashed), LineColor(kOrange + 7));
lxcpg.plotOn(frame, LineColor(kBlue));



//integrating functions and calculating MPV
t.setRange("integral", -10, 50);
RooAbsReal* Int = lxcpg.createIntegral(t, Range("integral"));
RooAbsReal* lxcInt = lxc.createIntegral(t, Range("integral"));
RooAbsReal* gaussInt = gauss.createIntegral(t, Range("integral"));
double LxCInt = lxcInt->getVal();
double GaussInt = gaussInt->getVal();
double FullInt = Int->getVal();
cout << "Integral of LxC = " << LxCInt << endl;
cout << "Integral of gauss = " << GaussInt << endl;
cout << "Full Integral = " << FullInt << endl;

RooAbsReal* dldt = lxc.derivative(t) ;
Double_t tpeak ;
Bool_t ok = RooBrentRootFinder(RooRealBinding(*dldt, t)).findRoot(tpeak, 0, 100, 0) ;
cout << "MPV = " << tpeak << endl;



//editing canvas and drawing it
TCanvas* c1 = new TCanvas("gausscrystal", "gausscrystal", 1000, 1000);
gPad->SetLeftMargin(0.15);
frame->GetYaxis()->SetTitleOffset(1.2);
frame->GetYaxis()->SetTitle("Frequency");
frame->GetXaxis()->SetTitle("ADC Count");
std::stringstream CanvasName;
CanvasName << HistName << " :: Frequency by ADC Count";
frame->SetTitle(CanvasName.str().c_str());
frame->Draw();

//adding text to show MPV
std::stringstream MPVText;
MPVText << "MPV = " << tpeak;
c1->Update();
double x = gPad->GetUxmax() - gPad->GetUxmin();
double y = gPad->GetUymax() - gPad->GetUymin();
TPaveText *text = new TPaveText(x * .8 + gPad->GetUxmin(),y * .9 + gPad->GetUymin(),x + gPad->GetUxmin(),y + gPad->GetUymin());
text->AddText(MPVText.str().c_str());
text->SetBorderSize(1);
text->SetFillColor(19);
text->Draw();

//printing to a PDF file
std::stringstream PlotName;
PlotName << "Plots/" << FileName << "_" << HistName << ".pdf";
c1->Print(PlotName.str().c_str());

}[/code]

Thanks,

4

Hi,

You should use the Componets(gauss) when plotting the Gaussian component of your pdf and not use plotOn on the single component pdf, otherwise the normalisation will not be correct.

See root.cern.ch/doc/master/rf201__ … te_8C.html

Best Regards

Lorenzo