Skip to content

Instantly share code, notes, and snippets.

@giorgiopizz

giorgiopizz/LatinosUtils.py

Created May 11, 2023 14:35
Show Gist options
  • Save giorgiopizz/f05ef9a6c4b3f3996307fc15ab30d1d7 to your computer and use it in GitHub Desktop.
Save giorgiopizz/f05ef9a6c4b3f3996307fc15ab30d1d7 to your computer and use it in GitHub Desktop.
Download ZIP
mkPlot.py and mkDatacards.py port to mkShapesRDF
Raw
LatinosUtils.py
def flatten_samples(samples):
## flatten the subsamples (create full samples named sample_subsample)
subsamplesmap = []
for sname in list(samples.keys()):
sample = samples[sname]
if 'subsamples' not in sample:
continue
subsamplesmap.append((sname, []))
for sub in sample['subsamples']:
samples['%s_%s' % (sname, sub)] = sample
subsamplesmap[-1][1].append('%s_%s' % (sname, sub))
sample.pop('subsamples')
samples.pop(sname)
return subsamplesmap
def flatten_cuts(cuts):
## flatten the categories (create full cuts named cut_category)
categoriesmap = []
for cname in list(cuts.keys()):
cut = cuts[cname]
if 'categories' not in cut:
continue
categoriesmap.append((cname, []))
for cat in cut['categories']:
cuts['%s_%s' % (cname, cat)] = cut
categoriesmap[-1][1].append('%s_%s' % (cname, cat))
cut.pop('categories')
cuts.pop(cname)
return categoriesmap
def update_variables_with_categories(variables, categoriesmap):
## variables can have "cuts" specifications
for variable in variables.items():
if 'cuts' not in variable:
continue
cutspec = variable['cuts']
for cname, categories in categoriesmap:
if cname not in cutspec:
continue
# original cut is in the spec
cutspec.remove(cname)
# if a category (subcut) is also in the spec, we won't touch this variable
if len(set(cutspec) & set(categories)) != 0:
continue
# otherwise we replace the cut with all the categories
cutspec.extend(categories)
def update_nuisances_with_subsamples(nuisances, subsamplesmap):
for nuisance in nuisances.items():
if 'samples' not in nuisance:
continue
samplespec = nuisance['samples']
for sname, subsamples in subsamplesmap:
if sname not in samplespec:
continue
# original sample is in the spec
sconfig = samplespec.pop(sname)
# if a subsample is also in the spec, we won't toucn this any more
if len(set(samplespec.keys()) & set(subsamples)) != 0:
continue
# otherwise we replace the sample with all the subsamples
samplespec.update((subsample, sconfig) for subsample in subsamples)
def update_nuisances_with_categories(nuisances, categoriesmap):
for nuisance in nuisances.values():
if 'cuts' not in nuisance:
continue
cutspec = nuisance['cuts']
for cname, categories in categoriesmap:
if cname not in cutspec:
continue
# original cut is in the spec
cutspec.remove(cname)
# if a category (subcut) is also in the spec, we won't touch this nuisance
if len(set(cutspec) & set(categories)) != 0:
continue
# otherwise we replace the cut with all the categories
cutspec.extend(categories)
Raw
mkDatacards.py
#!/usr/bin/env python
import json
import sys
argv = sys.argv
sys.argv = argv[:1]
import ROOT
ROOT.gROOT.SetBatch(True)
import optparse
#import LatinoAnalysis.Gardener.hwwtools as hwwtools
#import logging
import collections
import os.path
import shutil
# Common Tools & batch
#from LatinoAnalysis.Tools.commonTools import *
# ----------------------------------------------------- DatacardFactory --------------------------------------
class DatacardFactory:
#_logger = logging.getLogger('DatacardFactory')
# _____________________________________________________________________________
def __init__(self):
self._fileIn = None
self._skipMissingNuisance = False
# _____________________________________________________________________________
# a datacard for each "cut" and each "variable" will be produced, in separate sub-folders, names after "cut/variable"
# _____________________________________________________________________________
def makeDatacards( self, inputFile, outputDirDatacard, variables, cuts, samples, structureFile, nuisances):
print("=======================")
print("==== makeDatacards ====")
print("=======================")
if os.path.isdir(inputFile):
# ONLY COMPATIBLE WITH OUTPUTS MERGED TO SAMPLE LEVEL!!
self._fileIn = {}
for sampleName in samples:
self._fileIn[sampleName] = ROOT.TFile.Open(inputFile+'/plots_%s_ALL_%s.root' % (self._tag, sampleName))
if not self._fileIn[sampleName]:
raise RuntimeError('Input file for sample ' + sampleName + ' missing')
else:
self._fileIn = ROOT.TFile(inputFile, "READ")
ROOT.gROOT.GetListOfFiles().Remove(self._fileIn) #Speed up
# categorize the sample names
signal_ids = collections.OrderedDict() # id map of alternative_signals + cumulative_signals
alternative_signals = ['']
cumulative_signals = []
data = []
background_ids = collections.OrderedDict()
# divide the list of samples among signal, background and data
isig = 0
ibkg = 1
for sampleName in samples:
if structureFile[sampleName]['isSignal'] != 0:
signal_ids[sampleName] = isig
isig -= 1
if structureFile[sampleName]['isSignal'] == 2 :
alternative_signals.append(sampleName)
if structureFile[sampleName]['isSignal'] == 1 :
cumulative_signals.append(sampleName)
if structureFile[sampleName]['isData'] == 1 :
data.append(sampleName)
if structureFile[sampleName]['isSignal'] == 0 and structureFile[sampleName]['isData'] == 0:
background_ids[sampleName] = ibkg
ibkg += 1
print("Number of Signals: " + str(len(signal_ids)))
print("Number of Alternative Signals: " + str(len(alternative_signals) - 1))
print("Number of Cumulative Signals: " + str(len(cumulative_signals)))
print("Number of Backgrounds: " + str(len(background_ids)))
if not os.path.isdir(outputDirDatacard + "/") :
os.mkdir(outputDirDatacard + "/")
# loop over cuts. One directory per cut will be created
for cutName in cuts:
print("cut = ", cutName)
try:
shutil.rmtree(outputDirDatacard + "/" + cutName)
except OSError:
pass
os.mkdir(outputDirDatacard + "/" + cutName)
#
# prepare the signals and background list of samples
# after removing the ones not to be used in this specific phase space
#
cut_signals = []
cut_backgrounds = []
for snameList, idmap in [(cut_signals, signal_ids), (cut_backgrounds, background_ids)]:
for sampleName in idmap:
if 'removeFromCuts' in structureFile[sampleName] and cutName in structureFile[sampleName]['removeFromCuts']:
# remove from the list
print(' remove ', sampleName, ' from ', cutName)
else:
snameList.append(sampleName)
print(" sampleNames = ", cut_signals + cut_backgrounds)
# loop over variables
for variableName, variable in variables.items():
if 'cuts' in variable and cutName not in variable['cuts']:
continue
print(" variableName = ", variableName)
tagNameToAppearInDatacard = cutName
# e.g. hww2l2v_13TeV_of0j
# to be defined in cuts.py
os.mkdir(outputDirDatacard + "/" + cutName + "/" + variableName)
os.mkdir(outputDirDatacard + "/" + cutName + "/" + variableName + "/shapes/") # and the folder for the root files
self._outFile = ROOT.TFile.Open(outputDirDatacard + "/" + cutName + "/" + variableName + "/shapes/histos_" + tagNameToAppearInDatacard + ".root", 'recreate')
ROOT.gROOT.GetListOfFiles().Remove(self._outFile)
# prepare yields
yieldsSig = {}
yieldsBkg = {}
yieldsData = {}
# Bin to be killed
#killBinSig = {}
#killBinBkg = {}
for snameList, yields in [(cut_signals, yieldsSig), (cut_backgrounds, yieldsBkg), (data, yieldsData)]:
for sampleName in snameList:
histo = self._getHisto(cutName, variableName, sampleName)
histo.SetDirectory(self._outFile)
# get the integral == rate from histogram
if structureFile[sampleName]['isData'] == 1:
yields['data'] = histo.Integral()
histo.SetName("histo_Data")
else:
# if 'removeNegNomVal' in structureFile[sampleName] and structureFile[sampleName]['removeNegNomVal'] :
# for iBin in range(0,histo.GetNbinsX()+2) :
# BinContent = histo.GetBinContent(iBin)
# if BinContent < 0.1 and not BinContent == 0:
# print('Killing Bin :' , sampleName , iBin , BinContent)
# if not sampleName in killBinSig : killBinSig[sampleName] = []
# killBinSig[sampleName].append(iBin)
# histo.SetBinContent(iBin,0.)
if 'scaleSampleForDatacard' in structure[sampleName]:
scaleFactor = 1.
if type(structure[sampleName]['scaleSampleForDatacard']) is dict:
try:
scaleFactor = structure[sampleName]['scaleSampleForDatacard'][cutName]
except:
pass
if type(structure[sampleName]['scaleSampleForDatacard']) is int or type(structure[sampleName]['scaleSampleForDatacard']) is float:
scaleFactor = structure[sampleName]['scaleSampleForDatacard']
histo.Scale(scaleFactor)
yields[sampleName] = histo.Integral()
#
# Remove statistical uncertainty from the histogram
# This may be used to suppress the effect on "autoMCstat" of a specific
# sample, by means of putting its uncertainty to 0
#
if 'removeStatUnc' in structureFile[sampleName] and structureFile[sampleName]['removeStatUnc']:
self._removeStatUncertainty (histo)
self._outFile.cd()
histo.Write()
# Loop over alternative signal samples. One card per signal (there is always at least one entry (""))
for signalName in alternative_signals:
if signalName == '':
signals = [name for name in cumulative_signals if name in cut_signals]
else:
if signalName not in cut_signals:
continue
signals = [signalName] + [name for name in cumulative_signals if name in cut_signals]
processes = signals + cut_backgrounds
print(" Alternative signal: " + str(signalName))
alternativeSignalName = str(signalName)
alternativeSignalTitle = ""
if alternativeSignalName != "":
alternativeSignalTitle = "_" + str(signalName)
# start creating the datacard
cardPath = outputDirDatacard + "/" + cutName + "/" + variableName + "/datacard" + alternativeSignalTitle + ".txt"
print(' Writing to ' + cardPath )
card = open( cardPath ,"w")
print(' headers..')
card.write('## Shape input card\n')
card.write('imax 1 number of channels\n')
card.write('jmax * number of background\n')
card.write('kmax * number of nuisance parameters\n')
card.write('-'*100+'\n')
card.write('bin %s' % tagNameToAppearInDatacard+'\n')
if len(data) == 0:
self._logger.warning( 'no data, no fun! ')
#raise RuntimeError('No Data found!')
yieldsData['data'] = 0
card.write('observation %.0f\n' % yieldsData['data'])
card.write('shapes * * '+
'shapes/histos_' + tagNameToAppearInDatacard + ".root" +
' histo_$PROCESS histo_$PROCESS_$SYSTEMATIC' + '\n')
card.write('shapes data_obs * '+
'shapes/histos_' + tagNameToAppearInDatacard + ".root" +
' histo_Data' + '\n')
# shapes * * shapes/hww-19.36fb.mH125.of_vh2j_shape_mll.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC
# shapes data_obs * shapes/hww-19.36fb.mH125.of_vh2j_shape_mll.root histo_Data
columndef = 30
# adapt column length to long bin names
if len(tagNameToAppearInDatacard) >= (columndef -5) :
columndef = len(tagNameToAppearInDatacard) + 7
for name in signals :
if len(name)>= (columndef -5) :
columndef = len(name) + 7
print(' processes and rates..')
card.write('bin'.ljust(80))
card.write(''.join([tagNameToAppearInDatacard.ljust(columndef)] * (len(signals) + len(cut_backgrounds)))+'\n')
card.write('process'.ljust(80))
card.write(''.join(name.ljust(columndef) for name in signals))
card.write(''.join(name.ljust(columndef) for name in cut_backgrounds))
card.write('\n')
card.write('process'.ljust(80))
card.write(''.join(('%d' % signal_ids[name]).ljust(columndef) for name in signals))
card.write(''.join(('%d' % background_ids[name]).ljust(columndef) for name in cut_backgrounds))
card.write('\n')
card.write('rate'.ljust(80))
card.write(''.join(('%-.4f' % yieldsSig[name]).ljust(columndef) for name in signals))
card.write(''.join(('%-.4f' % yieldsBkg[name]).ljust(columndef) for name in cut_backgrounds))
card.write('\n')
card.write('-'*100+'\n')
print(' nuisances..')
nuisanceGroups = collections.defaultdict(list)
# add normalization and shape nuisances
for nuisanceName, nuisance in nuisances.items():
if 'type' not in nuisance:
raise RuntimeError('Nuisance ' + nuisanceName + ' is missing the type specification')
if nuisanceName == 'stat' or nuisance['type'] == 'rateParam':
# will deal with these later
continue
# check if a nuisance can be skipped because not in this particular cut
if 'cuts' in nuisance and cutName not in nuisance['cuts']:
continue
if nuisance['type'] in ['lnN', 'lnU']:
# why is adding CMS_ not the default for lnN/lnU? (Y.I. 2019.11.06)
entryName = nuisance['name']
if 'perRecoBin' in nuisance.keys() and nuisance['perRecoBin'] == True:
entryName += "_"+cutName
card.write(entryName.ljust(80-20))
if 'AsShape' in nuisance and float(nuisance['AsShape']) >= 1.:
print(">>>>>", nuisance['name'], " was derived as a lnN uncertainty but is being treated as a shape")
card.write(('shape').ljust(20))
for sampleName in processes:
if 'cuts_samples' in nuisance and sampleName in nuisance['cuts_samples'] and cutName not in nuisance['cuts_samples'][sampleName]:
# If the cuts_samples options is there and the sample is inserted in the dictionary
# check if the current cutName is included. Excluded the cuts not in the list
print("Removing nuisance ", nuisanceName, " for sample ", sampleName, " from cut ", cutName)
card.write(('-').ljust(columndef))
elif ('all' in nuisance and nuisance['all'] == 1) or \
('samples' in nuisance and sampleName in nuisance['samples']):
histo = self._getHisto(cutName, variableName, sampleName)
histoUp = histo.Clone('%s_%sUp' % (histo.GetName(), nuisance['name']))
histoDown = histo.Clone('%s_%sDown' % (histo.GetName(), nuisance['name']))
histoUp.SetDirectory(self._outFile)
histoDown.SetDirectory(self._outFile)
if 'scaleSampleForDatacard' in structure[sampleName]:
scaleFactor = 1.
if type(structure[sampleName]['scaleSampleForDatacard']) is dict:
try:
scaleFactor = structure[sampleName]['scaleSampleForDatacard'][cutName]
except:
pass
if type(structure[sampleName]['scaleSampleForDatacard']) is int or type(structure[sampleName]['scaleSampleForDatacard']) is float:
scaleFactor = structure[sampleName]['scaleSampleForDatacard']
histoUp.Scale(scaleFactor)
histoDown.Scale(scaleFactor)
if '/' in nuisance['samples'][sampleName]:
up, down = nuisance['samples'][sampleName].split('/')
histoUp.Scale(float(up))
histoDown.Scale(float(down))
else:
histoUp.Scale(float(nuisance['samples'][sampleName]))
histoDown.Scale(1. / float(nuisance['samples'][sampleName]))
self._outFile.cd()
histoUp.Write()
histoDown.Write()
card.write('1.000'.ljust(columndef))
else:
card.write(('-').ljust(columndef))
else:
card.write((nuisance['type']).ljust(20))
if 'all' in nuisance and nuisance ['all'] == 1: # for all samples
card.write(''.join(('%-.4f' % nuisance['value']).ljust(columndef) for _ in processes))
else:
# apply only to selected samples
for sampleName in processes:
if sampleName in nuisance['samples']:
# in this case nuisance['samples'] is a dict mapping sample name to nuisance values in string
card.write(('%s' % nuisance['samples'][sampleName]).ljust(columndef))
else:
card.write(('-').ljust(columndef))
elif nuisance['type'] == 'shape':
#
# 'skipCMS' is a flag not to introduce "CMS" automatically in the name
# of the nuisance.
# It may be needed for combination purposes with ATLAS
# and agreed upon for all CMS analyses.
# For example: theoretical uncertainties on ggH with migration scheme 2017
#
if 'skipCMS' in nuisance and nuisance['skipCMS'] == 1:
entryName = nuisance['name']
else:
entryName = 'CMS_' + nuisance['name']
if 'perRecoBin' in nuisance.keys() and nuisance['perRecoBin'] == True:
entryName += "_"+cutName
card.write(entryName.ljust(80-20))
if 'AsLnN' in nuisance and float(nuisance['AsLnN']) >= 1.:
print(">>>>>", nuisance['name'], " was derived as a shape uncertainty but is being treated as a lnN")
card.write(('lnN').ljust(20))
for sampleName in processes:
if 'cuts_samples' in nuisance and sampleName in nuisance['cuts_samples'] and cutName not in nuisance['cuts_samples'][sampleName]:
# If the cuts_samples options is there and the sample is inserted in the dictionary
# check if the current cutName is included. Excluded the cuts not in the list
print("Removing nuisance ", nuisanceName, " for sample ", sampleName, " from cut ", cutName)
card.write(('-').ljust(columndef))
elif ('all' in nuisance and nuisance['all'] == 1) or \
('samples' in nuisance and sampleName in nuisance['samples']):
histo = self._getHisto(cutName, variableName, sampleName)
histoUp = self._getHisto(cutName, variableName, sampleName, '_' + nuisance['name'] + 'Up')
histoDown = self._getHisto(cutName, variableName, sampleName, '_' + nuisance['name'] + 'Down')
if (not histoUp or not histoDown):
print('Histogram for nuisance', nuisance['name'], 'on sample', sampleName, 'missing')
if self._skipMissingNuisance:
card.write(('-').ljust(columndef))
continue
histoIntegral = histo.Integral()
histoUpIntegral = histoUp.Integral()
histoDownIntegral = histoDown.Integral()
if histoIntegral > 0. and histoUpIntegral > 0.:
diffUp = (histoUpIntegral - histoIntegral)/histoIntegral/float(nuisance['AsLnN'])
else:
diffUp = 0.
if histoIntegral > 0. and histoDownIntegral > 0.:
diffDo = (histoDownIntegral - histoIntegral)/histoIntegral/float(nuisance['AsLnN'])
else:
diffDo = 0.
if 'symmetrize' in nuisance and nuisance['symmetrize']:
diff = (diffUp - diffDo) * 0.5
if diff >= 1.:
# can't symmetrize
diffUp = diff * 2. - 0.999
diffDo = -0.999
elif diff <= -1.:
diffUp = -0.999
diffDo = -diff * 2. - 0.999
else:
diffUp = diff
diffDo = -diff
lnNUp = 1. + diffUp
lnNDo = 1. + diffDo
#
# avoid 0.00/XXX and XXX/0.00 in the lnN --> ill defined nuisance
# if this happens put 0.00 ---> 1.00, meaning *no* effect of this nuisance
# Done like this because it is very likely what you are experiencing
# is a MC fluctuation in the varied sample, that is the up/down histogram has 0 entries!
#
# NB: with the requirement "histoUpIntegral > 0" and "histoDownIntegral > 0" this should never happen,
# except for some strange coincidence of "AsLnN" ...
# This fix is left, just for sake of safety
#
if lnNUp==0: lnNUp = 1.
if lnNDo==0: lnNDo = 1.
if abs(lnNUp - 1.) < 5.e-4 and abs(lnNDo - 1.) < 5.e-4:
card.write(('-').ljust(columndef))
else:
card.write((('%-.4f' % lnNUp)+"/"+('%-.4f' % lnNDo)).ljust(columndef))
else:
card.write(('-').ljust(columndef))
else:
card.write('shape'.ljust(20))
for sampleName in processes:
if 'cuts_samples' in nuisance and sampleName in nuisance['cuts_samples'] and cutName not in nuisance['cuts_samples'][sampleName]:
# If the cuts_samples options is there and the sample is inserted in the dictionary
# check if the current cutName is included. Excluded the cuts not in the list
print("Removing nuisance ", nuisanceName, " for sample ", sampleName, " from cut ", cutName)
card.write(('-').ljust(columndef))
elif ('all' in nuisance and nuisance ['all'] == 1) or \
('samples' in nuisance and sampleName in nuisance['samples']):
# save the nuisance histograms in the root file
if ('skipCMS' in nuisance.keys()) and nuisance['skipCMS'] == 1:
suffixOut = None
else:
suffixOut = '_CMS_' + nuisance['name']
if 'perRecoBin' in nuisance.keys() and nuisance['perRecoBin'] == True:
if ('skipCMS' in nuisance.keys()) and nuisance['skipCMS'] == 1:
suffixOut = "_"+nuisance['name']
else:
suffixOut = '_CMS_' + nuisance['name']
suffixOut += "_"+cutName
symmetrize = 'symmetrize' in nuisance and nuisance['symmetrize']
saved = self._saveNuisanceHistos(cutName, variableName, sampleName, '_' + nuisance['name'], suffixOut, symmetrize)
if saved:
card.write('1.000'.ljust(columndef))
else:
# saved can be false if skipMissingNuisance is true and histogram cannot be found
card.write(('-').ljust(columndef))
else:
card.write(('-').ljust(columndef))
# closing block if type == lnN or shape
card.write('\n')
if 'group' in nuisance:
nuisanceGroups[nuisance['group']].append(entryName)
for group in sorted(nuisanceGroups.keys()):
card.write('%s group = %s\n' % (group, ' '.join(nuisanceGroups[group])))
# done with normalization and shape nuisances.
# now add the stat nuisances
if 'stat' in nuisances:
nuisance = nuisances['stat']
if nuisance['type'] == 'auto':
# use autoMCStats feature of combine
card.write('* autoMCStats ' + nuisance['maxPoiss'] + ' ' + nuisance['includeSignal'])
# nuisance ['maxPoiss'] = Number of threshold events for Poisson modelling
# nuisance ['includeSignal'] = Include MC stat nuisances on signal processes (1=True, 0=False)
card.write('\n')
else:
# otherwise we process stat nuisances sample-by-sample
for sampleName in processes:
if sampleName not in nuisance['samples']:
continue
sampleNuisance = nuisance['samples'][sampleName]
sampleIndex = processes.index(sampleName)
if sampleNuisance['typeStat'] == 'uni' : # unified approach
# save the nuisance histograms in the root file
saved = self._saveNuisanceHistos(cutName, variableName, sampleName, 'stat', '_CMS_' + tagNameToAppearInDatacard + "_" + sampleName + "_stat")
if saved:
card.write(('CMS_' + tagNameToAppearInDatacard + "_" + sampleName + "_stat").ljust(80-20))
card.write((nuisance['type']).ljust(20))
# write the line in datacard. Only the column for this sample gets 1.000
for idx in range(len(processes)):
if idx == sampleIndex:
card.write(('1.000').ljust(columndef))
else:
card.write(('-').ljust(columndef))
card.write('\n')
elif sampleNuisance['typeStat'] == 'bbb': # bin-by-bin
histoTemplate = self._getHisto(cutName, variableName, sampleName)
for iBin in range(1, histoTemplate.GetNbinsX()+1):
if 'correlate' in sampleNuisance:
#specify the sample that is source of the variation
tag = "_ibin"+sampleName+"_"
correlate = list(sampleNuisance["correlate"]) # list of sample names to correlate bin-by-bin with this sample
else:
tag = "_ibin_"
correlate = []
# save the nuisance histograms in the root file
saved = self._saveNuisanceHistos(cutName, variableName, sampleName, tag + str(iBin) + '_stat',
'_CMS_' + tagNameToAppearInDatacard + "_" + sampleName + '_ibin_' + str(iBin) + '_stat')
if not saved:
continue
for other in list(correlate):
saved = self._saveNuisanceHistos(cutName, variableName, other, tag + str(iBin) + '_stat',
'_CMS_' + tagNameToAppearInDatacard + "_" + sampleName + '_ibin_' + str(iBin) + '_stat')
if not saved:
correlate.remove(other)
card.write(('CMS_' + tagNameToAppearInDatacard + "_" + sampleName + "_ibin_" + str(iBin) + "_stat").ljust(100-20))
card.write((nuisance['type']).ljust(20))
# write line in datacard. One line per sample per bin
for idx in range(len(processes)):
if idx == sampleIndex or processes[idx] in correlate:
card.write(('1.000').ljust(columndef))
else:
card.write(('-').ljust(columndef))
card.write('\n')
# now add the "rateParam" for the normalization
# e.g.: z_norm rateParam htsearch zll 1
# see: https://twiki.cern.ch/twiki/bin/view/CMS/HiggsWG/SWGuideNonStandardCombineUses#Rate_Parameters
# 'rateParam' has a separate treatment -> it's just a line at the end of the datacard. It defines "free floating" samples
# I do it here and not before because I want the freee floating parameters at the end of the datacard
for nuisance in nuisances.values():
if nuisance['type'] != 'rateParam':
continue
# check if a nuisance can be skipped because not in this particular cut
if 'cuts' in nuisance and cutName not in nuisance['cuts']:
continue
card.write(nuisance['name'].ljust(80-20))
card.write('rateParam'.ljust(25))
card.write(tagNameToAppearInDatacard.ljust(columndef)) # the bin
# there can be only 1 sample per rateParam
if len(nuisance['samples']) != 1:
raise RuntimeError('Invalid rateParam: number of samples != 1')
sampleName, initialValue = list(nuisance['samples'].items())[0]
if sampleName not in samples:
raise RuntimeError('Invalid rateParam: unknown sample %s' % sampleName)
card.write(sampleName.ljust(25))
card.write(('%-.4f' % float(initialValue)).ljust(columndef))
card.write('\n')
# now add other nuisances
# Are there other kind of nuisances I forgot?
card.write('-'*100+'\n')
card.write('\n')
card.close()
print(' done.')
self._outFile.Close()
if type(self._fileIn) is dict:
for source in self._fileIn.values():
source.Close()
else:
self._fileIn.Close()
# _____________________________________________________________________________
def _saveNuisanceHistos(self, cutName, variableName, sampleName, suffixIn, suffixOut = None, symmetrize = False):
histoUp = self._getHisto(cutName, variableName, sampleName, suffixIn + 'Up')
histoDown = self._getHisto(cutName, variableName, sampleName, suffixIn + 'Down')
if not histoUp or not histoDown:
print('Up/down histogram for', cutName, variableName, sampleName, suffixIn, 'missing')
if self._skipMissingNuisance:
return False
# else let ROOT raise
if 'scaleSampleForDatacard' in structure[sampleName]:
scaleFactor = 1.
if type(structure[sampleName]['scaleSampleForDatacard']) is dict:
try:
scaleFactor = structure[sampleName]['scaleSampleForDatacard'][cutName]
except:
pass
if type(structure[sampleName]['scaleSampleForDatacard']) is int or type(structure[sampleName]['scaleSampleForDatacard']) is float:
scaleFactor = structure[sampleName]['scaleSampleForDatacard']
histoUp.Scale(scaleFactor)
histoDown.Scale(scaleFactor)
if symmetrize:
histoNom = self._getHisto(cutName, variableName, sampleName)
histoDiff = histoUp.Clone('histoDiff')
histoDiff.Add(histoDown, -1.)
histoDiff.Scale(0.5)
histoUp.Reset()
histoUp.Add(histoNom)
histoUp.Add(histoDiff)
histoDown.Reset()
histoDown.Add(histoNom)
histoDown.Add(histoDiff, -1.)
histoUp.SetDirectory(self._outFile)
histoDown.SetDirectory(self._outFile)
if suffixOut:
histoUp.SetName('histo_%s%sUp' % (sampleName, suffixOut))
histoDown.SetName('histo_%s%sDown' % (sampleName, suffixOut))
self._outFile.cd()
histoUp.Write()
histoDown.Write()
return True
# _____________________________________________________________________________
def _getHisto(self, cutName, variableName, sampleName, suffix = None):
shapeName = '%s/%s/histo_%s' % (cutName, variableName, sampleName)
if suffix:
shapeName += suffix
if type(self._fileIn) is dict:
# by-sample ROOT file
histo = self._fileIn[sampleName].Get(shapeName)
else:
# Merged single ROOT file
histo = self._fileIn.Get(shapeName)
if not histo:
print(shapeName, 'not found')
if 'met' in variableName.lower() and 'met' in suffix.lower():
print(variableName, 'does not contain', suffix)
sys.exit()
if suffix:
print('Getting the nominal instead of varied')
return self._getHisto(cutName, variableName, sampleName)
return histo
# _____________________________________________________________________________
def _removeStatUncertainty(self, histo):
for iBin in range(0,histo.GetNbinsX()+2) :
histo.SetBinError(iBin,0.)
class list_maker:
def __init__(self, var, sep=',', type=None ):
self._type= type
self._var = var
self._sep = sep
def __call__(self,option, opt_str, value, parser):
if not hasattr(parser.values,self._var):
setattr(parser.values,self._var,[])
try:
array = value.split(self._sep)
if self._type:
array = [ self._type(e) for e in array ]
setattr(parser.values, self._var, array)
except:
print('Malformed option (comma separated list expected):',value)
def main():
sys.argv = argv
print('''
--------------------------------------------------------------------------------------------------
__ \ | | \ | |
| | _` | __| _` | __| _` | __| _` | |\/ | _` | | / _ \ __|
| | ( | | ( | ( ( | | ( | | | ( | < __/ |
____/ \__,_| \__| \__,_| \___| \__,_| _| \__,_| _| _| \__,_| _|\_\ \___| _|
--------------------------------------------------------------------------------------------------
''')
usage = 'usage: %prog [options]'
parser = optparse.OptionParser(usage)
parser.add_option('--tag' , dest='tag' , help='Tag used for the shape file name' , default=None)
parser.add_option('--sigset' , dest='sigset' , help='Signal samples [SM]' , default='SM')
parser.add_option('--outputDirDatacard' , dest='outputDirDatacard' , help='output directory' , default='./datacards')
parser.add_option('--inputFile' , dest='inputFile' , help='input directory' , default='./input.root')
parser.add_option('--structureFile' , dest='structureFile' , help='file with datacard configurations' , default=None )
parser.add_option('--nuisancesFile' , dest='nuisancesFile' , help='file with nuisances configurations' , default=None )
parser.add_option('--cardList' , dest="cardList" , help="List of cuts to produce datacards" , default=[], type='string' , action='callback' , callback=list_maker('cardList',','))
parser.add_option('--skipMissingNuisance', dest='skipMissingNuisance', help="Don't write nuisance lines when histograms are missing", default=False, action='store_true')
# read default parsing options as well
# hwwtools.addOptions(parser)
# hwwtools.loadOptDefaults(parser)
(opt, args) = parser.parse_args()
#sys.argv.append( '-b' )
#print(" configuration file = ", opt.pycfg)
global cuts, plot
import glob
import json
from mkShapesRDF.shapeAnalysis.ConfigLib import ConfigLib
lFile = list(filter(lambda k: os.path.isfile(k), glob.glob('configs/*json')))
lFile.sort(key= lambda x: os.path.getmtime(x))
lFile = lFile[-1]
with open(lFile) as file:
d = json.load(file)
ConfigLib.loadDict(d, globals())
cuts = cuts['cuts']
groupPlot = plot['groupPlot']
legend = plot['legend']
plot = plot['plot']
inputFile = outputFolder + '/' + outputFile
# print(" inputFile = ", opt.inputFile)
# print(" outputDirDatacard = ", opt.outputDirDatacard)
# if not opt.debug:
# pass
# elif opt.debug == 2:
# print('Logging level set to DEBUG (%d)' % opt.debug)
# logging.basicConfig( level=logging.DEBUG )
# elif opt.debug == 1:
# print('Logging level set to INFO (%d)' % opt.debug)
# logging.basicConfig( level=logging.INFO )
# if opt.nuisancesFile == None :
# print(" Please provide the nuisances structure if you want to add nuisances ")
#
# if opt.structureFile == None :
# print(" Please provide the datacard structure ")
# exit ()
ROOT.TH1.SetDefaultSumw2(True)
factory = DatacardFactory()
factory._tag = tag
factory._skipMissingNuisance = opt.skipMissingNuisance
# ## load the samples
# ## Use flag to skip loading samples files, as currently done in mkPlot.py
# _samples_noload = True
# samples = {}
# if os.path.exists(opt.samplesFile) :
# handle = open(opt.samplesFile,'r')
# exec(handle)
# handle.close()
#
# ## load the cuts
# cuts = {}
# if os.path.exists(opt.cutsFile) :
# handle = open(opt.cutsFile,'r')
# exec(handle)
# handle.close()
#
# ## load the variables
# variables = {}
# if os.path.exists(opt.variablesFile) :
# handle = open(opt.variablesFile,'r')
# exec(handle)
# handle.close()
#
# ## load the nuisances
# nuisances = collections.OrderedDict()
# if os.path.exists(opt.nuisancesFile) :
# handle = open(opt.nuisancesFile,'r')
# exec(handle)
# handle.close()
import mkShapesRDF.shapeAnalysis.LatinosUtils as utils
subsamplesmap = utils.flatten_samples(samples)
categoriesmap = utils.flatten_cuts(cuts)
utils.update_variables_with_categories(variables, categoriesmap)
utils.update_nuisances_with_subsamples(nuisances, subsamplesmap)
utils.update_nuisances_with_categories(nuisances, categoriesmap)
#sys.exit()
# ## load the structure file (use flattened sample and cut names)
# structure = collections.OrderedDict()
# if os.path.exists(opt.structureFile) :
# handle = open(opt.structureFile,'r')
# exec(handle)
# handle.close()
## command-line cuts restrictions
if len(opt.cardList)>0:
try:
newCuts = []
for iCut in opt.cardList:
for iOptim in optim:
newCuts.append(iCut+'_'+iOptim)
opt.cardList = newCuts
print(opt.cardList)
except:
print("No optim dictionary")
cut2del = []
for iCut in cuts:
if not iCut in opt.cardList : cut2del.append(iCut)
for iCut in cut2del : del cuts[iCut]
factory.makeDatacards( inputFile , opt.outputDirDatacard, variables, cuts, samples, structure, nuisances)
if __name__ == '__main__':
main()
Raw
mkPlot.py
#!/usr/bin/env python
import sys
argv = sys.argv
sys.argv = argv[:1]
import optparse
#import LatinoAnalysis.Gardener.hwwtools as hwwtools
import os
import os.path
#import logging
from collections import OrderedDict
from multiprocessing import Process
#from multiprocessing import Process
#import os.path
# X. Janssen - 21 March 2018
# PlotFactory splitted from this file and moved to python to be able to use in other scripts (mkACPlot.py)
#
from mkShapesRDF.shapeAnalysis.PlotFactory import PlotFactory
if __name__ == '__main__':
main()
def main():
sys.argv = argv
usage = 'usage: %prog [options]'
parser = optparse.OptionParser(usage)
parser.add_option('--scaleToPlot' , dest='scaleToPlot' , help='scale of maxY to maxHistoY' , default=3.0 , type=float )
parser.add_option('--minLogC' , dest='minLogC' , help='min Y in log plots' , default=0.01 , type=float )
parser.add_option('--maxLogC' , dest='maxLogC' , help='max Y in log plots' , default=100 , type=float )
parser.add_option('--minLogCratio' , dest='minLogCratio' , help='min Y in log ratio plots' , default=0.001 , type=float )
parser.add_option('--maxLogCratio' , dest='maxLogCratio' , help='max Y in log ratio plots' , default=10 , type=float )
parser.add_option('--maxLinearScale' , dest='maxLinearScale' , help='scale factor for max Y in linear plots (1.45 magic number as default)' , default=1.45 , type=float )
parser.add_option('--outputDirPlots' , dest='outputDirPlots' , help='output directory' , default='./')
parser.add_option('--inputFile' , dest='inputFile' , help='input file with histograms' , default='input.root')
parser.add_option('--tag' , dest='tag' , help='Tag used for the shape file name. Used if inputFile is a directory', default=None)
parser.add_option('--nuisancesFile' , dest='nuisancesFile' , help='file with nuisances configurations' , default=None )
parser.add_option('--onlyVariable' , dest='onlyVariable' , help='draw only one variable (may be needed in post-fit plots)' , default=None)
parser.add_option('--onlyCut' , dest='onlyCut' , help='draw only one cut phase space (may be needed in post-fit plots)' , default=None)
parser.add_option('--onlyPlot' , dest='onlyPlot' , help='draw only specified plot type (comma-separated c, cratio, and/or cdifference)', default=None)
parser.add_option('--linearOnly' , dest='linearOnly' , help='Make linear plot only.', action='store_true', default=False)
parser.add_option('--logOnly' , dest='logOnly' , help='Make log plot only.', action='store_true', default=False)
parser.add_option('--fileFormats' , dest='fileFormats' , help='Output plot file formats (comma-separated png, pdf, root, C, and/or eps). Default "png,root"', default='png,root')
parser.add_option('--plotNormalizedIncludeData' , dest='plotNormalizedIncludeData' , help='plot also normalized distributions for data, for shape comparison purposes', default=None )
parser.add_option('--plotNormalizedDistributions' , dest='plotNormalizedDistributions' , help='plot also normalized distributions for optimization purposes' , action='store_true' , default=None )
parser.add_option('--plotNormalizedDistributionsTHstack' , dest='plotNormalizedDistributionsTHstack' , help='plot also normalized distributions for optimization purposes, with stacked sig and bkg' , action='store_true' , default=None )
parser.add_option('--showIntegralLegend' , dest='showIntegralLegend' , help='show the integral, the yields, in the legend' , default=0, type=float )
parser.add_option('--showRelativeRatio' , dest='showRelativeRatio' , help='draw instead of data-expected, (data-expected) / expected' , action='store_true', default=False)
parser.add_option('--showDataMinusBkgOnly', dest='showDataMinusBkgOnly', help='draw instead of data-expected, data-expected background only' , action='store_true', default=False)
parser.add_option('--removeWeight', dest='removeWeight', help='Remove weight S/B for PR plots, just do the sum' , action='store_true', default=False)
parser.add_option('--invertXY', dest='invertXY', help='Invert the weighting for X <-> Y. Instead of slices along Y, do slices along X' , action='store_true', default=False)
parser.add_option('--postFit', dest='postFit', help='Plot sum of post-fit backgrounds, and the data/post-fit ratio.' , default='n')
parser.add_option('--skipMissingNuisance', dest='skipMissingNuisance', help='Do not trigger errors if a nuisance is missing. To be used with absolute care!!!' , action='store_true', default=False)
parser.add_option('--removeMCStat', dest='removeMCStat', help='Do not plot the MC statistics contribution in the uncertainty band', action='store_true', default=False)
parser.add_option('--extraLegend' , dest='extraLegend' , help='User-specified additional legend' , default=None)
parser.add_option('--customize', dest='customizeKey', help="Optional parameters for the customizations script", default=None)
parser.add_option('--plotFancy', dest='plotFancy', help='Plot fancy data - bkg plot' , action='store_true', default=False)
parser.add_option('--NoPreliminary', dest='NoPreliminary', help='Remove preliminary status in plots' , action='store_true', default=False)
parser.add_option('--RemoveAllMC', dest='RemoveAllMC', help='Remove all MC in legend' , action='store_true', default=False)
parser.add_option('--parallelPlotting', dest='parallelPlotting', help='Plot each cut in parallel' , action='store_true', default=False)
# read default parsing options as well
# hwwtools.addOptions(parser)
# hwwtools.loadOptDefaults(parser)
(opt, args) = parser.parse_args()
# ROOT.gROOT.SetBatch()
print("")
#print(" configuration file =", opt.pycfg)
#print(" lumi =", opt.lumi)
print(" inputFile =", opt.inputFile)
print(" outputDirPlots =", opt.outputDirPlots)
print(" plotNormalizedDistributions =", opt.plotNormalizedDistributions)
print(" plotNormalizedIncludeData =", opt.plotNormalizedIncludeData )
print(" plotNormalizedDistributionsTHstack =", opt.plotNormalizedDistributionsTHstack)
print(" showIntegralLegend =", opt.showIntegralLegend)
print(" scaleToPlot =", opt.scaleToPlot)
print(" minLogC =", opt.minLogC)
print(" maxLogC =", opt.maxLogC)
print(" minLogCratio =", opt.minLogCratio)
print(" maxLogCratio =", opt.maxLogCratio)
print(" showRelativeRatio =", opt.showRelativeRatio)
print(" showDataMinusBkgOnly =", opt.showDataMinusBkgOnly)
print(" removeWeight =", opt.removeWeight)
print(" invertXY =", opt.invertXY )
print(" skipMissingNuisance =", opt.skipMissingNuisance)
print(" postFit =", opt.postFit)
print(" removeMCStat =", opt.removeMCStat)
print(" plotFancy =", opt.plotFancy)
print(" NoPreliminary =", opt.NoPreliminary )
print(" RemoveAllMC =", opt.RemoveAllMC )
print(" parallelPlotting =", opt.parallelPlotting)
print("")
opt.scaleToPlot = float(opt.scaleToPlot)
opt.minLogC = float(opt.minLogC)
opt.maxLogC = float(opt.maxLogC)
opt.minLogCratio = float(opt.minLogCratio)
opt.maxLogCratio = float(opt.maxLogCratio)
# if not opt.debug:
# pass
# elif opt.debug == 2:
# print 'Logging level set to DEBUG (%d)' % opt.debug
# logging.basicConfig( level=logging.DEBUG )
# elif opt.debug == 1:
# print 'Logging level set to INFO (%d)' % opt.debug
# logging.basicConfig( level=logging.INFO )
#samples = {}
#samples = OrderedDict()
# if opt.samplesFile == None :
# print " Please provide the samples structure (not strictly needed in mkPlot, since list of samples read from plot.py) "
# elif os.path.exists(opt.samplesFile) :
# # This line is needed for mkplot not to look for samples in eos.
# # Imagine the samples have been removed in eos, but the file with histograms
# # has been already generated, there is no need to check the existence of the samples on eos
# # NB: in samples.py the function "nanoGetSampleFiles" must handle this, if needed
# _samples_noload = True
# handle = open(opt.samplesFile,'r')
# exec(handle)
# handle.close()
# cuts = {}
# if os.path.exists(opt.cutsFile) :
# handle = open(opt.cutsFile,'r')
# exec(handle)
# handle.close()
#
# variables = {}
# if os.path.exists(opt.variablesFile) :
# handle = open(opt.variablesFile,'r')
# exec(handle)
# handle.close()
#
# nuisances = {}
# if opt.nuisancesFile == None :
# print " Please provide the nuisances structure if you want to add nuisances "
# elif os.path.exists(opt.nuisancesFile) :
# handle = open(opt.nuisancesFile,'r')
# exec(handle)
# handle.close()
import mkShapesRDF.shapeAnalysis.LatinosUtils as utils
# import glob
# import json
# lFile = list(filter(lambda k: os.path.isfile(k), glob.glob('configs/*json')))
# lFile.sort(key= lambda x: os.path.getmtime(x))
# lFile = lFile[-1]
# with open(lFile) as file:
# d = json.load(file)
from mkShapesRDF.shapeAnalysis.ConfigLib import ConfigLib
global plot
global cuts
ConfigLib.loadLatestJSON('configs', globals())
#ConfigLib.loadDict(d, globals())
print(dir())
print(globals().keys())
cuts = cuts['cuts']
groupPlot = plot['groupPlot']
legend = plot['legend']
plot = plot['plot']
inputFile = outputFolder + '/' + outputFile
subsamplesmap = utils.flatten_samples(samples)
categoriesmap = utils.flatten_cuts(cuts)
utils.update_variables_with_categories(variables, categoriesmap)
utils.update_nuisances_with_subsamples(nuisances, subsamplesmap)
utils.update_nuisances_with_categories(nuisances, categoriesmap)
# check if only one cut or only one variable
# is requested, and filter th elist of cuts and variables
# using this piece of information
if opt.onlyVariable != None :
list_to_remove = []
for variableName, variable in variables.iteritems():
if variableName != opt.onlyVariable :
list_to_remove.append(variableName)
for toRemove in list_to_remove:
del variables[toRemove]
print( " variables = ", variables)
if opt.onlyCut != None :
list_to_remove = []
for cutName, cutExtended in cuts.iteritems():
if cutName not in opt.onlyCut :
list_to_remove.append(cutName)
for toRemove in list_to_remove:
del cuts[toRemove]
print( " cuts = ", cuts)
# groupPlot = OrderedDict()
# plot = {}
# legend = {}
# if os.path.exists(opt.plotFile) :
# handle = open(opt.plotFile,'r')
# exec(handle)
# handle.close()
#energy = '13TeV'
#sys.exit()
#=====================
def launch_plot(inputFile, outputDirPlots, variables, cuts, samples, plot, nuisances, legend, groupPlot):
factory = PlotFactory()
factory._tag = tag
#factory._energy = opt.energy
factory._lumi = lumi
factory._plotNormalizedDistributions = opt.plotNormalizedDistributions
factory._plotNormalizedIncludeData = opt.plotNormalizedIncludeData
factory._plotNormalizedDistributionsTHstack = opt.plotNormalizedDistributionsTHstack
factory._showIntegralLegend = opt.showIntegralLegend
if opt.onlyPlot is not None:
factory._plotsToWrite = opt.onlyPlot.split(',')
factory._plotLinear = opt.linearOnly or not opt.logOnly
factory._plotLog = opt.logOnly or not opt.linearOnly
factory._scaleToPlot = opt.scaleToPlot
factory._minLogC = opt.minLogC
factory._maxLogC = opt.maxLogC
factory._minLogCratio = opt.minLogCratio
factory._maxLogCratio = opt.maxLogCratio
factory._maxLinearScale = opt.maxLinearScale
factory._minLogCdifference = opt.minLogCratio
factory._maxLogCdifference = opt.maxLogCratio
factory._showRelativeRatio = opt.showRelativeRatio
factory._showDataMinusBkgOnly = opt.showDataMinusBkgOnly
factory._removeWeight = opt.removeWeight
factory._invertXY = opt.invertXY
factory._fileFormats = opt.fileFormats.split(',')
factory._postFit = opt.postFit
factory._removeMCStat = opt.removeMCStat
factory._plotFancy = opt.plotFancy
factory._SkipMissingNuisance = opt.skipMissingNuisance
factory._extraLegend = opt.extraLegend
factory._preliminary = not opt.NoPreliminary
factory._removeAllMC = opt.RemoveAllMC
factory.makePlot(inputFile ,outputDirPlots, variables, cuts, samples, plot, nuisances, legend, groupPlot)
#===============================
if opt.parallelPlotting:
for cut in cuts:
p = Process(target=launch_plot( inputFile , plotPath, variables, [cut], samples, plot, nuisances,legend, groupPlot) )
p.start()
else:
launch_plot( inputFile , plotPath , variables, cuts, samples,
plot, nuisances, legend, groupPlot)
Raw
PlotFactory.py
#!/usr/bin/env python
import json
import sys
from sys import exit
import ROOT
ROOT.gROOT.SetBatch(True)
import optparse
#import LatinoAnalysis.Gardener.hwwtools as hwwtools
import os.path
import string
import logging
#import LatinoAnalysis.Gardener.odict as odict
import traceback
from array import array
from collections import OrderedDict
import math
import numpy as np
#import root_numpy as rnp
import re
def hist2array(h, overflow=False, copy=True):
shift = 0
if overflow:
shift = 1
d = []
for i in range(1 - shift, h.GetNbinsX() + 1 + shift):
d.append(h.GetBinContent(i))
d = np.array(d)
return d
def np_array(_TArrayD):
return np.array(_TArrayD)
# ----------------------------------------------------- PlotFactory --------------------------------------
class PlotFactory:
_logger = logging.getLogger('PlotFactory')
# _____________________________________________________________________________
def __init__(self):
self._tag = ''
variables = {}
self._variables = variables
cuts = {}
self._cuts = cuts
samples = OrderedDict()
self._samples = samples
self._plotsToWrite = ['c', 'cratio', 'cdifference']
self._plotLinear = True
self._plotLog = True
outputDirPlots = {}
self._outputDirPlots = outputDirPlots
self._showIntegralLegend = 1
# 0 is no
self._FigNamePF = ''
self._fileFormats = ['png', 'root']
self._preliminary = True
self._removeAllMC = False
# _____________________________________________________________________________
def makePlot(self, inputFile, outputDirPlots, variables, cuts, samples, plot, nuisances, legend, groupPlot):
print ("==================")
print ("==== makePlot ====")
print ("==================")
self.defineStyle()
self._variables = variables
self._samples = samples
self._cuts = cuts
self._outputDirPlots = outputDirPlots
os.system ("mkdir " + outputDirPlots + "/")
#
# prepare plotter.html
#
text_file_html = open(self._outputDirPlots + "/" + "plotter.html", "w")
text_file_html.write("<script type=\"text/javascript\" src=\"https://root.cern.ch/js/latest/scripts/JSRootCore.js?gui\"></script>\n")
text_file_html.write("<div id=\"simpleGUI\" path=\"./\" files=\"\n")
#tcanvas = ROOT.TCanvas( "cc", "cc" , 800, 600 )
#tcanvasRatio = ROOT.TCanvas( "ccRatio", "ccRatio", 800, 800 )
## weight_X_tcanvas = ROOT.TCanvas( "weight_X_tcanvas", "weight_X_tcanvas", 800, 800 )
#weight_X_tcanvasRatio = ROOT.TCanvas( "weight_X_tcanvasRatio", "weight_X_tcanvasRatio", 800, 800 )
ROOT.TH1.SetDefaultSumw2(True)
dataColor = 1
#thsData = ROOT.THStack ("thsData", "thsData")
#thsSignal = ROOT.THStack ("thsSignal", "thsSignal")
#thsBackground = ROOT.THStack ("thsBackground","thsBackground")
#thsSignal_grouped = ROOT.THStack ("thsSignal_grouped", "thsSignal_grouped")
#thsBackground_grouped = ROOT.THStack ("thsBackground_grouped","thsBackground_grouped")
ROOT.gROOT.cd()
list_thsData = {}
list_thsSignal = {}
list_thsBackground = {}
list_thsSignal_grouped = {}
list_thsBackground_grouped = {}
list_thsSignal_grouped_normalized = {}
list_thsBackground_grouped_normalized = {}
list_tcanvas = {}
list_tcanvasRatio = {}
list_weight_X_tcanvasRatio = {}
list_tcanvasDifference = {}
list_weight_X_tcanvasDifference = {}
list_tcanvasSigVsBkg = {}
list_tcanvasSigVsBkgTHstack = {}
# standalont data-bkg plot
list_tcanvasDifference_Fancy = {}
generalCounter = 0
if os.path.isdir(inputFile):
# ONLY COMPATIBLE WITH OUTPUTS MERGED TO SAMPLE LEVEL!!
fileIn = {}
allFiles = os.listdir(inputFile)
for sampleName in self._samples:
fileIn[sampleName] = ROOT.TFile.Open(inputFile+'/plots_%s_ALL_%s.root' % (self._tag, sampleName))
if not fileIn[sampleName]:
raise RuntimeError('Input file for sample ' + sampleName + ' missing')
if os.path.exists(inputFile+'/plots_total.root'):
fileIn['total'] = ROOT.TFile.Open(inputFile+'/plots_total.root')
else:
fileIn = ROOT.TFile(inputFile, "READ")
#---- save one TCanvas for every cut and every variable
for cutName in self._cuts :
print ("cut =", cutName)
for variableName, variable in self._variables.items():
if 'cuts' in variable and cutName not in variable['cuts']:
continue
if type(fileIn) is not dict and not fileIn.GetDirectory(cutName+"/"+variableName):
continue
print("variableName =", variableName)
if not "divideByBinWidth" in variable.keys():
variable["divideByBinWidth"] = 0
tcanvas = ROOT.TCanvas( "cc" + cutName + "_" + variableName, "cc" , 800, 600 )
tcanvasRatio = ROOT.TCanvas( "ccRatio" + cutName + "_" + variableName, "ccRatio", 800, 800 )
weight_X_tcanvasRatio = ROOT.TCanvas( "weight_X_tcanvasRatio" + cutName + "_" + variableName, "weight_X_tcanvasRatio", 800, 800 )
tcanvasDifference = ROOT.TCanvas( "ccDifference" + cutName + "_" + variableName, "ccDifference", 800, 800 )
weight_X_tcanvasDifference = ROOT.TCanvas( "weight_X_tcanvasDifference" + cutName + "_" + variableName, "weight_X_tcanvasDifference", 800, 800 )
if self._plotNormalizedDistributions :
tcanvasSigVsBkg = ROOT.TCanvas( "ccSigVsBkg" + cutName + "_" + variableName, "cc" , 800, 700 )
if self._plotNormalizedDistributionsTHstack :
tcanvasSigVsBkgTHstack = ROOT.TCanvas( "ccTHstackSigVsBkg" + cutName + "_" + variableName, "cc" , 800, 600 )
list_tcanvas [generalCounter] = tcanvas
list_tcanvasRatio [generalCounter] = tcanvasRatio
list_weight_X_tcanvasRatio [generalCounter] = weight_X_tcanvasRatio
list_tcanvasDifference [generalCounter] = tcanvasDifference
list_weight_X_tcanvasDifference [generalCounter] = weight_X_tcanvasDifference
if self._plotNormalizedDistributions :
list_tcanvasSigVsBkg [generalCounter] = tcanvasSigVsBkg
if self._plotNormalizedDistributionsTHstack :
list_tcanvasSigVsBkgTHstack [generalCounter] = tcanvasSigVsBkgTHstack
tcanvasDifference_Fancy = ROOT.TCanvas( "ccDifference_Fancy" + cutName + "_" + variableName, "ccDifference_Fancy", 800, 800 )
list_tcanvasDifference_Fancy [generalCounter] = tcanvasDifference_Fancy
histos = {}
histos_grouped = {}
canvasNameTemplateRatio = 'ccRatio_' + cutName + "_" + variableName
canvasNameTemplateDifference = 'ccDifference_' + cutName + "_" + variableName
#tcanvasRatio = ROOT.TCanvas( canvasNameTemplateRatio, variableName, 800, 800 )
canvasNameTemplate = 'c_' + cutName + "_" + variableName
#tcanvas = ROOT.TCanvas( canvasNameTemplate, variableName , 800, 600 )
tcanvas.cd()
#print(" and now this ...")
tgrData_vx = array('f')
tgrData_evx = array('f')
tgrData_vy = array('f')
tgrData_evy_up = array('f')
tgrData_evy_do = array('f')
# at least 1 "MC" should be around ... otherwise what are we plotting? Only data?
tgrMC_vx = array('f')
tgrMC_evx = array('f')
#these vectors are needed for nuisances accounting
nuisances_vy_up = {}
nuisances_vy_do = {}
ROOT.gROOT.cd()
thsData = ROOT.THStack ("thsData_" + cutName + "_" + variableName, "thsData_" + cutName + "_" + variableName)
#print('really before thstack ... one')
thsSignal = ROOT.THStack ("thsSignal_" + cutName + "_" + variableName, "thsSignal_" + cutName + "_" + variableName)
#print('really before thstack ... two')
thsBackground = ROOT.THStack ("thsBackground_" + cutName + "_" + variableName,"thsBackground_" + cutName + "_" + variableName)
#print('really before thstack ... three')
thsSignal_grouped = ROOT.THStack ("thsSignal_grouped_" + cutName + "_" + variableName, "thsSignal_grouped_" + cutName + "_" + variableName)
#print('really before thstack ... four')
thsBackground_grouped = ROOT.THStack ("thsBackground_grouped_" + cutName + "_" + variableName,"thsBackground_grouped_" + cutName + "_" + variableName)
list_thsData [generalCounter] = thsData
list_thsSignal [generalCounter] = thsSignal
list_thsBackground [generalCounter] = thsBackground
list_thsSignal_grouped [generalCounter] = thsSignal_grouped
list_thsBackground_grouped [generalCounter] = thsBackground_grouped
# for special case of plotting normalized
thsSignal_grouped_normalized = ROOT.THStack ("thsSignal_grouped_normalized_" + cutName + "_" + variableName, "thsSignal_grouped_normalized_" + cutName + "_" + variableName)
thsBackground_grouped_normalized = ROOT.THStack ("thsBackground_grouped_normalized_" + cutName + "_" + variableName,"thsBackground_grouped_normalized_" + cutName + "_" + variableName)
list_thsSignal_grouped_normalized [generalCounter] = thsSignal_grouped_normalized
list_thsBackground_grouped_normalized [generalCounter] = thsBackground_grouped_normalized
generalCounter += 1
#print('... after thstack ...')
sigSupList = []
sigSupList_grouped = []
# list of additional histograms to be used in the ratio plot
sigForAdditionalRatioList = {}
sigForAdditionalDifferenceList = {}
# enhanced list of nuisances, including bin-by-bin
mynuisances = {}
nexpected = 0
for sampleName, plotdef in plot.items():
if 'samples' in variable and sampleName not in variable['samples']:
continue
shapeName = cutName+"/"+variableName+'/histo_' + sampleName
print(' -> shapeName = ', shapeName)
if type(fileIn) is dict:
histo = fileIn[sampleName].Get(shapeName)
else:
histo = fileIn.Get(shapeName)
if not histo: continue
print(' --> ', histo)
print('new_histo_' + sampleName + '_' + cutName + '_' + variableName)
histos[sampleName] = histo.Clone('new_histo_' + sampleName + '_' + cutName + '_' + variableName)
#print(" -> sampleName = ", sampleName, " --> ", histos[sampleName].GetTitle(), " --> ", histos[sampleName].GetName(), " --> ", histos[sampleName].GetNbinsX())
#for iBinAmassiro in range(1, histos[sampleName].GetNbinsX()+1):
#print(" i = ", iBinAmassiro, " [" , sampleName, " ==> ", histos[sampleName].GetBinContent(iBinAmassiro))
# allow arbitrary scaling in MC (and DATA??), if needed
# for example to "see" a signal
if 'scale' in plotdef.keys() :
histos[sampleName].Scale(plotdef['scale'])
#print(" >> scale ", sampleName, " to ", plotdef['scale'])
# apply cut dependent scale factors
# for example when plotting different phase spaces
if 'cuts' in plotdef.keys() and cutName in plotdef['cuts']:
histos[sampleName].Scale( float( plotdef['cuts'][cutName] ) )
# data style
if plotdef['isData'] == 1 :
if variable['divideByBinWidth'] == 1:
histos[sampleName].Scale(1,"width")
#print(' plot[', sampleName, '][color] = ' , plotdef['color'])
histos[sampleName].SetMarkerColor(plotdef['color'])
histos[sampleName].SetMarkerSize(1)
histos[sampleName].SetMarkerStyle(20)
histos[sampleName].SetLineColor(self._getColor(plotdef['color']))
# blind data
if 'isBlind' in plotdef.keys() and plotdef['isBlind'] == 1:
histos[sampleName].Reset()
# Per variable blinding
if 'blind' in variable:
if cutName in variable['blind']:
blind_range = variable['blind'][cutName]
if blind_range == "full":
for iBin in range(1, histos[sampleName].GetNbinsX()+1):
histos[sampleName].SetBinContent(iBin, 0)
histos[sampleName].SetBinError (iBin, 0)
histos[sampleName].Reset()
elif type(blind_range) in [list,tuple] and len(blind_range)==2:
b0 = histos[sampleName].FindBin(blind_range[0])
b1 = histos[sampleName].FindBin(blind_range[1])
for iBin in range(1, histos[sampleName].GetNbinsX()+1):
if iBin >= b0 and iBin <= b1:
histos[sampleName].SetBinContent(iBin, 0)
histos[sampleName].SetBinError (iBin, 0)
# Allow to also pass arbitrary set of bin indexes to blind (e.g. for unrolled 2D histos)
elif type(blind_range) in [list,tuple] and len(blind_range)>2:
for iBin in blind_range:
histos[sampleName].SetBinContent(iBin, 0)
histos[sampleName].SetBinError (iBin, 0)
thsData.Add(histos[sampleName])
# first time fill vectors X axis
if len(tgrData_vx) == 0 :
dataColor = histos[sampleName].GetMarkerColor()
for iBin in range(1, histos[sampleName].GetNbinsX()+1):
tgrData_vx.append( histos[sampleName].GetBinCenter (iBin))
tgrData_evx.append( histos[sampleName].GetBinWidth (iBin) / 2.)
tgrData_vy.append( histos[sampleName].GetBinContent (iBin))
#print(" plot[", sampleName, "].keys() = ", plotdef.keys())
if ('isSignal' not in plotdef.keys() or plotdef['isSignal'] != 3) and not ('isBlind' in plotdef.keys() and plotdef['isBlind'] == 1) :
if variable['divideByBinWidth'] == 1:
tgrData_evy_up.append( self.GetPoissError(histos[sampleName].GetBinContent (iBin) * histos[sampleName].GetBinWidth (iBin) , 0, 1) / histos[sampleName].GetBinWidth (iBin) )
tgrData_evy_do.append( self.GetPoissError(histos[sampleName].GetBinContent (iBin) * histos[sampleName].GetBinWidth (iBin) , 1, 0) / histos[sampleName].GetBinWidth (iBin) )
else:
tgrData_evy_up.append( self.GetPoissError(histos[sampleName].GetBinContent (iBin) , 0, 1) )
tgrData_evy_do.append( self.GetPoissError(histos[sampleName].GetBinContent (iBin) , 1, 0) )
else :
tgrData_evy_up.append( 0 )
tgrData_evy_do.append( 0 )
else :
for iBin in range(1, histos[sampleName].GetNbinsX()+1):
tgrData_vx[iBin-1] = ( histos[sampleName].GetBinCenter (iBin))
tgrData_evx.append( histos[sampleName].GetBinWidth (iBin) / 2.)
tgrData_vy[iBin-1] += histos[sampleName].GetBinContent (iBin)
if 'isSignal' not in plotdef.keys() or plotdef['isSignal'] == 3 :
if variable['divideByBinWidth'] == 1:
tgrData_evy_up[iBin-1] = SumQ ( tgrData_evy_up[iBin-1], self.GetPoissError(histos[sampleName].GetBinContent (iBin) * histos[sampleName].GetBinWidth (iBin) , 0, 1) / histos[sampleName].GetBinWidth (iBin))
tgrData_evy_do[iBin-1] = SumQ ( tgrData_evy_do[iBin-1], self.GetPoissError(histos[sampleName].GetBinContent (iBin) * histos[sampleName].GetBinWidth (iBin) , 1, 0) / histos[sampleName].GetBinWidth (iBin))
else:
tgrData_evy_up[iBin-1] = SumQ ( tgrData_evy_up[iBin-1], self.GetPoissError(histos[sampleName].GetBinContent (iBin) , 0, 1) )
tgrData_evy_do[iBin-1] = SumQ ( tgrData_evy_do[iBin-1], self.GetPoissError(histos[sampleName].GetBinContent (iBin) , 1, 0) )
# if plotdef['isData'] == 1:
else:
# MC style
# only background "filled" histogram
if plotdef['isSignal'] == 0:
histos[sampleName].SetFillColorAlpha(self._getColor(plotdef['color']), 0.60)
histos[sampleName].SetLineColor(self._getColor(plotdef['color']))
if 'fill' in plotdef:
histos[sampleName].SetFillStype(plotdef['fill'])
#else:
# histos[sampleName].SetFillStyle(3001)
else :
histos[sampleName].SetFillStyle(0)
histos[sampleName].SetLineWidth(2)
histos[sampleName].SetLineColor(self._getColor(plotdef['color']))
# scale to luminosity if MC
#histos[sampleName].Scale(self._lumi) ---> NO! They are already scaled to luminosity in mkShape!
if plotdef['isSignal'] == 1 :
if variable['divideByBinWidth'] == 1:
histos[sampleName].Scale(1,"width")
thsSignal.Add(histos[sampleName])
elif plotdef['isSignal'] == 2 or plotdef['isSignal'] == 3 :
#print("SigSup histo: ", histos[sampleName])
if variable['divideByBinWidth'] == 1:
histos[sampleName].Scale(1,"width")
sigSupList.append(histos[sampleName])
if plotdef['isSignal'] == 3 :
#print("sigForAdditionalRatio histo: ", histos[sampleName])
sigForAdditionalRatioList[sampleName] = histos[sampleName]
sigForAdditionalDifferenceList[sampleName] = histos[sampleName]
else :
nexpected += histos[sampleName].Integral(1,histos[sampleName].GetNbinsX()) # it was (-1, -1) in the past, correct now. Overflow and underflow bins not taken into account
if variable['divideByBinWidth'] == 1:
histos[sampleName].Scale(1,"width")
# this is really background only, meaning if "isSignal" is 1,2,3 then it's not included here
thsBackground.Add(histos[sampleName])
#print(" adding to background: ", sampleName)
# handle 'stat' nuisance to create the bin-by-bin list of nuisances
# "massage" the list of nuisances accordingly
for nuisanceName, nuisance in nuisances.items():
if 'cuts' in nuisance and cutName not in nuisance['cuts']:
continue
# run only if this nuisance will affect the phase space defined in "cut"
#print(" nuisanceName = ", nuisanceName)
if nuisanceName == 'stat' : # 'stat' has a separate treatment, it's the MC/data statistics
#print(" nuisance = ", nuisance)
if 'samples' in nuisance.keys():
if sampleName in nuisance['samples'].keys() :
#print(" stat nuisances for ", sampleName)
if nuisance['samples'][sampleName]['typeStat'] == 'uni' : # unified approach
print('In principle nothing to be done here ... just wait')
if nuisance['samples'][sampleName]['typeStat'] == 'bbb' : # bin-by-bin
# add N ad hoc nuisances, one for each bin
for iBin in range(1, histos[sampleName].GetNbinsX()+1):
if ('ibin_' + str(iBin) + '_stat') not in mynuisances.keys() : # if new, add the new nuisance
# Name of the histogram: histo_" + sampleName + "_ibin_" + str(iBin) + "_statUp"
# Then the nuisance is "ibin_" + str(iBin) + "_stat"
mynuisances['ibin_' + str(iBin) + '_stat'] = {
'samples' : { sampleName : '1.00', },
}
else : # otherwise just add the new sample in the list of samples to be considered
mynuisances['ibin_' + str(iBin) + '_stat']['samples'][sampleName] = '1.00'
else :
if nuisanceName not in mynuisances.keys() :
if 'type' in nuisance.keys() and (nuisance['type'] == 'rateParam' or nuisance['type'] == 'lnU') :
pass
#print("skip this nuisance since 100 percent uncertainty :: ", nuisanceName)
else :
mynuisances[nuisanceName] = nuisances[nuisanceName]
nuisanceHistos = ({}, {})
for nuisanceName, nuisance in mynuisances.items():
# is this nuisance to be considered for this background?
if 'samples' in nuisance:
if sampleName not in nuisance['samples']:
continue
elif 'all' not in nuisance or nuisance['all'] != 1:
continue
if 'cuts' in nuisance and cutName not in nuisance['cuts']:
continue
if 'name' in nuisance:
shapeNameVars = tuple(cutName+"/"+variableName+'/histo_' + sampleName+"_"+nuisance['name']+var for var in ['Up', 'Down'])
else:
shapeNameVars = tuple(cutName+"/"+variableName+'/histo_' + sampleName+"_"+nuisanceName+var for var in ['Up', 'Down'])
if 'type' in nuisance and nuisance['type'] == 'lnN':
if 'samples' in nuisance:
values = nuisance['samples'][sampleName]
# example:
# 'samples' : {
# 'WW' : '1.00',
# 'ggH': '1.23/0.97'
# },
else: # 'all'
values = nuisance['value']
if '/' in values:
variations = map(float, values.split('/'))
else:
variations = (float(values), 2. - float(values))
# don't use histos[sampleName], or the second "scale" will fail!!!
for ivar, shapeNameVar in enumerate(shapeNameVars):
histoVar = histo.Clone(shapeNameVar.replace('/', '__'))
histoVar.Scale(variations[ivar])
nuisanceHistos[ivar][nuisanceName] = histoVar
else:
for ivar, shapeNameVar in enumerate(shapeNameVars):
if type(fileIn) is dict:
histoVar = fileIn[sampleName].Get(shapeNameVar)
else:
histoVar = fileIn.Get(shapeNameVar)
if histoVar != None :
nuisanceHistos[ivar][nuisanceName] = histoVar
elif not self._SkipMissingNuisance :
print(" This is bad, the nuisance ", nuisanceName, " is missing! You need to add it, maybe some jobs crashed?")
nuisanceHistos[ivar][nuisanceName]= fileIn.Get(cutName+"/"+variableName+'/histo_' + sampleName)
#nuisanceHistos[ivar][nuisanceName] = histoVar
else :
# if you had self._SkipMissingNuisance set to true, put the variation the same as the nominal
histoVar = histo.Clone(shapeNameVar.replace('/', '__'))
nuisanceHistos[ivar][nuisanceName] = histoVar
#
# I fill the nuisances_vy_up and nuisances_vy_do with the sum of all "up" and all "down" variations
#
for ivar, nuisances_vy in enumerate([nuisances_vy_up, nuisances_vy_do]):
for nuisanceName, nuisance in mynuisances.items():
#print(nuisanceName)
try:
histoVar = nuisanceHistos[ivar][nuisanceName]
#test = rnp.hist2array(histoVar, copy=False)
test = hist2array(histoVar, copy=False)
except KeyError:
# now, even if not considered this nuisance, I need to add it,
# so that in case is "empty" it will add the nominal value
# for this sample that is not affected by the nuisance
#print(nuisanceName, sampleName)
histoVar = histos[sampleName]
except TypeError:
histoVar = histos[sampleName]
else:
if 'scale' in plotdef:
histoVar.Scale(plotdef['scale'])
# apply cut dependent scale factors
# for example when plotting different phase spaces
if 'cuts' in plotdef and cutName in plotdef['cuts']:
histoVar.Scale(float(plotdef['cuts'][cutName]))
if variable["divideByBinWidth"] == 1:
histoVar.Scale(1., "width")
try:
vy = nuisances_vy[nuisanceName]
#print(sampleName, nuisanceName, vy)
except KeyError:
#vy = nuisances_vy[nuisanceName] = np.zeros_like(rnp.hist2array(histo, copy=False))
vy = nuisances_vy[nuisanceName] = np.zeros_like(hist2array(histo, copy=False))
# get the background sum
if plotdef['isSignal'] == 0: # ---> add the signal too????? See ~ 20 lines below
#vy += rnp.hist2array(histoVar, copy=False)
vy += hist2array(histoVar, copy=False)
# print(sampleName, nuisanceName, rnp.hist2array(histoVar, copy=False))
# create the group of histograms to plot
# this has to be done after the scaling of the previous lines
# andl also after all the rest, so that we inherit the style of the histograms
for sampleNameGroup, sampleConfiguration in groupPlot.items():
if sampleName in sampleConfiguration['samples']:
if sampleNameGroup in histos_grouped.keys() :
histos_grouped[sampleNameGroup].Add(histos[sampleName])
else :
histos_grouped[sampleNameGroup] = histos[sampleName].Clone('new_histo_group_' + sampleNameGroup + '_' + cutName + '_' + variableName)
# end sample loop
# set the colors for the groups of samples
for sampleNameGroup, sampleConfiguration in groupPlot.items():
if sampleNameGroup in histos_grouped.keys() :
histos_grouped[sampleNameGroup].SetLineColor(self._getColor(sampleConfiguration['color']))
if sampleConfiguration['isSignal'] == 0:
histos_grouped[sampleNameGroup].SetFillColorAlpha(self._getColor(sampleConfiguration['color']), 0.60)
histos_grouped[sampleNameGroup].SetLineColor(self._getColor(sampleConfiguration['color']))
if 'fill' in sampleConfiguration:
histos_grouped[sampleNameGroup].SetFillStyle(sampleConfiguration['fill'])
#else:
# histos_grouped[sampleNameGroup].SetFillStyle(3001)
else :
histos_grouped[sampleNameGroup].SetFillStyle(0)
histos_grouped[sampleNameGroup].SetLineWidth(2)
# fill the reference distribution with the background only distribution
# save the central values of the bkg sum for use for the nuisance band
#
# How could this be ==0 ?
# At least one MC sample should be defined ...
# but still, let's leave the possibility
#
if thsBackground.GetNhists() != 0:
last = thsBackground.GetStack().Last()
#tgrMC_vy = rnp.hist2array(last, copy=True)
tgrMC_vy = hist2array(last, copy=True)
for iBin in range(1,thsBackground.GetStack().Last().GetNbinsX()+1):
tgrMC_vx .append(thsBackground.GetStack().Last().GetBinCenter(iBin))
tgrMC_evx.append(thsBackground.GetStack().Last().GetBinWidth(iBin) / 2.)
# nuisances_err2_up = rnp.array(last.GetSumw2())[1:-1]
# nuisances_err2_do = rnp.array(last.GetSumw2())[1:-1]
nuisances_err2_up = np_array(last.GetSumw2())[1:-1]
nuisances_err2_do = np_array(last.GetSumw2())[1:-1]
if (self._removeMCStat):
nuisances_err2_up.fill(0)
nuisances_err2_do.fill(0)
else:
tgrMC_vy = np.zeros((0,))
nuisances_err2_up = np.zeros((0,))
nuisances_err2_do = np.zeros((0,))
#
# and now let's add the signal on top of the background stack
# It is important to do this after setting (without signal) tgrMC_vy
#
for sampleName, plotdef in plot.items():
if 'samples' in variable and sampleName not in variable['samples']:
continue
# MC style
if plotdef['isData'] == 0 :
if plotdef['isSignal'] == 1 :
thsBackground.Add(histos[sampleName])
#
# you need to add the signal as well, since the signal was considered in the nuisances vector
# otherwise you would introduce an uncertainty as big as the signal itself!!!
#
# IF these lines below are commented it means that the uncertainty band is drawn only on bkg-only stacked historgram,
# meaning that you will get the dashed band on bkg stacked, and signal (typicalli empty red) drawn on top
# without the uncertainty band.
# Is this what you want? I hope so ...
#
#if thsSignal.GetNhists() != 0:
#for iBin in range(1,thsSignal.GetStack().Last().GetNbinsX()+1):
#tgrMC_vy[iBin] += (thsSignal.GetStack().Last().GetBinContent(iBin))
#print(" nominal: ", iBin, " ===> ", thsBackground.GetStack().Last().GetBinContent(iBin))
#print(" tgrMC_vy = ", tgrMC_vy)
#else :
# for iBin in range(1, histos[sampleName].GetNbinsX()+1):
# tgrBkg_vx[iBin-1] = ( histos[sampleName].GetBinCenter (iBin))
# tgrBkg_evx.append( histos[sampleName].GetBinWidth (iBin) / 2.)
# tgrBkg_vy[iBin-1] += histos[sampleName].GetBinContent (iBin)
# tgrBkg_evy_up[iBin-1] = SumQ ( tgrBkg_evy_up[iBin-1], self.GetPoissError(histos[sampleName].GetBinContent (iBin) , 0, 1) )
# tgrBkg_evy_do[iBin-1] = SumQ ( tgrBkg_evy_do[iBin-1], self.GetPoissError(histos[sampleName].GetBinContent (iBin) , 1, 0) )
#print(">>>> Sample name: ", sampleName)
for nuisanceName in mynuisances.keys():
# now we need to tell wthether the variation is actually up or down ans sum in quadrature those with the same sign
up = nuisances_vy_up[nuisanceName]
do = nuisances_vy_do[nuisanceName]
#
# NB: the test underneath is to be read "for each bin check if up[bin] > tgrMC_vy[bin] and store in an array of True or False
# In other words, up_is_up is an array([False, False, False])
#
up_is_up = (up > tgrMC_vy)
#
# this one above is an array of booleans
#
dup2 = np.square(up - tgrMC_vy)
ddo2 = np.square(do - tgrMC_vy)
nuisances_err2_up += np.where(up_is_up, dup2, ddo2)
nuisances_err2_do += np.where(up_is_up, ddo2, dup2)
nuisances_err_up = np.sqrt(nuisances_err2_up)
nuisances_err_do = np.sqrt(nuisances_err2_do)
#
# NB: reminder: only background uncertainties have been considered in the nuisances, no impacts on the signal (isSignal = 1,2,3)
#
tgrData = ROOT.TGraphAsymmErrors(thsBackground.GetStack().Last().GetNbinsX())
for iBin in range(0, len(tgrData_vx)) :
tgrData.SetPoint (iBin, tgrData_vx[iBin], tgrData_vy[iBin])
tgrData.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], tgrData_evy_do[iBin], tgrData_evy_up[iBin])
tgrData.SetMarkerColor(dataColor)
tgrData.SetLineColor(dataColor)
## Default: --postFit 0 --> No additional line is drawn
## ----------------------------------------------------
## --postFit 1 --> line is prefit
if self._postFit == 'p':
tgrDataOverPF = tgrData.Clone("tgrDataOverPF") # use this for ratio with Post-Fit MC
histoPF = fileIn.Get(cutName+"/"+variableName+'/histo_total_prefit')
## --postFit 2 --> line is (S+B) postfit
if self._postFit == 's':
tgrDataOverPF = tgrData.Clone("tgrDataOverPF") # use this for ratio with Post-Fit MC
histoPF = fileIn.Get(cutName+"/"+variableName+'/histo_total_postfit_s')
## --postFit 3 --> line is B-only postfit
if self._postFit == 'b':
tgrDataOverPF = tgrData.Clone("tgrDataOverPF") # use this for ratio with Post-Fit MC
histoPF = fileIn.Get(cutName+"/"+variableName+'/histo_total_postfit_b')
# at this stage "thsBackground" and then "last" includes ALSO the signal
last = thsBackground.GetStack().Last()
tgrDataOverMC = tgrData.Clone("tgrDataOverMC")
tgrDataMinusMC = tgrData.Clone("tgrDataMinusMC")
tgrMCSigMinusMCBkg = tgrData.Clone("tgrMCSigMinusMCBkg") # is like saying "signal"
for iBin in range(0, len(tgrData_vx)) :
tgrDataOverMC.SetPoint (iBin, tgrData_vx[iBin], self.Ratio(tgrData_vy[iBin] , last.GetBinContent(iBin+1)) )
tgrDataOverMC.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], self.Ratio(tgrData_evy_do[iBin], last.GetBinContent(iBin+1)) , self.Ratio(tgrData_evy_up[iBin], last.GetBinContent(iBin+1)) )
if self._postFit == 'p' or self._postFit == 's' or self._postFit == 'b':
tgrDataOverPF.SetPoint(iBin, tgrData_vx[iBin], self.Ratio(tgrData_vy[iBin] , histoPF.GetBinContent(iBin+1)) )
tgrDataOverPF.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], self.Ratio(tgrData_evy_do[iBin], histoPF.GetBinContent(iBin+1)) , self.Ratio(tgrData_evy_up[iBin], last.GetBinContent(iBin+1)) )
print("Pre-fit ratio: " + str(self.Ratio(tgrData_vy[iBin] , histoPF.GetBinContent(iBin+1))))
print("Post-fit ratio: " + str(self.Ratio(tgrData_vy[iBin] , last.GetBinContent(iBin+1))))
print(iBin)
#
# data - MC :
# MC could be background only
# or it can include the signal.
# Default is background+signal (check isSignal = 1,2,3 options).
# tgrMC_vy is "background only" !!
# You can activate the data - "background only" by
# using the flag "showDataMinusBkgOnly".
# NB: this will change also the case of "(data - expected) / expected"
#
#
if self._showRelativeRatio :
if self._showDataMinusBkgOnly :
tgrDataMinusMC.SetPoint (iBin, tgrData_vx[iBin], self.Ratio( self.Difference(tgrData_vy[iBin] , tgrMC_vy[iBin]), tgrMC_vy[iBin] ) )
tgrDataMinusMC.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], self.Ratio(tgrData_evy_do[iBin], tgrMC_vy[iBin]) , self.Ratio(tgrData_evy_up[iBin], tgrMC_vy[iBin]) )
else :
tgrDataMinusMC.SetPoint (iBin, tgrData_vx[iBin], self.Ratio( self.Difference(tgrData_vy[iBin] , last.GetBinContent(iBin+1)), last.GetBinContent(iBin+1)) )
tgrDataMinusMC.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], self.Ratio(tgrData_evy_do[iBin], last.GetBinContent(iBin+1)) , self.Ratio(tgrData_evy_up[iBin], last.GetBinContent(iBin+1)) )
else :
if self._showDataMinusBkgOnly :
tgrDataMinusMC.SetPoint (iBin, tgrData_vx[iBin], self.Difference(tgrData_vy[iBin] , tgrMC_vy[iBin] ) )
tgrDataMinusMC.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], tgrData_evy_do[iBin] , tgrData_evy_up[iBin] )
tgrMCSigMinusMCBkg.SetPoint (iBin, tgrData_vx[iBin], self.Difference(last.GetBinContent(iBin+1) , tgrMC_vy[iBin] ) )
tgrMCSigMinusMCBkg.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], 0 , 0 ) # error set to 0 for sake of simplicity
else :
tgrDataMinusMC.SetPoint (iBin, tgrData_vx[iBin], self.Difference(tgrData_vy[iBin] , last.GetBinContent(iBin+1)) )
tgrDataMinusMC.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], tgrData_evy_do[iBin] , tgrData_evy_up[iBin] )
#
# if there is an histogram called 'histo_total'
# it means that post-fit plots are provided,
# and we can neglect about the nuisances as "quadrature sum" here
# and use directly the error bars (so far symmetric
# see https://hypernews.cern.ch/HyperNews/CMS/get/higgs-combination/995.html )
# from the histogram itself
#
# NB: the post-fit plot includes the signal!
#
#
special_shapeName = cutName+"/"+variableName+'/histo_total'
if type(fileIn) is dict:
if 'total' in fileIn:
histo_total = fileIn['total'].Get(special_shapeName)
else:
histo_total = None
else:
histo_total = fileIn.Get(special_shapeName)
if variable['divideByBinWidth'] == 1 and histo_total != None:
histo_total.Scale(1,"width")
print('--> histo_total = ', histo_total)
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMC = ROOT.TGraphAsymmErrors()
for iBin in range(0, len(tgrMC_vx)) :
tgrMC.SetPoint (iBin, tgrMC_vx[iBin], tgrMC_vy[iBin])
if histo_total:
#
# if there is histo_total, we plot the uncertainty band on top of sig+bkg, since histo_total IS sig+bkg and it's uncertainty accordingly
# This fix, in case of histo_total overrules the comment
# few lines above "Is this what you want? I hope so ..."
#
tgrMC.SetPoint(iBin, tgrMC_vx[iBin], histo_total.GetBinContent(iBin+1))
tgrMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], histo_total.GetBinError(iBin+1), histo_total.GetBinError(iBin+1))
else :
tgrMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], nuisances_err_do[iBin], nuisances_err_up[iBin])
tgrMCOverMC = tgrMC.Clone("tgrMCOverMC")
tgrMCMinusMC = tgrMC.Clone("tgrMCMinusMC")
for iBin in range(0, len(tgrMC_vx)) :
tgrMCOverMC.SetPoint (iBin, tgrMC_vx[iBin], 1.) # the ratio MC / MC is by construction 1
tgrMCMinusMC.SetPoint (iBin, tgrMC_vx[iBin], 0.) # the difference MC - MC is by construction 0
if histo_total:
# histo_total include also the signal
tgrMCOverMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], self.Ratio(histo_total.GetBinError(iBin+1), last.GetBinContent(iBin+1)), self.Ratio(histo_total.GetBinError(iBin+1), last.GetBinContent(iBin+1) ))
if self._showDataMinusBkgOnly :
if self._showRelativeRatio :
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], self.Ratio(histo_total.GetBinError(iBin+1), tgrMC_vy[iBin]), self.Ratio(histo_total.GetBinError(iBin+1), tgrMC_vy[iBin]))
else :
# ok, this should have been the error on background only, without signal, properly propagated ... in first approximation it's the uncertainty on sig+bkg = histo_total
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], histo_total.GetBinError(iBin+1), histo_total.GetBinError(iBin+1))
else :
if self._showRelativeRatio :
# not 100% sure if in the relative ratio I should put an uncertainty bar on the expected ... but I assume yes, as in the ratio plot, since expected-rate uncertainty is not propagated to "data"
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], self.Ratio(histo_total.GetBinError(iBin+1), last.GetBinContent(iBin+1)), self.Ratio(histo_total.GetBinError(iBin+1), last.GetBinContent(iBin+1)))
else :
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], histo_total.GetBinError(iBin+1), histo_total.GetBinError(iBin+1))
else :
# nuisances_err_do and nuisances_err_up do NOT include the signal
# thus in first approximation we say "uncertainty_(sig+bgk) / (sig+bkg) = uncertainty_(bkg) / bkg
# this is why everywhere here we have "tgrMC_vy[iBin]" instead of "last.GetBinContent(iBin+1)"
tgrMCOverMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], self.Ratio(nuisances_err_do[iBin], tgrMC_vy[iBin]), self.Ratio(nuisances_err_up[iBin], tgrMC_vy[iBin]))
if self._showDataMinusBkgOnly :
if self._showRelativeRatio :
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], self.Ratio(nuisances_err_do[iBin], tgrMC_vy[iBin]), self.Ratio(nuisances_err_up[iBin], tgrMC_vy[iBin]))
else :
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], nuisances_err_do[iBin], nuisances_err_up[iBin])
else :
if self._showRelativeRatio :
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], self.Ratio(nuisances_err_do[iBin], tgrMC_vy[iBin]), self.Ratio(nuisances_err_up[iBin], tgrMC_vy[iBin]))
else :
tgrMCMinusMC.SetPointError(iBin, tgrMC_evx[iBin], tgrMC_evx[iBin], nuisances_err_do[iBin], nuisances_err_up[iBin])
tgrRatioList = {}
for samplesToRatioName, samplesToRatio in sigForAdditionalRatioList.items() :
tgrDataOverMCTemp = tgrData.Clone("tgrDataOverMC"+samplesToRatioName)
for iBin in range(0, len(tgrData_vx)) :
tgrDataOverMCTemp.SetPoint (iBin, tgrData_vx[iBin], self.Ratio(tgrData_vy[iBin] , samplesToRatio.GetBinContent(iBin+1)) )
tgrDataOverMCTemp.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], self.Ratio(tgrData_evy_do[iBin], samplesToRatio.GetBinContent(iBin+1)) , self.Ratio(tgrData_evy_up[iBin], samplesToRatio.GetBinContent(iBin+1)) )
if variableName == 'events' :
print(' >> ratio[', cutName, '][', samplesToRatioName, '] = ', self.Ratio(tgrData_vy[0] , samplesToRatio.GetBinContent(0+1)) )
tgrDataOverMCTemp.SetLineColor(samplesToRatio.GetLineColor())
tgrDataOverMCTemp.SetMarkerColor(samplesToRatio.GetLineColor())
tgrDataOverMCTemp.SetMarkerSize(0.3)
tgrRatioList[samplesToRatioName] = tgrDataOverMCTemp
tgrDifferenceList = {}
for samplesToDifferenceName, samplesToDifference in sigForAdditionalDifferenceList.items() :
tgrDataMinusMCTemp = tgrData.Clone("tgrDataMinusMC"+samplesToDifferenceName)
for iBin in range(0, len(tgrData_vx)) :
tgrDataMinusMCTemp.SetPoint (iBin, tgrData_vx[iBin], self.Difference(tgrData_vy[iBin] , samplesToDifference.GetBinContent(iBin+1)) )
tgrDataMinusMCTemp.SetPointError(iBin, tgrData_evx[iBin], tgrData_evx[iBin], tgrData_evy_do[iBin] , tgrData_evy_up[iBin] )
if variableName == 'events' :
print(' >> difference[', cutName, '][', samplesToDifferenceName, '] = ', self.Difference(tgrData_vy[0] , samplesToDifference.GetBinContent(0+1)) )
tgrDataMinusMCTemp.SetLineColor(samplesToDifference.GetLineColor())
tgrDataMinusMCTemp.SetMarkerColor(samplesToDifference.GetLineColor())
tgrDataMinusMCTemp.SetMarkerSize(0.3)
tgrDifferenceList[samplesToDifferenceName] = tgrDataMinusMCTemp
groupFlag = False
#---- prepare the grouped histograms
for sampleNameGroup, sampleConfiguration in groupPlot.items():
if 'samples' in variable and len(set(sampleConfiguration['samples']) & set(variable['samples'])) == 0:
continue
if sampleConfiguration['isSignal'] == 1 :
print("############################################################## isSignal 1", sampleNameGroup)
#
# if, for some reason, you want to scale only the overlaid signal
# for example to show the shape of the signal, without affecting the actual stacked (true) distribution
#
if 'scaleMultiplicativeOverlaid' in sampleConfiguration.keys() :
# may this clone not mess up too much with "gDirectory", see TH1::Copy
temp_overlaid = histos_grouped[sampleNameGroup].Clone()
temp_overlaid.Scale(sampleConfiguration['scaleMultiplicativeOverlaid'])
thsSignal_grouped.Add(temp_overlaid)
else :
thsSignal_grouped.Add(histos_grouped[sampleNameGroup])
elif sampleConfiguration['isSignal'] == 2 :
print("############################################################## isSignal 2", sampleNameGroup)
groupFlag = True
sigSupList_grouped.append(histos_grouped[sampleNameGroup])
# the signal is added on top of the background
# the signal has to be the last one in the dictionary!
# make it sure in plot.py
if groupFlag == False:
thsBackground_grouped.Add(histos_grouped[sampleNameGroup])
#---- now plot
if thsBackground.GetNhists() != 0:
print(" MC = ", thsBackground.GetStack().Last().Integral())
for ihisto in range(thsBackground.GetNhists()) :
print(" - ",ihisto, " - ", ((thsBackground.GetHists().At(ihisto))).GetName(), " = ", ((thsBackground.GetHists().At(ihisto))).Integral() )
if thsData.GetNhists() != 0:
print(" DATA = ", thsData.GetStack().Last().Integral())
# - get axis range
minXused = 0.
maxXused = 1.
minYused = 1.
maxYused = 1.
for sampleName, plotdef in plot.items():
if 'samples' in variable and sampleName not in variable['samples']:
continue
if plotdef['isData'] == 1 :
histos[sampleName].Draw("p")
minXused = histos[sampleName].GetXaxis().GetBinLowEdge(1)
maxXused = histos[sampleName].GetXaxis().GetBinUpEdge(histos[sampleName].GetNbinsX())
maxY = self.GetMaximumIncludingErrors(histos[sampleName])
histos[sampleName].SetMaximum(self._scaleToPlot * maxY)
maxYused = self._scaleToPlot * maxY
minYused = self.GetMinimum(histos[sampleName])
if thsBackground.GetNhists() != 0:
thsBackground.Draw("hist")
maxY = thsBackground.GetMaximum ()
minXused = thsBackground.GetXaxis().GetBinLowEdge(1)
maxXused = thsBackground.GetXaxis().GetBinUpEdge(thsBackground.GetHistogram().GetNbinsX())
if (self._scaleToPlot * maxY) > maxYused :
maxYused = self._scaleToPlot * maxY
minY = thsBackground.GetMinimum ()
if (minY < minYused) :
minYused = minY
if thsSignal.GetNhists() != 0:
thsSignal.Draw("hist")
maxY = thsSignal.GetMaximum ()
minXused = thsSignal.GetXaxis().GetBinLowEdge(1)
maxXused = thsSignal.GetXaxis().GetBinUpEdge(thsSignal.GetHistogram().GetNbinsX())
if (self._scaleToPlot * maxY) > maxYused :
maxYused = self._scaleToPlot * maxY
minY = thsSignal.GetMinimum ()
if (minY < minYused) :
minYused = minY
#print(" X axis = ", minXused, " - ", maxXused)
frame = ROOT.TH1F
frame = tcanvas.DrawFrame(minXused, 0.0, maxXused, 1.0)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxis = frame.GetXaxis()
xAxis.SetNdivisions(6,5,0)
# setup axis names
# New proposal (following https://twiki.cern.ch/twiki/bin/viewauth/CMS/Internal/PubGuidelines)
# Set xaxis label if needed. Format should be variable name, followed by (units) or [units] if needed
if 'xaxis' in variable.keys() :
frame.GetXaxis().SetTitle(variable['xaxis'])
# If yaxis is set, we want to override normal conventions
if 'yaxis' in variable.keys():
frame.GetYaxis().SetTitle(variable['yaxis'])
else:
# Grab unit from x axis title -- capture part in () or [] brackets
xaxistitle = frame.GetXaxis().GetTitle()
unitpattern = '(?:\[|\()(\w+)(?:\]|\))'
unitsearch = re.search(unitpattern,xaxistitle)
unit = 'unit' if unitsearch is None else unitsearch.group(1)
# If dividing by bin width, yaxis should be "<Events / [unit]>"
if variable["divideByBinWidth"] == 1:
frame.GetYaxis().SetTitle("< Events / %s >"%unit)
#frame.GetYaxis().SetTitle("dN/d"+variable['xaxis'])
else:
# If using fixed bin width, yaxis should be "Events / bin size [unit]"
if len(variable['range']) == 3:
binsize = float(variable['range'][2] - variable['range'][1])/float(variable['range'][0])
frame.GetYaxis().SetTitle("Events / %g %s"%(binsize,unit))
# Otherwise, yaxis should be "Events / bin"
else:
frame.GetYaxis().SetTitle("Events / bin")
#
# - now draw
# - first the MC
# before the background+signal, and then only the signal alone
#
if len(groupPlot.keys()) == 0:
if thsBackground.GetNhists() != 0:
thsBackground.Draw("hist same")
if thsSignal.GetNhists() != 0:
#for ihisto in range(thsSignal.GetNhists()) :
#((thsSignal.GetHists().At(ihisto))).SetFillStyle(0)
#((thsSignal.GetHists().At(ihisto))).Draw("hist same")
thsSignal.Draw("hist same noclear")
else :
if thsBackground_grouped.GetNhists() != 0:
thsBackground_grouped.Draw("hist same")
if thsSignal_grouped.GetNhists() != 0:
thsSignal_grouped.Draw("hist same noclear")
if len(sigSupList_grouped) != 0:
for histo in sigSupList_grouped:
histo.Draw("hist same")
# if there is a systematic band draw it
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMC.SetLineColor(12)
tgrMC.SetFillColor(12)
tgrMC.SetLineWidth(2)
tgrMC.SetFillStyle(3004)
tgrMCOverMC.SetLineColor(12)
tgrMCOverMC.SetFillColor(12)
tgrMCOverMC.SetLineWidth(2)
tgrMCOverMC.SetFillStyle(3004)
tgrMC.Draw("2")
# - then the superimposed MC
if len(sigSupList) != 0 and groupFlag==False:
for hist in sigSupList:
hist.Draw("hist same")
# - then the DATA
if tgrData.GetN() != 0:
tgrData.Draw("P0")
else : # never happening if at least one data histogram is provided
for sampleName, plotdef in plot.items():
if 'samples' in variable and sampleName not in variable['samples']:
continue
if plotdef['isData'] == 1 :
histos[sampleName].Draw("p same")
#---- the Legend
tlegend = ROOT.TLegend(0.20, 0.65, 0.80, 0.88)
tlegend.SetFillColor(0)
tlegend.SetTextFont(42)
tlegend.SetTextSize(0.035)
tlegend.SetLineColor(0)
tlegend.SetShadowColor(0)
reversedSampleNames = list(self._samples)
reversedSampleNames.reverse()
if len(groupPlot.keys()) == 0:
for sampleName in reversedSampleNames:
try:
plotdef = plot[sampleName]
except KeyError:
continue
if plotdef['isData'] == 0 :
if 'nameHR' in plotdef.keys() :
if plotdef['nameHR'] != '' :
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos[sampleName], plotdef['nameHR'], "F")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX(),"width")
else:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX())
tlegend.AddEntry(histos[sampleName], plotdef['nameHR'] + " [" + str(round(nevents,1)) + "]", "F")
else :
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos[sampleName], sampleName, "F")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX(),"width")
else:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX())
tlegend.AddEntry(histos[sampleName], sampleName + " [" + str(round(nevents,1)) + "]", "F")
for sampleName in reversedSampleNames:
try:
plotdef = plot[sampleName]
except KeyError:
continue
if plotdef['isData'] == 1 :
if 'nameHR' in plotdef.keys() :
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos[sampleName], plotdef['nameHR'], "EPL")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX(),"width")
else:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX())
tlegend.AddEntry(histos[sampleName], plotdef['nameHR'] + " [" + str(round(nevents,1)) + "]", "EPL")
else :
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos[sampleName], sampleName, "EPL")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX(),"width")
else:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX())
tlegend.AddEntry(histos[sampleName], sampleName + " [" + str(round(nevents,1)) + "]", "EPL")
else :
for sampleNameGroup, sampleConfiguration in groupPlot.items():
if 'samples' in variable and len(set(sampleConfiguration['samples']) & set(variable['samples'])) == 0:
continue
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos_grouped[sampleNameGroup], sampleConfiguration['nameHR'], "F")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos_grouped[sampleNameGroup].Integral(1,histos_grouped[sampleNameGroup].GetNbinsX(),"width")
else:
nevents = histos_grouped[sampleNameGroup].Integral(1,histos_grouped[sampleNameGroup].GetNbinsX())
tlegend.AddEntry(histos_grouped[sampleNameGroup], sampleConfiguration['nameHR'] + " [" + str(round(nevents,1)) + "]" , "F")
for sampleName in reversedSampleNames:
if 'samples' in variable and sampleName not in variable['samples']:
continue
try:
plotdef = plot[sampleName]
except KeyError:
continue
if plotdef['isData'] == 1 :
if 'nameHR' in plotdef.keys() :
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos[sampleName], plotdef['nameHR'], "EPL")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX(),"width")
else:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX())
print(" nevents [", sampleName, "] = ", nevents)
tlegend.AddEntry(histos[sampleName], plotdef['nameHR'] + " [" + str(round(nevents,1)) + "]", "EPL")
else :
if self._showIntegralLegend == 0 :
tlegend.AddEntry(histos[sampleName], sampleName , "EPL")
else :
if variable["divideByBinWidth"] == 1:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX(),"width")
else:
nevents = histos[sampleName].Integral(1,histos[sampleName].GetNbinsX())
print(" nevents [", sampleName, "] = ", nevents)
tlegend.AddEntry(histos[sampleName], sampleName + " [" + str(round(nevents,1)) + "]", "EPL")
# add "All MC" in the legend
if not self._removeAllMC :
# if there is "histo_total" there is no need of explicit nuisances
if len(mynuisances.keys()) != 0 or histo_total!= None:
if self._showIntegralLegend == 0 :
tlegend.AddEntry(tgrMC, "Syst.", "F")
else :
print(" nexpected = ", nexpected)
tlegend.AddEntry(tgrMC, "Syst. [" + str(round(nexpected,1)) + "]", "F")
tlegend.SetNColumns(2)
tlegend.Draw()
#change the CMS_lumi variables (see CMS_lumi.py)
import mkShapesRDF.shapeAnalysis.CMS_lumi as CMS_lumi
CMS_lumi.lumi_7TeV = "4.8 fb^{-1}"
CMS_lumi.lumi_8TeV = "18.3 fb^{-1}"
CMS_lumi.lumi_13TeV = "100 fb^{-1}"
CMS_lumi.writeExtraText = 1
CMS_lumi.extraText = "Preliminary"
if not self._preliminary :
CMS_lumi.extraText = ""
CMS_lumi.relPosX = 0.14
CMS_lumi.lumi_sqrtS = "13 TeV" # used with iPeriod = 0, e.g. for simulation-only plots (default is an empty string)
if 'sqrt' in legend.keys() :
CMS_lumi.lumi_sqrtS = legend['sqrt']
if 'lumi' in legend.keys() :
CMS_lumi.lumi_13TeV = legend['lumi']
else:
CMS_lumi.lumi_13TeV = 'L = %.1f fb^{-1}' % self._lumi
# Simple example of macro: plot with CMS name and lumi text
# (this script does not pretend to work in all configurations)
# iPeriod = 1*(0/1 7 TeV) + 2*(0/1 8 TeV) + 4*(0/1 13 TeV)
# For instance:
# iPeriod = 3 means: 7 TeV + 8 TeV
# iPeriod = 7 means: 7 TeV + 8 TeV + 13 TeV
# iPeriod = 0 means: free form (uses lumi_sqrtS)
iPeriod = 4
iPos = 0
CMS_lumi.CMS_lumi(tcanvas, iPeriod, iPos)
print("- draw tlegend")
#---- the Legend (end)
tlegend.Draw()
frame.GetYaxis().SetRangeUser( 0, maxYused )
# draw back all the axes
#frame.Draw("AXIS")
tcanvas.RedrawAxis()
if 'c' in self._plotsToWrite:
if self._plotLinear:
self._saveCanvas(tcanvas, self._outputDirPlots + "/" + canvasNameTemplate + self._FigNamePF)
if self._plotLog:
# log Y axis
frame.GetYaxis().SetRangeUser( max(self._minLogC, minYused), self._maxLogC * maxYused ) # Jonatan
#frame.GetYaxis().SetRangeUser( min(self._minLogC, minYused), self._maxLogC * maxYused ) # Jonatan
tcanvas.SetLogy(True)
# if plotLinear is true, we have already saved root and C (if in the list of formats)
self._saveCanvas(tcanvas, self._outputDirPlots + "/log_" + canvasNameTemplate + self._FigNamePF, imageOnly=self._plotLinear)
tcanvas.SetLogy(False)
if 'root' in self._fileFormats:
text_file_html.write(canvasNameTemplate + self._FigNamePF + ".root;\n")
# ~~~~~~~~~~~~~~~~~~~~
# plot with ratio plot
print("- draw with ratio")
canvasRatioNameTemplate = 'cratio_' + cutName + "_" + variableName
tcanvasRatio.cd()
canvasPad1Name = 'pad1_' + cutName + "_" + variableName
pad1 = ROOT.TPad(canvasPad1Name,canvasPad1Name, 0, 1-0.72, 1, 1)
pad1.SetTopMargin(0.098)
pad1.SetBottomMargin(0.020)
pad1.Draw()
#pad1.cd().SetGrid()
pad1.cd()
#print(" pad1 = ", pad1)
canvasFrameDistroName = 'frame_distro_' + cutName + "_" + variableName
frameDistro = pad1.DrawFrame(minXused, 0.0, maxXused, 1.0, canvasFrameDistroName)
#print(" pad1 = ", pad1)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = frameDistro.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
xAxisDistro.SetLabelSize(0)
# setup axis names
# New proposal (following https://twiki.cern.ch/twiki/bin/viewauth/CMS/Internal/PubGuidelines)
# Set xaxis label if needed. Format should be variable name, followed by (units) or [units] if needed
if 'xaxis' in variable.keys() :
frameDistro.GetXaxis().SetTitle(variable['xaxis'])
# If yaxis is set, we want to override normal conventions
if 'yaxis' in variable.keys():
frameDistro.GetYaxis().SetTitle(variable['yaxis'])
else:
# Grab unit from x axis title -- capture part in () or [] brackets
xaxistitle = frameDistro.GetXaxis().GetTitle()
unitpattern = '(?:\[|\()(\w+)(?:\]|\))'
unitsearch = re.search(unitpattern,xaxistitle)
unit = 'unit' if unitsearch is None else unitsearch.group(1)
# If dividing by bin width, yaxis should be "<Events / [unit]>"
if variable["divideByBinWidth"] == 1:
frameDistro.GetYaxis().SetTitle("< Events / %s >"%unit)
#frameDistro.GetYaxis().SetTitle("dN/d"+variable['xaxis'])
else:
# If using fixed bin width, yaxis should be "Events / bin size [unit]"
if len(variable['range']) == 3:
binsize = float(variable['range'][2] - variable['range'][1])/float(variable['range'][0])
frameDistro.GetYaxis().SetTitle("Events / %g %s"%(binsize,unit))
# Otherwise, yaxis should be "Events / bin"
else:
frameDistro.GetYaxis().SetTitle("Events / bin")
#frameDistro.GetYaxis().SetRangeUser( 0, maxYused )
frameDistro.GetYaxis().SetRangeUser( min(0.001, minYused), maxYused )
if len(groupPlot.keys()) == 0:
if thsBackground.GetNhists() != 0:
thsBackground.Draw("hist same")
if thsSignal.GetNhists() != 0:
#for ihisto in range(thsSignal.GetNhists()) :
#((thsSignal.GetHists().At(ihisto))).SetFillStyle(0)
#((thsSignal.GetHists().At(ihisto))).Draw("hist same")
thsSignal.Draw("hist same noclear")
else :
if thsBackground_grouped.GetNhists() != 0:
thsBackground_grouped.Draw("hist same")
if thsSignal_grouped.GetNhists() != 0:
thsSignal_grouped.Draw("hist same noclear")
if len(sigSupList_grouped) != 0:
for histo in sigSupList_grouped:
histo.Draw("hist same")
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMC.Draw("2")
# - then the superimposed MC
if len(sigSupList) != 0 and groupFlag==False:
for hist in sigSupList:
hist.Draw("hist same")
# - then the DATA
if tgrData.GetN() != 0:
tgrData.Draw("P0")
tlegend.Draw()
#if 'lumi' in legend.keys() and 'sqrt' not in legend.keys():
#flag_lumi = ROOT.TLatex (minXused + (maxXused-minXused)*3./4., 0 + (maxYused-0)*3.9/4., legend['lumi'])
#flag_lumi.Draw()
#if 'sqrt' in legend.keys() and 'lumi' not in legend.keys():
#flag_sqrt = ROOT.TLatex (minXused + (maxXused-minXused)*3./4., 0 + (maxYused-0)*3.9/4., legend['sqrt'])
#flag_sqrt.Draw()
#if 'sqrt' in legend.keys() and 'lumi' in legend.keys():
#flag_lumi_sqrt = ROOT.TLatex (minXused + (maxXused-minXused)*2.5/4., 0 + (maxYused-0)*3.9/4., "#splitline{CMS preliminary}{#splitline{" + legend['lumi'] + "}{" + legend['sqrt'] + "} }")
#flag_lumi_sqrt.Draw()
if self._extraLegend is not None:
legExtra = ROOT.TLatex()
legExtra.SetTextSize(0.045)
legExtra.SetTextAngle(0)
legExtra.SetTextAlign(22)
legExtra.SetTextFont(62)
legExtra.DrawLatexNDC(0.85,0.8,self._extraLegend)
CMS_lumi.CMS_lumi(tcanvasRatio, iPeriod, iPos)
# draw back all the axes
#frameDistro.Draw("AXIS")
pad1.RedrawAxis()
tcanvasRatio.cd()
canvasPad2Name = 'pad2_' + cutName + "_" + variableName
pad2 = ROOT.TPad(canvasPad2Name,canvasPad2Name,0,0,1,1-0.72)
pad2.SetTopMargin(0.06)
pad2.SetBottomMargin(0.392)
pad2.Draw()
#pad2.cd().SetGrid()
pad2.cd()
#print(" pad1 = ", pad1)
#print(" pad2 = ", pad2, " minXused = ", minXused, " maxXused = ", maxXused)
canvasFrameRatioName = 'frame_ratio_' + cutName + "_" + variableName
#print(" canvasFrameRatioName = ", canvasFrameRatioName)
frameRatio = pad2.DrawFrame(minXused, 0.0, maxXused, 2.0, canvasFrameRatioName)
#print(" pad2 = ", pad2)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = frameRatio.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
if 'xaxis' in variable.keys() :
frameRatio.GetXaxis().SetTitle(variable['xaxis'])
else :
frameRatio.GetXaxis().SetTitle(variableName)
frameRatio.GetYaxis().SetTitle("Data/Expected")
#frameRatio.GetYaxis().SetTitle("Data/MC")
#frameRatio.GetYaxis().SetRangeUser( 0.0, 2.0 )
#frameRatio.GetYaxis().SetRangeUser( 0.8, 1.2 )
# changed to 0.5 1.5!
frameRatio.GetYaxis().SetRangeUser( 0.5, 1.5 )
# frameRatio.GetYaxis().SetRangeUser( 0.5, 1.5 )
self.Pad2TAxis(frameRatio)
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMCOverMC.Draw("2")
tgrDataOverMC.Draw("P0")
if self._postFit == 'p' or self._postFit == 's' or self._postFit == 'b':
#---- Ratio Legend
tlegendRatio = ROOT.TLegend(0.20, 0.40, 0.60, 0.55)
tlegendRatio.SetFillColor(0)
tlegendRatio.SetTextFont(42)
##tlegendRatio.SetTextSize(0.035)
tlegendRatio.SetLineColor(0)
tlegendRatio.SetShadowColor(0)
if self._postFit == 'p':
tlegendRatio.AddEntry(tgrDataOverMC, "post-fit", "PL")
tlegendRatio.AddEntry(tgrDataOverPF, "pre-fit", "PL")
if self._postFit == 's' or self._postFit == 'b':
tlegendRatio.AddEntry(tgrDataOverMC, "pre-fit", "PL")
tlegendRatio.AddEntry(tgrDataOverPF, "post-fit", "PL")
for sampleName, sample in self._samples.items():
##if sampleName.find('total') == 1:
## or sampleName == 'total_background_prefit' or sampleName == 'total_background_postfit_s' or sampleName == 'total_background_postfit_b':
if 'total' in sampleName:
tgrDataOverPF.SetMarkerColor(plot[sampleName]['color'])
tgrDataOverPF.SetLineColor(plot[sampleName]['color'])
# tgrDataOverPF.SetMarkerColor(2)
# tgrDataOverPF.SetLineColor(2)
tgrDataOverPF.Draw("PE,same")
tlegendRatio.Draw("same")
for samplesToRatioGrName, samplesGrToRatio in tgrRatioList.items() :
samplesGrToRatio.Draw("P")
oneLine2 = ROOT.TLine(frameRatio.GetXaxis().GetXmin(), 1, frameRatio.GetXaxis().GetXmax(), 1);
oneLine2.SetLineStyle(3)
oneLine2.SetLineWidth(3)
oneLine2.Draw("same")
# draw back all the axes
#frameRatio.Draw("AXIS")
pad2.RedrawAxis()
pad2.SetGrid()
if 'cratio' in self._plotsToWrite:
if self._plotLinear:
self._saveCanvas(tcanvasRatio, self._outputDirPlots + "/" + canvasRatioNameTemplate + self._FigNamePF)
if self._plotLog:
# log Y axis
#frameDistro.GetYaxis().SetRangeUser( max(self._minLogCratio, maxYused/1000), self._maxLogCratio * maxYused )
#frameDistro.GetYaxis().SetRangeUser( min(self._minLogCratio, maxYused/1000), self._maxLogCratio * maxYused )
if maxYused>0:
frameDistro.GetYaxis().SetRangeUser( min(self._minLogCratio, maxYused/1000), math.pow(10, math.log10(maxYused) + (math.log10(maxYused)-math.log10(self._minLogCratio)) * (self._maxLinearScale-1)) )
else:
frameDistro.GetYaxis().SetRangeUser( min(self._minLogCratio, maxYused/1000), self._maxLogCratio * maxYused )
pad1.SetLogy(True)
self._saveCanvas(tcanvasRatio, self._outputDirPlots + "/log_" + canvasRatioNameTemplate + self._FigNamePF, imageOnly=self._plotLinear)
pad1.SetLogy(False)
if 'root' in self._fileFormats:
text_file_html.write(canvasRatioNameTemplate + ".root;\n")
# ~~~~~~~~~~~~~~~~~~~~
# plot with difference plot
print("- draw with difference")
if self._showRelativeRatio :
canvasDifferenceNameTemplate = 'cdifference_relative_' + cutName + "_" + variableName
else :
canvasDifferenceNameTemplate = 'cdifference_' + cutName + "_" + variableName
tcanvasDifference.cd()
canvasPad1differenceName = 'pad1difference_' + cutName + "_" + variableName
pad1difference = ROOT.TPad(canvasPad1differenceName,canvasPad1differenceName, 0, 1-0.72, 1, 1)
pad1difference.SetTopMargin(0.098)
pad1difference.SetBottomMargin(0.000)
pad1difference.Draw()
#pad1difference.cd().SetGrid()
pad1difference.cd()
#print(" pad1difference = ", pad1difference)
canvasFrameDistroName = 'frame_distro_' + cutName + "_" + variableName
frameDistro = pad1difference.DrawFrame(minXused, 0.0, maxXused, 1.0, canvasFrameDistroName)
#print(" pad1difference = ", pad1difference)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = frameDistro.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
# setup axis names
# New proposal (following https://twiki.cern.ch/twiki/bin/viewauth/CMS/Internal/PubGuidelines)
# Set xaxis label if needed. Format should be variable name, followed by (units) or [units] if needed
if 'xaxis' in variable.keys() :
frameDistro.GetXaxis().SetTitle(variable['xaxis'])
# If yaxis is set, we want to override normal conventions
if 'yaxis' in variable.keys():
frameDistro.GetYaxis().SetTitle(variable['yaxis'])
else:
# Grab unit from x axis title -- capture part in () or [] brackets
xaxistitle = frameDistro.GetXaxis().GetTitle()
unitpattern = '(?:\[|\()(\w+)(?:\]|\))'
unitsearch = re.search(unitpattern,xaxistitle)
unit = 'unit' if unitsearch is None else unitsearch.group(1)
# If dividing by bin width, yaxis should be "<Events / [unit]>"
if variable["divideByBinWidth"] == 1:
frameDistro.GetYaxis().SetTitle("< Events / %s >"%unit)
#frameDistro.GetYaxis().SetTitle("dN/d"+variable['xaxis'])
else:
# If using fixed bin width, yaxis should be "Events / bin size [unit]"
if len(variable['range']) == 3:
binsize = float(variable['range'][2] - variable['range'][1])/float(variable['range'][0])
frameDistro.GetYaxis().SetTitle("Events / %g %s"%(binsize,unit))
# Otherwise, yaxis should be "Events / bin"
else:
frameDistro.GetYaxis().SetTitle("Events / bin")
#frameDistro.GetYaxis().SetRangeUser( 0, maxYused )
frameDistro.GetYaxis().SetRangeUser( min(0.001, minYused), maxYused )
if len(groupPlot.keys()) == 0:
if thsBackground.GetNhists() != 0:
thsBackground.Draw("hist same")
if thsSignal.GetNhists() != 0:
#for ihisto in range(thsSignal.GetNhists()) :
#((thsSignal.GetHists().At(ihisto))).SetFillStyle(0)
#((thsSignal.GetHists().At(ihisto))).Draw("hist same")
thsSignal.Draw("hist same noclear")
else :
if thsBackground_grouped.GetNhists() != 0:
thsBackground_grouped.Draw("hist same")
if thsSignal_grouped.GetNhists() != 0:
thsSignal_grouped.Draw("hist same noclear")
if len(sigSupList_grouped) != 0:
for histo in sigSupList_grouped:
histo.Draw("hist same")
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMC.Draw("2")
# - then the superimposed MC
if len(sigSupList) != 0 and groupFlag==False:
for hist in sigSupList:
hist.Draw("hist same")
# - then the DATA
if tgrData.GetN() != 0:
tgrData.Draw("P0")
tlegend.Draw()
#if 'lumi' in legend.keys() and 'sqrt' not in legend.keys():
#flag_lumi = ROOT.TLatex (minXused + (maxXused-minXused)*3./4., 0 + (maxYused-0)*3.9/4., legend['lumi'])
#flag_lumi.Draw()
#if 'sqrt' in legend.keys() and 'lumi' not in legend.keys():
#flag_sqrt = ROOT.TLatex (minXused + (maxXused-minXused)*3./4., 0 + (maxYused-0)*3.9/4., legend['sqrt'])
#flag_sqrt.Draw()
#if 'sqrt' in legend.keys() and 'lumi' in legend.keys():
#flag_lumi_sqrt = ROOT.TLatex (minXused + (maxXused-minXused)*2.5/4., 0 + (maxYused-0)*3.9/4., "#splitline{CMS preliminary}{#splitline{" + legend['lumi'] + "}{" + legend['sqrt'] + "} }")
#flag_lumi_sqrt.Draw()
CMS_lumi.CMS_lumi(tcanvasDifference, iPeriod, iPos)
# draw back all the axes
#frameDistro.Draw("AXIS")
pad1difference.RedrawAxis()
tcanvasDifference.cd()
canvasPad2differenceName = 'pad2difference_' + cutName + "_" + variableName
pad2difference = ROOT.TPad(canvasPad2differenceName,canvasPad2differenceName,0,0,1,1-0.72)
pad2difference.SetTopMargin(0.000)
pad2difference.SetBottomMargin(0.392)
pad2difference.Draw()
#pad2difference.cd().SetGrid()
pad2difference.cd()
canvasFrameDifferenceName = 'frame_difference_' + cutName + "_" + variableName
if self._showRelativeRatio :
frameDifference = pad2difference.DrawFrame(minXused, -1.0 , maxXused, 1.0 , canvasFrameDifferenceName)
else :
frameDifference = pad2difference.DrawFrame(minXused, int ( ROOT.TMath.MinElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) - 2 ), maxXused, int ( ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) + 2 ), canvasFrameDifferenceName)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = frameDifference.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
if 'xaxis' in variable.keys() :
frameDifference.GetXaxis().SetTitle(variable['xaxis'])
else :
frameDifference.GetXaxis().SetTitle(variableName)
if self._showRelativeRatio :
frameDifference.GetYaxis().SetTitle("#frac{Data - Expected}{Expected}")
frameDifference.GetYaxis().SetRangeUser( -1.0, 1.0 )
else :
frameDifference.GetYaxis().SetTitle("Data - Expected")
frameDifference.GetYaxis().SetRangeUser( int (ROOT.TMath.MinElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) - 2 ), int (ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) + 2 ) )
self.Pad2TAxis(frameDifference)
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMCMinusMC.SetLineColor(12)
tgrMCMinusMC.SetFillColor(12)
tgrMCMinusMC.SetLineWidth(2)
tgrMCMinusMC.SetFillStyle(3004)
tgrMCMinusMC.Draw("2")
tgrDataMinusMC.Draw("P0")
#print(" tgrDataMinusMC.GetMinimum()/Max = " , tgrDataMinusMC.GetMinimum(), " / " , tgrDataMinusMC.GetMaximum())
for samplesToDifferenceGrName, samplesGrToDifference in tgrDifferenceList.items() :
samplesGrToDifference.Draw("P")
oneLine2 = ROOT.TLine(frameDifference.GetXaxis().GetXmin(), 0, frameDifference.GetXaxis().GetXmax(), 0);
oneLine2.SetLineStyle(3)
oneLine2.SetLineWidth(3)
oneLine2.Draw("same")
# draw back all the axes
#frameDifference.Draw("AXIS")
pad2difference.RedrawAxis()
pad2difference.SetGrid()
if 'cdifference' in self._plotsToWrite:
if self._plotLinear:
self._saveCanvas(tcanvasDifference, self._outputDirPlots + "/" + canvasDifferenceNameTemplate + self._FigNamePF)
if self._plotLog:
# log Y axis
#frameDistro.GetYaxis().SetRangeUser( max(self._minLogCdifference, maxYused/1000), self._maxLogCdifference * maxYused )
frameDistro.GetYaxis().SetRangeUser( min(self._minLogC, maxYused/1000), self._maxLogCdifference * maxYused )
pad1difference.SetLogy(True)
self._saveCanvas(tcanvasDifference, self._outputDirPlots + "/log_" + canvasDifferenceNameTemplate + self._FigNamePF, imageOnly=self._plotLinear)
pad1difference.SetLogy(False)
if 'root' in self._fileFormats:
text_file_html.write(canvasDifferenceNameTemplate + ".root;\n")
# ~~~~~~~~~~~~~~~~~~~~
# plot with difference plot Fancy : data - bkg subtracted
#
# only IF we have the difference and not the ratio --> _showRelativeRatio
#
if self._plotFancy and not self._showRelativeRatio:
print("- draw with difference Fancy")
#blind data
if 'blind' in variable:
blind_range = variable['blind'][cutName]
b0 = histos[sampleName].FindBin(blind_range[0])
b1 = histos[sampleName].FindBin(blind_range[1])
for ip in range(tgrDataMinusMC.GetN(), tgrDataMinusMC.GetN()-(b1-b0+1),-1):
tgrDataMinusMC.RemovePoint(ip)
canvasDifferenceNameTemplate = 'cdifference_' + cutName + "_" + variableName + "_Fancy"
tcanvasDifference_Fancy.cd()
canvasPad1differenceName = 'pad1difference_' + cutName + "_" + variableName + "_Fancy"
pad1difference = ROOT.TPad(canvasPad1differenceName,canvasPad1differenceName, 0., 0., 1., 1.)
pad1difference.Draw()
pad1difference.cd()
canvasFrameDistroName = 'frame_fancy_' + cutName + "_" + variableName
frameDistro_Fancy = pad1difference.DrawFrame(minXused,
int (ROOT.TMath.MinElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) - int ( ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetEYlow ()) ) - 20 ),
maxXused,
int (ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) + int ( ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetEYhigh()) ) + 20 ),
canvasFrameDistroName)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = frameDistro_Fancy.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
frameDistro_Fancy.GetYaxis().SetMaxDigits(2)
xaxis = frameDistro_Fancy.GetXaxis()
xaxis.SetLabelFont ( 42)
xaxis.SetLabelOffset( 0.015)
xaxis.SetLabelSize ( 0.035)
xaxis.SetNdivisions ( 505)
xaxis.SetTitleFont ( 42)
xaxis.SetTitleOffset( 1.35)
xaxis.SetTitleSize ( 0.035)
yaxis = frameDistro_Fancy.GetYaxis()
yaxis.SetLabelFont ( 42)
yaxis.SetLabelOffset( 0.01)
yaxis.SetLabelSize ( 0.035)
yaxis.SetNdivisions ( 505)
yaxis.SetTitleFont ( 42)
yaxis.SetTitleOffset( 1.55)
yaxis.SetTitleSize ( 0.045)
if 'xaxis' in variable.keys() :
frameDistro_Fancy.GetXaxis().SetTitle(variable['xaxis'])
if variable["divideByBinWidth"] == 1:
if 'yaxis' in variable.keys() :
frameDistro_Fancy.GetYaxis().SetTitle("Data - Bkg " + variable['yaxis'])
else:
if "GeV" in variable['xaxis']:
### FIXME: it's maybe better to add a "yaxis" field in the variable to let the user choose the y axis name
frameDistro_Fancy.GetYaxis().SetTitle("Data - Bkg dN/d"+variable['xaxis'].replace("GeV","GeV^{-1}"))
else:
frameDistro_Fancy.GetYaxis().SetTitle("Data - Bkg dN/d"+variable['xaxis'])
else:
if 'yaxis' in variable.keys() :
frameDistro_Fancy.GetYaxis().SetTitle("Data - Bkg " + variable['yaxis'])
else :
frameDistro_Fancy.GetYaxis().SetTitle("Data - Bkg Events")
else :
frameDistro_Fancy.GetXaxis().SetTitle(variableName)
if variable["divideByBinWidth"] == 1:
frameDistro_Fancy.GetYaxis().SetTitle("Data - Bkg dN/d"+variableName)
else:
if 'yaxis' in variable.keys() :
frameDistro_Fancy.GetYaxis().SetTitle("Data - Expected " + variable['yaxis'])
else :
frameDistro_Fancy.GetYaxis().SetTitle("Data - Expected Events")
frameDistro_Fancy.GetYaxis().SetRangeUser( int (ROOT.TMath.MinElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) - int ( ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetEYlow ()) ) - 20 ),
int (ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetY()) + int ( ROOT.TMath.MaxElement(tgrDataMinusMC.GetN(),tgrDataMinusMC.GetEYhigh()) ) + 20 ) )
# if there is "histo_total" there is no need of explicit nuisances
if (not self._removeMCStat) or len(mynuisances.keys()) != 0 or histo_total!= None:
tgrMCMinusMC.SetLineColor(12)
tgrMCMinusMC.SetFillColor(12)
tgrMCMinusMC.SetLineWidth(2)
tgrMCMinusMC.SetFillStyle(3004)
tgrMCMinusMC.Draw("2")
#---- the Legend
special_tlegend = ROOT.TLegend(0.40, 0.75, 0.90, 0.90)
special_tlegend.SetFillColor(0)
special_tlegend.SetTextFont(42)
special_tlegend.SetTextSize(0.035)
special_tlegend.SetLineColor(0)
special_tlegend.SetShadowColor(0)
special_tlegend.AddEntry( tgrDataMinusMC , 'Data - Bkg', "PL")
special_tlegend.AddEntry( tgrMCMinusMC, "Tot uncertainty", "F")
if self._showDataMinusBkgOnly :
tgrMCSigMinusMCBkg.SetLineWidth(3)
tgrMCSigMinusMCBkg.SetLineColor(2) # red
tgrMCSigMinusMCBkg.SetMarkerColor(2) # red
tgrMCSigMinusMCBkg.SetMarkerSize(0)
tgrMCSigMinusMCBkg.Draw("P")
special_tlegend.AddEntry( tgrMCSigMinusMCBkg , 'Signal', "PL")
# draw the data - MC
# if blind remove the last points
# maxip = 0
# for ip in range(0, tgrDataMinusMC.GetN()):
# x = tgrDataMinusMC.GetPointX(ip)
# if x > 0.6:
# maxip = ip
tgrDataMinusMC.Draw("P")
oneLine2 = ROOT.TLine(frameDistro_Fancy.GetXaxis().GetXmin(), 0, frameDistro_Fancy.GetXaxis().GetXmax(), 0);
oneLine2.SetLineStyle(3)
oneLine2.SetLineWidth(3)
oneLine2.Draw("same")
special_tlegend.Draw()
CMS_lumi.CMS_lumi(tcanvasDifference_Fancy, iPeriod, iPos)
# draw back all the a xes
pad1difference.RedrawAxis()
pad1difference.SetGrid()
if 'cdifference' in self._plotsToWrite:
self._saveCanvas(tcanvasDifference_Fancy, self._outputDirPlots + "/" + canvasDifferenceNameTemplate + self._FigNamePF)
if 'root' in self._fileFormats:
text_file_html.write(canvasDifferenceNameTemplate + ".root;\n")
#
# draw weighted plot
#
if 'doWeight' in variable.keys() :
if variable['doWeight'] == 1 :
if 'binX' in variable.keys() and 'binY' in variable.keys() :
nbinX = variable['binX']
nbinY = variable['binY']
#
# Add weight 1D
# - sample the 1D histogram in 'nbinX' slices long 'nbinY'.
# - calculate the integral of S and B for each of them
# - add a weight on them S/B, for signal, background and data
# - add the weighted histogram to the final histogram, made of 'nbinY' bins
#
# trick to exchange X <-> Y
if self._invertXY :
temp_nbinX = nbinX
nbinX = nbinY
nbinY = temp_nbinX
#
# check if I have to remove the signal on the ackground stack here
# --> ok, it works and it is correct
#
if thsBackground.GetNhists() != 0 and thsSignal.GetNhists() != 0 :
weight_X_thsData = ROOT.THStack ("weight_X_thsData", "weight_X_thsData")
weight_X_thsSignal = ROOT.THStack ("weight_X_thsSignal", "weight_X_thsSignal")
weight_X_thsBackground = ROOT.THStack ("weight_X_thsBackground","weight_X_thsBackground")
#
# the final histogram should be scaled such that the integral
# is the expected signal + background yield
# --> or only background ?
#
if len(groupPlot.keys()) == 0:
totalBkg = thsBackground.GetStack().Last().Integral()
totalSig = thsSignal.GetStack().Last().Integral()
else :
totalBkg = thsBackground_grouped.GetStack().Last().Integral()
totalSig = thsSignal_grouped.GetStack().Last().Integral()
totalBkgSig = totalBkg + totalSig
totalWeightedIntegralBkg = 0.
totalWeightedIntegralSig = 0.
weight_X_list_Data = []
weight_X_list_Signal = []
weight_X_list_Background = []
weight_X_list_weights = []
for sliceX in range(nbinX) :
integral_bkg = 0.
#for ibin in range( thsBackground.GetStack().Last().GetNbinsX() )
for ibin in range( nbinY ) :
if len(groupPlot.keys()) == 0:
if self._invertXY :
integral_bkg += thsBackground.GetStack().Last().GetBinContent(ibin*nbinX +1 + sliceX)
else :
integral_bkg += thsBackground.GetStack().Last().GetBinContent(ibin+1 + sliceX * nbinY)
else :
if self._invertXY :
integral_bkg += thsBackground_grouped.GetStack().Last().GetBinContent(ibin*nbinX +1 + sliceX)
else :
integral_bkg += thsBackground_grouped.GetStack().Last().GetBinContent(ibin+1 + sliceX * nbinY)
integral_sig = 0.
if len(groupPlot.keys()) == 0:
for ibin in range( thsSignal.GetStack().Last().GetNbinsX() ) :
if self._invertXY :
integral_sig += thsSignal.GetStack().Last().GetBinContent(ibin*nbinX +1 + sliceX)
else :
integral_sig += thsSignal.GetStack().Last().GetBinContent(ibin+1 + sliceX * nbinY)
else :
for ibin in range( thsSignal_grouped.GetStack().Last().GetNbinsX() ) :
if self._invertXY :
integral_sig += thsSignal_grouped.GetStack().Last().GetBinContent(ibin*nbinX +1 + sliceX)
else :
integral_sig += thsSignal_grouped.GetStack().Last().GetBinContent(ibin+1 + sliceX * nbinY)
# this is because the signal was added into the background stack before
integral_bkg = integral_bkg - integral_sig
weight = 1
if integral_bkg != 0 :
weight = integral_sig / integral_bkg
else :
weight = 1
#
# remove weight: use just 1 for each line,
# meaning we are just adding the bins together
#
if self._removeWeight == True :
weight = 1
weight_X_list_weights.append(weight)
if len(groupPlot.keys()) == 0:
for ihisto in range(thsBackground.GetNhists()):
hentry = thsBackground.GetHists().At(ihisto)
histo = ROOT.TH1F('h_weigth_X_' + cutName + '_' + variableName + '_' + hentry.GetName() + '_slice_'+ str(sliceX), '-' , nbinY, 0, nbinY)
histo = self.FixBins (histo, tgrData_vx, tgrData_evx)
histo.SetFillColor( hentry.GetFillColor())
histo.SetFillStyle( hentry.GetFillStyle())
histo.SetLineColor( hentry.GetLineColor())
histo.SetLineWidth( hentry.GetLineWidth())
for ibin in range( nbinY ) :
if self._invertXY :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin*nbinX +1 + sliceX) ) )
else :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin+1 + sliceX * nbinY) ) )
if sliceX != 0:
weight_X_list_Background[ihisto].Add(histo)
else :
weight_X_list_Background.append(histo)
# the minus signal is because the signal was added into the background stack before
if ihisto < (thsBackground.GetNhists() - thsSignal.GetNhists()) :
totalWeightedIntegralBkg += histo.Integral()
for ihisto in range(thsSignal.GetNhists()):
hentry = thsSignal.GetHists().At(ihisto)
histo = ROOT.TH1F('h_weigth_X_' + cutName + '_' + variableName + '_' + (()).GetName() + '_slice_'+ str(sliceX), "-", nbinY, 0, nbinY)
histo = self.FixBins (histo, tgrData_vx, tgrData_evx)
histo.SetFillColor( hentry.GetFillColor())
histo.SetFillStyle( hentry.GetFillStyle())
histo.SetLineColor( hentry.GetLineColor())
histo.SetLineWidth( hentry.GetLineWidth())
for ibin in range( nbinY ) :
if self._invertXY :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin*nbinX +1 + sliceX) ) )
else :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin+1 + sliceX * nbinY) ) )
if sliceX != 0:
weight_X_list_Signal[ihisto].Add(histo)
else :
weight_X_list_Signal.append(histo)
totalWeightedIntegralSig += histo.Integral()
for ihisto in range(thsData.GetNhists()) :
hentry = thsData.GetHists().At(ihisto)
histo = ROOT.TH1F('h_weigth_X_' + cutName + '_' + variableName + '_' + hentry.GetName() + '_slice_'+ str(sliceX), "-", nbinY, 0, nbinY)
histo = self.FixBins (histo, tgrData_vx, tgrData_evx)
for ibin in range( nbinY ) :
if self._invertXY :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin*nbinX +1 + sliceX ) ) )
else :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin+1 + sliceX * nbinY) ) )
if sliceX != 0:
weight_X_list_Data[ihisto].Add(histo)
else :
weight_X_list_Data.append(histo)
#weight_X_list_Data.append(histo) ## aaaargh!
else :
for ihisto in range(thsBackground_grouped.GetNhists()):
hentry = thsBackground_grouped.GetHists().At(ihisto)
histo = ROOT.TH1F('h_weigth_X_' + cutName + '_' + variableName + '_' + hentry.GetName() + '_slice_'+ str(sliceX), '-' , nbinY, 0, nbinY)
histo = self.FixBins (histo, tgrData_vx, tgrData_evx)
histo.SetFillColor( hentry.GetFillColor())
histo.SetFillStyle( hentry.GetFillStyle())
histo.SetLineColor( hentry.GetLineColor())
histo.SetLineWidth( hentry.GetLineWidth())
for ibin in range( nbinY ) :
if self._invertXY :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin*nbinX +1 + sliceX) ) )
else :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin+1 + sliceX * nbinY) ) )
if sliceX != 0:
weight_X_list_Background[ihisto].Add(histo)
else :
weight_X_list_Background.append(histo)
# the minus signal is because the signal was added into the background stack before
if ihisto < (thsBackground_grouped.GetNhists() - thsSignal_grouped.GetNhists()) :
totalWeightedIntegralBkg += histo.Integral()
for ihisto in range(thsSignal_grouped.GetNhists()):
hentry = thsSignal_grouped.GetHists().At(ihisto)
histo = ROOT.TH1F('h_weigth_X_' + cutName + '_' + variableName + '_' + hentry.GetName() + '_slice_'+ str(sliceX), "-", nbinY, 0, nbinY)
histo = self.FixBins (histo, tgrData_vx, tgrData_evx)
histo.SetFillColor( hentry.GetFillColor())
histo.SetFillStyle( hentry.GetFillStyle())
histo.SetLineColor( hentry.GetLineColor())
histo.SetLineWidth( hentry.GetLineWidth())
for ibin in range( nbinY ) :
if self._invertXY :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin*nbinX +1 + sliceX) ) )
else :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin+1 + sliceX * nbinY) ) )
if sliceX != 0:
weight_X_list_Signal[ihisto].Add(histo)
else :
weight_X_list_Signal.append(histo)
totalWeightedIntegralSig += histo.Integral()
for ihisto in range(thsData.GetNhists()):
hentry = thsData.GetHists().At(ihisto)
histo = ROOT.TH1F('h_weigth_X_' + cutName + '_' + variableName + '_' + hentry.GetName() + '_slice_'+ str(sliceX), "-", nbinY, 0, nbinY)
histo = self.FixBins (histo, tgrData_vx, tgrData_evx)
for ibin in range( nbinY ) :
if self._invertXY :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin*nbinX +1 + sliceX) ) )
else :
histo.SetBinContent(ibin+1, weight * ( hentry.GetBinContent(ibin+1 + sliceX * nbinY) ) )
if sliceX != 0:
weight_X_list_Data[ihisto].Add(histo)
else :
weight_X_list_Data.append(histo)
#weight_X_list_Data.append(histo) ## aaaargh!
#
# gloabal scale factor so that the total number of events is such that
# the expected signal events is unchanged
#
#global_normalization = totalBkg / totalWeightedIntegralBkg
global_normalization = totalSig / totalWeightedIntegralSig if totalWeightedIntegralSig != 0 else 1.0
for histo in weight_X_list_Data:
histo.Scale(global_normalization)
for histo in weight_X_list_Background:
histo.Scale(global_normalization)
for histo in weight_X_list_Signal:
histo.Scale(global_normalization)
for histo in weight_X_list_Data:
weight_X_thsData.Add(histo)
for histo in weight_X_list_Background:
weight_X_thsBackground.Add(histo)
for histo in weight_X_list_Signal:
weight_X_thsSignal.Add(histo)
#print(" weight_X_list_weights = ", weight_X_list_weights)
# create the weighted data distribution
weight_X_tgrData = ROOT.TGraphAsymmErrors()
for sliceX in range( nbinX ) :
for ibin in range( nbinY ) :
x = tgrData_vx[ibin]
#print(" sliceX, ibin = ", sliceX, " , ", ibin, " -> ", tgrData_vx[ibin])
number_of_bin = ibin + sliceX * nbinY
if self._invertXY :
number_of_bin = ibin*nbinX + sliceX
y = weight_X_list_weights[sliceX] * global_normalization * tgrData_vy[number_of_bin]
exlow = tgrData_evx[number_of_bin]
exhigh = tgrData_evx[number_of_bin]
eylow = weight_X_list_weights[sliceX] * global_normalization * tgrData_evy_do[number_of_bin]
eyhigh = weight_X_list_weights[sliceX] * global_normalization * tgrData_evy_up[number_of_bin]
if sliceX != 0:
y += weight_X_tgrData.GetY()[ibin]
eylow = self.SumQ ( eylow, weight_X_tgrData.GetErrorYlow(ibin) )
eyhigh = self.SumQ ( eyhigh, weight_X_tgrData.GetErrorYhigh(ibin) )
#print(" eylow = ", eylow, " eyhigh = ", eyhigh, " y = ", y, " x = ", x )
weight_X_tgrData.SetPoint (ibin, x, y)
weight_X_tgrData.SetPointError (ibin, exlow, exhigh, eylow, eyhigh)
#if sliceX == ( nbinX - 1) :
#print(" ibin,x,y = ", ibin, ", ", x, ", ", y)
# create the weighted data distribution
weight_X_tgrMC = ROOT.TGraphAsymmErrors()
for sliceX in range(nbinX) :
for ibin in range( nbinY ) :
x = tgrData_vx[ibin]
number_of_bin = ibin + sliceX * nbinY
if self._invertXY :
number_of_bin = ibin*nbinX + sliceX
y = weight_X_list_weights[sliceX] * global_normalization * tgrMC_vy[number_of_bin]
exlow = tgrData_evx[number_of_bin]
exhigh = tgrData_evx[number_of_bin]
eylow = 0
eyhigh = 0
if len(nuisances_err_do) != 0:
if histo_total:
eylow = weight_X_list_weights[sliceX] * global_normalization * histo_total.GetBinError(number_of_bin + 1)
eyhigh = weight_X_list_weights[sliceX] * global_normalization * histo_total.GetBinError(number_of_bin + 1)
else :
eylow = weight_X_list_weights[sliceX] * global_normalization * nuisances_err_do[number_of_bin]
eyhigh = weight_X_list_weights[sliceX] * global_normalization * nuisances_err_up[number_of_bin]
if sliceX != 0 :
y += weight_X_tgrMC.GetY()[ibin]
eylow = self.SumQ ( eylow, weight_X_tgrMC.GetErrorYlow(ibin) )
eyhigh = self.SumQ ( eyhigh, weight_X_tgrMC.GetErrorYhigh(ibin) )
#print(" eylow = ", eylow, " eyhigh = ", eyhigh, " y = ", y, " x = ", x )
weight_X_tgrMC.SetPoint (ibin, x, y)
weight_X_tgrMC.SetPointError (ibin, exlow, exhigh, eylow, eyhigh)
#
# create the weighted data over MC distribution
#
weight_X_tgrDataOverMC = weight_X_tgrData.Clone("tgrDataOverMCweighted")
last = weight_X_thsBackground.GetStack().Last()
for ibin in range( nbinY ) :
x = weight_X_tgrDataOverMC.GetX()[ibin]
y = self.Ratio(weight_X_tgrData.GetY()[ibin] , last.GetBinContent(ibin+1) )
number_of_bin = ibin + sliceX * nbinY
if self._invertXY :
number_of_bin = ibin*nbinX + sliceX
exlow = tgrData_evx[number_of_bin]
exhigh = tgrData_evx[number_of_bin]
eylow = self.Ratio(weight_X_tgrData.GetErrorYlow(ibin), last.GetBinContent(ibin+1) )
eyhigh = self.Ratio(weight_X_tgrData.GetErrorYhigh(ibin), last.GetBinContent(ibin+1) )
weight_X_tgrDataOverMC.SetPoint (ibin, x, y)
weight_X_tgrDataOverMC.SetPointError (ibin, exlow, exhigh, eylow, eyhigh)
#print(" Ratio:: ibin,x,y = ", ibin, ", ", x, ", ", y, ", ", eylow, ", ", eyhigh, " <-- ", weight_X_tgrData.GetY()[ibin], " / ", last.GetBinContent(ibin+1) )
#
# create the weighted data minus MC distribution
#
weight_X_tgrDataMinusMC = weight_X_tgrData.Clone("tgrDataMinusMCweighted")
last = weight_X_thsBackground.GetStack().Last()
for ibin in range( nbinY ) :
x = weight_X_tgrDataMinusMC.GetX()[ibin]
if self._showRelativeRatio :
y = self.Ratio ( self.Difference(weight_X_tgrData.GetY()[ibin] , last.GetBinContent(ibin+1) ) , last.GetBinContent(ibin+1) )
else :
# if show only "data - bkg", subtract "data - (bkg+sig) + sig "
if self._showDataMinusBkgOnly :
y = self.Difference(weight_X_tgrData.GetY()[ibin] , last.GetBinContent(ibin+1) ) + weight_X_thsSignal.GetStack().Last().GetBinContent(ibin+1)
else :
y = self.Difference(weight_X_tgrData.GetY()[ibin] , last.GetBinContent(ibin+1) )
number_of_bin = ibin + sliceX * nbinY
if self._invertXY :
number_of_bin = ibin*nbinX + sliceX
exlow = tgrData_evx[number_of_bin]
exhigh = tgrData_evx[number_of_bin]
if self._showRelativeRatio :
eylow = self.Ratio(weight_X_tgrData.GetErrorYlow(ibin), last.GetBinContent(ibin+1) )
eyhigh = self.Ratio(weight_X_tgrData.GetErrorYhigh(ibin), last.GetBinContent(ibin+1) )
else :
eylow = weight_X_tgrData.GetErrorYlow(ibin)
eyhigh = weight_X_tgrData.GetErrorYhigh(ibin)
weight_X_tgrDataMinusMC.SetPoint (ibin, x, y)
weight_X_tgrDataMinusMC.SetPointError (ibin, exlow, exhigh, eylow, eyhigh)
#
# create the weighted MC over MC distribution
#
weight_X_tgrMCOverMC = weight_X_tgrData.Clone("tgrMCOverMCweighted")
last = weight_X_thsBackground.GetStack().Last()
for ibin in range( nbinY ) :
x = weight_X_tgrMCOverMC.GetX()[ibin]
y = 1
number_of_bin = ibin + sliceX * nbinY
if self._invertXY :
number_of_bin = ibin*nbinX + sliceX
exlow = tgrData_evx[number_of_bin]
exhigh = tgrData_evx[number_of_bin]
eylow = self.Ratio(weight_X_tgrMC.GetErrorYlow(ibin), last.GetBinContent(ibin+1) )
eyhigh = self.Ratio(weight_X_tgrMC.GetErrorYhigh(ibin), last.GetBinContent(ibin+1) )
weight_X_tgrMCOverMC.SetPoint (ibin, x, y)
weight_X_tgrMCOverMC.SetPointError (ibin, exlow, exhigh, eylow, eyhigh)
#
# create the weighted MC over MC distribution
#
weight_X_tgrMCMinusMC = weight_X_tgrData.Clone("tgrMCMinusMCweighted")
last = weight_X_thsBackground.GetStack().Last()
for ibin in range( nbinY ) :
x = weight_X_tgrMCMinusMC.GetX()[ibin]
y = 0
number_of_bin = ibin + sliceX * nbinY
if self._invertXY :
number_of_bin = ibin*nbinX + sliceX
exlow = tgrData_evx[number_of_bin]
exhigh = tgrData_evx[number_of_bin]
if self._showRelativeRatio :
eylow = self.Ratio(weight_X_tgrMC.GetErrorYlow(ibin), last.GetBinContent(ibin+1) )
eyhigh = self.Ratio(weight_X_tgrMC.GetErrorYhigh(ibin), last.GetBinContent(ibin+1) )
else :
eylow = weight_X_tgrMC.GetErrorYlow(ibin)
eyhigh = weight_X_tgrMC.GetErrorYhigh(ibin)
#print(" DIFF:: eylow = ", eylow, " eyhigh = ", eyhigh, " y = ", y, " x = ", x )
weight_X_tgrMCMinusMC.SetPoint (ibin, x, y)
weight_X_tgrMCMinusMC.SetPointError (ibin, exlow, exhigh, eylow, eyhigh)
#
# now plot
#
# - recalculate the maxY
# _maxLinearScale --> 1.45 in the past
maxYused = self._maxLinearScale * self.GetMaximumIncludingErrors(last)
# recalculate min-max X due to weighting rolling
minXused = weight_X_tgrMCMinusMC.GetX()[0] - tgrData_evx[0]
maxXused = weight_X_tgrMCMinusMC.GetX()[nbinY-1] + tgrData_evx[nbinY-1]
weight_X_canvasRatioNameTemplate = 'cratio_weight_X_' + cutName + '_' + variableName
weight_X_tcanvasRatio.cd()
canvasPad1Name = 'weight_X_pad1_' + cutName + "_" + variableName
weight_X_pad1 = ROOT.TPad(canvasPad1Name,canvasPad1Name, 0, 1-0.72, 1, 1)
weight_X_pad1.SetTopMargin(0.098)
weight_X_pad1.SetBottomMargin(0.000)
weight_X_pad1.Draw()
weight_X_pad1.cd()
weight_X_canvasFrameDistroName = 'weight_X_frame_distro_' + cutName + "_" + variableName
weight_X_frameDistro = weight_X_pad1.DrawFrame(minXused, 0.0, maxXused, 1.0, weight_X_canvasFrameDistroName)
#weight_X_frameDistro = weight_X_pad1.DrawFrame(minXused, 0.0, maxXused, 1.0, weight_X_canvasFrameDistroName)
#weight_X_frameDistro = weight_X_pad1.DrawFrame(0.0, 0.0, nbinY, 1.0, weight_X_canvasFrameDistroName)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = weight_X_frameDistro.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
if 'xaxis' in variable.keys() :
weight_X_frameDistro.GetXaxis().SetTitle(variable['xaxis'])
else :
weight_X_frameDistro.GetXaxis().SetTitle(variableName)
weight_X_frameDistro.GetYaxis().SetTitle("S/B weighted Events")
if self._removeWeight == True :
weight_X_frameDistro.GetYaxis().SetTitle("Events")
weight_X_frameDistro.GetYaxis().SetRangeUser( min(0.001, minYused), maxYused )
if weight_X_thsBackground.GetNhists() != 0:
weight_X_thsBackground.Draw("hist same")
if weight_X_thsSignal.GetNhists() != 0:
weight_X_thsSignal.Draw("hist same noclear")
# if there is "histo_total" there is no need of explicit nuisances
if len(mynuisances.keys()) != 0 or histo_total!= None:
weight_X_tgrMC.SetLineColor(12)
weight_X_tgrMC.SetFillColor(12)
weight_X_tgrMC.SetFillStyle(3004)
weight_X_tgrMC.Draw("2")
#weight_X_tgrMC.Draw("P0")
#print(" -------------------------> here ")
# - then the DATA
if weight_X_tgrData.GetN() != 0:
weight_X_tgrData.Draw("P0")
#print(" -------------------------> here data ")
tlegend.Draw()
CMS_lumi.CMS_lumi(weight_X_tcanvasRatio, iPeriod, iPos)
# draw back all the axes
#weight_X_frameDistro.Draw("AXIS")
weight_X_pad1.RedrawAxis()
weight_X_tcanvasRatio.cd()
canvasPad2Name = 'weight_X_weight_X_pad2_' + cutName + "_" + variableName
weight_X_pad2 = ROOT.TPad(canvasPad2Name,canvasPad2Name,0,0,1,1-0.72)
weight_X_pad2.SetTopMargin(0.000)
weight_X_pad2.SetBottomMargin(0.392)
weight_X_pad2.Draw()
#weight_X_pad2.cd().SetGrid()
weight_X_pad2.cd()
weight_X_canvasFrameRatioName = 'weight_X_frame_ratio_' + cutName + "_" + variableName
#weight_X_frameRatio = weight_X_pad2.DrawFrame(minXused, 0.0, nbinY, 2.0, weight_X_canvasFrameRatioName)
weight_X_frameRatio = weight_X_pad2.DrawFrame(minXused, 0.0, maxXused, 2.0, weight_X_canvasFrameRatioName)
#print(" minXused = " , minXused)
#print(" maxXused = " , maxXused)
#print(" nbinY = " , nbinY)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = weight_X_frameRatio.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
if 'xaxis' in variable.keys() :
weight_X_frameRatio.GetXaxis().SetTitle(variable['xaxis'])
else :
weight_X_frameRatio.GetXaxis().SetTitle(variableName)
weight_X_frameRatio.GetYaxis().SetTitle("Data/Expected")
weight_X_frameRatio.GetYaxis().SetRangeUser( 0.5, 1.5 )
self.Pad2TAxis(weight_X_frameRatio)
# if there is "histo_total" there is no need of explicit nuisances
if len(mynuisances.keys()) != 0 or histo_total!= None:
weight_X_tgrMCOverMC.SetLineColor(12)
weight_X_tgrMCOverMC.SetFillColor(12)
weight_X_tgrMCOverMC.SetFillStyle(3004)
weight_X_tgrMCOverMC.Draw("2")
weight_X_tgrDataOverMC.Draw("P0")
oneLine2 = ROOT.TLine(weight_X_frameRatio.GetXaxis().GetXmin(), 1, weight_X_frameRatio.GetXaxis().GetXmax(), 1);
oneLine2.SetLineStyle(3)
oneLine2.SetLineWidth(3)
oneLine2.Draw("same")
# draw back all the axes
#weight_X_frameRatio.Draw("AXIS")
weight_X_pad2.RedrawAxis()
self._saveCanvas(weight_X_tcanvasRatio, self._outputDirPlots + "/" + weight_X_canvasRatioNameTemplate + self._FigNamePF)
if 'root' in self._fileFormats:
text_file_html.write(weight_X_canvasRatioNameTemplate + ".root;\n")
# save also all the TH1F separately for later combination
temp_file = ROOT.TFile (self._outputDirPlots + "/" + weight_X_canvasRatioNameTemplate + self._FigNamePF + ".root", "UPDATE")
histo_global_normalization = ROOT.TH1F("histo_global_normalization", "", 1, 0, 1)
histo_global_normalization.Fill(0.5, global_normalization)
histo_global_normalization.Write()
weight_X_tgrMCOverMC.Write()
weight_X_tgrDataOverMC.Write()
# if there is "histo_total" there is no need of explicit nuisances
if len(mynuisances.keys()) != 0 or histo_total!= None:
weight_X_tgrMC.Write("weight_X_tgrMC")
if weight_X_tgrData.GetN() != 0:
weight_X_tgrData.Write("weight_X_tgrData")
if weight_X_thsBackground.GetNhists() != 0:
weight_X_thsBackground.Write()
if weight_X_thsSignal.GetNhists() != 0:
weight_X_thsSignal.Write()
for histo in weight_X_list_Data:
histo.Write()
for histo in weight_X_list_Background:
histo.Write()
for histo in weight_X_list_Signal:
histo.Write()
temp_file.Close()
# log Y axis
if self._plotLog:
weight_X_frameDistro.GetYaxis().SetRangeUser( min(0.001, maxYused/1000), 10 * maxYused )
weight_X_pad1.SetLogy(True)
self._saveCanvas(weight_X_tcanvasRatio, self._outputDirPlots + "/log_" + weight_X_canvasRatioNameTemplate, imageOnly=self._plotLinear)
weight_X_pad1.SetLogy(False)
#
# Now plot difference
#
weight_X_pad2.cd()
if self._showRelativeRatio :
weight_X_frameRatio = weight_X_pad2.DrawFrame(minXused, -1.0, maxXused, 1.0, weight_X_canvasFrameRatioName)
else :
weight_X_frameRatio = weight_X_pad2.DrawFrame(minXused, int( ROOT.TMath.MinElement(weight_X_tgrDataMinusMC.GetN(),weight_X_tgrDataMinusMC.GetY()) - 2 ), maxXused, int ( ROOT.TMath.MaxElement(weight_X_tgrDataMinusMC.GetN(),weight_X_tgrDataMinusMC.GetY()) + 2 ), weight_X_canvasFrameRatioName)
# style from https://ghm.web.cern.ch/ghm/plots/MacroExample/myMacro.py
xAxisDistro = weight_X_frameRatio.GetXaxis()
xAxisDistro.SetNdivisions(6,5,0)
if 'xaxis' in variable.keys() :
weight_X_frameRatio.GetXaxis().SetTitle(variable['xaxis'])
else :
weight_X_frameRatio.GetXaxis().SetTitle(variableName)
if self._showRelativeRatio :
weight_X_frameRatio.GetYaxis().SetRangeUser( -1.0 , 1.0 )
weight_X_frameRatio.GetYaxis().SetTitle("#frac{Data - Expected}{Expected}")
else :
weight_X_frameRatio.GetYaxis().SetRangeUser( int( ROOT.TMath.MinElement(weight_X_tgrDataMinusMC.GetN(),weight_X_tgrDataMinusMC.GetY()) - 2 ), int ( ROOT.TMath.MaxElement(weight_X_tgrDataMinusMC.GetN(),weight_X_tgrDataMinusMC.GetY()) + 2 ) )
weight_X_frameRatio.GetYaxis().SetTitle("Data - Expected")
self.Pad2TAxis(weight_X_frameRatio)
# if there is "histo_total" there is no need of explicit nuisances
if len(mynuisances.keys()) != 0 or histo_total!= None:
weight_X_tgrMCMinusMC.SetLineColor(12)
weight_X_tgrMCMinusMC.SetFillColor(12)
weight_X_tgrMCMinusMC.SetFillStyle(3004)
weight_X_tgrMCMinusMC.Draw("2")
weight_X_tgrDataMinusMC.Draw("P0")
#print(" BINS = " , weight_X_tgrDataMinusMC.GetN())
oneLine2 = ROOT.TLine(weight_X_frameRatio.GetXaxis().GetXmin(), 0, weight_X_frameRatio.GetXaxis().GetXmax(), 0);
oneLine2.SetLineStyle(3)
oneLine2.SetLineWidth(3)
oneLine2.Draw("same")
weight_X_pad2.RedrawAxis()
if self._showRelativeRatio :
weight_X_canvasDifferenceNameTemplate = 'cdifference_relative_weight_X_' + cutName + '_' + variableName
else :
weight_X_canvasDifferenceNameTemplate = 'cdifference_weight_X_' + cutName + '_' + variableName
self._saveCanvas(weight_X_tcanvasRatio, self._outputDirPlots + "/" + weight_X_canvasDifferenceNameTemplate)
if 'root' in self._fileFormats:
text_file_html.write(weight_X_canvasDifferenceNameTemplate + ".root;\n")
#
# This is performed at the end because it will change the "FillStyle" of the histograms
# and you don't want to change it in the previous plots!
# All histograms will become "transparent" as far as fill style is concerned
#
if self._plotNormalizedDistributions :
# ~~~~~~~~~~~~~~~~~~~~
# plot signal vs background normalized
tcanvasSigVsBkg.cd()
frameNorm = ROOT.TH1F
frameNorm = tcanvasSigVsBkg.DrawFrame(minXused, 0.0, maxXused, 1.0)
frameNorm.GetYaxis().SetRangeUser( 0, 1.5 )
# setup axis names
if 'xaxis' in variable.keys() :
frameNorm.GetXaxis().SetTitle(variable['xaxis'])
frameNorm.GetYaxis().SetTitle("a.u.")
tcanvasSigVsBkg.RedrawAxis()
maxY_normalized=0.0
for hentry in thsBackground_grouped.GetHists():
num_bins = hentry.GetNbinsX()
if hentry.Integral() > 0.:
y_normalized = hentry.GetBinContent(hentry.GetMaximumBin())/hentry.Integral()
if y_normalized > maxY_normalized:
maxY_normalized = y_normalized
for ibin in range( num_bins ) :
hentry.SetBinError(ibin+1, 0.000001)
hentry.SetFillStyle(0)
hentry.SetLineWidth(3)
for sampleNameGroup, sampleConfiguration in groupPlot.items():
if ('new_histo_group_' + sampleNameGroup + '_' ) in hentry.GetName() :
if 'line' in sampleConfiguration.keys():
hentry.SetLineStyle(self._getLine(sampleConfiguration['line']))
hentry.DrawNormalized("hist,same")
if thsSignal_grouped.GetNhists() != 0:
for hentry in thsSignal_grouped.GetHists():
num_bins = hentry.GetNbinsX()
if hentry.Integral() > 0.:
y_normalized = hentry.GetBinContent(hentry.GetMaximumBin())/hentry.Integral()
if y_normalized > maxY_normalized:
maxY_normalized = y_normalized
for ibin in range( num_bins ) :
hentry.SetBinError(ibin+1, 0.000001)
hentry.SetFillStyle(0)
hentry.SetLineWidth(3)
for sampleNameGroup, sampleConfiguration in groupPlot.items():
if ('new_histo_group_' + sampleNameGroup + '_' ) in hentry.GetName() :
if 'line' in sampleConfiguration.keys():
hentry.SetLineStyle(self._getLine(sampleConfiguration['line']))
hentry.DrawNormalized("hist,same")
# ~~~~~~~~~~~~~~~~~~~~
# include data only if required
if self._plotNormalizedIncludeData :
for sampleName, plotdef in plot.items():
if plotdef['isData'] == 1 :
if 'line' in plotdef.keys():
histos[sampleName].SetLineStyle(self._getLine(plotdef['line']))
histos[sampleName].DrawNormalized("p, same")
frameNorm.GetYaxis().SetRangeUser(0, 1.8*maxY_normalized)
CMS_lumi.CMS_lumi(tcanvasSigVsBkg, iPeriod, iPos)
tlegend.Draw()
tcanvasSigVsBkg.RedrawAxis()
self._saveCanvas(tcanvasSigVsBkg, self._outputDirPlots + "/" + 'cSigVsBkg_' + cutName + "_" + variableName + self._FigNamePF, imageOnly=True)
if self._plotLog:
# log Y axis
#frameDistro.GetYaxis().SetRangeUser( max(self._minLogCdifference, maxYused/1000), self._maxLogCdifference * maxYused )
frameNorm.GetYaxis().SetRangeUser( min(self._minLogC, maxY_normalized/1000), self._maxLogC * maxY_normalized )
tcanvasSigVsBkg.SetLogy(True)
self._saveCanvas(tcanvasSigVsBkg, self._outputDirPlots + "/log_cSigVsBkg_" + cutName + "_" + variableName + self._FigNamePF, imageOnly=self._plotLinear)
tcanvasSigVsBkg.SetLogy(False)
if self._plotNormalizedDistributionsTHstack :
# ~~~~~~~~~~~~~~~~~~~~
#
# Plot signal vs background normalized
# All the backgrounds or signals will be shown as stacked
# All contributions will be shown as well as in the normal stack distribution
# keeping though the integral of background and signal set to 1
#
tcanvasSigVsBkgTHstack.cd()
frameNormTHstack = ROOT.TH1F
frameNormTHstack = tcanvasSigVsBkgTHstack.DrawFrame(minXused, 0.0, maxXused, 1.0)
frameNormTHstack.GetYaxis().SetRangeUser( 0, 1.5 )
# setup axis names
if 'xaxis' in variable.keys() :
frameNormTHstack.GetXaxis().SetTitle(variable['xaxis'])
tcanvasSigVsBkgTHstack.RedrawAxis()
maxY_normalized=0.0
h_sum_of_backgrounds = thsBackground_grouped.GetStack().Last()
h_sum_of_signals = thsSignal_grouped.GetStack().Last()
normalization_factor_background = 1. / h_sum_of_backgrounds.Integral()
normalization_factor_signal = 1. / h_sum_of_signals.Integral()
if h_sum_of_backgrounds.Integral() > 0.:
maxY_normalized = h_sum_of_backgrounds.GetBinContent(h_sum_of_backgrounds.GetMaximumBin())/h_sum_of_backgrounds.Integral()
if h_sum_of_signals.Integral() > 0.:
temp_maxY_normalized = h_sum_of_signals.GetBinContent(h_sum_of_signals.GetMaximumBin())/h_sum_of_signals.Integral()
if (temp_maxY_normalized > maxY_normalized) :
maxY_normalized = temp_maxY_normalized
for hentry in thsBackground_grouped.GetHists():
if (thsSignal_grouped.GetNhists()==0) or (hentry not in thsSignal_grouped.GetHists()) : # since signal is part of the "background" for plotting reason
num_bins = hentry.GetNbinsX()
for ibin in range( num_bins ) :
hentry.SetBinError(ibin+1, 0.000001)
hentry.SetFillStyle(0)
hentry.SetLineWidth(3)
hentry.Scale(normalization_factor_background)
thsBackground_grouped_normalized.Add(hentry)
if thsSignal_grouped.GetNhists() != 0:
for hentry in thsSignal_grouped.GetHists():
num_bins = hentry.GetNbinsX()
for ibin in range( num_bins ) :
hentry.SetBinError(ibin+1, 0.000001)
hentry.SetFillStyle(0)
hentry.SetLineWidth(3)
hentry.Scale(normalization_factor_signal)
thsSignal_grouped_normalized.Add(hentry)
thsSignal_grouped_normalized.Draw("hist same noclear")
thsBackground_grouped_normalized.Draw("hist same noclear")
frameNormTHstack.GetYaxis().SetRangeUser(0, 1.8*maxY_normalized)
tlegend.Draw()
self._saveCanvas(tcanvasSigVsBkgTHstack, self._outputDirPlots + "/" + 'ccTHstackSigVsBkg_' + cutName + "_" + variableName + self._FigNamePF, imageOnly=True)
# some cleaning
#print(" cleaning ...")
#thsData.Delete()
#print(" cleaning ...")
#thsSignal.Delete()
#print(" cleaning ...")
#thsBackground.Delete()
#print(" cleaning ...")
#thsSignal_grouped.Delete()
#print(" cleaning ...")
#thsBackground_grouped.Delete()
print(" >> end:", variableName)
print(" >> all end")
print(" >> all but really all ")
#
# close plotter.html
#
text_file_html.write(" \"></div> \"\n")
text_file_html.close()
#sys.exit(0)
#quit()
#raise SystemExit()
os._exit(0)
#exit()
# ... or it will remain hanging forever ...
# _____________________________________________________________________________
# --- squared sum
def Pad2TAxis(self, hist):
xaxis = hist.GetXaxis()
xaxis.SetLabelFont ( 42)
xaxis.SetLabelOffset( 0.025)
xaxis.SetLabelSize ( 0.1)
xaxis.SetNdivisions ( 505)
xaxis.SetTitleFont ( 42)
xaxis.SetTitleOffset( 1.35)
xaxis.SetTitleSize ( 0.11)
yaxis = hist.GetYaxis()
yaxis.CenterTitle ( )
yaxis.SetLabelFont ( 42)
yaxis.SetLabelOffset( 0.02)
yaxis.SetLabelSize ( 0.1)
yaxis.SetNdivisions ( 505)
yaxis.SetTitleFont ( 42)
yaxis.SetTitleOffset( .65)
yaxis.SetTitleSize ( 0.11)
# _____________________________________________________________________________
# --- squared sum
def SumQ(self, A, B):
return math.sqrt(A*A + B*B)
# _____________________________________________________________________________
# --- Ratio: if denominator is zero, then put 0!
def Ratio(self, A, B):
if B == 0:
#print("divide by 0")
return 0.
else :
return A / B
# _____________________________________________________________________________
# --- Difference
def Difference(self, A, B):
return A - B
# _____________________________________________________________________________
# --- poissonian error bayesian 1sigma band
# 1/0 1/0
def GetPoissError(self, numberEvents, down, up):
alpha = (1-0.6827)
L = 0
if numberEvents!=0 :
L = ROOT.Math.gamma_quantile (alpha/2,numberEvents,1.)
U = 0
if numberEvents==0 :
#
# Poisson error agreed in the CMS statistics committee
# see: https://hypernews.cern.ch/HyperNews/CMS/get/statistics/263.html
# and https://hypernews.cern.ch/HyperNews/CMS/get/HIG-16-042/32/1/1/1/1/1.html
# and https://twiki.cern.ch/twiki/bin/viewauth/CMS/PoissonErrorBars
# to avoid flip-flop.
# The commented version would have created 1.147 for 0 observed events
# while now we get 1.84 in the case of 0 observed events
#
U = ROOT.Math.gamma_quantile_c (alpha/2,numberEvents+1,1.)
#U = ROOT.Math.gamma_quantile_c (alpha,numberEvents+1,1.)
#print("u = ", U)
else :
U = ROOT.Math.gamma_quantile_c (alpha/2,numberEvents+1,1.)
# the error
L = numberEvents - L
if numberEvents > 0 :
U = U - numberEvents
#else :
#U = 1.14 # --> bayesian interval Poisson with 0 events observed
#1.14790758039 from 10 lines above
if up and not down :
return U
if down and not up :
return L
if up and down :
return (L,U)
# _____________________________________________________________________________
def GetMaximumIncludingErrors(self, histo):
maxWithErrors = 0.
for iBin in range(1, histo.GetNbinsX()+1):
binHeight = histo.GetBinContent (iBin) + histo.GetBinError (iBin)
if binHeight > maxWithErrors :
maxWithErrors = binHeight
return maxWithErrors;
# _____________________________________________________________________________
def GetMinimum(self, histo):
minimum = -1.
for iBin in range(1, histo.GetNbinsX()+1):
binHeight = histo.GetBinContent (iBin)
if binHeight < minimum or minimum<0:
minimum = binHeight
return minimum;
# _____________________________________________________________________________
def defineStyle(self):
print("==================")
import mkShapesRDF.shapeAnalysis.tdrStyle as tdrStyle
tdrStyle.setTDRStyle()
ROOT.TGaxis.SetExponentOffset(-0.08, 0.00,"y")
# _____________________________________________________________________________
# --- fix binning
#
def FixBins(self, histo, reference_x, reference_x_err):
#
# variable bin width
#
nbins = len (reference_x)
#print(" nbins = ", nbins)
binning = [ reference_x[ibin]-reference_x_err[ibin] for ibin in range (0, nbins) ]
#binning = [ reference_histo.GetXaxis().GetBinLowEdge(ibin) for ibin in reference_histo.GetNbinsX()+1 ]
binning.append (reference_x[nbins-1]+reference_x_err[nbins-1])
#print(" >>> histo.GetName() ::", histo.GetName(), " ::> " , binning)
hnew = ROOT.TH1F("new_" + histo.GetName(),"", len(binning)-1, array('d', binning ))
for ibin in range (0, nbins+1) :
y = histo.GetBinContent(ibin)
x = histo.GetXaxis().GetBinCenter(ibin)
hnew.SetBinContent(ibin,y)
hnew.SetFillColor(histo.GetFillColor())
hnew.SetLineColor(histo.GetLineColor())
hnew.SetFillStyle(histo.GetFillStyle())
return hnew
def _saveCanvas(self, tcanvas, nameBase, imageOnly=False):
if 'png' in self._fileFormats:
tcanvas.SaveAs(nameBase + ".png")
if 'pdf' in self._fileFormats:
tcanvas.SaveAs(nameBase + ".pdf")
if 'eps' in self._fileFormats:
tcanvas.SaveAs(nameBase + ".eps")
if not imageOnly:
if 'root' in self._fileFormats:
tcanvas.SaveAs(nameBase + ".root")
if 'C' in self._fileFormats:
tcanvas.SaveAs(nameBase + ".C")
def _getColor(self, color):
if type(color) == int:
return color
elif type(color) == tuple:
# RGB
return ROOT.TColor.GetColor(*color)
elif type(color) == str:
# hex string
return ROOT.TColor.GetColor(color)
def _getLine(self, line):
if type(line) == int:
return line
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment