[RF] Range not considered when computing integral of RooParamHistFunc

[lmoneta]

The range is not used in RooParamHIstFunc::analyticalIntegralWN, see https://root.cern/doc/master/RooParamHistFunc_8cxx_source.html#l00247
as reported in https://root-forum.cern.ch/t/createintegral-giving-wrong-results/43508.

Simple code to reproduce:

void testParamHistFunc() {                                                                                                                                
                                                                                                                                                          
   RooRealVar x("x","x",0,10);                                                                                                                            
   TH1D h1("h1","h1",20,0,10);                                                                                                                            
   TF1 f1("f1","gaus"); f1.SetParameters(1,5,1);                                                                                                          
   h1.FillRandom("f1",1000);                                                                                                                              
                                                                                                                                                          
   RooDataHist dh("dh","dh",x,&h1);                                                                                                                       
   RooParamHistFunc phf("phf","",x,dh);                                                                                                                   
   phf.Print("t");                                                                                                                                        
   x.setRange("R1",0,5);                                                                                                                                  
                                                                                                                                                                                                                
   auto ig1n = phf.createIntegral(x,x);                                                                                                                   
   std::cout << ig1n->getVal() << std::endl;                                                                                                              
   auto igr1 = phf.createIntegral(x,x,"R1");                                                                                                              
   std::cout << igr1->getVal() << std::endl;                                                                                                              
                                                                                                            
}        

[lmoneta]

The same problem exists also for the RooHIstFunc.
It has been reported already here
https://sft.its.cern.ch/jira/browse/ROOT-7413?filter=-2
using createRunningIntegral with a possible solution that needs to be investigated

[guitargeek]

While working on this issue, I noticed that also RooHistPdf is also not bug free. Here is what you can do to get strange results with it:

// g++ $(root-config --cflags) -o testHistPdf testHistPdf.cc $(root-config --libs) -lRooFitCore -lRooFit -g

#include "RooRealVar.h"
#include "RooHistPdf.h"
#include "RooDataHist.h"
#include "TH2D.h"
#include "TF2.h"

using namespace RooFit;

int main(int argc, char const *argv[]) {

   RooRealVar x("x","x",0, 10);
   RooRealVar y("y","y",0.05);
   TH2D h2("h2","h2",20,0,10, 30, 0, 10);
   TF2 f2("f2","y < 0.1");
   h2.FillRandom("f2",1000);

   RooArgSet argSet{x, y};

   RooDataHist dh("dh","dh",argSet,&h2);
   RooHistPdf phf("phf","",argSet,dh);
   x.setRange("R1",0,5);
   y.setRange("R1",0,10);
   auto int1 = phf.createIntegral(argSet,argSet);
   std::cout << int1->getVal() << std::endl;
   auto int2 = phf.createIntegral(argSet,argSet, "R1");
   std::cout << int2->getVal() << std::endl;

   auto int3 = phf.createIntegral(x,x);
   std::cout << int3->getVal() << std::endl;
   auto int4 = phf.createIntegral(x,x,"R1");
   std::cout << int4->getVal() << std::endl;

  return 0;
}

Obviously the integral of a constant function if you only take half of the x range should be half of the value you get for the full range, even if there is an additional variable in the slice set. However, this is what the program above outputs in ROOT master:

1
0.493
1
0.164333

The first 3 integrals are correct, but the final one (custom range for x and slice of y) gives the wrong result (0.493 expected).

I thought about how the integration should be done correctly for all of the RooHistPdf, RooHistFunc and RooParamHistFunc family. Some updates to the code to fix all the issues can be found in PR #7271.

[guitargeek]

The RooParamHistFunc story continues.

It was noticed that the range is still not considered in the integral when you clone the integral just like in this example:

void testRooParamHistFunc2() {
   using namespace RooFit;

   RooRealVar x("x","x",0,10);
   x.setRange("R1",0,5);
   TF1 f1("f1","1");

   TH1D h1("h1","h1",10,0,10);
   h1.FillRandom("f1",50);
   TH1D h2("h2","h2",10,0,10);
   h2.FillRandom("f1",50);

   RooDataHist dh1("dh1","dh1",x,&h1);
   RooDataHist dh2("dh2","dh2",x,&h2);

   RooParamHistFunc ph("ph","",x,dh1);

   RooUniform uni("uni", "uni", RooArgList(x));
   RooRealVar frac("frac", "frac", 0.5, 0.0, 1.0);
   RooRealSumPdf model{"model", "model", ph, uni, frac};

   auto fitResult = ph.fitTo(dh2, PrintLevel(0), Save());

   auto integral = model.createIntegral(x,x,"R1");
   auto integralClone = static_cast<RooAbsReal*>(integral->cloneTree());

   cout << integral->getVal(x) << endl;
   cout << integralClone->getVal(x) << endl;
}

So this issue can't be closed yet.

[guitargeek]

Fixed in master by #7198, #7271, and #7478.

Fixed in 6.24 by the #7198, #7448, and #7479. The fist PR is the same as for master because that was before the 6.24 branch was created, and the other two PRs are backports.