commit 9747ccd9f09d9053e211124fb57edf50268ef415
Author: Bakar Chargeishvili <bakar@bcharge.de>
Date: Tue, 29 Mar 2022 17:44:53 +0200
Import the project
Diffstat:
6 files changed, 448 insertions(+), 0 deletions(-)
diff --git a/LICENSE b/LICENSE
@@ -0,0 +1,29 @@
+BSD 3-Clause License
+
+Copyright (c) 2021, Bakar Chargeishvili
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+ list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+ this list of conditions and the following disclaimer in the documentation
+ and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its
+ contributors may be used to endorse or promote products derived from
+ this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/README.md b/README.md
@@ -0,0 +1,10 @@
+# pyDipole - python module to automatize dipole subtraction term calculations
+
+Instructions will follow...
+
+
+## References:
+
+- Massless - [hep-ph/0201036](https://arxiv.org/abs/hep-ph/9605323)
+- Massive - [hep-ph/9605323](https://arxiv.org/abs/hep-ph/0201036)
+- AutoDipole - [0911.4371 [hep-ph]](https://arxiv.org/abs/0911.4371)
diff --git a/calculate.py b/calculate.py
@@ -0,0 +1,314 @@
+import subprocess
+import logging
+import sympy as sy
+import re
+
+from pyDipole import Process
+
+logging.getLogger().setLevel(logging.INFO)
+
+
+proc = Process(['g','g'],['t','tbar'],['t'])
+
+expressions = ''
+par1_rules = ''
+par1_defs = 'S '
+nu_rules = ''
+gamma6_rules = ''
+
+
+def eq_c27(id1, id2):
+ global expressions
+ expressions += '''\
+L I{0}{1} = -Alfas/(2*pi)*InvGamma(1-eps)*Denom(Col({0},{0}))*Nu0({0})*Col({0},{1})*(4*pi*mu^2*s({0},{1})^(-1))^eps;\n
+'''.format(id1,id2)
+
+def eq_616(id1,id2):
+ global expressions
+ expressions += '''\
+L I{0}{1} = -Alfas*(4*pi)**eps/(2*pi)*InvGamma(1-eps)*Denom(Col({0},{0}))*Col({0},{1})*(Col({0},{0})*(mu^2*s({0},{1})^(-1))^eps*
+(Nu616({0},{1})-pi^2/3)+Gamma616({0})+Gamm({0})*log(mu^2*s({0},{1})^(-1))+Gamm({0})+K({0}));
+'''.format(id1,id2)
+
+def eq_652(id1,id2):
+ global expressions
+ expressions += '''\
+L I{0}{1} = -Alfas*(4*pi)**eps/(2*pi)*InvGamma(1-eps)*Denom(Col({0},{0}))*Col({0},{1})*(Col({0},{0})*(mu^2*s({0},{1})^(-1))^eps*
+(Nu616({0},{1})-pi^2/3)+Gamma616({0})+Gamm({0})*log(mu^2*s({0},{1})^(-1))+Gamm({0})+K({0}));
+'''.format(id1,id2)
+
+for par1 in proc.all_particles:
+ if not par1.isQCD:
+ continue
+
+ par1_defs += 'm{0}, p{0}, '.format(par1.id)
+
+ if par1.name == 'g':
+ par1_rules += '''\
+id Col({0},{0}) = CA;
+id Gamm({0}) = CA*11/6-Tr*Nfl*2/3;
+id K({0}) = (67/18-pi^2/6)*CA-Tr*Nfl*10/9;\n
+'''.format(par1.id)
+ else:
+ par1_rules += '''\
+id Col({0},{0}) = CF;
+id Gamm({0}) = CF*3/2;
+id K({0}) = (7/2-pi^2/6)*CF;
+'''.format(par1.id)
+
+ for par2 in proc.all_particles:
+ if par1.id == par2.id:
+ continue
+ if not par2.isQCD:
+ continue
+
+ logging.info('Evaluating I_{0}{1}...'.format(par1.id,par2.id))
+
+ #final fermion
+ if (par1.name != 'g' and not par1.initial):
+ if not par2.initial:
+ if not par1.isMassive and not par2.isMassive:
+ #C.27
+ logging.info('(i,k)=(f,k), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ else:
+ #6.16
+ logging.info('(i,k)=(f,k), massive, using eq 6.16')
+ eq_616(par1.id, par2.id)
+ else:
+ if not par1.isMassive:
+ #C.27
+ logging.info('(i,k)=(f,b), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ else:
+ #6.52
+ logging.info('(i,k)=(f,b), massive, using eq 6.52')
+ eq_652(par1.id, par2.id)
+ #final gluon
+ elif (par1.name == 'g' and not par1.initial):
+ if not par2.initial:
+ if not par1.isMassive and not proc.mF_list:
+ #C.27
+ logging.info('(i,k)=(g,k), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ else:
+ #6.16
+ logging.info('(i,k)=(g,k), massive, using eq 6.16')
+ eq_616(par1.id, par2.id)
+ else:
+ if not proc.mF_list:
+ #C.27
+ logging.info('(i,k)=(g,b), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ else:
+ #6.52
+ logging.info('(i,k)=(g,b), massive, using eq 6.52')
+ eq_652(par1.id, par2.id)
+ #initial fermion
+ elif (par1.name != 'g' and par1.initial):
+ if not par2.initial:
+ if not par2.isMassive:
+ #C.27
+ logging.info('(a,k)=(f,k), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ else:
+ #6.52
+ logging.info('(a,k)=(f,k), massive, using eq 6.52')
+ eq_652(par1.id, par2.id)
+ else:
+ #C.27
+ logging.info('(a,b)=(f,b), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ #initial gluon
+ elif (par1.name == 'g' and par1.initial):
+ if not par2.initial:
+ if not par2.isMassive:
+ #C.27
+ logging.info('(a,k)=(g,k), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+ else:
+ #6.52
+ logging.info('(a,k)=(g,k), massive, using eq 6.52')
+ eq_652(par1.id, par2.id)
+ else:
+ #C.27
+ logging.info('(a,b)=(g,b), massless, using eq C.27')
+ eq_c27(par1.id, par2.id)
+
+ #Nu replacement rules
+ if (par1.isMassive and par2.isMassive):
+ nu_rules += '''\
+id Nu616S({0},{1}) = Denom(v({0},{1}))*(eps^(-1)*log(rho({0},{1}))-
+1/4*log(rhon({0},{1},{0})^2)^2-1/4*log(rhon({0},{1},{1})^2)^2-pi^2/6+
+log(rho({0},{1}))*log(Qij({0},{1})^2/s({0},{1})));
+'''.format(par1.id, par2.id)
+ elif (par1.isMassive and not par2.isMassive):
+ nu_rules += '''\
+id Nu616S({0},{1}) = eps^-2/2+eps^(-1)/2*log(m({0})^2/s({0},{1}))-
+1/4*log(m({0})^2/s({0},{1}))-pi^2/12-1/2*log(m({0})^2/s({0},{1}))*log(s({0},{1})*Qij({0},{1})^-2)-
+1/2*log(m({0})^2*Qij({0},{1})^-2)*log(s({0},{1})*Qij({0},{1})^-2);
+'''.format(par1.id, par2.id)
+ elif (not par1.isMassive and par2.isMassive):
+ nu_rules += '''\
+id Nu616S({1},{0}) = eps^-2/2+eps^(-1)/2*log(m({0})^2/s({0},{1}))-
+1/4*log(m({0})^2/s({0},{1}))-pi^2/12-1/2*log(m({0})^2/s({0},{1}))*log(s({0},{1})*Qij({0},{1})^-2)-
+1/2*log(m({0})^2*Qij({0},{1})^-2)*log(s({0},{1})*Qij({0},{1})^-2);
+'''.format(par2.id, par1.id)
+ elif (not par1.isMassive and not par2.isMassive):
+ nu_rules += '''\
+id Nu616S({0},{1}) = eps^-2;
+'''.format(par1.id, par2.id)
+
+ #gamma616 replacement rules (singular part only)
+ if not par1.isMassive:
+ gamma6_rules += 'id Gamma616({0}) = eps^(-1)*Gamm({0});\n'.format(par1.id)
+ else:
+ gamma6_rules += 'id Gamma616({0}) = eps^(-1)*CF;\n'.format(par1.id)
+
+print(expressions)
+
+
+code = '''\
+#-
+off stats;
+
+S Alfas, pi, eps, CF, CA, Tr, Nfl, mu;
+CF Denom;
+CF I,InvGamma, Col, Nu0, Nu616 Gamm, K, Gamma616;
+CF Nu616S, Nu616Ns;
+CF Mu0, Mu616, log;
+CF Qij, s, rho, rhon, v, m;
+CF NuSmm, NuSmz, NuSzm;
+CF NuNSmm, NuNSmz, NuNSzm;
+CF sqrt, mun;
+AutoDeclare I i;
+'''
+
+par1_defs = par1_defs[:-2]+';\n\n'
+code += par1_defs
+
+code += expressions+'\n\n'
+
+code += '''\
+id Nu0(i1?) = Col(i1,i1)*(eps^(-2)-pi^2/3)+Gamm(i1)*eps^(-1)+Gamm(i1)+K(i1);
+id Nu616(?k) = Nu616Ns(?k) + Nu616S(?k);
+\n\n
+'''
+
+code += nu_rules
+
+rep_rules = '''\
+id rho(i1?,i2?) = sqrt((1-v(i1,i2))*Denom(1+v(i1,i2)));
+id rhon(i1?,i2?,i3?) = sqrt((1-v(i1,i2)+2*mun(i1,i2,i3)^2*Denom(1-mun(i1,i2,i1)^2-mun(i1,i2,i2)^2))*
+Denom(1+v(i1,i2)+2*mun(i1,i2,i3)^2*Denom(1-mun(i1,i2,i1)^2-mun(i1,i2,i2)^2)));
+
+id v(i1?,i2?) = (s(i1,i2))^(-1)*sqrt((s(i1,i2))^(2)-(4*m(i1)^2*m(i2)^2) );
+id mun(i1?,i2?,i3?) = m(i3)*Denom(Qij(i1,i2));
+id Qij(i1?,i2?) = sqrt(s(i1,i2) + m(i1)^2 + m(i2)^2);
+
+'''
+
+code += '''\
+argument log;
+{0}
+endargument;
+
+argument Denom;
+{0}
+endargument;
+
+argument sqrt;
+{0}
+endargument;
+
+argument log;
+argument sqrt;
+{0}
+endargument;
+endargument;
+
+argument log;
+argument Denom;
+{0}
+endargument;
+endargument;
+
+argument sqrt;
+argument Denom;
+{0}
+endargument;
+endargument;
+
+argument log;
+argument sqrt;
+argument Denom;
+{0}
+endargument;
+endargument;
+endargument;
+
+{0}
+\n\n
+'''.format(rep_rules)
+
+code += gamma6_rules+'\n\n'
+
+code += par1_rules
+
+code += '''\
+\nargument Denom;
+{0}
+endargument;
+id Denom(?k) = exp_(?k,-1);
+.sort
+Print;
+B InvGamma,Alfas,pi,Mu0;
+.end
+'''.format(par1_rules)
+
+form_file_name = '/tmp/FORM_Code.frm'
+with open(form_file_name,'w') as f:
+ f.write(code)
+
+"""
+Run FORM, read the result back in python and execute SymPy simplifications
+"""
+
+form_res = subprocess.run(["form","-f",form_file_name],capture_output=True,text=True).stdout
+form_res=form_res.replace('\n','')
+form_res=form_res.replace(' ','')
+form_res = form_res.split(';')[:-1]
+
+#Prepare SymPy variables
+InvGamma = sy.Function('InvGamma')
+MuFactor = sy.Function('MuFactor')
+Col = sy.Function('Col')
+Denom = sy.Function('Denom')
+m = sy.Function('m')
+eps = sy.Symbol('eps')
+a,b = map(sy.Wild, 'ab')
+alphas = sy.Symbol('alphas')
+ep,pi,mu = sy.symbols('ep pi mu')
+mom = dict()
+masses = dict()
+for p in proc.all_particles:
+ mom[p.id]=sy.Symbol('p{0}'.format(p.id))
+ if p.isQCD and p.isMassive:
+ masses[p.id] = sy.Symbol(p.mass)
+
+
+eq_dict = dict()
+for eq in form_res:
+ eq_name = re.search(r"^(.*)=",eq).group(0)[:-1]
+ eq_dict[eq_name] = []
+
+ eq = sy.together(sy.sympify(eq.replace(eq_name+'=','')))
+ eq = eq.replace(InvGamma(a), 1/sy.gamma(a))
+ eq = eq.replace(Denom(a), 1/a)
+ eq = eq.replace(m(a), lambda a: masses[a])
+
+ series = sy.series(eq,eps,n=0)
+
+ eq_dict[eq_name].append(eq)
+ eq_dict[eq_name].append(series.coeff(eps**-1).simplify())
+ eq_dict[eq_name].append(series.coeff(eps**-2).simplify())
diff --git a/pyDipole/__init__.py b/pyDipole/__init__.py
@@ -0,0 +1,2 @@
+from .particles import Particle, particles
+from .process import Process
diff --git a/pyDipole/particles.py b/pyDipole/particles.py
@@ -0,0 +1,62 @@
+__all__ = [
+ 'Particle',
+ 'particles',
+ ]
+"""
+Style: name: [mass(string),isQCD(bool)]
+"""
+
+particle_table = {
+ 'u' : ['0' , True],
+ 'ubar' : ['0' , True],
+ 't' : ['mt', True],
+ 'tbar' : ['mt', True],
+ 'e' : ['0' , False],
+ 'ebar' : ['0' , False],
+ 'g' : ['0' , True]
+ }
+
+par_id = 1
+
+
+class Particle(object):
+ def __init__(self,name):
+ self.name = name
+ try:
+ self.mass = particle_table[name][0]
+ self.isQCD = particle_table[name][1]
+ self.isMassive = False if self.mass == '0' else True
+ except KeyError:
+ raise KeyError("particle {} is not defined".format(name))
+ global par_id
+ self.id = par_id
+ par_id = par_id + 1
+
+ @property
+ def name(self):
+ return self.__name;
+
+ @name.setter
+ def name(self, name):
+ if not isinstance(name,str):
+ raise TypeError("Particle name should be a string")
+ else:
+ self.__name = name
+
+ def __str__(self):
+ return '{0}({1})'.format(self.__class__.__name__, self.__dict__)
+
+ def __setattr__(self,name,value):
+ super(Particle,self).__setattr__(name,value)
+
+def particles(name):
+ if not isinstance(name, str):
+ raise TypeError("Particle names should be given as a space separated string")
+ else:
+ particles = name.split()
+ output = tuple()
+ for p in particles:
+ particle=Particle(p)
+ output = output + (particle,)
+
+ return output
diff --git a/pyDipole/process.py b/pyDipole/process.py
@@ -0,0 +1,31 @@
+from .particles import Particle
+
+__all__ = [
+ 'Process',
+ ]
+
+class Process(object):
+ def __init__(self, initial, final, mF_list=[]):
+ """
+ initial and final should be a list of initial and final state particle
+ names for the Born process.
+ For the NLO correction all possible gluon emissions will be assumed.
+
+ particle name should be given as a string
+ """
+ self.all_particles = []
+ self.mF_list = []
+ for par in initial:
+ p = Particle(par)
+ setattr(p,'initial', True)
+ self.all_particles.append(p)
+
+ for par in final:
+ p = Particle(par)
+ setattr(p,'initial', False)
+ self.all_particles.append(p)
+
+ for par in final:
+ p = Particle(par)
+ setattr(p,'initial', False)
+ self.mF_list.append(p)
