Fitting a 2-D tilted gaussian with RooFit

[pinchun]

Dear experts,

PART 1
I would like to fit a 2-D gaussian, whose axes are not parallel to the x and y axes.
For example:

So, I tried to use RooFormulaVar to define x1 and y1 as:

RooRealVar theta_rot("theta_rot","theta_rot",3.1415926/5);
RooFormulaVar x1("x1","x1","cos(theta_rot)*(x-5.28)+sin(theta_rot)*y+5.28",RooArgList(theta_rot,x,y));
RooFormulaVar y1("y1","y1","-sin(theta_rot)*(x-5.28)+cos(theta_rot)*y",RooArgList(theta_rot,x,y));

And then I constructed a PDF and generated the toy data.
RooFitRotate_testsig.cc (2.3 KB)

The histogram of the generated toy data (left) and the PDF (right) look consistent:

However, in the X projection (left), I found that the red line does not fit the data points well:

Why does this happen? How to fix it?

---

PART 2
The situation becomes worse when I try to float the PDF parameters like this:

RooRealVar meanx("mean1","mean of gaussian x1",5.28,5.27,5.29) ;
RooRealVar meany("mean2","mean of gaussian y1",0,-0.1,0.1) ;
RooRealVar sigmax("sigmax","width of gaussian x1",0.0025,0.0001,0.01) ;
RooRealVar sigmay("sigmay","width of gaussian y1",0.01,0.001,0.1) ;

And then it continues to print lots of warning messages like:

[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 6.93939e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 6.3406e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 8.49857e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 9.52489e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 9.07818e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 9.17086e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 6.69686e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 4.6005e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 8.2697e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 7.62193e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 6.55296e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 1.02203e+09
[#0] WARNING:Eval -- RooAbsPdf::getLogVal(sigpdf) WARNING: large likelihood value: 3.38952e+08

How to solve this? Thanks!

[StephanH]

Hi @pinchun,

Part 1

@pinchun

However, in the X projection (left), I found that the red line does not fit the data points well:

This looks just like a plotting issue. You can see that the core fits perfectly, but apparently the sampling is just a bit inaccurate in the tail regions. How did you create the projection? Is it a “real” integral over the PDF or is it a sampled one? It could be that you need more sampling points.

A fast way to get the projection accurate is to add
ProjWData(*toysig, true) to plotOn.

Part 2

I don’t know yet. It’s probably the numeric integral that’s going wrong. One way is to switch to the MonteCarlo integrator. Changing integrators is explained in the tutorial rf901:
https://root.cern/doc/master/rf901__numintconfig_8C.html

[pinchun]

Hi @StephanH ,

Thank you for your prompt reply.

For Part 1, after adding ProjWData(*toysig, true) to plotOn for x-frame, the plotting issue is solved:

For Part 2, I tried adding:

RooAbsReal::defaultIntegratorConfig()->method2D().setLabel("RooMCIntegrator");
RooAbsReal::defaultIntegratorConfig()->methodND().setLabel("RooMCIntegrator");

But the problem still remains…
I even tried to modify nIntPerDim, nRefinePerDim, nRefineIter, genType, and samplingMode to different values, but it doesn’t work either.

[StephanH]

Ok, the integrator is apparently not the problem. Usually when this happens, it is a problem of a very small Gaussian vs. a wide fit range. I tried to play a bit with fit ranges and sigmas, but it didn’t change anything. Maybe you want to give that another try?

Can you also test if switching off the rotation by setting theta to zero or even not using the transformed x1, y1 gets rid of the problem?

[pinchun]

After setting theta to zero, the problem still remains.
If I fit x,y rather than x1,y1, then it would work, but what I would like to fit is the tilted gaussian (theta =/= 0) rather than the untilted one (theta = 0). Or are there any other method to fit such a gaussian with theta =/= 0?

[StephanH]

Ok, so the problem is apparently to get the integrals right when the variable transformation is in use. One way I see to get around this is to use a RooGenericPdf instead of the RooProdPdf, and to implement the multiplication of the 2D Gaussians yourself.
Normally, it’s better to use built-in PDFs because the integrals are faster, but since you need the transformation formulae, you can’t use analytic integrals anyway. If you do this manually with a GenericPdf, numeric integrals will be used throughout all steps, but these should get the normalisation right.

[pinchun]

Yes, now it works using RooGenericPdf:

char x1[51] =  "cos(theta_rot)*(x-5.28)+sin(theta_rot)*y+5.28";
char y1[51] = "-sin(theta_rot)*(x-5.28)+cos(theta_rot)*y";

char gaussx[301]; 
sprintf(gaussx,"exp(-0.5*((%s-%s)/%s)*((%s-%s)/%s))/(%s*sqrt(2*3.14159265359))",x1,"meanx","sigmax",x1,"meanx","sigmax","sigmax");

char gaussy[301]; 
sprintf(gaussy,"exp(-0.5*((%s-%s)/%s)*((%s-%s)/%s))/(%s*sqrt(2*3.14159265359))",y1,"meany","sigmay",y1,"meany","sigmay","sigmay");

char gaussxy[1001]; 
sprintf(gaussxy,"%s*%s",gaussx,gaussy);

RooGenericPdf sigpdf("sigpdf","gaussx*gaussy",gaussxy,RooArgList(theta_rot,x,y,meanx,sigmax,meany,sigmay));

and

 **********
 **   14 **MINOS        2500           5
 **********
 FCN=-775711 FROM MINOS     STATUS=SUCCESSFUL     40 CALLS         307 TOTAL
                     EDM=2.14674e-06    STRATEGY= 1      ERROR MATRIX ACCURATE 
  EXT PARAMETER                  PARABOLIC         MINOS ERRORS        
  NO.   NAME      VALUE            ERROR      NEGATIVE      POSITIVE   
   1  meanx        5.28000e+00   7.92335e-06  -7.92895e-06   7.91788e-06
   2  meany        6.36402e-06   3.16054e-05  -3.16050e-05   3.16059e-05
   3  sigmax       2.50558e-03   5.60267e-06  -5.58683e-06   5.61856e-06
   4  sigmay       9.99436e-03   2.23517e-05  -2.23289e-05   2.23745e-05
   5  theta_rot    6.27452e-01   8.45993e-04  -8.46055e-04   8.45945e-04
                               ERR DEF= 0.5

Thanks for your help!