4422.sagews

metric = Matrix([[1,0,0,0],[0,-1,0,0],[0,0,-1,0],[0,0,0,-1]])
class FourVector :
    vector = Matrix()
    def __init__ (self, vector_in) :
        self.vector = vector_in
    def __str__ (self) :
        return str(self.vector)
    def p4prod(self, p4other) :
        return FourVector(self.vector.transpose()*metric*p4other.vector)
    def __add__ (self, p4other) :
        return FourVector(self.vector+p4other.vector)
    def __multiply__ (self, scalar) :
        return FourVector(self.vector*scalar)
    def E(self) :
        return self.vector[0][0]
    def P(self) :
        return sqrt(self.Px()^2+self.Py()^2+self.Pz()^2)
    def Px(self) :
        return self.vector[1][0]
    def Py(self) :
        return self.vector[2][0]
    def Pz(self) :
        return self.vector[3][0]
    def Pt(self) :
        return sqrt(self.Px()*self.Px()+self.Py()*self.Py())
    def Phi(self) :
        return atan(self.Py()/self.Px())
    def Theta(self) :
        return atan(self.Pt()/self.Pz())
    def Eta(self) :
        return -log(tan(self.Theta()/2))
    def M(self) :
        return sqrt(self.E()^2-self.P()^2)

from sympy.utilities.codegen import codegen
Eb, pxb, pyb, pzb = var('Eb pxb pyb pzb')
El, pxl, pyl, pzl = var('El pxl pyl pzl')
p4b = FourVector(Matrix([[Eb], [pxb], [pyb], [pzb]]))
p4l = FourVector(Matrix([[El], [pxl], [pyl], [pzl]]))
p4bl = p4b+p4l
mass = p4bl.M()
ptb = var('ptb')
phib = var('phib')
etab = var('etab')
mb = var('mb')
print mass
mass = mass(Eb=sqrt(pxb^2+pyb^2+pzb^2+mb^2))
print mass
mass = mass(pxb=ptb*cos(phib), pyb=ptb*sin(phib), pzb = ptb*sinh(etab))
print mass

sqrt(-(pzb + pzl)^2 - (pyb + pyl)^2 - (pxb + pxl)^2 + (Eb + El)^2) sqrt(-(pzb + pzl)^2 - (pyb + pyl)^2 - (pxb + pxl)^2 + (El + sqrt(mb^2 + pxb^2 + pyb^2 + pzb^2))^2) sqrt((El + sqrt(ptb^2*sin(phib)^2 + ptb^2*cos(phib)^2 + ptb^2*sinh(etab)^2 + mb^2))^2 - (ptb*sinh(etab) + pzl)^2 - (ptb*cos(phib) + pxl)^2 - (ptb*sin(phib) + pyl)^2)

dpt = diff(mass, ptb)
dphi = diff(mass, phib)
deta = diff(mass, etab)
print dpt

-((ptb*sinh(etab) + pzl)*sinh(etab) + (ptb*cos(phib) + pxl)*cos(phib) + (ptb*sin(phib) + pyl)*sin(phib) - (El + sqrt(ptb^2*sin(phib)^2 + ptb^2*cos(phib)^2 + ptb^2*sinh(etab)^2 + mb^2))*(ptb*sin(phib)^2 + ptb*cos(phib)^2 + ptb*sinh(etab)^2)/sqrt(ptb^2*sin(phib)^2 + ptb^2*cos(phib)^2 + ptb^2*sinh(etab)^2 + mb^2))/sqrt((El + sqrt(ptb^2*sin(phib)^2 + ptb^2*cos(phib)^2 + ptb^2*sinh(etab)^2 + mb^2))^2 - (ptb*sinh(etab) + pzl)^2 - (ptb*cos(phib) + pxl)^2 - (ptb*sin(phib) + pyl)^2)

codegen(('dphi', dphi.simplify()), 'C', 'test')

