Commit cdecd505 authored by janlukas.bosse's avatar janlukas.bosse
Browse files

ran autopep8 to get somewhat more acceptable code

parent d70da821
from numpy import pi, log
NC = 3.
TR = 1./2.
TR = 1. / 2.
CA = NC
CF = (NC*NC-1.)/(2.*NC)
CF = (NC * NC - 1.) / (2. * NC)
class AlphaS:
def __init__(self,mz,asmz,order=1,mb=4.75,mc=1.3):
def __init__(self, mz, asmz, order=1, mb=4.75, mc=1.3):
self.order = order
self.mc2 = mc*mc
self.mb2 = mb*mb
self.mz2 = mz*mz
self.mc2 = mc * mc
self.mb2 = mb * mb
self.mz2 = mz * mz
self.asmz = asmz
self.asmb = self(self.mb2)
self.asmc = self(self.mc2)
#print("\\alpha_s({0}) = {1}".format(mz,self(self.mz2)))
def beta0(self,nf):
return 11./6.*CA-2./3.*TR*nf
def beta0(self, nf):
return 11. / 6. * CA - 2. / 3. * TR * nf
def beta1(self,nf):
return 17./6.*CA*CA-(5./3.*CA+CF)*TR*nf
def beta1(self, nf):
return 17. / 6. * CA * CA - (5. / 3. * CA + CF) * TR * nf
def as0(self,t):
def as0(self, t):
if t >= self.mb2:
tref = self.mz2
asref = self.asmz
b0 = self.beta0(5)/(2.*pi)
b0 = self.beta0(5) / (2. * pi)
elif t >= self.mc2:
tref = self.mb2
asref = self.asmb
b0 = self.beta0(4)/(2.*pi)
b0 = self.beta0(4) / (2. * pi)
else:
tref = self.mc2
asref = self.asmc
b0 = self.beta0(3)/(2.*pi)
return 1./(1./asref+b0*log(t/tref))
b0 = self.beta0(3) / (2. * pi)
return 1. / (1. / asref + b0 * log(t / tref))
def as1(self,t):
def as1(self, t):
if t >= self.mb2:
tref = self.mz2
asref = self.asmz
b0 = self.beta0(5)/(2.*pi)
b1 = self.beta1(5)/pow(2.*pi,2)
b0 = self.beta0(5) / (2. * pi)
b1 = self.beta1(5) / pow(2. * pi, 2)
elif t >= self.mc2:
tref = self.mb2
asref = self.asmb
b0 = self.beta0(4)/(2.*pi)
b1 = self.beta1(4)/pow(2.*pi,2)
b0 = self.beta0(4) / (2. * pi)
b1 = self.beta1(4) / pow(2. * pi, 2)
else:
tref = self.mc2
asref = self.asmc
b0 = self.beta0(3)/(2.*pi)
b1 = self.beta1(3)/pow(2.*pi,2)
w = 1.+b0*asref*log(t/tref)
return asref/w*(1.-b1/b0*asref*log(w)/w)
b0 = self.beta0(3) / (2. * pi)
b1 = self.beta1(3) / pow(2. * pi, 2)
w = 1. + b0 * asref * log(t / tref)
return asref / w * (1. - b1 / b0 * asref * log(w) / w)
def __call__(self,t):
if self.order == 0: return self.as0(t)
def __call__(self, t):
if self.order == 0:
return self.as0(t)
return self.as1(t)
......@@ -5,13 +5,15 @@ from utils.vector import Vec4
from utils.particle import Particle
from utils.histogram import Histo1D, Scatter2D
class Algorithm:
def Yij(self,p,q):
pq = p.px*q.px+p.py*q.py+p.pz*q.pz
return 2.*pow(min(p.E,q.E),2)*(1.0-min(max(pq/m.sqrt(p.P2()*q.P2()),-1.0),1.0))/self.ecm2
def Yij(self, p, q):
pq = p.px * q.px + p.py * q.py + p.pz * q.pz
return 2. * pow(min(p.E, q.E), 2) * (1.0 - min(max(pq /
m.sqrt(p.P2() * q.P2()), -1.0), 1.0)) / self.ecm2
def Cluster(self,event):
def Cluster(self, event):
# TODO: implement the clustering algorithm here, and return a list of
# splitting scales Yij, ordered from smallest Yij to largest Yij
......@@ -24,28 +26,30 @@ class Algorithm:
return []
class Analysis:
def __init__(self, algorithm):
self.n = 0.
self.ynm = [ Histo1D(100,-4.3,-0.3,'/LL_JetRates/log10_y_{0}{1}'.format(i+2,i+3))
for i in range(4) ]
self.ynm_integ = [ Scatter2D(100,-4.3,-0.3,'/LL_JetRates/integ_log10_y_{0}'.format(i+2)) for i in range(5) ]
self.ynm = [Histo1D(100, -4.3, -0.3, '/LL_JetRates/log10_y_{0}{1}'.format(i + 2, i + 3))
for i in range(4)]
self.ynm_integ = [Scatter2D(
100, -4.3, -0.3, '/LL_JetRates/integ_log10_y_{0}'.format(i + 2)) for i in range(5)]
self.duralg = algorithm
def Analyze(self,event,w):
def Analyze(self, event, w):
self.n += 1.
# fill differential j -> (j+1) splitting scale distributions
kt2 = self.duralg.Cluster(event)
for j in range(len(self.ynm)):
self.ynm[j].Fill(m.log10(kt2[-1-j]) if len(kt2)>j else -5.,w)
self.ynm[j].Fill(m.log10(kt2[-1 - j]) if len(kt2) > j else -5., w)
# fill integrated j-jet rates
previous_logy = 1e20
for j in range(len(self.ynm_integ)-1):
for j in range(len(self.ynm_integ) - 1):
s = self.ynm_integ[j]
logy = m.log10(kt2[-1-j]) if len(kt2)>j else -5.0
logy = m.log10(kt2[-1 - j]) if len(kt2) > j else -5.0
for p in s.points:
if p.x > logy and p.x < previous_logy:
p.y += w
......@@ -54,11 +58,13 @@ class Analysis:
if p.x < previous_logy:
p.y += w
def Finalize(self,name):
for h in self.ynm: h.ScaleW(1./self.n)
for s in self.ynm_integ: s.ScaleY(1./self.n)
file = open(name+".yoda","w")
file.write("\n\n".join([ str(h) for h in self.ynm ]))
def Finalize(self, name):
for h in self.ynm:
h.ScaleW(1. / self.n)
for s in self.ynm_integ:
s.ScaleY(1. / self.n)
file = open(name + ".yoda", "w")
file.write("\n\n".join([str(h) for h in self.ynm]))
file.write("\n\n")
file.write("\n\n".join([ str(s) for s in self.ynm_integ ]))
file.write("\n\n".join([str(s) for s in self.ynm_integ]))
file.close()
import sys
class Bin1D:
def __init__(self,xmin,xmax):
def __init__(self, xmin, xmax):
self.xmin = xmin
self.xmax = xmax
self.w = 0.
......@@ -16,57 +17,58 @@ class Bin1D:
return str(self)
def __str__(self):
return "{0:10.6e}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5:10.6e}\t{6}".format\
(self.xmin,self.xmax,self.w,self.w2,self.wx,self.wx2,int(self.n))
return "{0:10.6e}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5:10.6e}\t{6}".format(
self.xmin, self.xmax, self.w, self.w2, self.wx, self.wx2, int(self.n))
def Format(self,tag):
return "{0}\t{0}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5}".format\
(tag,self.w,self.w2,self.wx,self.wx2,int(self.n))
def Format(self, tag):
return "{0}\t{0}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5}".format(
tag, self.w, self.w2, self.wx, self.wx2, int(self.n))
def Fill(self,x,w):
def Fill(self, x, w):
self.w += w
self.w2 += w*w
self.wx += w*x
self.wx2 += w*w*x
self.w2 += w * w
self.wx += w * x
self.wx2 += w * w * x
self.n += 1.
def ScaleW(self,scale):
def ScaleW(self, scale):
self.w *= scale
self.w2 *= scale*scale
self.w2 *= scale * scale
self.wx *= scale
self.wx2 *= scale*scale
self.wx2 *= scale * scale
class Point2D:
def __init__(self, xmin, xmax):
self.x = (xmax + xmin)/2
self.xerr = [(xmax - xmin)/2]*2
self.x = (xmax + xmin) / 2
self.xerr = [(xmax - xmin) / 2] * 2
self.y = 0
self.yerr = [0]*2
self.yerr = [0] * 2
def __repr__(self):
return str(self)
def __str__(self):
return "{0:10.6e}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5:10.6e}".format\
(self.x,self.xerr[0],self.xerr[1],self.y,self.yerr[0],self.yerr[1])
return "{0:10.6e}\t{1:10.6e}\t{2:10.6e}\t{3:10.6e}\t{4:10.6e}\t{5:10.6e}".format(
self.x, self.xerr[0], self.xerr[1], self.y, self.yerr[0], self.yerr[1])
def ScaleY(self, scale):
self.y *= scale
self.yerr = [err * scale for err in self.yerr]
class Histo1D:
def __init__(self,nbin,xmin,xmax,name="/MC/untitled"):
def __init__(self, nbin, xmin, xmax, name="/MC/untitled"):
self.name = name
self.bins = []
self.uflow = Bin1D(-sys.float_info.max,xmin)
width = (xmax-xmin)/nbin
self.uflow = Bin1D(-sys.float_info.max, xmin)
width = (xmax - xmin) / nbin
for i in range(nbin):
self.bins.append(Bin1D(xmin+i*width,xmin+(i+1)*width))
self.oflow = Bin1D(xmax,sys.float_info.max)
self.total = Bin1D(-sys.float_info.max,sys.float_info.max)
self.bins.append(Bin1D(xmin + i * width, xmin + (i + 1) * width))
self.oflow = Bin1D(xmax, sys.float_info.max)
self.total = Bin1D(-sys.float_info.max, sys.float_info.max)
self.scale = 1.
def __repr__(self):
......@@ -78,34 +80,36 @@ class Histo1D:
s += "ScaledBy={0}\n".format(self.scale)
s += "Title=\nType=Histo1D\n"
s += "# ID\tID\tsumw\tsumw2\tsumwx\tsumwx2\tnumEntries\n"
s += self.total.Format("Total")+"\n"
s += self.uflow.Format("Underflow")+"\n"
s += self.oflow.Format("Overflow")+"\n"
s += self.total.Format("Total") + "\n"
s += self.uflow.Format("Underflow") + "\n"
s += self.oflow.Format("Overflow") + "\n"
s += "# xlow\txhigh\tsumw\tsumw2\tsumwx\tsumwx2\tnumEntries\n"
for i in range(len(self.bins)):
s += "{0}\n".format(self.bins[i])
s += "# END YODA_HISTO1D\n"
return s
def Fill(self,x,w):
def Fill(self, x, w):
l = 0
r = len(self.bins)-1
c = (l+r)//2
r = len(self.bins) - 1
c = (l + r) // 2
a = self.bins[c].xmin
while r-l > 1:
if x < a: r = c
else: l = c
c = (l+r)//2
while r - l > 1:
if x < a:
r = c
else:
l = c
c = (l + r) // 2
a = self.bins[c].xmin
if x > self.bins[r].xmin:
self.oflow.Fill(x,w)
self.oflow.Fill(x, w)
elif x < self.bins[l].xmin:
self.uflow.Fill(x,w)
self.uflow.Fill(x, w)
else:
self.bins[l].Fill(x,w)
self.total.Fill(x,w)
self.bins[l].Fill(x, w)
self.total.Fill(x, w)
def ScaleW(self,scale):
def ScaleW(self, scale):
self.total.ScaleW(scale)
self.uflow.ScaleW(scale)
self.oflow.ScaleW(scale)
......@@ -121,14 +125,16 @@ class Histo1D:
heights.append(heights[-1])
plt.step(lefts, heights, where='post')
class Scatter2D:
def __init__(self,npoints,xmin,xmax,name="/MC/untitled"):
def __init__(self, npoints, xmin, xmax, name="/MC/untitled"):
self.name = name
self.points = []
width = (xmax-xmin)/npoints
width = (xmax - xmin) / npoints
for i in range(npoints):
self.points.append(Point2D(xmin+i*width,xmin+(i+1)*width))
self.points.append(
Point2D(xmin + i * width, xmin + (i + 1) * width))
self.scale = 1.
def __repr__(self):
......@@ -144,7 +150,7 @@ class Scatter2D:
s += "# END YODA_SCATTER2D\n"
return s
def ScaleY(self,scale):
def ScaleY(self, scale):
for i in range(len(self.points)):
self.points[i].ScaleY(scale)
self.scale = scale
......@@ -3,39 +3,43 @@ import math as m
from utils.vector import Vec4
class Particle:
def __init__(self,pdgid,momentum,col=[0,0]):
self.Set(pdgid,momentum,col)
def __init__(self, pdgid, momentum, col=[0, 0]):
self.Set(pdgid, momentum, col)
def __repr__(self):
return "{0} {1} {2}".format(self.pid,self.mom,self.col)
return "{0} {1} {2}".format(self.pid, self.mom, self.col)
def __str__(self):
return "{0} {1} {2}".format(self.pid,self.mom,self.col)
return "{0} {1} {2}".format(self.pid, self.mom, self.col)
def Set(self,pdgid,momentum,col=[0,0]):
def Set(self, pdgid, momentum, col=[0, 0]):
self.pid = pdgid
self.mom = momentum
self.col = col
def ColorConnected(self,p):
def ColorConnected(self, p):
return (self.col[0] > 0 and self.col[0] == p.col[1]) or \
(self.col[1] > 0 and self.col[1] == p.col[0])
def CheckEvent(event):
psum = Vec4()
csum = {}
for p in event:
psum += p.mom
if p.col[0] > 0:
csum[p.col[0]] = csum.get(p.col[0],0) + 1
if csum[p.col[0]] == 0: del csum[p.col[0]]
if p.col[0] > 0:
csum[p.col[0]] = csum.get(p.col[0], 0) + 1
if csum[p.col[0]] == 0:
del csum[p.col[0]]
if p.col[1] > 0:
csum[p.col[1]] = csum.get(p.col[1],0) - 1
if csum[p.col[1]] == 0: del csum[p.col[1]]
return (m.fabs(psum.E)<1.e-12 and \
m.fabs(psum.px)<1.e-12 and \
m.fabs(psum.py)<1.e-12 and \
m.fabs(psum.pz)<1.e-12 and \
csum[p.col[1]] = csum.get(p.col[1], 0) - 1
if csum[p.col[1]] == 0:
del csum[p.col[1]]
return (m.fabs(psum.E) < 1.e-12 and
m.fabs(psum.px) < 1.e-12 and
m.fabs(psum.py) < 1.e-12 and
m.fabs(psum.pz) < 1.e-12 and
len(csum) == 0)
......@@ -5,136 +5,149 @@ from utils.vector import Vec4
from utils.particle import Particle, CheckEvent
NC = 3.
TR = 1./2.
TR = 1. / 2.
CA = NC
CF = (NC*NC-1.)/(2.*NC)
CF = (NC * NC - 1.) / (2. * NC)
class Kernel:
def __init__(self,flavs):
def __init__(self, flavs):
self.flavs = flavs
class Pqq (Kernel):
def Value(self,z,y):
return CF*(2./(1.-z*(1.-y))-(1.+z))
def Value(self, z, y):
return CF * (2. / (1. - z * (1. - y)) - (1. + z))
def Estimate(self,z):
return CF*2./(1.-z)
def Estimate(self, z):
return CF * 2. / (1. - z)
def Integral(self,zm,zp):
return CF*2.*log((1.-zm)/(1.-zp))
def Integral(self, zm, zp):
return CF * 2. * log((1. - zm) / (1. - zp))
def GenerateZ(self, zm, zp):
return 1. + (zp - 1.) * pow((1. - zm) / (1. - zp), random.random())
def GenerateZ(self,zm,zp):
return 1.+(zp-1.)*pow((1.-zm)/(1.-zp),random.random())
class Pgg (Kernel):
def Value(self,z,y):
return CA/2.*(2./(1.-z*(1.-y))-2.+z*(1.-z))
def Value(self, z, y):
return CA / 2. * (2. / (1. - z * (1. - y)) - 2. + z * (1. - z))
def Estimate(self, z):
return CA / (1. - z)
def Estimate(self,z):
return CA/(1.-z)
def Integral(self, zm, zp):
return CA * log((1. - zm) / (1. - zp))
def Integral(self,zm,zp):
return CA*log((1.-zm)/(1.-zp))
def GenerateZ(self, zm, zp):
return 1. + (zp - 1.) * pow((1. - zm) / (1. - zp), random.random())
def GenerateZ(self,zm,zp):
return 1.+(zp-1.)*pow((1.-zm)/(1.-zp),random.random())
class Pgq (Kernel):
def Value(self,z,y):
return TR/2.*(1.-2.*z*(1.-z))
def Value(self, z, y):
return TR / 2. * (1. - 2. * z * (1. - z))
def Estimate(self,z):
return TR/2.
def Estimate(self, z):
return TR / 2.
def Integral(self,zm,zp):
return TR/2.*(zp-zm)
def Integral(self, zm, zp):
return TR / 2. * (zp - zm)
def GenerateZ(self, zm, zp):
return zm + (zp - zm) * random.random()
def GenerateZ(self,zm,zp):
return zm+(zp-zm)*random.random()
class Shower:
def __init__(self,alpha,t0=1.0):
def __init__(self, alpha, t0=1.0):
self.t0 = t0
self.alpha = alpha
self.alphamax = alpha(self.t0)
self.kernels = [ Pqq([fl,fl,21]) for fl in [-5,-4,-3,-2,-1,1,2,3,4,5] ]
self.kernels += [ Pgq([21,fl,-fl]) for fl in [1,2,3,4,5] ]
self.kernels += [ Pgg([21,21,21]) ]
def MakeKinematics(self,z,y,phi,pijt,pkt):
Q = pijt+pkt
rkt = sqrt(Q.M2()*y*z*(1.-z))
self.kernels = [Pqq([fl, fl, 21])
for fl in [-5, -4, -3, -2, -1, 1, 2, 3, 4, 5]]
self.kernels += [Pgq([21, fl, -fl]) for fl in [1, 2, 3, 4, 5]]
self.kernels += [Pgg([21, 21, 21])]
def MakeKinematics(self, z, y, phi, pijt, pkt):
Q = pijt + pkt
rkt = sqrt(Q.M2() * y * z * (1. - z))
kt1 = pijt.Cross(pkt)
if kt1.P() < 1.e-6: kt1 = pijt.Cross(Vec4(0.,1.,0.,0.))
kt1 *= rkt*cos(phi)/kt1.P()
if kt1.P() < 1.e-6:
kt1 = pijt.Cross(Vec4(0., 1., 0., 0.))
kt1 *= rkt * cos(phi) / kt1.P()
kt2cms = Q.Boost(pijt).Cross(kt1)
kt2cms *= rkt*sin(phi)/kt2cms.P()
kt2cms *= rkt * sin(phi) / kt2cms.P()
kt2 = Q.BoostBack(kt2cms)
pi = z*pijt + (1.-z)*y*pkt + kt1 + kt2
pj = (1.-z)*pijt + z*y*pkt - kt1 - kt2
pk = (1.-y)*pkt
return [pi,pj,pk]
pi = z * pijt + (1. - z) * y * pkt + kt1 + kt2
pj = (1. - z) * pijt + z * y * pkt - kt1 - kt2
pk = (1. - y) * pkt
return [pi, pj, pk]
def MakeColors(self,flavs,colij,colk):
def MakeColors(self, flavs, colij, colk):
self.c += 1
if flavs[0] != 21:
if flavs[0] > 0:
return [ [self.c,0], [colij[0],self.c] ]
return [[self.c, 0], [colij[0], self.c]]
else:
return [ [0,self.c], [self.c,colij[1]] ]
return [[0, self.c], [self.c, colij[1]]]
else:
if flavs[1] == 21:
if colij[0] == colk[1]:
if colij[1] == colk[0] and random.random()>0.5:
return [ [colij[0],self.c], [self.c,colij[1]] ]
return [ [self.c,colij[1]], [colij[0],self.c] ]
if colij[1] == colk[0] and random.random() > 0.5:
return [[colij[0], self.c], [self.c, colij[1]]]
return [[self.c, colij[1]], [colij[0], self.c]]
else:
return [ [colij[0],self.c], [self.c,colij[1]] ]
return [[colij[0], self.c], [self.c, colij[1]]]
else:
if flavs[1] > 0:
return [ [colij[0],0], [0,colij[1]] ]
return [[colij[0], 0], [0, colij[1]]]
else:
return [ [0,colij[1]], [colij[0],0] ]
return [[0, colij[1]], [colij[0], 0]]
def GeneratePoint(self,event):