Allow multiple ranges in HFitInterface::FillData

[belmakier]

This Pull request:

Allows multiple (discontinuous) ranges for HFitInterface::FillData, see here

@lmoneta - one thing I didn't implement was checking if the different ranges overlapped at all. If several of the ranges are overlapping (i.e. [0,50], [25, 75]), then at present the bins in both regions will be added twice, not sure what sort of behavior this would cause. I think the best way to get around that is to check when the range are added, I note that one such check is already in place in DataRange::CleanRangeSet, this could be extended. Wanted to flag with you first before I changed too much though.

Example for 1D use:

  TF1 *back = new TF1("back", "[0] + x*[1]", -10, 10);
  TF1 *comb = new TF1("comb", "[0] + x*[1] + [2]*exp(-((x-[3])**2)/[4])", -10, 10);
  
  TH1F *th1 = new TH1F("th1", "th1", 100, -10, 10);  

  comb->SetParameter(0, 10.);
  comb->SetParameter(1, -0.5);
  comb->SetParameter(2, 15.0);
  comb->SetParameter(3, 1.0);
  comb->SetParameter(4, 1.5);

  TCanvas *c1 = new TCanvas();
  th1->FillRandom("comb", 5000);
  th1->Draw("colz");

  ROOT::Fit::DataRange range;
  //x ranges to miss the peak
  range.AddRange(0, -10, -2.5);
  range.AddRange(0, 4.5, 10);
  
  ROOT::Fit::DataOptions opt;
  ROOT::Fit::BinData data(opt, range);
  ROOT::Fit::FillData(data, th1);
  
  back->SetParameter(0, 10.);
  back->SetParameter(1, -0.5);
  
  ROOT::Math::WrappedMultiTF1 fitFunc(*back, back->GetNdim());
  ROOT::Fit::Fitter fitter;
  fitter.SetFunction(fitFunc, false);
  fitter.LikelihoodFit(data);
  fitter.Result().Print(std::cout);

  th1->Draw("hist");
  back->Draw("same");

or 2D use:

  TF2 *back = new TF2("back", "[0] + (x-[1])**2/[2]**2 + (y-[3])**2/[4]**2", -10, 10, -10, 10);
  TF2 *comb = new TF2("comb", "[0] + (x-[1])**2/[2]**2 + (y-[3])**2/[4]**2 + [5]*TMath::Gaus(x,[6],[7])*TMath::Gaus(y,[8],[9])", -10, 10, -10, 10);
  
  TH2F *th2 = new TH2F("th2", "th2", 100, -10, 10, 100, -10, 10);  

  comb->SetParameter(0, 3.);
  comb->SetParameter(1, 2.5);
  comb->SetParameter(2, 6.);
  comb->SetParameter(3, -1.0);
  comb->SetParameter(4, 3.5);

  comb->SetParameter(5, 50);
  comb->SetParameter(6, -0.5);
  comb->SetParameter(7, 0.6);
  comb->SetParameter(8, 0.2);
  comb->SetParameter(9, 0.5);

  //TF2 *comb = new TF2("comb", "back+gaus", -10, 10, -10, 10);

  TCanvas *c1 = new TCanvas();
  th2->FillRandom("comb", 1000000);
  th2->Draw("colz");

  ROOT::Fit::DataRange range;
  //x ranges to miss the peak
  range.AddRange(0, -10, -2.5); 
  range.AddRange(0, 1.5, 10);
  
  range.AddRange(1, -10, 10);
  
  ROOT::Fit::DataOptions opt;
  ROOT::Fit::BinData data(opt, range);
  ROOT::Fit::FillData(data, th2);
  
  back->SetParameter(0, 2.);
  back->SetParameter(1, 1.5);
  back->SetParameter(2, 8.);
  back->SetParameter(3, -0.5);
  back->SetParameter(4, 4.0);
  
  ROOT::Math::WrappedMultiTF1 fitFunc(*back, back->GetNdim());
  ROOT::Fit::Fitter fitter;
  fitter.SetFunction(fitFunc, false);
  fitter.LikelihoodFit(data);
  fitter.Result().Print(std::cout);

  TCanvas *c2 = new TCanvas();
  back->Draw("colz");
  back->SetMinimum(0);
  back->SetMaximum(th2->GetMaximum());

Checklist:

• tested changes locally
• updated the docs (if necessary)

This PR fixes #

[phsft-bot]

Can one of the admins verify this patch?

[lmoneta]

Thank you for this nice contribution.
I think it is important adding the case of overlapping ranges.
You might extend the Range functionality to do this. When adding an overlapping range CleanRangeSet should remove un-needed range and return the minimal set.
I think now it is working deleting a previous existing range if it is smaller than the one provided.

[lmoneta]

