Convoluting Crystal balll + Gaussian

SimN · February 25, 2019, 2:54pm

Hello,

I’m trying to fit a histogram with a convoluted (crystal ball + Gaussian ) for the signal part and Chebychev polynomial for the background part.

The idea is to try to fit the signal part in the data with the convoluted Monte carlo shape (crystal ball) and a gaussian.

The fit is failing. you can see the code down below, could you please let me know what I’m doing wrong here?

  RooRealVar x("x","MM_{(#pi #pi)} GeV/c^{2}",3.04,3.17) ;
  RooDataHist dh("dh","dataset",x,h1);
  RooPlot* frame = x.frame(Title("DATA : convoluted Sig(crystal ball + gaussian) + back(pol) ")) ;
  dh.plotOn(frame,Name("dh")) ;

  //Gaussian 
 RooRealVar mean1("#mu","mean of gaussians",3.0968,3.0985);
  RooRealVar sigma1("#sigma","width of gaussians",0.0025,0.0055);
  RooGaussian sig1("sig1", "sig1", x, mean1, sigma1);
//Crystal Ball with fixed parameters from MC 
  RooRealVar mean2("#mu","mean of gaussians",3.096401);
  RooRealVar sigma2("#sigma","width of gaussians",0.005126);
  RooRealVar alpha2("#alpha","alpha",1.780);
  RooRealVar n2("n2","n",3.27);
  RooCBShape sig2("sig2", "sig2", x, mean2, sigma2, alpha2, n2);

//Convolution 
  RooFFTConvPdf sigg("sigg","gauss",x,sig1,sig2) ;
  RooRealVar a0("a0","a0",-2.,2.) ;

//Background
  RooChebychev bkg("bkg","background p.d.f.",x,RooArgList(a0));
  RooRealVar nsig("N_{SIG}","signal events",0,1000000);
  RooRealVar nbkg("N_{BKG}","signal background even0ts",0,100000000);
  
 //Total
  RooAddPdf all("all","model",RooArgList(sigg,bkg),RooArgList(nsig,nbkg));
  RooFitResult* r = all.fitTo(dh,Extended(kTRUE),Save()) ;
  all.paramOn(frame,Layout(0.5,0.90,0.55));
  all.plotOn(frame,Components(bkg),LineStyle(kDashed));
  all.plotOn(frame,Name("all"));
  r->Print();
  frame->GetXaxis()->SetTitle("");
  frame->GetXaxis()->SetLabelSize(1);
  frame->GetYaxis()->SetLabelSize(0.05);
  frame->SetMinimum(20000);
  frame->Draw();

StephanH · February 26, 2019, 8:45am

It’s always helpful to have the output as well. I can’t reproduce anything from your code because there’s no data.

SimN · February 26, 2019, 9:15am

I attached the data root file.

4230_original.root (4.4 KB)

and this is the full code :

#include "RooFFTConvPdf.h"
#include "RooAbsCachedPdf.h"
#include "TH1.h"
#include "TMath.h"
#include "TF1.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TChain.h"
#include <TROOT.h>
#include <TStyle.h>
#include <TLegend.h>
#include "TROOT.h"
#include "TClassRef.h"
#include "TVirtualPad.h"
#include "TBrowser.h"
#include "TPad.h"
#include "TString.h"
#include "TMath.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TGaxis.h"
#include "TLine.h"
#include "TVirtualFitter.h"
#include "TFitResult.h"
#include "THStack.h"


using namespace std;
using namespace RooFit;
using namespace RooAbsPdf;


TCanvas *c1 = new TCanvas("c1","Jpsi",1200,1200);
TCanvas *c2 = new TCanvas("c2","Jpsi",1200,1200);

void conv_4230(){

  TFile *f = new TFile("4230_original.root");
  TH1F * h1 = new TH1F("h1","hist4230",150,3.04,3.17);
  h1 = (TH1F*)f->Get("hist4230");

  RooRealVar x("x","MM_{(#pi #pi)} GeV/c^{2}",3.04,3.17) ;
  RooDataHist dh("dh","dataset",x,h1);
  RooPlot* frame = x.frame(Title("DATA : convoluted Sig(crystal ball + gaussian) + back(pol) ")) ;
  dh.plotOn(frame,Name("dh")) ;
  RooRealVar mean1("#mu","mean of gaussians",3.0968,3.0985);
  RooRealVar sigma1("#sigma","width of gaussians",0.0025,0.0055);

  RooGaussian sig1("sig1", "sig1", x, mean1, sigma1);

  RooRealVar mean2("#mu","mean of gaussians",3.096401);
  RooRealVar sigma2("#sigma","width of gaussians",0.005126);
  RooRealVar alpha2("#alpha","alpha",1.780);
  RooRealVar n2("n2","n",3.27);


  RooCBShape sig2("sig2", "sig2", x, mean2, sigma2, alpha2, n2);

  RooFFTConvPdf sigg("sigg","gauss",x,sig1,sig2) ;
  RooRealVar a0("a0","a0",-2.,2.) ;

  RooChebychev bkg("bkg","background p.d.f.",x,RooArgList(a0));
  RooRealVar nsig("N_{SIG}","signal events",0,1000000);
  RooRealVar nbkg("N_{BKG}","signal background even0ts",0,100000000);
  RooAddPdf all("all","model",RooArgList(sigg,bkg),RooArgList(nsig,nbkg));
  RooFitResult* r = all.fitTo(dh,Extended(kTRUE),Save()) ;
  all.paramOn(frame,Layout(0.5,0.90,0.55));
  all.plotOn(frame,Components(bkg),LineStyle(kDashed));
  all.plotOn(frame,Name("all"));
  r->Print();
  frame->GetXaxis()->SetTitle("");
  frame->GetXaxis()->SetLabelSize(1);
  frame->GetYaxis()->SetLabelSize(0.05);
  frame->SetMinimum(20000);
  frame->Draw();

}

