python/physicsutils/QGLikelihoodCalculator.py

0001 from __future__ import print_function
0002 from builtins import range
0003 import ROOT
0004 import math
0005
0006
0007
0008 def getBinNumber( bins, value ) :
0009
0010   if value<bins[0] or value>bins[len(bins)-1]: return False
0011   for i in range(0,len(bins)-1):
0012     if value>bins[i] and value<bins[i+1]:
0013       return i
0014
0015   return -1
0016
0017
0018
0019
0020 class QGLikelihoodCalculator:
0021
0022   def __init__(self, filename) :
0023     self.pdfs = {}
0024     self.etaBins = []
0025     self.ptBins  = []
0026     self.rhoBins = []
0027     self.varNames = { 0:'mult', 1:'ptD', 2:'axis2' }
0028     self.qgNames = { 0:'quark', 1:'gluon' }
0029     self.init(filename)
0030
0031   def init(self, filename) :
0032
0033     print("[QGLikelihoodCalculator]: Initializing from file: " + filename)
0034
0035     f = ROOT.TFile.Open(filename)
0036     if f.IsZombie() : return False
0037     try :
0038       self.etaBins = f.Get("etaBins")
0039       self.ptBins  = f.Get("ptBins")
0040       self.rhoBins = f.Get("rhoBins")
0041     except :
0042       return False
0043
0044
0045
0046     self.pdfs = {}
0047
0048     for it in range(0,len(self.qgNames)):
0049       self.pdfs[self.qgNames[it]] = {}
0050       for iv in range(0,len(self.varNames)):
0051         self.pdfs[self.qgNames[it]][self.varNames[iv]] = {}
0052         for ie in range(0,len(self.etaBins)):
0053           self.pdfs[self.qgNames[it]][self.varNames[iv]][ie] = {}
0054           for ip in range(0,len(self.ptBins)):
0055             self.pdfs[self.qgNames[it]][self.varNames[iv]][ie][ip] = {}
0056             for ir in range(0,len(self.rhoBins)):
0057               self.pdfs[self.qgNames[it]][self.varNames[iv]][ie][ip][ir] = 0
0058
0059
0060
0061     print("[QGLikelihoodCalculator]: Initialized binning of pdfs...")
0062
0063     keys = f.GetListOfKeys()
0064     for key in keys :
0065       if key.IsFolder() == False: continue
0066       hists = key.ReadObj().GetListOfKeys()
0067       for hist in hists :
0068         pieces = hist.GetName().split("_")
0069         varName = pieces[0]
0070         qgType  = pieces[1]
0071         etaStr  = pieces[2]
0072         ptStr   = pieces[3]
0073         rhoStr  = pieces[4]
0074         etaBin  = int(etaStr.split("eta")[1])
0075         ptBin   = int( ptStr.split("pt") [1])
0076         rhoBin  = int(rhoStr.split("rho")[1])
0077         histogram = hist.ReadObj()
0078         histogram.SetDirectory(0)
0079         self.pdfs[qgType][varName][etaBin][ptBin][rhoBin] = histogram
0080
0081     print("[QGLikelihoodCalculator]: pdfs initialized...")
0082
0083
0084     return True
0085
0086
0087
0088
0089   def isValidRange( self, pt, rho, eta ) :
0090
0091     if pt < self.ptBins[0]: return False
0092     if pt > self.ptBins[len(self.ptBins)-1]: return False
0093     if rho < self.rhoBins[0]: return False
0094     if rho > self.rhoBins[len(self.rhoBins)-1]: return False
0095     if math.fabs(eta) < self.etaBins[0]: return False
0096     if math.fabs(eta) > self.etaBins[len(self.etaBins)-1]: return False
0097     return True
0098
0099
0100
0101
0102   def findEntry( self, eta, pt, rho, qgType, varName ) :
0103
0104     etaBin = getBinNumber( self.etaBins, math.fabs(eta))
0105     if etaBin==-1: return None
0106     ptBin = getBinNumber( self.ptBins, pt )
0107     if ptBin==-1 : return None
0108
0109     rhoBin = getBinNumber( self.rhoBins, rho )
0110     if rhoBin==-1 : return None
0111
0112     return self.pdfs[qgType][varName][etaBin][ptBin][rhoBin]
0113
0114
0115
0116
0117
0118   def computeQGLikelihood( self, jet, rho ):
0119
0120     if self.isValidRange(jet.pt(), rho, jet.eta())==False:  return -1
0121
0122     # careful!!! this needs to be in the same order of self.varNames
0123     vars = {0:jet.mult, 1:jet.ptd, 2:jet.axis2}
0124
0125     Q=1.
0126     G=1.
0127
0128     #print "----------------------"
0129     #if jet.partonFlavour()==21 :
0130     #  print "this jet is a GLUON"
0131     #elif math.fabs(jet.partonFlavour())<4 and jet.partonFlavour()!=0:
0132     #  print "this jet is a QUARK"
0133     #print "pt: " + str(jet.pt()) + " eta: " + str(jet.eta()) + " rho: " + str(rho)
0134     #print "multi: " + str(jet.mult) + " ptd: " + str(jet.ptd) + " axis2: " + str(jet.axis2)
0135
0136     for i in vars :
0137
0138       #print self.varNames[i] + ": " + str(vars[i])
0139       # quarks
0140       qgEntry = self.findEntry(jet.eta(), jet.pt(), rho, "quark", self.varNames[i])
0141
0142       if qgEntry == None:  return -1
0143       Qi = qgEntry.GetBinContent(qgEntry.FindBin(vars[i]))
0144       mQ = qgEntry.GetMean()
0145       #print "Qi: " + str(Qi)
0146
0147       # gluons
0148       qgEntry = self.findEntry(jet.eta(), jet.pt(), rho, "gluon", self.varNames[i])
0149
0150       if qgEntry == None: return -1
0151       Gi = qgEntry.GetBinContent(qgEntry.FindBin(vars[i]))
0152       mG = qgEntry.GetMean()
0153       #print "Gi: " + str(Gi)
0154
0155       epsilon=0.
0156       delta=0.000001
0157       if Qi <= epsilon and Gi <= epsilon :
0158         if mQ>mG :
0159           if vars[i] > mQ :
0160              Qi = 1.-delta
0161              Gi = delta
0162           elif vars[i] < mG :
0163              Qi = delta
0164              Gi = 1.-delta
0165         else :
0166           if vars[i]<mQ :
0167              Qi = 1.-delta
0168              Gi = delta
0169           elif vars[i]>mG :
0170              Qi = delta
0171              Gi = 1.-delta
0172
0173       Q*=Qi
0174       G*=Gi
0175
0176
0177     #print "Q: " + str(Q)
0178     #print "G: " + str(G)
0179     if Q==0. : return 0.
0180     else : return Q/(Q+G)
0181