I think you don't need to call ExamineRange multiple times. You can call it later when looping on the histogram bins and then get the number of points for each range. In that case you don't call BinData::Initialize(n, dim) but you call BinData::Append(newPoints, dim) for each examined range.

[belmakier]

Ok, that makes sense. I've updated the flow - is that closer to what you had in mind?

[lmoneta]

Yes, this is exactly what I asked for !

[belmakier]

I've updated CleanRangeSet() to handle cases where the new range partially overlaps with an existing one. You can test with the following:

  ROOT::Fit::DataRange range;
  std::cout << "Start with [0,5]" << std::endl;
  range.AddRange(0,5);
  std::cout << "Size = 1 = " << range.Size() << std::endl;
  std::cout << "range [0,5] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;
                               
  std::cout << "Add [2,3]" << std::endl;
  range.AddRange(2,3);
  
  assert(range.Size() == 1);
  std::cout << "Size = 1 = " << range.Size() << std::endl;

  assert((range(0,0).first == 0));
  assert(range(0,0).second == 5);

  std::cout << "range [0,5] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;

  std::cout << "Add [-1,6]" << std::endl;
  range.AddRange(-1,6);
  assert(range.Size() == 1);
  std::cout << "Size = 1 = " << range.Size() << std::endl;

  assert(range(0,0).first == -1);
  assert(range(0,0).second == 6);
  std::cout << "range [-1,6] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;

  std::cout << "Add [-2,4]" << std::endl;
  range.AddRange(-2,4);
  assert(range.Size() == 1);
  std::cout << "Size = 1 = " << range.Size() << std::endl;

  assert(range(0,0).first == -2);
  assert(range(0,0).second == 6);
  std::cout << "range [-2,6] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;

  std::cout << "Add [5,7]" << std::endl;
  range.AddRange(5,7);
  assert(range.Size() == 1);
  std::cout << "Size = 1 = " << range.Size() << std::endl;

  assert(range(0,0).first == -2);
  assert(range(0,0).second == 7);
  std::cout << "range [-2,7] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;

  std::cout << "Add [20,25]" << std::endl;
  range.AddRange(20,25);
  assert(range.Size() == 2);
  std::cout << "Size = 2 = " << range.Size() << std::endl;

  assert(range(0,0).first == -2);
  assert(range(0,0).second == 7);
  std::cout << "range [-2,7] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;


  assert(range(0,1).first == 20);
  assert(range(0,1).second == 25);
  std::cout << "range [20,25] = [" << range(0,1).first << "," << range(0,1).second << "]" << std::endl;

  std::cout << "Add [24,26]" << std::endl;
  range.AddRange(24,26);
  assert(range.Size() == 2);
  std::cout << "Size = 2 = " << range.Size() << std::endl;
  assert(range(0,1).first == 20);
  assert(range(0,1).second == 26);
  std::cout << "range [20,26] = [" << range(0,1).first << "," << range(0,1).second << "]" << std::endl;

  std::cout << "Add [19,20]" << std::endl;
  range.AddRange(19,20);
  assert(range(0,1).first == 19);
  assert(range(0,1).second == 26);
  std::cout << "range [16,26] = [" << range(0,1).first << "," << range(0,1).second << "]" << std::endl;

  std::cout << "Add [6,20]" << std::endl;
  range.AddRange(6,20);
  assert(range.Size() == 1);
  std::cout << "Size = 1 = " << range.Size() << std::endl;
  assert(range(0,0).first == -2);
  assert(range(0,0).second == 26);
  std::cout << "range [-2,26] = [" << range(0,0).first << "," << range(0,0).second << "]" << std::endl;

This shows several behaviors: adding a range fully encapsulated by an existing one, adding a range which fully encapsulates and existing one, extending an existing range by adding an overlapping one, merging two existing ranges by adding one in the middle.

Tim

[lmoneta]

It is great you add the support for overlapping range in DataRange.
Can you also please convert the test you have above in a simple Google test that can be added in math/mathcore/test ?
You can just replace the assert with EXPECT_EQ.
A simple example to look is this test:
https://github.com/root-project/root/blob/master/hist/hist/test/test_TH1.cxx

Thank you for this contribution!

[belmakier]

Hi @lmoneta, is this ok now?

[guitargeek]

Thank you very much @belmakier! This is a very nice PR that should be merged, as these multi-range fits are quite common, e.g. for excluding the signal region.

Unfortunately, our CI can't deal with PRs where the PR branch is also called master. So I have re-created a PR with your changes:

• [hist][math] Allow multiple ranges in HFitInterface::FillData #22861

Other changes I did besides rebasing:

• Add more tests
• Avoid indentation changes for now, for easier review
• Some fixes for corner cases, like empty data handling etc