"""
All plots inherit from the :py:class:`~MPF.plotStore.PlotStore`.
.. inheritance-diagram:: MPF.plotStore.PlotStore
MPF.plot.Plot
MPF.dataMCRatioPlot.DataMCRatioPlot
MPF.bgContributionPlot.BGContributionPlot
MPF.significanceScanPlot.SignificanceScanPlot
MPF.signalRatioPlot.SignalRatioPlot
MPF.efficiencyPlot.EfficiencyPlot
To make a plot, create an instance of one of the :ref:`plot classes
<plotClassesList>`::
from MPF.plot import Plot
plot = Plot()
Most options concerning the plot and its style are set when the plot is instanciated::
plot = Plot(logy=True, yMin=1e-5, yMax=1e5)
See :py:meth:`~MPF.plotStore.PlotStore` for a list of all options.
Then use :py:meth:`~MPF.plotStore.PlotStore.registerHist` to register
`ROOT histograms <https://root.cern.ch/doc/master/classTH1.html>`_ to
the plot::
plot.registerHist(hist1, style="background", process="Bkg1")
The basic "styles" are `background`, `signal` and `data`. Background
histograms will be added to a `ROOT THStack
<https://root.cern.ch/doc/master/classTHStack.html>`_, each signal is
overlaid as a separate histogram and data is drawn as black dots with
error bars. If multiple histograms are added for the same process name
they are added up.
Finally the plot can be drawn using :py:meth:`~MPF.plotStore.PlotStore.saveAs`::
plot.saveAs(outputFilename)
The output format is determined by the filename extension.
"""
import math
import itertools
import ROOT
from .commonHelpers.logger import logger
# from MPF import logger
logger = logger.getChild(__name__)
from .histograms import getHM
from .histograms import Histogram, HistogramStack
from .canvas import Canvas
from .legend import Legend
from .arrow import Arrow
from .errorBands import AE, getPoissonRatio
from . import pyrootHelpers as PH
from .commonHelpers import interaction as IH
from . import globalStyle as gst
from .line import Line
from .atlasStyle import setAtlasStyle
from .labels import LumiLabel, CMELabel, ATLASLabel, ProcessLabel, InfoLabel
[docs]class missingHistogramException(Exception):
pass
[docs]class PlotStore(object):
'''
Store histograms and provide basic utilities to be implemented by plot classes
.. _plotClassesList:
Some plots (with example code):
* :py:meth:`~MPF.plot.Plot`
* :py:meth:`~MPF.dataMCRatioPlot.DataMCRatioPlot`
* :py:meth:`~MPF.bgContributionPlot.BGContributionPlot`
* :py:meth:`~MPF.significanceScanPlot.SignificanceScanPlot`
* :py:meth:`~MPF.signalRatioPlot.SignalRatioPlot`
* :py:meth:`~MPF.efficiencyPlot.EfficiencyPlot`
'''
def __init__(self,
splitting=None,
atlasLabel='Internal',
infoLabel=None,
processLabel=None,
targetLumi=None,
centerOfMassEnergy=13,
xTitle=None,
yTitle=None,
binLabels=None,
unit=None,
debugText=None,
customTexts=None,
lowerCut=None,
upperCut=None,
ratioUp=1.75,
ratioDown=0.25,
ignoreNumErrors=False,
ignoreDenErrors=False,
ratioMode="pois",
drawTotalBGError=True,
drawTotalBGErrorLegend=False,
drawRatioArrows=True,
drawRatioBkgErrorBand=True,
drawRatioLine=True,
addOverflowToLastBin=False,
addOverflowBinNote=False,
yMin=None,
yMax=None,
xAxisLabelsOption = None,
logy=False,
insideTopMargin=0.25,
normStack=False,
sortStackByIntegral=True,
raptorFunction=None):
""":param splitting: will be passed to :py:meth:`~MPF.canvas.Canvas`
:param atlasLabel: Text next to "ATLAS" (set None to deactivate)
:param infoLabel: Text underneath the lumiLabel/CMELabel (set None to deactivate)
:param processLabel: If not set to None a label with this process will be created
:param targetLumi: If not set to None a label with this integrated luminosity will be created
:param centerOfMassEnergy: If not set to None a label with this center of mass energy will be created
:param xTitle: will be assigned to all histograms added by :py:meth:`~MPF.plotStore.PlotStore.registerHist`
:param yTitle: y-Axis title for the main pad
:param binLabels: will be assigned to all histograms added by :py:meth:`~MPF.plotStore.PlotStore.registerHist`
:param unit: will be assigned to all histograms added by :py:meth:`~MPF.plotStore.PlotStore.registerHist`
:param debugText: Take list of lines and adds them to the main pad (might be deprecated soon)
:param customTexts: Take list of tuples with (xpos, ypos, text) and adds them as TLatex labels to the plot
:param lowerCut: print a line for the lower cut
:param upperCut: print a line for the upper cut
:param ratioUp/ratioDown: y-axis limits in ratio plots
:param ignoreNumErrors: ignore numerator errors in ratio plots
:param ignoreDenErrors: ignore denominator errors in ratio plots
:param ratioMode: how to calculate the ratio and the errors (see
:py:meth:`~MPF.plotStore.PlotStore.getRatioHistogram` for possible options)
:param drawTotalBGError: draw the error of the total bkg
:param drawTotalBGErrorLegend: draw the error of the total bkg in the legend
:param drawRatioArrows: draw arrows in ratio plots
:param drawRatioBkgErrorBand: draw error band in ratio plots that
indicates the denominator errors
:param drawRatioLine: draw a line at 1 in ratio plots
:param addOverflowToLastBin: Merge the overflow bin and the last bin which is plotted
:param addOverflowBinNote: Add a text that informs if the overflow was added to the last bin
:param yMin: Minimum on the mainpad y-Axis
:param yMax: Maximum on the mainpad y-Axis
:param xAxisLabelsOption: Pass option to the x-Axis labels (e.g. "v" for vertical labels) cf. `ROOT TAxis::LabelsOption <https://root.cern.ch/doc/master/classTAxis.html#a05dd3c5b4c3a1e32213544e35a33597c>`_
:param logy: Set Log scale on the mainpad y-Axis
:param insideTopMargin: Inside top margin in the mainpad
(between histograms maximum and the top - relativ to canvas
size)
:param sortStackByIntegral: If set to false the histograms will be stacked in the order they are registered
:param normStack: Normalise background stack to total Integral
:param raptorFunction: (Experimental) - Execute this function before saving the canvas (gets
plotStore as parameter). You should return all the stuff that should not
get out of scope or garbage collected
For setting options for the legend and labels see:
* :py:meth:`~MPF.plotStore.PlotStore.setLegendOptions`
* :py:meth:`~MPF.plotStore.PlotStore.setATLASLabelOptions`
* :py:meth:`~MPF.plotStore.PlotStore.setCMELabelOptions`
* :py:meth:`~MPF.plotStore.PlotStore.setLumiLabelOptions`
* :py:meth:`~MPF.plotStore.PlotStore.setProcessLabelOptions`
* :py:meth:`~MPF.plotStore.PlotStore.setInfoLabelOptions`
"""
logger.debug('create plotstore')
setAtlasStyle()
self.atlasLabel = atlasLabel
self.targetLumi = targetLumi
self.centerOfMassEnergy = centerOfMassEnergy
self.processLabel = processLabel
self.infoLabel = infoLabel
self.xTitle = xTitle
self.yTitle = yTitle
self.binLabels = binLabels
self.unit = unit
self.debugText = debugText
self.customTexts = customTexts
self.ratioUp = ratioUp
self.ratioDown = ratioDown
self.ignoreNumErrors = ignoreNumErrors
self.ignoreDenErrors = ignoreDenErrors
self.ratioMode = ratioMode
self.drawTotalBGError = drawTotalBGError
self.drawTotalBGErrorLegend = drawTotalBGErrorLegend
self.drawRatioArrows = drawRatioArrows
self.drawRatioBkgErrorBand = drawRatioBkgErrorBand
self.drawRatioLine = drawRatioLine
self.addOverflowToLastBin = addOverflowToLastBin
self.addOverflowBinNote = addOverflowBinNote
self.yMin = yMin
self.yMax = yMax
self.xAxisLabelsOption = xAxisLabelsOption
self.logy = logy
self.insideTopMargin = insideTopMargin
self.sortStackByIntegral = sortStackByIntegral
self.normStack = normStack
self.objects = []
self.sysHists = {}
self.stackOnTopHists = []
self.canvas = Canvas(splitting)
self.totalBG = None
self.totalBGVariations = None
self.upperCut = upperCut
self.lowerCut = lowerCut
# old stuff
self.totalData = None
self.colorScheme = {}
if not gst.randomColors:
self.colorGen = PH.yieldColors()
else:
self.colorGen = self.yieldRandomColors()
self.totalBGError = None
self.totalBGErrorLegend = "MC stat"
self.legendOptions = ((), {})
self.ATLASLabelOptions = ((), dict(text=self.atlasLabel))
if not gst.mergeCMEIntoLumiLabel:
self.lumiLabelOptions = ((), dict(lumi=self.targetLumi))
else:
self.lumiLabelOptions = ((), dict(lumi=self.targetLumi, cme=self.centerOfMassEnergy))
self.CMELabelOptions = ((), dict(cme=self.centerOfMassEnergy))
self.processLabelOptions = ((), dict(text=self.processLabel))
self.infoLabelOptions = ((), dict(text=self.infoLabel))
# for adding multiple ratios where a process has another process as denominator
self.ratioDenominators = {}
self.raptorFunction = raptorFunction
self.markerStyles = itertools.count(20)
[docs] def registerHist(self, hist,
style="signal",
drawString="",
process="",
legendTitle=None,
drawLegend=True,
hide=False,
lineStyle=None,
lineWidth=None,
fillStyle=None,
markerStyle=None,
color=None,
lineColor=None,
fillColor=None,
markerColor=None,
drawErrorBand=False,
maskBins=None,
xTitle=None,
unit=None,
ratioDenominatorProcess=None,
stackOnTop=False):
'''
Add histogram to plot
:param hist: `ROOT histogram <https://root.cern.ch/doc/master/classTH1.html>`_
:param style: Possible values: "signal", "background", "data"
:param drawString: Custom drawing option for this histogram
(see `ROOT THistPainter
<https://root.cern.ch/doc/master/classTHistPainter.html>`_)
:param process: Label this histogram with a process name
:param legendTitle: Custom legend title (if not given, process label is used)
:param drawLegend: Set to false if no legend should be drawn for this process
:param hide: Don't explicitely draw this histogram and don't include it in the legend
:param lineStyle: `ROOT line style <https://root.cern.ch/doc/master/classTAttLine.html#L3>`_ for this histogram
:param lineWidth: `ROOT line width <https://root.cern.ch/doc/master/classTAttLine.html#L2>`_ for this histogram
:param fillStyle: if given use this fill style instead of the
predifined one (see `ROOT::TAttFill <https://root.cern.ch/doc/master/classTAttFill.html>`_)
:param markerStyle: if given use this marker style (for data) instead of the
predifined one (see `ROOT::TAttMarker <https://root.cern.ch/doc/master/classTAttMarker.html>`_)
:param color: can be given as a `ROOT TColor enum
<https://root.cern.ch/doc/master/classTColor.html>`_ or as a
string (see :py:meth:`~MPF.pyrootHelpers.rootStyleColor`)
:param lineColor: set this if you want a different color for the line
:param fillColor: set this if you want a different color for the filled area
:param markerColor: set this if you want a different color for the marker
:param drawErrorBand: draw an errorband in the style of the usual total background error for this process
:param maskBins: truncate these bins (by bin numbers)
:type maskBins: list
:param xTitle: x-Axis title (overwriting what is in the ROOT histogram)
:param unit: if given it will be shown both on x and y Axis
:param ratioDenominatorProcess: if multiple ratios are to be drawn use
the histogram with the given process name as denominator
:param stackOnTop: stack this histogram on top of the total background (even if it is a signal)
'''
if hist.ClassName()[:3] == 'TH1':
hm = getHM(hist.Clone(next(PH.tempNames)))
hm.SetDirectory(0)
if hm.GetEntries() == 0:
logger.warning("The Histogram {} you are registering for process \"{}\" is empty".format(hm, process))
if self.addOverflowToLastBin:
logger.debug("Adding overflow to last bin for hist {} ({})".format(hist, process))
hm.addOverflowToLastBin()
# Set standard options
hm.drawString += drawString
hm.process = process
if legendTitle is not None:
hm.legend = legendTitle
else:
hm.legend = process
hm.hide = hide
hm.drawLegend = drawLegend
hm.lineStyle = lineStyle
hm.lineWidth = lineWidth
hm.fillStyle = fillStyle
hm.markerStyle = markerStyle
hm.lineColor = lineColor
hm.markerColor = markerColor
hm.fillColor = fillColor
hm.drawErrorBand = drawErrorBand
if markerStyle is None and style == "data":
hm.markerStyle = next(self.markerStyles)
self.determineColor(hm, color=color)
logger.debug('registering {} histogram {} with content {}, color {}, legend {}'
.format(style, hist.GetName(), hist.Integral(), hm.color, hm.process))
if maskBins is not None:
hm.maskedBins += maskBins
# probably not really useful to be able to set this here
if xTitle is not None:
hm.xTitle = xTitle
if unit is not None:
hm.unit = unit
# overwrite if defined in plot
if self.xTitle is not None:
hm.xTitle = self.xTitle
if self.unit is not None:
hm.unit = self.unit
# Set bin labels if necessary
if self.binLabels is not None:
if len(self.binLabels) != hm.GetNbinsX():
logger.error("length of binlabels ({}) does not match number of bins".format(len(self.binLabels), hm.GetNbinsX()))
for ibin in range(1, hm.GetNbinsX()+1):
hm.GetXaxis().SetBinLabel(ibin, self.binLabels[ibin-1])
# Set options from input
hm.style = style
if style == 'signal':
hm.zIndex = 2
elif style == 'data':
hm.zIndex = 1
if color is None:
hm.color = "kBlack"
elif style == 'background':
hm.zIndex = -1
else:
logger.error('histogram {} does not have a proper style'.format(hm))
self.objects.append(hm)
elif hist.ClassName()[:3] == 'TH2':
self.twoDHists.append(hist)
else:
logger.warning('could not add {} to plot {} not supported (only TH1/2)'.format(hist.GetName(), hist.ClassName()))
if ratioDenominatorProcess is not None:
self.setRatioDenominator(hm, ratioDenominatorProcess)
if stackOnTop:
self.stackOnTopHists.append(hm)
[docs] def setRatioDenominator(self, numeratorHist, targetProcessName):
for hist in self.getHists():
if hist.process == targetProcessName:
self.ratioDenominators[numeratorHist] = hist
break
else:
raise KeyError("Can't find process {}".format(targetProcessName))
[docs] def registerSysHist(self, nomProcess, sysName, hist):
hm = getHM(hist.Clone(next(PH.tempNames)))
if self.addOverflowToLastBin:
logger.debug("Adding overflow to sys hist {} for {}".format(sysName, nomProcess))
hm.addOverflowToLastBin()
if not sysName in self.sysHists:
self.sysHists[sysName] = {}
self.sysHists[sysName][nomProcess] = hm
self.totalBGErrorLegend = "sys #oplus MC stat"
[docs] def setLegendOptions(self, *args, **kwargs):
"""Set the options to be passed to :py:meth:`~MPF.legend.Legend`"""
self.legendOptions = (args, kwargs)
[docs] def setLumiLabelOptions(self, *args, **kwargs):
"""Set the options to be passed to :py:meth:`~MPF.labels.LumiLabel`"""
self.lumiLabelOptions = (args, kwargs)
[docs] def setCMELabelOptions(self, *args, **kwargs):
"""Set the options to be passed to :py:meth:`~MPF.labels.CMELabel`"""
self.CMELabelOptions = (args, kwargs)
[docs] def setATLASLabelOptions(self, *args, **kwargs):
"""Set the options to be passed to :py:meth:`~MPF.labels.ATLASLabel`"""
self.ATLASLabelOptions = (args, kwargs)
[docs] def setProcessLabelOptions(self, *args, **kwargs):
"""Set the options to be passed to :py:meth:`~MPF.labels.ProcessLabel`"""
self.processLabelOptions = (args, kwargs)
[docs] def setInfoLabelOptions(self, *args, **kwargs):
"""Set the options to be passed to :py:meth:`~MPF.labels.InfoLabel`"""
self.infoLabelOptions = (args, kwargs)
[docs] def getHists(self, *selection):
"""
return hists for the selections:
all
backgrounds
signal
data
noData
"""
if not self.objects:
raise missingHistogramException('no objects: {}'.format(self.objects))
for hist in self.objects:
# if type(histogram) is HistogramStack:
# continue
if not selection:
yield hist
continue
if 'nodata' in selection:
selection += ['signal', 'backgrounds']
if hist.style == 'signal':
if 'signal' in selection:
yield hist
elif hist.style == 'data':
if 'data' in selection:
yield hist
else:
if 'backgrounds' in selection:
yield hist
selection = []
[docs] def getTotalBG(self):
if self.totalBG is not None:
return self.totalBG
# create the nominal hist
totalBGNom = None
for hist in self.getHists("backgrounds"):
if totalBGNom is None:
totalBGNom = hist.clone()
else:
totalBGNom.add(hist)
totalBGVariations = []
# create the variations
for sysName, histsDict in self.sysHists.items():
totalBG = None
for nomHist in self.getHists("backgrounds"):
if nomHist.process in histsDict:
hist = histsDict[nomHist.process]
else:
hist = nomHist
if totalBG is None:
totalBG = hist.clone()
else:
totalBG.add(hist)
# for debugging
totalBG.process = sysName
totalBGVariations.append(totalBG)
self.totalBGVariations = totalBGVariations
if totalBGNom is not None:
totalBGNom.addSystematicError(*totalBGVariations)
self.totalBG = totalBGNom
if gst.noLinesForBkg and self.totalBG is not None:
# in this case the totalBG hist should have a line
self.totalBG.color = 0
self.totalBG.lineColor = 0
self.totalBG.fillStyle = 0
self.totalBG.zIndex = 0
self.totalBG.lineWidth = gst.totalBGLineWidth
return self.totalBG
[docs] def saveAs(self, path, **kwargs):
"""
Save the canvas. Arguments are passed to :py:meth:`MPF.canvas.Canvas.saveAs`
"""
if self.xAxisLabelsOption is not None:
self.canvas.xAxisLabelsOption = self.xAxisLabelsOption
if self.raptorFunction is not None:
kwargs.pop("noSave", "")
self.canvas.saveAs(path, **dict(noSave=True, **kwargs))
stuff = self.raptorFunction(self)
self.canvas.SaveAs(path)
self.canvas.Close()
return path
return self.canvas.saveAs(path, **kwargs)
def __repr__(self):
lines = ['####', 'PlotStore content:']
lines.append('Background:')
for hist in self.getHists('backgrounds'):
lines.append(" {}, {}".format(hist, hist.process))
lines.append('Signal:')
for hist in self.getHists('signal'):
lines.append(" {}, {}".format(hist, hist.process))
for i in range(hist.GetNbinsX()+1):
lines.append("{}: {}".format(i, hist.GetBinContent(i)))
lines.append('Data:')
for hist in self.getHists('data'):
lines.append(" {}, {}".format(hist, hist.process))
lines.append('####')
return '\n'.join(lines)
# Helper functions
[docs] def buildMainPad(self, **kwargs):
"""
Builds the main plot with stacked histograms, signals and data
"""
self.getCanvas()
try:
self.stackHistograms()
except ValueError:
pass
# First draw the main Pad (should always be defined)
mainpad = self.canvas.pads['main']
if self.atlasLabel is not None:
atlasArgs, atlasKwargs = self.ATLASLabelOptions
mainpad.drawables.append(ATLASLabel(*atlasArgs, **atlasKwargs))
if self.targetLumi:
lumiArgs, lumiKwargs = self.lumiLabelOptions
mainpad.drawables.append(LumiLabel(*lumiArgs, **lumiKwargs))
if self.centerOfMassEnergy:
cmeArgs, cmeKwargs = self.CMELabelOptions
if not self.targetLumi:
# shift to the left if no lumi label
xOffset = cmeKwargs.get("xOffset", 0) - (gst.CMELabelX-gst.lumiLabelX)
cmeKwargs = dict(cmeKwargs, xOffset=xOffset)
mainpad.drawables.append(CMELabel(*cmeArgs, **cmeKwargs))
if self.infoLabel:
infoArgs, infoKwargs = self.infoLabelOptions
mainpad.drawables.append(InfoLabel(*infoArgs, **infoKwargs))
if self.debugText is not None:
mainpad.debugText = self.debugText
if self.customTexts is not None:
mainpad.customTexts = self.customTexts
mainpad.targetLumi = self.targetLumi
mainpad.centerOfMassEnergy = self.centerOfMassEnergy
mpAdd = mainpad.drawables.append
if gst.rainbowMode:
colors = self.yieldRainbowColors(len(list(self.getHists("signal"))))
for histogram in self.getHists():
if (not histogram.style == 'background'
and not histogram.hide
and not histogram in self.stackOnTopHists):
if gst.rainbowMode:
try:
histogram.color = next(colors)
except StopIteration:
pass
mpAdd(histogram)
try:
mpAdd(self.hStack)
logger.debug('added stack')
except AttributeError:
pass
if self.drawTotalBGError:
totalBGError = AE(self.getTotalBG())
totalBGError.SetFillColor(gst.totalBGErrorColor)
totalBGError.SetFillStyle(gst.totalBGFillStyle)
mpAdd(totalBGError)
logger.debug("Added total BG error")
self.totalBGError = totalBGError
if gst.noLinesForBkg and self.getTotalBG():
mpAdd(self.getTotalBG())
if self.processLabel:
processArgs, processKwargs = self.processLabelOptions
mainpad.drawables.append(ProcessLabel(*processArgs, **processKwargs))
self.setupLegend(mainpad)
if self.upperCut is not None:
mainpad.addVertLine(self.upperCut, direction="left")
if self.lowerCut is not None:
mainpad.addVertLine(self.lowerCut, direction="right")
if self.yMin is not None:
mainpad.yMinimum = self.yMin
if self.yMax is not None:
mainpad.yMaximum = self.yMax
if self.yTitle is not None:
mainpad.yTitle = self.yTitle
mainpad.logy = self.logy
if self.insideTopMargin is not None:
mainpad.insideTopMargin = self.insideTopMargin
if self.addOverflowToLastBin and self.addOverflowBinNote:
mainpad.customTexts.append((0.2, 0.78, "Overflow added to last bin"))
[docs] def determineColor(self, hm, color):
if hm.process in self.colorScheme:
hm.color = self.colorScheme[hm.process]
elif color is None:
hm.color = next(self.colorGen)
else:
hm.color = PH.rootStyleColor(color)
return
[docs] def setupLegend(self, pad):
if pad.legend is not None:
logger.error("legend for pad {} already defined!".format(pad))
else:
legArgs, legKwargs = self.legendOptions
pad.legend = Legend(*legArgs, **legKwargs)
bgHists = [hist for hist in self.getHists('backgrounds')]
if self.sortStackByIntegral:
bgHists.sort(key=lambda hist: hist.Integral(), reverse=True)
for hist in bgHists:
pad.legend.addEntry(hist, 'F')
for hist in self.getHists('signal'):
if hist.fillStyle is None:
pad.legend.addEntry(hist, 'l')
else:
pad.legend.addEntry(hist, 'F')
for hist in self.getHists('data'):
pad.legend.addEntry(hist, 'P')
if self.totalBGError is not None and self.drawTotalBGErrorLegend:
self.totalBGError.legend = self.totalBGErrorLegend
logger.debug("add {} to legend".format(self.totalBGError))
self.totalBGError.hide = False
pad.legend.addEntry(self.totalBGError, "f")
[docs] def yieldRainbowColors(self, nCol):
ROOT.gStyle.SetPalette(ROOT.kRainBow)
for i in range(nCol):
yield ROOT.TColor.GetColorPalette(int(i*255./float(nCol)))
[docs] def yieldRandomColors(self):
while(True):
r = int(256*ROOT.gRandom.Rndm())
g = int(256*ROOT.gRandom.Rndm())
b = int(256*ROOT.gRandom.Rndm())
yield ROOT.TColor.GetColor(r, b, g)
[docs] def stackHistograms(self):
""" stacks background histograms with labels from histograms[0] """
bgHists = [hist for hist in self.getHists('backgrounds')]
integrals = [hist.Integral() for hist in bgHists]
if sum(integrals) == 0.0:
raise ValueError('no backgrund')
if self.sortStackByIntegral:
bgHists.sort(key=lambda hist: hist.Integral())
self.hStack = getHM(ROOT.THStack('hStack{}'.format(next(PH.tempNames)), ''))
if gst.rainbowMode:
colors = self.yieldRainbowColors(len(bgHists))
for histo in bgHists:
histo.setColorAndStyle()
if self.normStack:
histo.Scale(1./sum(integrals))
if not histo in self.stackOnTopHists:
if gst.rainbowMode:
histo.color = next(colors)
histo.setColorAndStyle()
self.hStack.add(histo)
for histo in self.stackOnTopHists:
histo.setColorAndStyle()
self.hStack.add(histo)
if self.normStack:
# we scaled the hists, so in case totalBG was already
# calcuated, need to recalculate when used next
self.totalBG = None
[docs] def getCanvas(self):
try:
return self.canvas
except AttributeError:
# self.canvas = Canvas(self.splitting)
self.canvas = Canvas('')
return self.canvas
[docs] def getYields(self, totalBG=True):
""" yield dictionaries with keys:
process
yield
error
fraction
"""
d = ROOT.Double()
totBG = 0.
if not self.getTotalBG():
totalBG = False
if totalBG:
totBG = self.getTotalBG().IntegralAndError(self.getTotalBG().GetXaxis().GetFirst(), self.getTotalBG().GetXaxis().GetLast(), d)
yield {"process": "bg tot", "yield": "{0:.2f}".format(totBG), "error": "{0:.2f}".format(d), "fraction": "{0:.2%}".format(1)}
for histo in self.getHists():
if type(histo) is HistogramStack:
continue
d = ROOT.Double()
integral = histo.IntegralAndError(histo.GetXaxis().GetFirst(), histo.GetXaxis().GetLast(), d)
try:
fraction = integral/totBG
except ZeroDivisionError:
fraction = -1
yield {"process": histo.process, "yield": "{0:.2f}".format(integral), "error": "{0:.2f}".format(d), "fraction": "{0:.2%}".format(fraction)}
[docs] def yieldTable(self):
""" returns table with yields"""
table = [["process", "yield", "stat. error"]]
for y in self.getYields():
table.append([y["process"],y["yield"], y["error"]])
return table
[docs] def dumpYields(self):
print (IH.printTable(self.yieldTable()))
[docs] def BGContributionHisto(self):
bgHists = [hist for hist in self.getHists("backgrounds")]
if self.sortStackByIntegral:
bgHists.sort(key=lambda hist: hist.Integral())
if len(bgHists) < 1:
logger.error("No backgrounds found")
# stackHists = []
hStack = getHM(ROOT.THStack('hStack{}'.format(next(PH.tempNames)), ''))
for hist in bgHists:
h = hist.clone()
h.Reset()
for i in range(1, hist.GetNbinsX()+1):
try:
h.SetBinContent(i, hist.GetBinContent(i)/self.getTotalBG().GetBinContent(i))
except ZeroDivisionError:
pass
hStack.add(h)
return hStack
[docs] def addDataMCRatio(self, pad, **kwargs):
ratioMode = kwargs.pop("ratioMode", self.ratioMode)
drawRatioLine = kwargs.pop("drawRatioLine", self.drawRatioLine)
if kwargs:
raise KeyError("got unexpected kwargs: {}".format(kwargs))
dataHists = [hist for hist in self.getHists("data")]
if len(dataHists) != 1:
logger.error("Sorry, we want exactly one data histogram, not {}".format(dataHists))
numerator = dataHists[0]
pad.addReferenceHist(numerator)
denominator = PH.getMergedHist([hist for hist in self.getHists("backgrounds")])
pad.rangeDown = numerator.GetXaxis().GetBinLowEdge(1)
pad.rangeUp = numerator.GetXaxis().GetBinUpEdge(numerator.GetNbinsX())
logger.debug("Pad range: {} - {}".format(pad.rangeDown, pad.rangeUp))
pad.yMinimum = self.ratioDown
pad.yMaximum = self.ratioUp
ratio = self.getRatioHistogram(numerator, denominator,
mode=ratioMode,
ignoreDenErrors=self.ignoreDenErrors,
ignoreNumErrors=self.ignoreNumErrors)
pad.drawables.append(ratio)
if self.drawRatioBkgErrorBand:
ratioErrorBand = AE(self.getTotalBG(), relErr=True)
ratioErrorBand.SetFillColor(gst.ratioErrorBandColor)
ratioErrorBand.SetFillStyle(gst.ratioErrorBandFillStyle)
ratioErrorBand.zIndex = 0
pad.drawables.append(ratioErrorBand)
if drawRatioLine:
line = Line(pad.rangeDown, 1, pad.rangeUp, 1,
lineStyle=1,
lineWidth=1,
lineColor=ROOT.kRed)
pad.drawables.append(line)
pad.yMinimum = self.ratioDown
pad.yMaximum = self.ratioUp
[docs] def addSignalRatios(self, pad, style="signal", addLegend=False, register=False, applyRatioLimits=True):
"""
Create multiple ratio histograms/graphs and add them to the given pad.
Uses :py:meth:`~MPF.plotStore.PlotStore.getRatioHistogram` (options are set
via :py:meth:`~MPF.plotStore.PlotStore`)
:param pad: Add the ratios to this pad.
:param style: Create ratios for all histograms of this style
:param addLegend: add legend entries for the ratios
:param register: register the ratio histograms to the plot (useful for creating data/mc ratio)
:param applyRatioLimits: apply the ratioUp, ratioDown parameters to the given pad
"""
sigHists = [hist for hist in self.getHists(style)]
if len(sigHists) < 2:
logger.error("Sorry, we want more than two histograms to create a ratio, not {}".format(sigHists))
denominator = sigHists[0]
pad.addReferenceHist(denominator)
numerators = sigHists[1:]
pad.rangeDown = denominator.GetXaxis().GetBinLowEdge(1)
pad.rangeUp = denominator.GetXaxis().GetBinUpEdge(denominator.GetNbinsX())
logger.debug("Pad range: {} - {}".format(pad.rangeDown, pad.rangeUp))
firstHistRegistered = False
for numerator in numerators:
if numerator in self.ratioDenominators.values():
# this hist is a denominator!
continue
if numerator in self.ratioDenominators:
denominator = self.ratioDenominators[numerator]
logger.debug("Numerator ({}): {}".format(style, numerator))
logger.debug("Denominator ({}): {}".format(style, denominator))
ratio = self.getRatioHistogram(numerator, denominator,
mode=self.ratioMode,
ignoreDenErrors=self.ignoreDenErrors,
ignoreNumErrors=self.ignoreNumErrors)
ratio.SetLineColor(PH.rootStyleColor(numerator.color))
ratio.SetMarkerColor(PH.rootStyleColor(numerator.color))
pad.drawables.append(ratio)
if register and not firstHistRegistered:
if issubclass(type(ratio), ROOT.TGraphAsymmErrors):
registerHist = getHM(PH.histFromGraph(ratio))
else:
registerHist = getHM(ratio)
Histogram.cloneAttributes(registerHist, numerator)
if style == "data":
registerHist.style = "data"
else:
registerHist.style = "background"
self.objects.append(registerHist)
# in this case the numerator and denominator hists are not used anymore
self.objects.remove(numerator)
firstHistRegistered = True
if addLegend:
legHist = ratio
if numerator.legend is not None:
legHist.legend = numerator.legend
else:
legHist.legend = numerator.process
pad.legend.addEntry(legHist, "p")
if register:
# now we also don't use the denominator anymore
self.objects.remove(denominator)
for sigHist in sigHists:
if sigHist.drawErrorBand and (sigHist is denominator):
ratioErrorBand = AE(sigHist, relErr=True)
ratioErrorBand.SetFillColor(gst.ratioErrorBandColor)
ratioErrorBand.SetFillStyle(gst.ratioErrorBandFillStyle)
ratioErrorBand.zIndex = 0
pad.drawables.append(ratioErrorBand)
if self.drawRatioLine:
line = Line(pad.rangeDown, 1, pad.rangeUp, 1,
lineStyle=1,
lineWidth=1,
lineColor=ROOT.kRed)
pad.drawables.append(line)
if applyRatioLimits:
pad.yMinimum = self.ratioDown
pad.yMaximum = self.ratioUp
[docs] def addSignificanceScan(self, pad, processesOfInterest=None, scanDirection="forward", **kwargs):
"""Adds a significance scan for each signal in the given pad. The
kwargs are passed to getSignificanceHistogram"""
if processesOfInterest is None:
signalHists = self.getHists("signal")
else:
signalHists = []
for hist in self.objects:
if hist.process in processesOfInterest:
signalHists.append(hist)
if not scanDirection in ["forward", "backwards", "both"]:
raise TypeError("Unexpected arg for scanDirection: {}".format(scanDirection))
for signal in signalHists:
if scanDirection == "forward":
pad.drawables.append(self.getSignificanceHistogram(signal, self.getTotalBG(), **kwargs))
elif scanDirection == "backwards":
pad.drawables.append(self.getSignificanceHistogram(signal, self.getTotalBG(), forward=False, **kwargs))
elif scanDirection == "both":
pad.drawables.append(self.getSignificanceHistogram(signal, self.getTotalBG(), **kwargs))
pad.drawables.append(self.getSignificanceHistogram(signal, self.getTotalBG(), forward=False, overrideLineStyle=3, **kwargs))
if self.ratioUp is not None:
pad.yMaximum = self.ratioUp
if self.ratioDown is not None:
pad.yMinimum = self.ratioDown
[docs] def addCumulativeStack(self, pad, **kwargs):
"""Adds cumulative distributions for all backgrounds and signals to
the given pad. The kwargs are passed to getCumulativeStack"""
for signal in self.getHists("signal"):
cSignal = self.getCumulativeStack(signal, **kwargs)
pad.drawables.append(cSignal)
backgroundStack = self.getCumulativeStack(*self.getHists("backgrounds"), **kwargs)
pad.drawables.append(backgroundStack)
[docs] def getRatioHistogram(self, num, den,
drawArrows=True, mode="rawpois",
ignoreDenErrors=True, ignoreNumErrors=False):
"""create a ratio histogram.
:param num: numerator histogram
:param den: denominator histogram
:param drawArrows: draw Arrows when the ratio points reach out of the pad
:param mode: Mode for calculating the ratio (see options listed below)
:param ignoreDenErrors: Ignore the denominator errors
:param ignoreNumErrors: Ignore the numerator errors
The following modes are available:
* "rawpois" (default) : Use the raw yield to calculate poisson errors
* "pois" : Use `ROOT TGraphAsymmErrors::Divide <https://root.cern.ch/doc/master/classTGraphAsymmErrors.html#a37a202762b286cf4c7f5d34046be8c0b>`_
with option pois
* "hist" : Use `ROOT TH1::Divide <https://root.cern.ch/doc/master/classTH1.html#a4ef2329285beb6091515b836e160bd9f>`_
"""
arrows = []
if ignoreNumErrors:
num = getHM(num.Clone(next(PH.tempNames)))
# ToDo: maybe reimplement the ratio for this case (less hacky)
num.truncateErrors(0.0001)
logger.warn("ignoring numerator errors - this is not implemented perfectly - might "
"lead to strange results")
if ignoreDenErrors:
den = getHM(den.Clone(next(PH.tempNames)))
# ToDo: maybe reimplement the ratio for this case (less hacky)
den.truncateErrors(0.0001)
logger.debug("ignoring denominator errors - this is not implemented perfectly - "
"might lead to strange results")
if mode == "rawpois":
ae = getPoissonRatio(num, den,
ignoreDenErrors=ignoreDenErrors,
ignoreNumErrors=ignoreNumErrors)
elif mode == "pois":
ae = AE()
ae.Divide(num, den, "pois")
elif mode == "hist":
# not an ae in this case but who cares ;)
ae = num.clone()
ae.Divide(den)
# this is not working in that case
drawArrows = False
elif mode == "binomial":
ae = num.clone()
ae.Divide(num, den, 1, 1, "B")
ae.style = "data"
drawArrows = False
elif mode == "bayes":
ae = AE()
ae.Divide(num, den)
drawArrows = False
else:
raise ValueError("Mode {} not supported".format(mode))
#ae.xTitle = num.getXTitle()
ae.data = True
if drawArrows:
ae.addRatioArrows(self.ratioUp, self.ratioDown)
return ae
[docs] @staticmethod
def getSignificanceHistogram(signal, background,
mode="BinomialExpZ",
customFunction=None,
forward=True,
overrideLineStyle=None):
"""Create histogram which shows for each bin the significance when
cutting on the value on the x-axis.
:param mode: mode for signficance calculation:
:param customFunction: if given, use this function for
significance calculation (takes s, b, db and returns z)
:param forward: integrate to the right? (default: True)
Possible significance modes (in both cases the total bkg uncertainty is taken from the histogram):
* *BinomialExpZ*: Binomial Significance using `RooStats::NumberCountingUtils::BinomialExpZ <https://root.cern.ch/doc/master/namespaceRooStats_1_1NumberCountingUtils.html#acc62261fa2a9671cf80a70285133ff40>`_
* *PoissonAsimov*: using the significance based on asymptotic profile likelihood with poisson constraint (see `<http://www.pp.rhul.ac.uk/~cowan/stat/medsig/medsigNote.pdf>`_)
* *PoissonAsimovCLs*: same assumptions as PoissonAsimov, but applied for CLs (exclusion) - see :py:meth:`~MPF.pyrootHelpers.getExpCLsAsimov`
"""
if customFunction is not None:
significance = customFunction
elif mode == "BinomialExpZ":
significance = lambda s, b, db : ROOT.RooStats.NumberCountingUtils.BinomialExpZ(s, b, db/b)
elif mode == "PoissonAsimov":
significance = lambda s, b, db : math.sqrt(2*((s+b)*math.log(((s+b)*(b+db**2))/(b**2+(s+b)*db**2))-(b**2)/(db**2)*math.log(1+(db**2*s)/(b*(b+db**2)))))
elif mode == "PoissonAsimovCLs":
significance = lambda s, b, db : ROOT.RooStats.PValueToSignificance(PH.getExpCLsAsimov(s, b, db))
else:
raise ValueError("Significance mode {} unknown!".format(mode))
zHist = signal.clone()
zHist.Reset()
for i in range(1, background.GetNbinsX()+1):
if forward:
imin = i
imax = background.GetNbinsX()+1
else:
imin = 0
imax = i
s = signal.Integral(imin, imax)
db = ROOT.Double()
b = background.IntegralAndError(imin, imax, db)
db = float(db)
try:
z = significance(s, b, db)
except (ZeroDivisionError, ValueError) as e:
logger.warning("Division by zero for significance at bin {} (s, b, db = {}, {}, {}) - setting to 0".format(zHist.GetBinLowEdge(i), s, b, db))
z = 0
if abs(z) == float('inf'):
logger.warning("Got infinite significance at bin {} (s, b, db = {}, {}, {}) - setting to 0".format(zHist.GetBinLowEdge(i), s, b, db))
z = 0
zHist.SetBinContent(i, z)
zHist.style = 'signal'
if overrideLineStyle is not None:
zHist.lineStyle = overrideLineStyle
zHist.SetFillStyle(0) # To follow up why I have to do this here...
return zHist
[docs] @staticmethod
def getCumulativeStack(*hists, **kwargs):
"""Create a stack of cumulative histograms. If only one hist
is given, no stack is created (one cumulative histogram
instead)
:param hists: Histograms to be stacked
:param forward: Integrate to the right? (default: True)
"""
forward = kwargs.pop("forward", True)
if kwargs:
raise KeyError("got unexpected kwargs: {}".format(kwargs))
if len(hists) > 1:
cStack = getHM(ROOT.THStack('cStack{}'.format(next(PH.tempNames)), ''))
for hist in hists:
cHist = getHM(hist.GetCumulative(forward, next(PH.tempNames)))
Histogram.cloneAttributes(cHist, hist)
if len(hists) < 2:
cHist = getHM(cHist)
Histogram.cloneAttributes(cHist, hist)
return cHist
else:
cStack.add(cHist)
return cStack