Strange behavior trying to use Conditional Observables?

Roofit and RooStats
1

I’m attempting to implement a simultaneous fit of the mass and lifetime of j/psis as a precursor to fitting the j/psi polarization accounting for the non-prompt fraction in the polarization distributions. I use various parts of the (mass,c-tau) plane to fit regions enriched in prompt, nonprompt and background simultaneously to get the mass/lifetime distribution of prompt j/psis.

To describe the lifetime distribution I’m using a gaussian with per-event scaling of the standard deviation, requiring that the scaling factor be a conditional parameter for the prompt, non-prompt and background lifetime descriptions.

For some reason, when performing the convolution in the RooDecay and even for the RooGaussModel itself, RooFit attempts to integrate over the c-tau error. As far as I know this is wrong behavior and I have tried many workarounds with limited or no success. At best the fit converges only sometimes, and right now I have it setup such that the RooGaussModel forces a flat integral over the conditional parameter, but this isn’t the right way to do it either. No matter what I try, on some level RooFit always tries to integrate over the conditional observable, and I am completely stumped.

You can find the code here: cmssw.cvs.cern.ch/cgi-bin/cmssw. … iew=markup

Thanks in advance!

2

Hi,

I had a look at your code, and at first sight I don’t see anything wrong with it. However,
since I cannot run it (no data), I can’t redirectly reproduce your problem.

Instead I have written a small test program that builds a pdf of similar complexity
(although w/o simultaneous fit) to see if I can reproduce what you and and I cannot.

I insert my test code here:

void jpsi_test()
{
RooWorkspace w(“w”,1) ;

// Basic model with conditional observable
w.factory(“Decay::decay(t[0,20],tau[1.547,1,2],RooGaussModel::rmodel(t,bias[0],resol[1,0.9,1.2],terr[0.1,20]),RooDecay::SingleSided)”) ;
w.factory(“Uniform::dummy(terr)”) ;
w.factory(“PROD::model(dummy,decay|terr)”) ;

// Now add dummy flat term and construct extended pdf with sum
w.factory(“Decay::decay2(t,tau2[20],RooGaussModel::rmodel2(t,bias2[0],resol2[1,0.9,1.2],terr),RooDecay::SingleSided)”) ;
w.factory(“Uniform::dummy2(terr)”) ;
w.factory(“PROD::model2(dummy2,decay2|terr)”) ;

w.factory(“SUM::model3(N1[100,0,1000]*model,N2[100,0,1000]*model2)”) ;

RooDataSet* d = w::model3.generate(RooArgSet(w::t,w::terr),1000) ;

w.Print(“t”) ;

//RooMsgService::instance().addStream(DEBUG,Topic(Eval),ClassName(“RooProdPdf”)) ;

cout << “now fitting model” << endl ;
w::model3.fitTo(*d) ;
cout << “now done” << endl ;

}

To see what happens in terms of normalization you should uncomment the line ‘RooMsgService::instance()’ which will then show for each call to RooProdPdf what normalization it assumes for the call to each of its components. In my case it just lists ‘t’, which is correct for the conditional construction made here.

What I suggest is that you

  • Run this macro and see if it works OK for you
  • Run your own code with the ‘RooMsgService…’ line prepended and see what output that produces.

Also if my macro works OK for you, could you extend following your own macro to the point where it starts failing and then send it back to me? That would allow me to have a closer look.

Wouter

3

Hi Wouter,

Thanks for your reply, I’ve taken your code, run it and experimented with my own. I can now get the conditional observables to be handled correctly.

I also found the reason for the instability of my fits. If I fit the mass and lifetime distributions separately & simultaneously the fit converges without fail. This should be equivalent to what I was doing before, i.e. fitting the product of the mass and lifetime distributions, since in the log this forms two separate sums over each PDF. In looking into why I get a convergent fit when I manually separate the mass and lifetime pdfs and then fit them simultaneously, what I understand is that when RooNLLVar calls getLogVal() on a RooProdPdf it actually evaluates the log of the product instead of the sum of the logs of factorizable terms in the product. This can be highly numerically unstable and is a property shared with all classes derived from RooAbsPdf since they all use the same getLogVal() function. In addition, this leads to the rather interesting case that when you do getLogVal() for a RooGaussian, you evaluate log(norm*e^(-(x-mu)^2/(2s^2)) instead of the much more numerically accurate and less computationally expensive -(x-mu)^2/(2s^2)*log(norm). If this hasn’t been addressed already in newer versions of RooFit, I think it may be wise to spend some time developing less expensive and more accurate implementations of getLogVal() for each derived Pdf class.

Please let me know what you think.

Thanks,
Lindsey

4

Hi,

I believe that I am suffering from a very similar problem to this one. I am performing a 2D mass-lifetime fit where the decay time resolution width is conditional on the per-event error (multiplied by a common scale factor). I use a RooHistPdf to take the distribution of the error from data. For the moment I take the same distribution in signal and background. As with the previous poster, the PDF appears to be integrating over the per-event error. The fit usually converges, but the value returned for the resolution scale factor is definitely wrong. Looking at it’s NLL profile, it is exponential rather than something closer to a parabola.

It appears that I can prevent the integration over the per-event error if I use the ConditionalObservables(*TAUERR) flag to the fitTo() method, but the fit result still looks wrong. Can anyone suggest what I need in order to correctly account for the per-event error. I can provide code examples if necessary.

Best regards,
Greig