Intersection of two TF1s

[Erdrick]

---

ROOT Version: root/6.14.06
Platform: Red Hat Enterprise Linux Server release 7.7 (Maipo)
Compiler: gcc/8.3.0

---

I am trying to find the intersection of two TF1s. I am following the example from How to get intersection coordinates. The minimal example below shows my problem.

The attached file minimal.root (8.3 KB) contains a single histogram that has been fit with three TF1s. One TF1 is the sum of two Gaussians (black dotted line in the canvas below), while the other two are single Gaussians (red dotted line and green dotted line in the canvas below).

Now I want to find the intersection of the two Gaussians (the red and the green in the canvas above). I use this script (attached as minimal.C (1.4 KB)):

#include <iostream>

#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include "TCanvas.h"

void minimal() {

  TFile* file = TFile::Open("minimal.root", "READ");

  TH1F* hist = (TH1F*) file->Get("run1396_hippo_elastic_scattering_beam_left_hit");
  std::cout << hist->GetName() << std::endl;
  std::cout << hist->GetListOfFunctions()->At(1)->GetName() << std::endl;
  std::cout << hist->GetListOfFunctions()->At(2)->GetName() << std::endl;
  std::cout << hist->GetListOfFunctions()->At(3)->GetName() << std::endl;

  TCanvas* canvas = new TCanvas("canvas", "canvas");
  gPad->SetLogy();

  hist->Draw();

  ((TF1*) hist->GetListOfFunctions()->At(1))->SetLineColor(1);
  ((TF1*) hist->GetListOfFunctions()->At(2))->SetLineColor(2);
  ((TF1*) hist->GetListOfFunctions()->At(3))->SetLineColor(3);

  ((TF1*) hist->GetListOfFunctions()->At(1))->SetLineStyle(2);
  ((TF1*) hist->GetListOfFunctions()->At(2))->SetLineStyle(2);
  ((TF1*) hist->GetListOfFunctions()->At(3))->SetLineStyle(2);

  ((TF1*) hist->GetListOfFunctions()->At(1))->Draw("same");
  ((TF1*) hist->GetListOfFunctions()->At(2))->Draw("same");
  ((TF1*) hist->GetListOfFunctions()->At(3))->Draw("same");

  TF1* component_A = (TF1*) hist->GetListOfFunctions()->At(2);
  TF1* component_B = (TF1*) hist->GetListOfFunctions()->At(3);
  TF1* difference = new TF1("difference",
                            Form("abs(%s-%s)", component_A->GetName(), component_B->GetName()),
                            0, 10000);

}

I then compile with ACLiC and get this error:

root [0] .L minimal.C++g
Info in <TUnixSystem::ACLiC>: creating shared library /path/to/file/./minimal_C.so
Warning in cling::IncrementalParser::CheckABICompatibility():
  Possible C++ standard library mismatch, compiled with __GLIBCXX__ '20170502'
  Extraction of runtime standard library version was: '20190222'
