Make conditional histograms with ROOT

example_cond.png

make_conditional.py

import ROOT


def make_conditional(histo):
    """
    return a "conditional" histogram, where the entries have
    been normalized by columns
    """
    result = histo.Clone(histo.GetName() + "_condX")
    proj = histo.ProjectionX()

    for xbin in xrange(histo.GetNbinsX() + 2):
        norm = float(proj.GetBinContent(xbin))
        if norm == 0:
            continue
        for ybin in xrange(histo.GetNbinsY() + 2):
            value = histo.GetBinContent(xbin, ybin)
            result.SetBinContent(xbin, ybin, value / norm)
    return result


if __name__ == "__main__":
    ROOT.gStyle.SetPalette(55)
    h = ROOT.TH2F("h", "P(x, y) = exp(#lambda x) #times N(y|3 + x, 3);x;y", 100, 0, 20, 100, 0, 20)

    import random
    NTOYS = 1000000
    for itoy in xrange(NTOYS):
        x = random.expovariate(0.3)
        y = random.normalvariate(3 + x, 3)
        h.Fill(x, y)
    h.Scale(1. / NTOYS)

    canvas = ROOT.TCanvas()
    canvas.Divide(2, 1)
    canvas.cd(1)
    h.Draw("colz")
    h.SetContour(100)
    canvas.cd(2)
    h_cond = make_conditional(h)
    h_cond.SetTitle("P(y|x) = P(x, y) / #scale[0.5]{#int} P(x, y) dy")
    h_cond.Draw("colz")
    h_cond.SetMaximum(0.04)
    h_cond.SetContour(100)

    canvas.SaveAs("example_cond.png")
    raw_input()

make_conditional_example.C

#include <string>


/**
 * return a "conditional" histogram, where the entries have
 * been normalized by columns. User must delete the returned
 * histogram
 **/

TH2* make_conditional(const TH2& histo)
{
    auto result = dynamic_cast<TH2*>(histo.Clone((std::string(histo.GetName()) + "_condX").c_str()));
    const auto proj = histo.ProjectionX();
    
    for (xbin = 0; xbin != histo.GetNbinsX(); ++xbin) {
	double norm = proj->GetBinContent(xbin);
	if (norm == 0) continue;
	for (ybin = 0; ybin != histo.GetNbinsY(); ++ybin) {
	    const auto value = histo.GetBinContent(xbin, ybin);
	    result->SetBinContent(xbin, ybin, value / norm);
	}
    }
    return result;
}


make_conditional_example()
{
    gStyle->SetPalette(55);

    TH2F h("h", "P(x, y) = exp(#lambda x) #times N(y|3 + x, 3);x;y", 100, 0, 20, 100, 0, 20);

    const int NTOYS = 100000;
    for (int itoy = 0; itoy < NTOYS; ++itoy) {
	const auto x = gRandom->Exp(1);
	const auto y = gRandom->Gaus(x + 3, 3);
	h.Fill(x, y);
    }
    h.Scale(1. / NTOYS);
    
    TCanvas canvas;
    canvas.Divide(2, 1);
    canvas.cd(1);
    h.Draw("colz");
    h.SetContour(100);
    canvas.cd(2);
    auto h_cond = make_conditional(h);
    h_cond->SetTitle("P(y|x) = P(x, y) / #scale[0.5]{#int} P(x, y) dy");
    h_cond->Draw("colz");
    h_cond->SetMaximum(0.04);
    h_cond->SetContour(100);
    canvas.SaveAs("example_cond_cpp.png");
    
    delete h_cond;
}