Issue interpolating histogram shapes with RooIntegralMorph

Hi all,

I am attempting to interpolate between two known histograms to produce a shape between them. I am running into a problem where the output histogram has the correct shape but the height and integral are less than either of the input histograms. The problem is illustrated in the attached photo, you will see the black histogram which was produced by interpolating between the red and green histograms.

I am also attaching the code that made this with the caveat that it is not a completely working example as attached. This is part of a much larger code and reads in MC Signal files which I am unable to share at this time. If this is not enough to solve the issue, let me know and I can try to create a short version of the code.

The interpolation is on lines 74-81, which I will copy below:

    RHL = ROOT.RooDataHist("HL_".format(un), ";DiCluster Mass [GeV];Events/GeV", ROOT.RooArgList(self._x), HL)
    RHLR = ROOT.RooHistPdf("HL_AbsReal_{}".format(un), "", ROOT.RooArgSet(self._x), RHL)
    RHH = ROOT.RooDataHist("HH_{}".format(un), ";DiCluster Mass [GeV];Events/GeV", ROOT.RooArgList(self._x), HH)
    RHHR = ROOT.RooHistPdf("HH_AbsReal_{}".format(un), "", ROOT.RooArgSet(self._x), RHH)

    RHIM = ROOT.RooIntegralMorph("Hmorph_{}".format(un), "", RHHR, RHLR, self._x, rmass)
    self.xframe = self._x.frame(ROOT.RooFit.Title(";DiCluster Mass [GeV];Events/GeV"), ROOT.RooFit.Range(0, 10000))
    RHI = RHIM.createHistogram("Hinterpo_{}".format(un), self._x)

The lines after this (85 - 107) plot the output of this and produce the attached plot. As you can see, I am exactly plotting the two input and output histograms and for some reason the output does not seem to make sense, it would seem impossible that the interpolation would produce a histogram that is lower in height than either of the two input histograms.

Thank you in advance for you help, let me know if there are other questions I can answer

The forum will not let me upload a link since I am a new user, so I am putting the whole code snippet here:

import ROOT
import numpy
import os
import math
import sys

def computeBoundingIndices(M, anchors):
  lowI, highI = 0,0

  minlDiff=9999
  minhDiff=9999
  for ind, N in enumerate(anchors):
    if(N < M):
      if (M - N < minlDiff): 
        minlDiff = M - N
        lowI = ind

    elif(N > M):
      if (N - M < minhDiff): 
        minhDiff = N - M
        highI = ind

  return lowI, highI

def integralInterpo(Min, INTS, M, log=True):
  MSS = [float(k) for k in Min]

  SPLINE = 0
  if log:
    IL, IH = computeBoundingIndices(M, MSS)
    TF = ROOT.TF1("tempF", "[1]*TMath::Exp((x-[0])*TMath::Log([3]/[1])/([2]-[0]))", M-50, M+50)
    TF.SetParameter(0, MSS[IL])
    TF.SetParameter(1, INTS[IL])
    TF.SetParameter(2, MSS[IH])
    TF.SetParameter(3, INTS[IH])
    SPLINE = TF.Eval(M)

  else:
    LINTS = [numpy.log(INTS) for sl in range(len(INTS.keys()))]
    TG = ROOT.TGraph(len(MSS), numpy.array(MSS), numpy.array(LINTS[sl]))
    TG.SetBit(ROOT.TGraph.kIsSortedX)
    SPLINE = TG.Eval(M)

  return SPLINE


