RooFormula bad numerical expression

weishi · December 18, 2020, 12:40am

Dear experts,

I was modeling my background shape pdf called HighMassFit2018_m1:

RooRealVar m1_high_mass("m1_high_mass", "m_{#mu#mu_{1}}", m_SR3_min, m_SR3_max, "GeV");
RooGenericPdf HighMassFit2018_m1("HighMassFit2018_m1", "2018 m1 estimated BKG shape at SR3", "Limit_cfg::My_BKGShapem1(m1_high_mass)", RooArgSet(m1_high_mass) );

The “Limit_cfg::My_BKGShapem1(m1_high_mass)” is a single-valued function/formula that uses simple interpolation between many mass points. The interpolation is coded this way:

double MCBinCenterMass[14] = {12.75, 16.25, 19.75, 23.25, 26.75, 30.25, 33.75, 37.25, 40.75, 44.25, 47.75, 51.25, 54.75, 58.25};
double MCBinContentm1[14] = {0.0260007, 0.070103, 0.457687, 0.161783, 0.157428, 0.62272, 0.599198, 0.887516, 1.56432, 2.07059, 0.707749, 1.12067, 2.29358, 0.954618};

namespace Limit_cfg {
double My_BKGShapem1(double m1) {
double mybkg = 0.0;

if ( m1 >= 11.000              && m1 < MCBinCenterMass[0] )   mybkg = MCBinContentm1[0];//use flat profile of bin content for the first half bin
if ( m1 >= MCBinCenterMass[0]  && m1 < MCBinCenterMass[1] )   mybkg = MCBinContentm1[0]  + ( m1 - MCBinCenterMass[0] )*( MCBinContentm1[1] - MCBinContentm1[0] )/( MCBinCenterMass[1] - MCBinCenterMass[0] );
if ( m1 >= MCBinCenterMass[1]  && m1 < MCBinCenterMass[2] )   mybkg = MCBinContentm1[1]  + ( m1 - MCBinCenterMass[1] )*( MCBinContentm1[2] - MCBinContentm1[1] )/( MCBinCenterMass[2] - MCBinCenterMass[1] );
if ( m1 >= MCBinCenterMass[2]  && m1 < MCBinCenterMass[3] )   mybkg = MCBinContentm1[2]  + ( m1 - MCBinCenterMass[2] )*( MCBinContentm1[3] - MCBinContentm1[2] )/( MCBinCenterMass[3] - MCBinCenterMass[2] );
if ( m1 >= MCBinCenterMass[3]  && m1 < MCBinCenterMass[4] )   mybkg = MCBinContentm1[3]  + ( m1 - MCBinCenterMass[3] )*( MCBinContentm1[4] - MCBinContentm1[3] )/( MCBinCenterMass[4] - MCBinCenterMass[3] );
if ( m1 >= MCBinCenterMass[4]  && m1 < MCBinCenterMass[5] )   mybkg = MCBinContentm1[4]  + ( m1 - MCBinCenterMass[4] )*( MCBinContentm1[5] - MCBinContentm1[4] )/( MCBinCenterMass[5] - MCBinCenterMass[4] );
if ( m1 >= MCBinCenterMass[5]  && m1 < MCBinCenterMass[6] )   mybkg = MCBinContentm1[5]  + ( m1 - MCBinCenterMass[5] )*( MCBinContentm1[6] - MCBinContentm1[5] )/( MCBinCenterMass[6] - MCBinCenterMass[5] );
if ( m1 >= MCBinCenterMass[6]  && m1 < MCBinCenterMass[7] )   mybkg = MCBinContentm1[6]  + ( m1 - MCBinCenterMass[6] )*( MCBinContentm1[7] - MCBinContentm1[6] )/( MCBinCenterMass[7] - MCBinCenterMass[6] );
if ( m1 >= MCBinCenterMass[7]  && m1 < MCBinCenterMass[8] )   mybkg = MCBinContentm1[7]  + ( m1 - MCBinCenterMass[7] )*( MCBinContentm1[8] - MCBinContentm1[7] )/( MCBinCenterMass[8] - MCBinCenterMass[7] );
if ( m1 >= MCBinCenterMass[8]  && m1 < MCBinCenterMass[9] )   mybkg = MCBinContentm1[8]  + ( m1 - MCBinCenterMass[8] )*( MCBinContentm1[9] - MCBinContentm1[8] )/( MCBinCenterMass[9] - MCBinCenterMass[8] );
if ( m1 >= MCBinCenterMass[9]  && m1 < MCBinCenterMass[10]  ) mybkg = MCBinContentm1[9]  + ( m1 - MCBinCenterMass[9] )*( MCBinContentm1[10] - MCBinContentm1[9] )/( MCBinCenterMass[10] - MCBinCenterMass[9] );
if ( m1 >= MCBinCenterMass[10]  && m1 < MCBinCenterMass[11] ) mybkg = MCBinContentm1[10] + ( m1 - MCBinCenterMass[10] )*( MCBinContentm1[11] - MCBinContentm1[10] )/( MCBinCenterMass[11] - MCBinCenterMass[10] );
if ( m1 >= MCBinCenterMass[11]  && m1 < MCBinCenterMass[12] ) mybkg = MCBinContentm1[11] + ( m1 - MCBinCenterMass[11] )*( MCBinContentm1[12] - MCBinContentm1[11] )/( MCBinCenterMass[12] - MCBinCenterMass[11] );
if ( m1 >= MCBinCenterMass[12]  && m1 < MCBinCenterMass[13] ) mybkg = MCBinContentm1[12] + ( m1 - MCBinCenterMass[12] )*( MCBinContentm1[13] - MCBinContentm1[12] )/( MCBinCenterMass[13] - MCBinCenterMass[12] );
if ( m1 >= MCBinCenterMass[13]  && m1 <= 60.000             ) mybkg = MCBinContentm1[13];//use flat profile of bin content for the last half bin

return mybkg;}
}

