Fitting TH1 to TH1

Newbie
1

ROOT Experts,

I am looking to decompose a TH1 histogram into a linear combination of other TH1 histograms. I am familiar with fitting histograms to functions, but instead I am looking to do something of the flavor h3->Fit([0]*h1+[1]*h2) where [0] and [1] are the weights required to fit the histogram. What is the simplest way to do this?

A simple example to explain what I am trying to do is below.

`
void RootTalk() {

// Run as > root.exe RootTalk.C

TH1F *h1 = new TH1F(“h1”, “h1 title”, 5, 0, 5);
TH1F *h2 = new TH1F(“h2”, “h2 title”, 5, 0, 5);
TH1F *h3 = new TH1F(“h3”, “h3 title”, 5, 0, 5);

h1->SetBinContent(1,1);
h1->SetBinContent(2,2);

h2->SetBinContent(2,1);

h3->SetBinContent(1,3);
h3->SetBinContent(2,4);

// Insert code to fit some amount of h1 and some amount of h2 to equal h3.
// The “correct” answer is the h3 = 2*h1 + h2. Can I make a fit to do this?

}
`

2

Hello,

Welcome to the ROOT community!
Perhaps by encancapsulating the TH1s in a single function, and use that in the fit. Alternatively, you can use RooFit, and leverage a sum of RooHistPdf instances.

I hope this helps!

Best,
Danilo

3

Danilo,

My updated dummy example is below. I converted my TH1’s which I want to fit with to TGraphs. I am unsure how to convert them further to TF1 so that they can fit “h3”. Is this the sort of thing you have in mind with “encancapsulating the TH1s in a single function, and use that in the fit”?

I am open to RooFit and RooHistPdf instances, but I cannot make sense of the documentation. Are there examples I can look at? Or is it possible to demonstrate with my toy example below?

void RootTalk() {

  // Run as > root.exe RootTalk.C

  TH1F *h1 = new TH1F("h1", "h1 title", 5, 0, 5);
  TH1F *h2 = new TH1F("h2", "h2 title", 5, 0, 5);
  TH1F *h3 = new TH1F("h3", "h3 title", 5, 0, 5);

  TGraph *gr1 = new TGraph(h1->GetNbinsX());
  TGraph *gr2 = new TGraph(h2->GetNbinsX());

  h1->SetBinContent(1,1);
  h1->SetBinContent(2,2);

  h2->SetBinContent(2,1);

  h3->SetBinContent(1,3);
  h3->SetBinContent(2,4);

  for(int i=1; i < h1->GetNbinsX(); ++i) {

    gr1->SetPoint(i,h1->GetBinCenter(i),h1->GetBinContent(i));
    gr2->SetPoint(i,h2->GetBinCenter(i),h2->GetBinContent(i));

  }

//  TF1 *f = new TF1("f", "[0]*gr1+[1]*gr2"); // bad syntax
//  h3->Fit("f");
// If successful, the "correct" answer is h3 = 2*h1 + h2. Can I make a fit to do this?

}
4

Hello Will,

Thanks for sharing the code. You can have a look to this documentation, that shows how to wrap a generic C++ function in a TF1. In your case, the C++ code would foresee the “weighted evaluation” of the two TGraphs interpolations.

I hope this helps.

Cheers,
D

5

Hi Will,
below you can find the code for fitting an histogram as sum of two other histograms.

I changed a bit your code.
The first change concerns how build the graph, nothing wrong with your method the final result is basically the same, but with my method the final function resemble an histogram.

The TGraphs are now declared as global variable otherwise I was not able to use it inside the function definition.

Other difference, I changed the values you set in h3.
With your value the correct answer was 3 * h1 - 2 * h2 and the program was able to found it.
WIth the current values in h3 the expected answer is 2 * h1+1 * h2

TGraph *gr1 = new TGraph();
TGraph *gr2 = new TGraph();