class HC:
  def __init__(self, histArr, massArr):
    self._massArr = massArr
    self._histArr = histArr
    self._x  = ROOT.RooRealVar("x","x",histArr[0].GetXaxis().GetXmin(),histArr[0].GetXaxis().GetXmax())
    self._x.setBins(histArr[0].GetNbinsX())
    self._histInts = [h.Integral() for h in histArr]
    self._inxhists = []
    self._cutEff = []

  def morph(self, MM, var, nameM, N):
    self._lowI, self._hiI = computeBoundingIndices(MM, self._massArr)
    print(self._massArr[self._lowI], self._massArr[self._hiI])

    HL = self._histArr[self._lowI].Clone(self._histArr[self._lowI].GetName() + "HL")
    HH = self._histArr[self._hiI].Clone(self._histArr[self._hiI].GetName() + "HH")

    inxhists = [HL, HH]
    un = nameM

    alpha = (float(MM) - float(self._massArr[self._lowI]))/(float(self._massArr[self._hiI]) - float(self._massArr[self._lowI]))
    tlM, thM = self._massArr[self._lowI], self._massArr[self._hiI]
    lM = min(tlM, thM)
    hM = max(tlM, thM)
    wpoint = float(MM - lM) / float(hM - lM)
    rmass = ROOT.RooRealVar("rm_{}".format(un), "rmass", wpoint, 0., 1.)

    RHL = ROOT.RooDataHist("HL_".format(un), ";DiCluster Mass [GeV];Events/GeV", ROOT.RooArgList(self._x), HL)
    RHLR = ROOT.RooHistPdf("HL_AbsReal_{}".format(un), "", ROOT.RooArgSet(self._x), RHL)
    RHH = ROOT.RooDataHist("HH_{}".format(un), ";DiCluster Mass [GeV];Events/GeV", ROOT.RooArgList(self._x), HH)
    RHHR = ROOT.RooHistPdf("HH_AbsReal_{}".format(un), "", ROOT.RooArgSet(self._x), RHH)

    RHIM = ROOT.RooIntegralMorph("Hmorph_{}".format(un), "", RHHR, RHLR, self._x, rmass)
    self.xframe = self._x.frame(ROOT.RooFit.Title(";DiCluster Mass [GeV];Events/GeV"), ROOT.RooFit.Range(0, 10000))
    RHI = RHIM.createHistogram("Hinterpo_{}".format(un), self._x)

    ##
    ##
    c1 = ROOT.TCanvas()
    c1.cd()
    ll = ROOT.TLegend(0.6,0.5,0.8,0.75)
    ll.SetBorderSize(0)
    HH.SetTitle("HH")
    HH.SetLineColor(ROOT.kRed)
    HL.SetLineColor(ROOT.kGreen)
    HL.SetLineWidth(2)
    HL.SetTitle("HL")

    ll.AddEntry(HL, "In_Low")
    ll.AddEntry(HH, "In_High")
    ##
    rr = RHI.Clone(un+N)
    rr.SetTitle("OUT")
    rr.SetLineColor(ROOT.kBlack)
    ll.AddEntry(rr, "OUT")
    FindAndSetMax([HH, HL, rr])
    HH.Draw("hist")
    HL.Draw("histsame")
    rr.Draw("histsame")
    ll.Draw("same")
    c1.Print("tc3.png")
    ##

    return RHI.Clone(un+N), inxhists

Hi @SClarkPhysics; I am sure @couet can help you with this.

Cheers,
J.

I think @moneta will know better.

I have learned a bit more about this problem and found a workaround, although I don’t fully understand why it works and wouldn’t call it a “solution”

It seems as if the shortened height of the output histogram is coming as a result of the long tails on the two input histograms. If I take the following line from init function in the code above:

self._x  = ROOT.RooRealVar("x","x",histArr[0].GetXaxis().GetXmin(),histArr[0].GetXaxis().GetXmax())

and I replace it with

self._x  = ROOT.RooRealVar("x","x",histArr[0].GetXaxis().GetXmin(), 2100.)

then I get the output attached here, which is what I would expect. Now I don’t know exactly why this work and what the role of this RooRealVar is. I chose 2100 as the upper limit as it looks like where the histograms go to zero from the initial plots, so I don’t know why exactly the output is changing. Perhaps bins with 0 bin count are treated differently by the RooIntegralMorph function? Any insight?

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.