python/physicsutils/QGLikelihoodCalculator.py

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