root [1] minimal()
run1396_hippo_elastic_scattering_beam_left_hit
fitting_function_run1396_hippo_elastic_scattering_beam_left_hit
component_A_run1396_hippo_elastic_scattering_beam_left_hit
component_B_run1396_hippo_elastic_scattering_beam_left_hit
input_line_65:2:126: error: expected ')'
Double_t TFormula____id13879319982141213843(){ return TMath::Abs({component_A_run1396_hippo_elastic_scattering_beam_left_hit}-{component_B_run1396_hippo_elastic_scattering_beam_...
                                                                                                                             ^
input_line_65:2:65: note: to match this '('
Double_t TFormula____id13879319982141213843(){ return TMath::Abs({component_A_run1396_hippo_elastic_scattering_beam_left_hit}-{component_B_run1396_hippo_elastic_scattering_beam_...
                                                                ^
Error in <TFormula::PrepareEvalMethod>: Can't compile function TFormula____id13879319982141213843 prototype with arguments 
Error in <TFormula::InputFormulaIntoCling>: Error compiling formula expression in Cling
Error in <TFormula::ProcessFormula>: "component_A_run1396_hippo_elastic_scattering_beam_left_hit" has not been matched in the formula expression
Error in <TFormula::ProcessFormula>: "component_B_run1396_hippo_elastic_scattering_beam_left_hit" has not been matched in the formula expression
Error in <TFormula::ProcessFormula>: Formula "abs(component_A_run1396_hippo_elastic_scattering_beam_left_hit-component_B_run1396_hippo_elastic_scattering_beam_left_hit)" is invalid !

I do not understand the error messages. Any help would be appreciated.

Thanks

[couet]

Yes the example cross.C is very similar to what you are doing and is working.
I tried to see what are the definitions of the functions by printing their titles:

  printf("%s\n",component_A->GetTitle());
  printf("%s\n",component_B->GetTitle());

It gives me:

gaus(0)
gaus(0)

That does not seems correct …

If you try to define a function that way you get similar errors messages:

root [0] TF1 *a=new TF1("a","gauss(0)",40,1000);
input_line_10:2:55: error: use of undeclared identifier 'gauss'
Double_t TFormula____id11207996195491148296(){ return gauss(0) ; }
                                                      ^
input_line_11:2:55: error: use of undeclared identifier 'gauss'
Double_t TFormula____id11207996195491148296(){ return gauss(0) ; }
                                                      ^
Error in <prepareMethod>: Can't compile function TFormula____id11207996195491148296 prototype with arguments 
Error in <TFormula::InputFormulaIntoCling>: Error compiling formula expression in Cling
Error in <TFormula::ProcessFormula>: Formula "gauss(0)" is invalid !

[Erdrick]

Hi couet,

Thanks for the reply. However, the way I have defined the function is absolutely correct. I obviously have defined my TF1s correctly if I have used them to fit the histogram…

The way I have defined my functions can be seen in the ROOT documentation. You can see all of the examples yourself here: https://root.cern.ch/doc/master/classTF1.html

With the code snippet you posted, you spelled “gaus” wrong. You spelled it with two s’s. Notice it works perfectly fine with one s:

root [0] TF1* spelled_incorrectly = new TF1("spelled_incorrectly", "gauss(0)", 0, 10)
input_line_11:2:55: error: use of undeclared identifier 'gauss'
Double_t TFormula____id17246998206739916000(){ return gauss(0) ; }
                                                      ^
Error in <TFormula::PrepareEvalMethod>: Can't compile function TFormula____id17246998206739916000 prototype with arguments 
Error in <TFormula::InputFormulaIntoCling>: Error compiling formula expression in Cling
Error in <TFormula::ProcessFormula>: Formula "gauss(0)" is invalid !
(TF1 *) 0x1f19520
root [1] 
root [1] TF1* spelled_correctly = new TF1("spelled_correctly", "gaus(0)", 0, 10)
(TF1 *) 0x1494e30

Can anyone else please help me?

[couet]

you are right … my bad … I put two “s” at gaus !! …
so

root [0] TF1 *a=new TF1("a","gaus(0)",40,1000);

is ok…
So back to your initial question I am not sure why you get these error messages …
May be @moneta has an idea.

[Wile_E_Coyote]

Try:

  component_A->SetParNames("Constant_A", "Mean_A", "Sigma_A");
  gROOT->GetListOfFunctions()->Add(component_A);
  component_B->SetParNames("Constant_B", "Mean_B", "Sigma_B");
  gROOT->GetListOfFunctions()->Add(component_B);
  TF1 *difference = new TF1("difference", Form("abs(%s(x)-%s(x))",
                            component_A->GetName(), component_B->GetName()),
                            0., 10000.);
  difference->SetNpx(10000);
  difference->Print("V");

BTW. You have two gaussians so you can derive an analytical formula for their intersection, which makes much more sense than trying to do it numerically.

[Erdrick]

Hi Wile_E_Coyote,

Thank you very much for your help! Yes, I realize I can derive an analytical formula for the intersection. But this is a minimal, simple example. With more complicated functions, I would like to do it numerically.

I still have some questions, however. The script that incorporates your help (attached as minimal.C (1.8 KB)) is:

#include <iostream>

#include "TFile.h"
#include "TH1F.h"
#include "TF1.h"
#include "TCanvas.h"

void minimal() {

  TFile* file = TFile::Open("minimal.root", "READ");

  TH1F* hist = (TH1F*) file->Get("run1396_hippo_elastic_scattering_beam_left_hit");
  std::cout << hist->GetName() << std::endl;
  std::cout << hist->GetListOfFunctions()->At(1)->GetName() << std::endl;
  std::cout << hist->GetListOfFunctions()->At(2)->GetName() << std::endl;
  std::cout << hist->GetListOfFunctions()->At(3)->GetName() << std::endl;
  
  TCanvas* canvas = new TCanvas("canvas", "canvas");
  gPad->SetLogy();

  hist->Draw();

  ((TF1*) hist->GetListOfFunctions()->At(1))->SetLineColor(1);
  ((TF1*) hist->GetListOfFunctions()->At(2))->SetLineColor(2);
  ((TF1*) hist->GetListOfFunctions()->At(3))->SetLineColor(3);

  ((TF1*) hist->GetListOfFunctions()->At(1))->SetLineStyle(2);
  ((TF1*) hist->GetListOfFunctions()->At(2))->SetLineStyle(2);
  ((TF1*) hist->GetListOfFunctions()->At(3))->SetLineStyle(2);

  ((TF1*) hist->GetListOfFunctions()->At(1))->Draw("same");
  ((TF1*) hist->GetListOfFunctions()->At(2))->Draw("same");
  ((TF1*) hist->GetListOfFunctions()->At(3))->Draw("same");

  TF1* component_A = (TF1*) hist->GetListOfFunctions()->At(2);
  TF1* component_B = (TF1*) hist->GetListOfFunctions()->At(3);

  gROOT->GetListOfFunctions()->Add(component_A);
  gROOT->GetListOfFunctions()->Add(component_B);

  component_A->SetParNames("Constant_A", "Mean_A", "Sigma_A"); // Why are these two lines necessary?!
  component_B->SetParNames("Constant_B", "Mean_B", "Sigma_B"); // Why are these two lines necessary?!

  TF1* difference = new TF1("difference",
			    Form("abs(%s-%s)", component_A->GetName(), component_B->GetName()),
			    0, 10000);

  difference->SetNpx(10000);
  difference->Draw("same");
  std::cout << difference->GetMinimumX() << std::endl;

}

My questions are:

1. So I need to “register my TF1s with ROOT” (via the gROOT::GetListOfFunctions::Add command) in order to use them in the TF1 that is the difference of the two? I assume/know the answer is yes based on what makes the script actually run, it just seems a little odd.
2. Why are the two lines of TF1::SetParNames necessary? If I comment out those two lines, the difference TF1 is not drawn and the result from std::cout << difference->GetMinimumX() << std::endl; is 1.0001. To me, having to set the parameter names to achieve the desired result is very odd.
3. Finally, if I run the script exactly as it is attached, the drawing of the difference TF1 looks correct. However, the result from std::cout << difference->GetMinimumX() << std::endl; is 2275.46. This value is not correct. Why?

Thanks again for all of your help.

[Wile_E_Coyote]

  difference->Print("V");
  Double_t left = TMath::Min(component_A->GetParameter(1), component_B->GetParameter(1)); // 1 = "Mean"
  Double_t right = TMath::Max(component_A->GetParameter(1), component_B->GetParameter(1)); // 1 = "Mean"
  std::cout << difference->GetMinimumX(left, right) << std::endl;