However, ROOT seems to be reporting error for this formula:

Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem1(m1_high_mass)"

I used similar interpolation functions in other places, those don’t seem to cause any issue. Could you help me indentify what happens here? I can provide whole script if you want.

Cheers,
Wei

_ROOT Version: 6.12/07
_Platform: CentOS Linux release 7.8.2003 (Core)
_Compiler: gcc version 7.3.1 20180127 (GCC)

etejedor · December 18, 2020, 11:50am

@moneta can perhaps give a hand? Thanks!

moneta · December 18, 2020, 2:17pm

Hi,
Difficult to help with this limited information. We would need all the code to reproduce this.
But the function Limit_cfg::My_BKGShapem1(m1_high_mass) is available when used it ?
It has ben loaded so it is recognised by Cling ?

Lorenzo

weishi · December 18, 2020, 3:29pm

Hi Lorenzo,

I think the function is available when used. But I pushed the code here on GitHub: https://github.com/weishi10141993/MuJetAnalysis_LimitSetting/blob/test/macros/makeWorkSpace_H2A4Mu.C#L363-L364

The interpolation function is defined in the Config.h in the main directory: https://github.com/weishi10141993/MuJetAnalysis_LimitSetting/blob/test/Config.h#L163-L189

It does not report any error when I simply run the “makeWorkSpace_H2A4Mu.C” script to make workspace to store the interpolated background shapes for a mass point, things seem to be ok:

Processing makeWorkSpace_H2A4Mu.C(53.0000)...
...
[#1] INFO:ObjectHandling -- RooWorkspace::import(w_H2A4Mu) importing RooGenericPdf::HighMassFit2018_m1
[#1] INFO:ObjectHandling -- RooWorkspace::import(w_H2A4Mu) importing RooGenericPdf::HighMassFit2018_m2
...
p.d.f.s
-------
RooProdPdf::HighMassBKG[ HighMassFit2018_m1 * HighMassFit2018_m2 ] = 0.610897
RooGenericPdf::HighMassFit2018_m1[ actualVars=(m1_high_mass) formula="Limit_cfg::My_BKGShapem1(m1_high_mass)" ] = 0.743357
RooGenericPdf::HighMassFit2018_m2[ actualVars=(m2_high_mass) formula="Limit_cfg::My_BKGShapem2(m2_high_mass)" ] = 0.821809
...

The compile error appears when I use a combine command to evaluate the shape (the combine command is used for upper limit calculation, so maybe not too relevant here, I think it evaluates the stored background shape in the workspace):

combine -n .H2A4Mu_mA_53.0000_GeV_LHC_T30000_0 -m 125 -M HybridNew --saveHybridResult --expectedFromGrid 0.500 --rule CLs --testStat LHC --cl 0.95 --rAbsAcc 0.01 --rRelAcc 0.001 -s -1 -T 30000 --fork 50 Datacards/2018/datacard_H2A4Mu_mA_53.0000_GeV.txt > macros/sh/2018/output/output_0.500_53.0000_T30000_0.txt
Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem1(m1_high_mass)"
Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem2(m2_high_mass)"

I hope this is helpful for you?

Cheers,
Wei

moneta · December 18, 2020, 5:38pm

HI,
Thank you for the explanations. I think the problem is that when using in combine you don’t have the definition of that function available.
If you would be within root, I would tell you to just include the file, in combine, since it is a separate application, I am not sure what is the best way to do it. Probably you would need to add a line like
gROOT->ProcessLine(".L Config.h") in the combine code before using the function.

Lorenzo

weishi · December 18, 2020, 6:09pm

Hi Lorenzo,

Thanks for the advice. I’ll try it.
But this is also where I am a bit confused as I also use another similar interpolation function called “Limit_cfg::My_MassWindow(m1, m2)” in exactly the same config file and same combine command:

github.com

weishi10141993/MuJetAnalysis_LimitSetting/blob/test/macros/makeWorkSpace_H2A4Mu.C#L138-L140


RooGenericPdf dia1( "dia1", "generic PDF for diagonal region at SR1", "fabs(m1_below_Jpsi - m2_below_Jpsi) < Limit_cfg::My_MassWindow(m1_below_Jpsi, m2_below_Jpsi)", RooArgSet(m1_below_Jpsi,m2_below_Jpsi) );
RooGenericPdf dia2( "dia2", "generic PDF for diaginal region at SR2", "fabs(m1_above_Jpsi - m2_above_Jpsi) < Limit_cfg::My_MassWindow(m1_above_Jpsi, m2_above_Jpsi)", RooArgSet(m1_above_Jpsi,m2_above_Jpsi) );
RooGenericPdf dia3( "dia3", "generic PDF for diaginal region at SR3", "fabs(m1_high_mass  - m2_high_mass)  < Limit_cfg::My_MassWindow(m1_high_mass, m2_high_mass)", RooArgSet(m1_high_mass,m2_high_mass) );

github.com

weishi10141993/MuJetAnalysis_LimitSetting/blob/test/macros/makeWorkSpace_H2A4Mu.C#L322


  w_H2A4Mu->import(signal2016_m2);
  w_H2A4Mu->factory("PROD::signal1(signal2016_m1, signal2016_m2)");
  w_H2A4Mu->factory("PROD::signal2(signal2016_m1, signal2016_m2)");
}
TFile* file = new TFile(inputFile1); // defined in ../Config.h
RooWorkspace *w = (RooWorkspace*) file->Get("w");
if(year==2017 || year==2018){
  w_H2A4Mu->import( *w->pdf("template1D_m1_below_Jpsi") );
  w_H2A4Mu->import( *w->pdf("template1D_m2_below_Jpsi") );
  w_H2A4Mu->factory("PROD::BBbar_below_Jpsi_2D(template1D_m1_below_Jpsi, template1D_m2_below_Jpsi)*dia1");
  w_H2A4Mu->import( *w->pdf("template1D_m1_above_Jpsi") );
  w_H2A4Mu->import( *w->pdf("template1D_m2_above_Jpsi") );
  w_H2A4Mu->factory("PROD::BBbar_above_Jpsi_2D(template1D_m1_above_Jpsi, template1D_m2_above_Jpsi)*dia2");
}
if(year==2016){
  w_H2A4Mu->import( *w->pdf("template1D_m1") );
  w_H2A4Mu->import( *w->pdf("template1D_m2") );
  w_H2A4Mu->factory("PROD::BBbar_below_Jpsi_2D(template1D_m1, template1D_m2)*dia2016");
  w_H2A4Mu->factory("PROD::BBbar_above_Jpsi_2D(template1D_m1, template1D_m2)*dia2016");//keep same name as 2017 and 2018
}

I don’t see any error reported on this function. I don’t get what’s the difference between these two functions.

Wei

weishi · December 18, 2020, 7:06pm

I just tried using the function Limit_cfg::My_MassWindow(m1, m2) to replace the Limit_cfg::My_BKGShapem1(m1_high_mass) in combine, it also reports similar compile error:

Error in <RooFormula::Compile>: Bad numerical expression : "Limit_cfg::My_MassWindow(m1_high_mass,m1_high_mass)"
[#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m1): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_MassWindow(m1_high_mass, m1_high_mass)
Error in <RooFormula::Compile>: Bad numerical expression : "Limit_cfg::My_MassWindow(m2_high_mass,m2_high_mass)"
[#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m2): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m2): Formula doesn't compile: Limit_cfg::My_MassWindow(m2_high_mass, m2_high_mass)

So this is not related to the function itself. Apparently it’s complaining this formula expression in RooGenericPdf:

"Limit_cfg::My_MassWindow(m1_high_mass,m1_high_mass)",

but this expression is recogonized in RooGenericPdf:

fabs(m1_high_mass  - m2_high_mass)  < Limit_cfg::My_MassWindow(m1_high_mass, m2_high_mass)

How should I resolve this then if I want to just use “Limit_cfg::My_BKGShapem1(m1_high_mass)” in RooGenericPdf?

Wei

moneta · December 20, 2020, 2:04pm

Hi,
Have you tried to include the function definition before using it, as I suggested using

gROOT->ProcessLine(".L Config.h");

(sorry not gSystem->ProcessLine(…) as I said in the previous post that I have corrected now)

Lorenzo

weishi · December 20, 2020, 3:58pm

Hi Lorenzo,

I did try your suggestion.

Since the combine command line is issued outside ROOT, so I did this:

(Parrot) [ws13@terra3 MuJetAnalysis_LimitSetting]$ echo 'gROOT->ProcessLine(".L Constants.h"); gROOT->ProcessLine(".L Config.h")' | root -l -b
(long) 0
(Parrot) [ws13@terra3 MuJetAnalysis_LimitSetting]$ combine -n .H2A4Mu_mA_53.0000_GeV_LHC_T30000_0 -m 125 -M HybridNew --saveHybridResult --expectedFromGrid 0.500 --rule CLs --testStat LHC --cl 0.95 --rAbsAcc 0.01 --rRelAcc 0.001 -s -1 -T 30000 --fork 50 Datacards/2018/datacard_H2A4Mu_mA_53.0000_GeV.txt > macros/sh/2018/output/output_0.500_53.0000_T30000_0.txt

Error in <RooFormula::Compile>: Bad numerical expression : "Limit_cfg::My_BKGShapem1(m1_high_mass)"
Error in <RooFormula::Compile>: Bad numerical expression : "Limit_cfg::My_BKGShapem2(m2_high_mass)"

It doesn’t seem to work.
As you see, when I run gROOT->ProcessLine(".L Config.h"), it returned (long) 0. Is this what you expected if it ran successfully?

Should I hack the combine application somehow to add the function in the combine source code?
Also, do you have any idea why it works well for the other similar function defined in a same fashion I reported in the previous thread?

RooGenericPdf dia3( "dia3", "generic PDF for diaginal region at SR3", "fabs(m1_high_mass - m2_high_mass) < Limit_cfg::My_MassWindow(m1_high_mass, m2_high_mass)", RooArgSet(m1_high_mass,m2_high_mass) );

Thanks,
Wei

moneta · December 21, 2020, 10:16am

Hi,

There is actually another better way to have your custom function code available in combine.
When you create and write your workspace in a file, you can add this line of code:

w_H2A4Mu->importClassCode("Limit_cfg");

This will import the code for the class (our namespace Limit_cfg) in the workspace and make it available for further use, e.g. when reading the workspace in combine

Best regards

Lorenzo

weishi · December 21, 2020, 4:47pm

Hi Lorenzo,

I tried importClassCode:

w_H2A4Mu->importClassCode("Limit_cfg");
RooRealVar testvar("testvar", "testvar", Limit_cfg::My_BKGShapem1(mA_GeV));
cout<<"testvar:" << Limit_cfg::My_BKGShapem1(mA_GeV) <<endl;
w_H2A4Mu->import(testvar);

Printing the value of the function is fine. Still, things don’t work in combine. I added debug mode for combine, here is what it returns:

...
Will work with unbinned datasets
Observables: ['m1_high_mass,m2_high_mass']
Will use category 'CMS_channel' to identify the 1 channels
Creating pdfs for individual modes (1): Creating RooAddPdf pdf_binA with 2 elements
. done.
Importing combined pdf model_s
Importing combined pdf model_b

RooWorkspace(w) w contents

variables
---------(CMS_H2A4Mu_BKG_norm,CMS_H2A4Mu_BKG_norm_In,CMS_H2A4Mu_BKG_syst,CMS_H2A4Mu_BKG_syst_In,CMS_H2A4Mu_effdimu_mass,CMS_H2A4Mu_effdimu_mass_In,CMS_H2A4Mu_nnlo_pt,CMS_H2A4Mu_nnlo_pt_In,CMS_channel,CMS_eff_mu_hlt,CMS_eff_mu_hlt_In,CMS_eff_mu_id,CMS_eff_mu_id_In,CMS_eff_mu_iso,CMS_eff_mu_iso_In,CMS_eff_mu_pu_eff,CMS_eff_mu_pu_eff_In,MH,QCDscale_ggH,QCDscale_ggH_In,Xsec_BR_decay,Xsec_BR_decay_In,lumi_13TeV,lumi_13TeV_In,m1_high_mass,m2_high_mass,r,signal3_alpha,signal3_mA,signal3_n,signal3_sigma)

p.d.f.s
-------
Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem1(m1_high_mass)"
RooSimultaneousOpt::model_b[ indexCat=CMS_channel A=pdf_binA_bonly extraConstraints=() channelMasks=() ] = [#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m1): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem2(m2_high_mass)"
[#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m2): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m2): Formula doesn't compile: Limit_cfg::My_BKGShapem2(m2_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
0
...

I also tried remove the prefix class name Limit_cfg::,

RooGenericPdf HighMassFit2018_m1("HighMassFit2018_m1", "2018 m1 estimated BKG shape at SR3", "My_BKGShapem1(m1_high_mass)", RooArgSet(m1_high_mass) );

This gives me compile error even when the workspace file was created before used in combine:

Error in <RooFormula::Compile>: Bad numerical expression : "My_BKGShapem1(m1_high_mass)"
RooProdPdf::HighMassBKG[ HighMassFit2018_m1 * HighMassFit2018_m2 ] = [#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m1): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: My_BKGShapem1(m1_high_mass)
Error in <RooFormula::Compile>: Bad numerical expression : "My_BKGShapem2(m2_high_mass)"
[#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m2): compile error

Wei

weishi · December 21, 2020, 8:17pm

Hi Lorenzo,

I did a check on the workspace file outside combine by opening it and printout the contents of the workspace:

root -l ws_H2A4Mu_mA_57.0000_GeV.root
w_H2A4Mu->Print()

It seems the function is not correctly written to the workspace in the first place even before its use in combine:

p.d.f.s
-------
Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem1(m1_high_mass)"
RooProdPdf::HighMassBKG[ HighMassFit2018_m1 * HighMassFit2018_m2 ] = [#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m1): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
Error in <RooFormula::Compile>:  Bad numerical expression : "Limit_cfg::My_BKGShapem2(m2_high_mass)"
[#0] ERROR:InputArguments -- RooFormula::RooFormula(HighMassFit2018_m2): compile error
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m2): Formula doesn't compile: Limit_cfg::My_BKGShapem2(m2_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
0
RooGenericPdf::HighMassFit2018_m1[ actualVars=(m1_high_mass) formula="Limit_cfg::My_BKGShapem1(m1_high_mass)" ] = [#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m1): Formula doesn't compile: Limit_cfg::My_BKGShapem1(m1_high_mass)
0
RooGenericPdf::HighMassFit2018_m2[ actualVars=(m2_high_mass) formula="Limit_cfg::My_BKGShapem2(m2_high_mass)" ] = [#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m2): Formula doesn't compile: Limit_cfg::My_BKGShapem2(m2_high_mass)
[#0] ERROR:Eval -- RooFormula::eval(HighMassFit2018_m2): Formula doesn't compile: Limit_cfg::My_BKGShapem2(m2_high_mass)
0

So I wonder is there an example of how to write a custom function for its use as a pdf and correctly write it to the workspace? I saw some discussions on using a custom pdf class, but didn’t seem to find a complete example.

Wei

moneta · December 22, 2020, 11:50am

Hi,

I think then the only solution is to make the class or function definition available in combine. You would need to hack to add the gROOT->ProcessLine function.
I don’t know if combine can offer some other possibility to load custom code. It is maybe a question to ask to the combine authors

Best regards

Lorenzo

system · January 5, 2021, 11:50am

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