void RootTalk() {

  // Run as > root.exe RootTalk.C

  TH1F *h1 = new TH1F("h1", "h1 title", 5, 0, 5);
  TH1F *h2 = new TH1F("h2", "h2 title", 5, 0, 5);
  TH1F *h3 = new TH1F("h3", "h3 title", 5, 0, 5);

  
  h1->SetBinContent(1,1);
 
  h3->SetBinContent(1,2);//h3->SetBinContent(1,3) this was your previous setting
  //2*1+ 1*0=2
  
  h1->SetBinContent(2,2);
  h2->SetBinContent(2,1);
  h3->SetBinContent(2,5);//h3->SetBinContent(2,4) this was your previous setting
  //2*2 + 1*1 = 5

  

  for(int i=1; i <=h1->GetNbinsX(); ++i) {



    gr1->SetPoint(gr1->GetN(),h1->GetBinLowEdge(i)+0.001*h1->GetBinWidth(i),h1->GetBinContent(i));
    gr1->SetPoint(gr1->GetN(),h1->GetBinLowEdge(i+1)-0.001*h1->GetBinWidth(i),h1->GetBinContent(i));

    gr2->SetPoint(gr2->GetN(),h2->GetBinLowEdge(i)+0.001*h2->GetBinWidth(i),h2->GetBinContent(i));
    gr2->SetPoint(gr2->GetN(),h2->GetBinLowEdge(i+1)-0.001*h2->GetBinWidth(i),h2->GetBinContent(i));


   //gr1->SetPoint(gr1->GetN(),h1->GetBinCenter(i),h1->GetBinContent(i));
  //gr2->SetPoint(gr1->GetN(),h2->GetBinCenter(i),h2->GetBinContent(i));

  }

  TF1 * f_sum=new TF1("f_sum",[&](double* x, double* p){return p[0]*gr1->Eval(x[0]) + p[1]*gr2->Eval(x[0]);},0.,5,2);

f_sum->SetNpx(1000);
	
	

h3->Draw("");

	
 h3->Fit("f_sum");



//  TF1 *f = new TF1("f", "[0]*gr1+[1]*gr2"); // bad syntax
//  h3->Fit("f");
// If successful, the "correct" answer is h3 = 2 * h1 + h2. Can I make a fit to do this?

}

Best,
Stefano

6

Stefano,

Thank you! Thank you!

Two things:
-Your post had left off the declaration of the TGraphs. I put this back in for easy reference of others.
-Though the fit gives the correct answer, I get a seg fault style error (below). This persists whether I run as “root.exe RootTalk.C” or go into root and run as “.L RootTalk.C” then “RootTalk()”. I am running root 6.30/04 on lxplus.

Regardless, you have solved my main issue and I’m deeply grateful!

root [0] 
Processing RootTalk.C...
Info in <TCanvas::MakeDefCanvas>:  created default TCanvas with name c1
****************************************
Minimizer is Minuit2 / Migrad
Chi2                      =  9.12785e-26
NDf                       =            0
Edm                       =  9.13279e-26
NCalls                    =           29
p0                        =            2   +/-   1.41421     
p1                        =            1   +/-   3.60555     
terminate called after throwing an instance of 'cling::InvalidDerefException'
  what():  Trying to access a pointer that points to an invalid memory address.
Aborted (core dumped)

