Commit 65334088 authored by bothmann's avatar bothmann
Browse files

Add reference data and python libraries

parent da1ea1c4
__pycache__
This diff is collapsed.
from numpy import pi, log
NC = 3.
TR = 1./2.
CA = NC
CF = (NC*NC-1.)/(2.*NC)
class AlphaS:
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.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 beta1(self,nf):
return 17./6.*CA*CA-(5./3.*CA+CF)*TR*nf
def as0(self,t):
if t >= self.mb2:
tref = self.mz2
asref = self.asmz
b0 = self.beta0(5)/(2.*pi)
elif t >= self.mc2:
tref = self.mb2
asref = self.asmb
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))
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)
elif t >= self.mc2:
tref = self.mb2
asref = self.asmb
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)
def __call__(self,t):
if self.order == 0: return self.as0(t)
return self.as1(t)
import math as m
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 Cluster(self,event):
# TODO: implement the clustering algorithm here, and return a list of splitting scales Yij, ordered from smallest Yij to largest Yij
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.duralg = algorithm
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)
# fill integrated j-jet rates
previous_logy = 1e20
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
for p in s.points:
if p.x > logy and p.x < previous_logy:
p.y += w
previous_logy = logy
for p in self.ynm_integ[-1].points:
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 ]))
file.write("\n\n")
file.write("\n\n".join([ str(s) for s in self.ynm_integ ]))
file.close()
import sys
class Bin1D:
def __init__(self,xmin,xmax):
self.xmin = xmin
self.xmax = xmax
self.w = 0.
self.w2 = 0.
self.wx = 0.
self.wx2 = 0.
self.n = 0.
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}\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 Fill(self,x,w):
self.w += w
self.w2 += w*w
self.wx += w*x
self.wx2 += w*w*x
self.n += 1.
def ScaleW(self,scale):
self.w *= scale
self.w2 *= scale*scale
self.wx *= scale
self.wx2 *= scale*scale
class Point2D:
def __init__(self, xmin, xmax):
self.x = (xmax + xmin)/2
self.xerr = [(xmax - xmin)/2]*2
self.y = 0
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])
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"):
self.name = name
self.bins = []
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.scale = 1.
def __repr__(self):
return str(self)
def __str__(self):
s = "# BEGIN YODA_HISTO1D {0}\n".format(self.name)
s += "Path={0}\n".format(self.name)
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 += "# 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):
l = 0
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
a = self.bins[c].xmin
if x > self.bins[r].xmin:
self.oflow.Fill(x,w)
elif x < self.bins[l].xmin:
self.uflow.Fill(x,w)
else:
self.bins[l].Fill(x,w)
self.total.Fill(x,w)
def ScaleW(self,scale):
self.total.ScaleW(scale)
self.uflow.ScaleW(scale)
self.oflow.ScaleW(scale)
for i in range(len(self.bins)):
self.bins[i].ScaleW(scale)
self.scale = scale
def Plot(self):
import matplotlib.pyplot as plt
lefts = [b.xmin for b in self.bins]
lefts.append(self.bins[-1].xmax)
heights = [b.w / (lefts[1] - lefts[0]) for b in self.bins]
heights.append(heights[-1])
plt.step(lefts, heights, where='post')
class Scatter2D:
def __init__(self,npoints,xmin,xmax,name="/MC/untitled"):
self.name = name
self.points = []
width = (xmax-xmin)/npoints
for i in range(npoints):
self.points.append(Point2D(xmin+i*width,xmin+(i+1)*width))
self.scale = 1.
def __repr__(self):
return str(self)
def __str__(self):
s = "# BEGIN YODA_SCATTER2D {0}\n".format(self.name)
s += "Path={0}\n".format(self.name)
s += "Title=\nType=Histo1D\n"
s += "# xval\txerr-\txerr-\tyval\tyerr-\tyerr+\n"
for i in range(len(self.points)):
s += "{0}\n".format(self.points[i])
s += "# END YODA_SCATTER2D\n"
return s
def ScaleY(self,scale):
for i in range(len(self.points)):
self.points[i].ScaleY(scale)
self.scale = scale
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 __repr__(self):
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)
def Set(self,pdgid,momentum,col=[0,0]):
self.pid = pdgid
self.mom = momentum
self.col = col
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[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 \
len(csum) == 0)
import random
from numpy import pi, sqrt, tan, arctan, log, cos, sin
from utils.vector import Vec4
from utils.particle import Particle, CheckEvent
NC = 3.
TR = 1./2.
CA = NC
CF = (NC*NC-1.)/(2.*NC)
class Kernel:
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 Estimate(self,z):
return CF*2./(1.-z)
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())
class Pgg (Kernel):
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 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())
class Pgq (Kernel):
def Value(self,z,y):
return TR/2.*(1.-2.*z*(1.-z))
def Estimate(self,z):
return TR/2.
def Integral(self,zm,zp):
return TR/2.*(zp-zm)
def GenerateZ(self,zm,zp):
return zm+(zp-zm)*random.random()
class Shower:
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))
kt1 = pijt.Cross(pkt)
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()
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]
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] ]
else:
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] ]
else:
return [ [colij[0],self.c], [self.c,colij[1]] ]
else:
if flavs[1] > 0:
return [ [colij[0],0], [0,colij[1]] ]
else:
return [ [0,colij[1]], [colij[0],0] ]
def GeneratePoint(self,event):
while self.t > self.t0:
t = self.t0
for split in event[2:]:
for spect in event[2:]:
if spect == split: continue
if not split.ColorConnected(spect): continue
for sf in self.kernels:
if sf.flavs[0] != split.pid: continue
m2 = (split.mom+spect.mom).M2()
if m2 < 4.*self.t0: continue
zp = .5*(1.+sqrt(1.-4.*self.t0/m2))
g = self.alphamax/(2.*pi)*sf.Integral(1.-zp,zp)
tt = self.t*pow(random.random(),1./g)
if tt > t:
t = tt
s = [ split, spect, sf, m2, zp ]
self.t = t
if t > self.t0:
z = s[2].GenerateZ(1.-s[4],s[4])
y = t/s[3]/z/(1.-z)
if y < 1.:
f = (1.-y)*self.alpha(t)*s[2].Value(z,y)
g = self.alphamax*s[2].Estimate(z)
if f/g > random.random():
phi = 2.*pi*random.random()
moms = self.MakeKinematics(z,y,phi,s[0].mom,s[1].mom)
cols = self.MakeColors(s[2].flavs,s[0].col,s[1].col)
event.append(Particle(s[2].flavs[2],moms[1],cols[1]))
s[0].Set(s[2].flavs[1],moms[0],cols[0])
s[1].mom = moms[2]
return
def Run(self,event,t):
self.c = 1
self.t = t
while self.t > self.t0:
self.GeneratePoint(event)
import math as m
class Vec4:
def __init__(self,E=0,px=0,py=0,pz=0):
self.E = E
self.px = px
self.py = py
self.pz = pz
def __getitem__(self,i):
if i == 0: return self.E
if i == 1: return self.px
if i == 2: return self.py
if i == 3: return self.pz
raise Exception('Vec4D')
def __repr__(self):
return '({0},{1},{2},{3})'.format(self.E,self.px,self.py,self.pz)
def __str__(self):
return '({0},{1},{2},{3})'.format(self.E,self.px,self.py,self.pz)
def __add__(self,v):
return Vec4(self.E+v.E,self.px+v.px,self.py+v.py,self.pz+v.pz)
def __sub__(self,v):
return Vec4(self.E-v.E,self.px-v.px,self.py-v.py,self.pz-v.pz)
def __neg__(self):
return Vec4(-self.E,-self.px,-self.py,-self.pz)
def __mul__(self,v):
if isinstance(v,Vec4):
return self.E*v.E-self.px*v.px-self.py*v.py-self.pz*v.pz
return Vec4(self.E*v,self.px*v,self.py*v,self.pz*v)
def __rmul__(self,v):
if isinstance(v,Vec4):
return self.E*v.E-self.px*v.px-self.py*v.py-self.pz*v.pz
return Vec4(self.E*v,self.px*v,self.py*v,self.pz*v)
def __div__(self,v):
return Vec4(self.E/v,self.px/v,self.py/v,self.pz/v)
def M2(self):
return self*self
def M(self):
return m.sqrt(self.M2())
def P2(self):
return self.px*self.px+self.py*self.py+self.pz*self.pz
def P(self):
return m.sqrt(self.P2())
def PT2(self):
return self.px*self.px+self.py*self.py
def PT(self):
return m.sqrt(self.PT2())
def Theta(self) :
return m.acos(self.pz/self.P())
def Phi(self) :
if self.px==0 and self.py==0:
return 0.0
else:
return m.atan2(self.py,self.px)
def Cross(self,v):
return Vec4(0.0,
self.py*v.pz-self.pz*v.py,
self.pz*v.px-self.px*v.pz,
self.px*v.py-self.py*v.px)
def Boost(self,v):
rsq = self.M()
v0 = (self.E*v.E-self.px*v.px-self.py*v.py-self.pz*v.pz)/rsq;
c1 = (v.E+v0)/(rsq+self.E)
return Vec4(v0,
v.px-c1*self.px,
v.py-c1*self.py,
v.pz-c1*self.pz)
def BoostBack(self,v):
rsq = self.M()
v0 = (self.E*v.E+self.px*v.px+self.py*v.py+self.pz*v.pz)/rsq;
c1 = (v.E+v0)/(rsq+self.E)
return Vec4(v0,
v.px+c1*self.px,
v.py+c1*self.py,
v.pz+c1*self.pz)
#!/usr/bin/env python3
import matplotlib.pyplot as plt
import numpy as np
lw = 1.0
colors = ['red', 'blue']
single_plot = False
def data_objects(filenames, yodatype="HISTO1D"):
names = set()
objects = {}
for filename in filenames:
with open(filename) as file:
objects[filename] = {}
current = None
for line in file:
words = line.split()
if len(words) == 0:
continue
if len(words) > 3 and words[2] == "YODA_" + yodatype:
objects[filename][words[3]] = {"scale_factor": 1.0, "bins": []}
current = objects[filename][words[3]]
names.add(words[3])
continue
if len(words) > 1 and "END" == words[1]:
current = None
continue