[('test.c', '/******************************************************************************\n * Code generated with sympy 0.7.1 *\n * *\n * See http://www.sympy.org/ for more information. *\n * *\n * This file is part of \'project\' *\n ******************************************************************************/\n#include "test.h"\n#include <math.h>\n\ndouble dphi(double El, double etab, double mb, double phib, double ptb, double pxl, double pyl, double pzl) {\n\n return (-ptb*(ptb*sin(phib) + pyl)*cos(phib) + ptb*(ptb*cos(phib) + pxl)*sin(phib))/sqrt(pow(El + sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)), 2) - pow(ptb*sin(phib) + pyl, 2) - pow(ptb*cos(phib) + pxl, 2) - pow(ptb*sinh(etab) + pzl, 2));\n\n}\n'), ('test.h', "/******************************************************************************\n * Code generated with sympy 0.7.1 *\n * *\n * See http://www.sympy.org/ for more information. *\n * *\n * This file is part of 'project' *\n ******************************************************************************/\n\n\n#ifndef PROJECT__TEST__H\n#define PROJECT__TEST__H\n\ndouble dphi(double El, double etab, double mb, double phib, double ptb, double pxl, double pyl, double pzl);\n\n#endif\n\n")]

codegen(('deta', deta.simplify()), 'C', 'test')

[('test.c', '/******************************************************************************\n * Code generated with sympy 0.7.1 *\n * *\n * See http://www.sympy.org/ for more information. *\n * *\n * This file is part of \'project\' *\n ******************************************************************************/\n#include "test.h"\n#include <math.h>\n\ndouble deta(double El, double etab, double mb, double phib, double ptb, double pxl, double pyl, double pzl) {\n\n return (pow(ptb, 2)*(El + sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)))*sinh(etab)*cosh(etab)/sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)) - ptb*(ptb*sinh(etab) + pzl)*cosh(etab))/sqrt(pow(El + sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)), 2) - pow(ptb*sin(phib) + pyl, 2) - pow(ptb*cos(phib) + pxl, 2) - pow(ptb*sinh(etab) + pzl, 2));\n\n}\n'), ('test.h', "/******************************************************************************\n * Code generated with sympy 0.7.1 *\n * *\n * See http://www.sympy.org/ for more information. *\n * *\n * This file is part of 'project' *\n ******************************************************************************/\n\n\n#ifndef PROJECT__TEST__H\n#define PROJECT__TEST__H\n\ndouble deta(double El, double etab, double mb, double phib, double ptb, double pxl, double pyl, double pzl);\n\n#endif\n\n")]

codegen(('dpt', dpt.simplify()), 'C', 'test')

[('test.c', '/******************************************************************************\n * Code generated with sympy 0.7.1 *\n * *\n * See http://www.sympy.org/ for more information. *\n * *\n * This file is part of \'project\' *\n ******************************************************************************/\n#include "test.h"\n#include <math.h>\n\ndouble dpt(double El, double etab, double mb, double phib, double ptb, double pxl, double pyl, double pzl) {\n\n return -(-(El + sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)))*(ptb*pow(sin(phib), 2) + ptb*pow(cos(phib), 2) + ptb*pow(sinh(etab), 2))/sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)) + (ptb*sin(phib) + pyl)*sin(phib) + (ptb*cos(phib) + pxl)*cos(phib) + (ptb*sinh(etab) + pzl)*sinh(etab))/sqrt(pow(El + sqrt(pow(mb, 2) + pow(ptb, 2)*pow(sin(phib), 2) + pow(ptb, 2)*pow(cos(phib), 2) + pow(ptb, 2)*pow(sinh(etab), 2)), 2) - pow(ptb*sin(phib) + pyl, 2) - pow(ptb*cos(phib) + pxl, 2) - pow(ptb*sinh(etab) + pzl, 2));\n\n}\n'), ('test.h', "/******************************************************************************\n * Code generated with sympy 0.7.1 *\n * *\n * See http://www.sympy.org/ for more information. *\n * *\n * This file is part of 'project' *\n ******************************************************************************/\n\n\n#ifndef PROJECT__TEST__H\n#define PROJECT__TEST__H\n\ndouble dpt(double El, double etab, double mb, double phib, double ptb, double pxl, double pyl, double pzl);\n\n#endif\n\n")]

print dpt( El=47.582720710921379,
pxl=-12.889458145985095,
pyl=-45.28588817992447,
pzl=-6.8677151012416644,
ptb=66.54874415133537,
phib=-0.84628181284034631,
etab=-1.4691815611915151,
mb=4.1999999999990649)
print mass( El=47.582720710921379,
pxl=-12.889458145985095,
pyl=-45.28588817992447,
pzl=-6.8677151012416644,
ptb=66.54874415133537,
phib=-0.84628181284034631,
etab=-1.4691815611915151,
pzb=-136.94278362460605,
mb=4.1999999999990649)

0.720573472840344 96.2044411395947