Ratio plot of two normalized histograms

waheede · March 4, 2016, 9:37am

Hi,
i have trouble getting the ratio plot of two normalized histograms. The error bars i get are very large though the mean is calculated right. This question was asked before as well but never addressed.

couet · March 4, 2016, 9:53am

We plan to develop something in that area this summer.
Meanwhile do you have a small reproducer of your problem ?

waheede · March 6, 2016, 7:05am

Thank you for your reply. This is part of the code of taking the ratio part only and plotting it.
The output i get is attached with large error bars.

[code]
TCanvas *c = new TCanvas(“c”,“Efficiency”,900,400);
c->Print(“ratiotest.pdf[”);

    hw->SetStats(0);
    h_w1_up->SetStats(0);

    hw->SetLineColor(kBlue);
    h_w1_up->SetLineColor(kRed);

    h_w1_up->Scale(1/h_w1_up->Integral());
    hw->Scale(1/hw->Integral());
    h_w1_up->SetMaximum(h_w1_up->GetMaximum()*1.5);
    h_w1_up->SetMinimum(0);
    hw->SetMinimum(0);
    h_w1_up->Draw("hist");
    hw->Draw("same");

    TLegend  *leg1g = new TLegend(0.75,0.85,0.96,0.95);
    leg1g->AddEntry(h_w1_up,"theory1","l");
    leg1g->AddEntry(hw,"Signal MC","l");
    leg1g->SetTextSize(0.025);
    leg1g->SetHeader("");
    leg1g->Draw();
    c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

double ymax =1.0;

    h_w_r1->GetYaxis()->SetTitleOffset(1.2);
    h_w_r1->SetStats(0);

    h_w_r1->SetMarkerStyle(21);
    h_w_r1->Divide(hw, h_w1_up,1,1);
    h_w_r1->SetMaximum(1.1);
    h_w_r1->SetMinimum(0.9);
    h_w_r1->Draw("ep");
    TLine *line2 =new TLine(1,ymax,1.504,ymax);
    line2->SetLineColor(kRed);
    line2->Draw();

    c->Print("ratiotest.pdf",  "Title: "  "w Theory ");

    c->Print("ratiotest.pdf]");
   [/code]

ratiotest.pdf (16 KB)

yus · March 6, 2016, 1:56pm

Hi,

I think you need to correct the stat. errors for your histograms, dividing them by a corresponding integral. Try something like this:

float red_histo_integral = hw->Integral();
float blue_histo_integral = h_w1_up->Integral();
float rel_uncert_overall = 0;
for (int i=1; i<=(h_w_r1->GetNbinsX()); i++)
{
    rel_uncert_from_blue_histo = (hw->GetBinContent(i) == 0)?0.:hw->GetBinError(i)/blue_histo_integral/hw->GetBinContent(i);
    rel_uncert_from_red_histo = (h_w1_up->GetBinContent(i) == 0)?0.:h_w1_up->GetBinError(i)/red_histo_integral/h_w1_up->GetBinContent(i);
    rel_uncert_overall = sqrt(pow(rel_uncert_from_blue_histo,2)+pow(rel_uncert_from_red_histo,2));
    h_w_r1->SetBinError(i, rel_uncert_overall*h_w_r1->GetBinContent(i));
}

waheede · March 7, 2016, 12:20am

[quote=“yus”]Hi,

I think you need to correct the stat. errors for your histograms, dividing them by a corresponding integral. Try something like this:

float red_histo_integral = hw->Integral(); float blue_histo_integral = h_w1_up->Integral(); float rel_uncert_overall = 0; for (int i=1; i<=(h_w_r1->GetNbinsX()); i++) { rel_uncert_from_blue_histo = (hw->GetBinContent(i) == 0)?0.:hw->GetBinError(i)/blue_histo_integral/hw->GetBinContent(i); rel_uncert_from_red_histo = (h_w1_up->GetBinContent(i) == 0)?0.:h_w1_up->GetBinError(i)/red_histo_integral/h_w1_up->GetBinContent(i); rel_uncert_overall = sqrt(pow(rel_uncert_from_blue_histo,2)+pow(rel_uncert_from_red_histo,2)); h_w_r1->SetBinError(i, rel_uncert_overall*h_w_r1->GetBinContent(i)); }[/quote]

Hi, Thank you so much for your reply. I have also tried to calculate the errors by hand as well.
But following your procedure. I get the errors =0 . And i think this is because of the line " h_w_r1->SetBinError(i, rel_uncert_overall*h_w_r1->GetBinContent(i));" I don’t understand that why do we have to multiply the bin content of h_w_r1 histogram which is empty before that is why it is giving zero.
For me h_w_r1 was the histogram that calculates the errors.
Waiting for your response .

yus · March 7, 2016, 1:36pm

Hi,

on the other hand, there is an easier solution. Here’s my standalone script that demonstrates it all, along with a resulting pdf with histograms. If you run the script, you’ll see that uncertainties are calculated correctly.

P.S. I don’t have your histograms, so I came up with my own ones.

void test_07Mar2016()
{
	TCanvas *c = new TCanvas("c","Efficiency",900,400);

	TH1F* hw = new TH1F("hw", "hw;X;Y", 10, 0., 10.);
	hw->SetStats(0);
	hw->SetLineColor(kBlue);

	TH1F* h_w1_up = new TH1F("h_w1_up", "h_w1_up;X;Y", 10, 0., 10.);
	h_w1_up->SetStats(0);
	h_w1_up->SetLineColor(kRed);

	int i;

//filling both histograms bin-by-bin with some numbers
	for (i=1; i<=hw->GetNbinsX(); i++)
	{
		hw->SetBinContent(i, i);
		h_w1_up->SetBinContent(i, 11-i);
	}

//looking at stat. errors of both histograms (again, bin-by-bin)
	for (i=1; i<=hw->GetNbinsX(); i++)
		printf("Stat. errors for non-scaled histograms: bin #%02i: red = %07.3f%%, blue = %07.3f%%\n", i, h_w1_up->GetBinError(i)*100/h_w1_up->GetBinContent(i), hw->GetBinError(i)*100/hw->GetBinContent(i));

	h_w1_up->SetMaximum(h_w1_up->GetMaximum()*1.5);

//getting integrals before scaling histograms:
	float blue_histo_integral = hw->Integral();
	float red_histo_integral = h_w1_up->Integral();

	hw->Sumw2();
	h_w1_up->Sumw2();

	printf("Now scaling histograms...\n");

	hw->Scale(1/hw->Integral());
	h_w1_up->Scale(1/h_w1_up->Integral());

	for (i=1; i<=hw->GetNbinsX(); i++)
	{
		printf("Stat. errors for scaled histograms: bin #%02i: red = %07.3f%%, blue = %07.3f%%\n", i, h_w1_up->GetBinError(i)*100/h_w1_up->GetBinContent(i), hw->GetBinError(i)*100/hw->GetBinContent(i));
	}

	h_w1_up->Draw();
	hw->Draw("same");

	TLegend *leg1g = new TLegend(0.75,0.85,0.96,0.95);
	leg1g->AddEntry(h_w1_up,"theory1","l");
	leg1g->AddEntry(hw,"Signal MC","l");
	leg1g->SetTextSize(0.025);
	leg1g->SetHeader("");
	leg1g->Draw();
	c->Update();
	
	c->Print("ratiotest2.pdf(");







	TCanvas *c1 = new TCanvas("c1","Efficiency",900,400);
//now dividing histograms:
	TH1F h_w_r1 = (*hw)/(*h_w1_up);
	h_w_r1.GetYaxis()->SetTitleOffset(1.2);
	h_w_r1.SetStats(0);

	h_w_r1.SetMarkerStyle(21);
//	h_w_r1.SetMaximum(1.1);
//	h_w_r1.SetMinimum(0.9);
	h_w_r1.Draw("ep");
	c1->Update();

//now looking at the stat. errors of this ratio histogram. They equal sqrt(rel_unc_of_histo_in_numerator^2 + rel_unc_of_histo_in_denominator^2)

	for (i=1; i<=h_w_r1.GetNbinsX(); i++)
		printf("Stat. errors for a ratio histograms: bin #%02i: %07.3f%%\n", i, h_w_r1.GetBinError(i)*100/h_w_r1.GetBinContent(i));

	double ymax =1.0;
	TLine *line2 =new TLine(0, ymax, 10, ymax);
	line2->SetLineColor(kRed);

	line2->Draw();

	c1->Print("ratiotest2.pdf)");

}

ratiotest2.pdf (16.6 KB)

moneta · March 7, 2016, 3:54pm

Hi,

If you have bins with low statistics you should use the correct Poisson description for the bins.
The correct errors are then compute in ROOT using the TGraphAsymErrors::Divide function using the “pois” option, that indicates “Poisson ratio”.
I have attached an example where the ratio plot for two normalised histograms is computed using TH1::Divide (normal approximation) or
TGraphAsymErrors::Divide.

In your first example, I guess the original histogram do not have the right errors. Remember you should call
TH1::Sumw2() before scaling the histogram

Best Regards

LorenzoexampleRatioPlot.C (1.14 KB)

waheede · March 8, 2016, 10:04am

Hi,
Thank you all of you guys for helping but despite of trying every way mentioned i couldn’t get the desired output.
I have tried all the ways and attached my whole code along with my root file containing the histogram and the output pdf.
I have no clue whats going on.
mc_test.root (8.38 KB)
ratiotest.pdf (27 KB)
ratiotest.C (8.33 KB)