void RootTalk() {

  // Run as >   root.exe
  //        [0] .L RootTalk.C
  //        [1] RootTalk() 

  TH1F *h1 = new TH1F("h1", "h1 title", 5, 0, 5);
  TH1F *h2 = new TH1F("h2", "h2 title", 5, 0, 5);
  TH1F *h3 = new TH1F("h3", "h3 title", 5, 0, 5);

  TGraph *gr1 = new TGraph(h1->GetNbinsX());
  TGraph *gr2 = new TGraph(h2->GetNbinsX());

  h1->SetBinContent(1,1);

  h3->SetBinContent(1,2);//h3->SetBinContent(1,3) this was your previous setting
  //2*1+ 1*0=2

  h1->SetBinContent(2,2);
  h2->SetBinContent(2,1);
  h3->SetBinContent(2,5);//h3->SetBinContent(2,4) this was your previous setting
  //2*2 + 1*1 = 5

  for(int i=1; i <=h1->GetNbinsX(); ++i) {



    gr1->SetPoint(gr1->GetN(),h1->GetBinLowEdge(i)+0.001*h1->GetBinWidth(i),h1->GetBinContent(i));
    gr1->SetPoint(gr1->GetN(),h1->GetBinLowEdge(i+1)-0.001*h1->GetBinWidth(i),h1->GetBinContent(i));

    gr2->SetPoint(gr2->GetN(),h2->GetBinLowEdge(i)+0.001*h2->GetBinWidth(i),h2->GetBinContent(i));
    gr2->SetPoint(gr2->GetN(),h2->GetBinLowEdge(i+1)-0.001*h2->GetBinWidth(i),h2->GetBinContent(i));


    //gr1->SetPoint(gr1->GetN(),h1->GetBinCenter(i),h1->GetBinContent(i));
    //gr2->SetPoint(gr1->GetN(),h2->GetBinCenter(i),h2->GetBinContent(i));

  }

  TF1 * f_sum=new TF1("f_sum",[&](double* x, double* p){return p[0]*gr1->Eval(x[0]) + p[1]*gr2->Eval(x[0]);},0.,5,2);
  f_sum->SetNpx(1000);

  h3->Draw("");
  h3->Fit("f_sum");
}
7

Making the TGraphs global solved my seg fault issue.

Thanks again!
Will

#include <TH1.h>
#include <TGraph.h>
#include <TF1.h>
#include <TCanvas.h>

TGraph *gr1;
TGraph *gr2;

TCanvas* RootTalk() {

  // Run as >   root.exe
  //        [0] .L RootTalk.C+
  //        [1] RootTalk() 

  TH1F *h1 = new TH1F("h1", "h1 title", 5, 0, 5);
  TH1F *h2 = new TH1F("h2", "h2 title", 5, 0, 5);
  TH1F *h3 = new TH1F("h3", "h3 title", 5, 0, 5);

  gr1 = new TGraph(h1->GetNbinsX());
  gr2 = new TGraph(h2->GetNbinsX());

  h1->SetBinContent(1,1);

  h3->SetBinContent(1,2);//h3->SetBinContent(1,3) this was your previous setting
  //2*1+ 1*0=2

  h1->SetBinContent(2,2);
  h2->SetBinContent(2,1);
  h3->SetBinContent(2,5);//h3->SetBinContent(2,4) this was your previous setting
  //2*2 + 1*1 = 5

  for(int i=1; i <=h1->GetNbinsX(); ++i) {

/*
    gr1->SetPoint(gr1->GetN(),h1->GetBinLowEdge(i)+0.001*h1->GetBinWidth(i),h1->GetBinContent(i));
    gr1->SetPoint(gr1->GetN(),h1->GetBinLowEdge(i+1)-0.001*h1->GetBinWidth(i),h1->GetBinContent(i));

    gr2->SetPoint(gr2->GetN(),h2->GetBinLowEdge(i)+0.001*h2->GetBinWidth(i),h2->GetBinContent(i));
    gr2->SetPoint(gr2->GetN(),h2->GetBinLowEdge(i+1)-0.001*h2->GetBinWidth(i),h2->GetBinContent(i));
*/

    gr1->SetPoint(gr1->GetN(),h1->GetBinCenter(i),h1->GetBinContent(i));
    gr2->SetPoint(gr1->GetN(),h2->GetBinCenter(i),h2->GetBinContent(i));

  }

  TF1 * f_sum=new TF1("f_sum",[&](double* x, double* p){return p[0]*gr1->Eval(x[0]) + p[1]*gr2->Eval(x[0]);},0.,5,2);
  f_sum->SetNpx(1000);

  h3->Fit("f_sum");
  TCanvas* c=new TCanvas("c1","c1",800,800);
  h3->Draw("");
  return c;

}
8

Yes as I told you in my previous post, the TGraphs should be declared as global otherwise you will have an error from the TF1 trying to access the TGraphs.

9

It’s sufficient to change “[&]” into “[=]”.