StephanH · March 1, 2019, 9:47am

Do you remember that I asked for the output? When I run your example, this is what I get:

[#0] WARNING:Minization -- RooMinimizerFcn: Minimized function has error status.
Returning maximum FCN so far (-6.8699e+07) to force MIGRAD to back out of this region. Error log follows
Parameter values: #mu=3.09841, #sigma=0.00250845, N_{BKG}=2.95444e+07, N_{SIG}=608058, a0=1.83447
RooChebychev::bkg[ x=x coefList=(a0) ]
     p.d.f value is less than zero (-0.822243), forcing value to zero @ x=x=3.04043, coefList=(a0 = 1.83447)
[...]

The problem is that your background is negative. I fiddled a bit with your code to get you started with a model that is always positive. From here, you have to try to get it working like you wanted to have it.

void conv_4230(){

  TFile *f = new TFile("4230_original.root");
  TH1F * h1 = new TH1F("h1","hist4230",150,3.04,3.17);
  h1 = (TH1F*)f->Get("hist4230");

  RooRealVar x("x","MM_{(#pi #pi)} GeV/c^{2}",3.04,3.17) ;
  RooDataHist dh("dh","dataset",x,h1);
  RooPlot* frame = x.frame(Title("DATA : convoluted Sig(crystal ball + gaussian) + back(pol) ")) ;
  
  RooRealVar mean1("mu","mean of gaussians",3.0968,3.0985);
  RooRealVar sigma1("sigma","width of gaussians",0.0025,0.0055);

  RooGaussian sig1("sig1", "sig1", x, mean1, sigma1);

  RooRealVar mean2("mu2","mean of gaussians",3.096401);
  RooRealVar sigma2("sigma2","width of gaussians",0.005126);
  RooRealVar alpha2("alpha2","alpha",1.780);
  RooRealVar n2("n2","n",3.27);


  RooCBShape sig2("sig2", "sig2", x, mean2, sigma2, alpha2, n2);

  RooFFTConvPdf sigg("sigg","gauss",x,sig1,sig2) ;
  RooRealVar a0("a0","a0",-10000, 10000);
  RooRealVar a1("a1","a1",-100, 100);

  RooPolynomial bkg("bkg","background p.d.f.",x,RooArgList(a0, a1));
  RooRealVar nsig("N_{SIG}","signal events",0,1000000);
  RooRealVar nbkg("N_{BKG}","signal background even0ts",1E5,1E7);
  RooAddPdf all("all","model",RooArgList(sig1,bkg),RooArgList(nsig,nbkg));
  RooFitResult* r = all.fitTo(dh,Extended(kTRUE),Save()) ;
  dh.plotOn(frame,Name("dh")) ;
  all.paramOn(frame,Layout(0.5,0.90,0.55));
  all.plotOn(frame,Components(bkg),LineStyle(kDashed));
  all.plotOn(frame,Name("all"));
  r->Print("v");
  frame->GetXaxis()->SetTitle("");
  frame->GetXaxis()->SetLabelSize(1);
  frame->GetYaxis()->SetLabelSize(0.05);
  frame->SetMinimum(20000);
  frame->Draw();

}

Watch out for the following things:

Keep the PDFs positive. A PDF cannot be negative. RooFit will tell you if you run into this problem.
Check that the coefficients are not at the limit. Look at the output of Minuit to check if this problem occurs, and change the ranges of the coefficients if necessary.
Don’t set extreme ranges for coefficients (i.e. don’t set -1.E9 to 1.E9). Minuit will get into trouble because it uses the range to compute a reasonable step size.

SimN · March 4, 2019, 10:12am

Hi @StephanH,
Thank you very much for your help.
just one more question :
in this line, why do you add to the model sig1 and not sigg ?

RooAddPdf all("all","model",RooArgList(sig1,bkg),RooArgList(nsig,nbkg));

StephanH · March 4, 2019, 11:06am

Because convolutions complicate the whole thing. I needed something to get it to work fast, and a Gaussian distribution seemed to be fine for the signal distribution.
If you really need the convolution, you will have to try finding good parameters.

system · March 18, 2019, 11:06am

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.