Skip to content

Instantly share code, notes, and snippets.

@jyoshida-sci

jyoshida-sci/1README.md

Last active January 30, 2017 06:50
Show Gist options
  • Save jyoshida-sci/b73149cf2bee393809c0708926b7a996 to your computer and use it in GitHub Desktop.
Save jyoshida-sci/b73149cf2bee393809c0708926b7a996 to your computer and use it in GitHub Desktop.
Download ZIP
scientific calculation with python
Raw
1README.md

scientific calculation with python

##linear_gaussian_fitting.py linear fitting and gaussian fitting

##SciPy_Constants.py Physics constants

##solving_polyeq.py intersection point of 2 pol3-curvews.

##rollingmin_2darrayoperation.py How to find max and min on a moving window

##affine.py class for affine transformation

##decaymodes.py !not-validated! counting possible decay modes of Proton+Neutrons+Lambda-system.

##Kinematic fit tbdkinematicfit.py: for three-body-decay kinematicfit_2body.py: for two-body-decay, momthetaphi.txt

##logAnalyzer.py logAnalyzer_app.py and log.txt

Raw
affine.py
import numpy as np
import unittest
class Affine:
lorgx = []
lorgy = []
ldstx = []
ldsty = []
abp = np.array([[ 1.], [0.], [ 0.]])
cdq = np.array([[ 0.], [1.], [ 0.]])
iabp = np.array([[ 1.], [0.], [ 0.]])
icdq = np.array([[ 0.], [1.], [ 0.]])
def __init__(self):
self.lorgx = []
self.lorgy = []
self.ldstx = []
self.ldsty = []
self.abp = np.array([[ 1.], [0.], [ 0.]])
self.cdq = np.array([[ 0.], [1.], [ 0.]])
self.iabp = np.array([[ 1.], [0.], [ 0.]])
self.icdq = np.array([[ 0.], [1.], [ 0.]])
def add(self, orgx, orgy, dstx, dsty):
self.lorgx.append(orgx)
self.lorgy.append(orgy)
self.ldstx.append(dstx)
self.ldsty.append(dsty)
def solveEq(self, ox, oy, v):
#minimum least square method
m = np.zeros((3,3))
b = np.zeros((3,1))
for i in range(len(ox)):
m[0][0] += ox[i] * ox[i]
m[0][1] += ox[i] * oy[i]
m[0][2] += ox[i]
m[1][0] += ox[i] * oy[i]
m[1][1] += oy[i] * oy[i]
m[1][2] += oy[i]
m[2][0] += ox[i]
m[2][1] += oy[i]
m[2][2] += 1.0
b[0][0] += v[i] * ox[i]
b[1][0] += v[i] * oy[i]
b[2][0] += v[i]
x = np.linalg.solve(m, b)
return x
def SolveHelmert(self,lx, ly, lX, lY):
x=0.0
X=0.0
y=0.0
Y=0.0
xx=0.0
yy=0.0
xX=0.0
yY=0.0
yX=0.0
xY=0.0
N=0.0
for i in range(len(lx)):
x += lx[i]
X += lX[i]
y += ly[i]
Y += lY[i]
xx += lx[i] * lx[i]
yy += ly[i] * ly[i]
xX += lx[i] * lX[i]
yY += ly[i] * lY[i]
yX += ly[i] * lX[i]
xY += lx[i] * lY[i]
N += 1.0
prms = []
d = x*x + y*y - N*(xx+yy)
prms.append( (x*X + y*Y - N*(xX+yY)) / d )
prms.append( (y*X - x*Y - N*(yX-xY)) / d )
prms.append(-prms[1])
prms.append( prms[0])
prms.append( (X - prms[0] * x - prms[1] * y)/N )
prms.append( (Y - prms[0] * y + prms[1] * x)/N )
return prms;
def calc(self):
if len(self.lorgx) < 2:
return False
elif len(self.lorgx) == 2:
prms = self.SolveHelmert(self.lorgx, self.lorgy, self.ldstx, self.ldsty);
self.abp[0][0] = prms[0]
self.abp[1][0] = prms[1]
self.abp[2][0] = prms[4]
self.cdq[0][0] = prms[2]
self.cdq[1][0] = prms[3]
self.cdq[2][0] = prms[5]
iprms = self.SolveHelmert(self.ldstx, self.ldsty, self.lorgx, self.lorgy);
self.iabp[0][0] = iprms[0]
self.iabp[1][0] = iprms[1]
self.iabp[2][0] = iprms[4]
self.icdq[0][0] = iprms[2]
self.icdq[1][0] = iprms[3]
self.icdq[2][0] = iprms[5]
else:
self.abp = self.solveEq(self.lorgx, self.lorgy, self.ldstx)
self.cdq = self.solveEq(self.lorgx, self.lorgy, self.ldsty)
self.iabp = self.solveEq(self.ldstx, self.ldsty, self.lorgx)
self.icdq = self.solveEq(self.ldstx, self.ldsty, self.lorgy)
return True
def printParams(self):
print "x y |--> X Y) {0:.8f} {1:.8f} {2:.8f} {3:.8f} {4:.1f} {5:.1f}".format(
self.abp[0][0], self.abp[1][0], self.cdq[0][0], self.cdq[1][0], self.abp[2][0], self.cdq[2][0])
print "X Y |--> x y) {0:.8f} {1:.8f} {2:.8f} {3:.8f} {4:.1f} {5:.1f}".format(
self.iabp[0][0], self.iabp[1][0], self.icdq[0][0], self.icdq[1][0], self.iabp[2][0], self.icdq[2][0])
def affTrans(self, x, y):
X = x*self.abp[0][0] + y*self.abp[1][0] + self.abp[2][0]
Y = x*self.cdq[0][0] + y*self.cdq[1][0] + self.cdq[2][0]
return (X, Y)
def invTrans(self, X, Y):
x = X*self.iabp[0][0] + Y*self.iabp[1][0] + self.iabp[2][0]
y = X*self.icdq[0][0] + Y*self.icdq[1][0] + self.icdq[2][0]
return (x, y)
def affTransTuple(self, p):
return self.affTrans(p[0], p[1])
def invTransTuple(self, P):
return self.invTrans(P[0], P[1])
def a(self):
return self.abp[0][0]
def b(self):
return self.abp[1][0]
def p(self):
return self.abp[2][0]
def c(self):
return self.cdq[0][0]
def d(self):
return self.cdq[1][0]
def q(self):
return self.cdq[2][0]
def inva(self):
return self.iabp[0][0]
def invb(self):
return self.iabp[1][0]
def invp(self):
return self.iabp[2][0]
def invc(self):
return self.icdq[0][0]
def invd(self):
return self.icdq[1][0]
def invq(self):
return self.icdq[2][0]
class TestAffine(unittest.TestCase):
def setUp(self):
self.af = Affine()
def test_case1(self):
# self.assertEqual(0.0011, 0.001)#case: failed
self.assertEqual(self.af.a(), 1.0)
self.assertEqual(self.af.b(), 0.0)
self.assertEqual(self.af.c(), 0.0)
self.assertEqual(self.af.d(), 1.0)
self.assertEqual(self.af.p(), 0.0)
self.assertEqual(self.af.q(), 0.0)
self.assertEqual(self.af.inva(), 1.0)
self.assertEqual(self.af.invb(), 0.0)
self.assertEqual(self.af.invc(), 0.0)
self.assertEqual(self.af.invd(), 1.0)
self.assertEqual(self.af.invp(), 0.0)
self.assertEqual(self.af.invq(), 0.0)
ret = self.af.calc()
self.assertEqual(ret, False)
self.af.add(0, 0, 0, 0)
self.af.add(1, 0, 0, 2)
self.af.add(0, 1, -2, 0)
ret = self.af.calc()
self.assertEqual(ret, True)
self.assertEqual(self.af.a(), 0.0)
self.assertEqual(self.af.b(), -2.0)
self.assertEqual(self.af.c(), 2.0)
self.assertEqual(self.af.d(), 0.0)
self.assertEqual(self.af.p(), 0.0)
self.assertEqual(self.af.q(), 0.0)
self.assertEqual(self.af.inva(), 0.0)
self.assertEqual(self.af.invb(), 0.5)
self.assertEqual(self.af.invc(), -0.5)
self.assertEqual(self.af.invd(), 0.0)
self.assertEqual(self.af.invp(), 0.0)
self.assertEqual(self.af.invq(), 0.0)
myxy = (2, 4.2)
myXY = self.af.affTransTuple(myxy)
myxy2 = self.af.invTransTuple(myXY)
self.assertEqual(myxy, myxy2)
myxy = (-1, 999)
myXY = self.af.affTransTuple(myxy)
myxy2 = self.af.invTransTuple(myXY)
self.assertEqual(myxy, myxy2)
self.af.__init__()
self.assertEqual(self.af.a(), 1.0)
self.assertEqual(self.af.b(), 0.0)
self.assertEqual(self.af.c(), 0.0)
self.assertEqual(self.af.d(), 1.0)
self.assertEqual(self.af.p(), 0.0)
self.assertEqual(self.af.q(), 0.0)
self.assertEqual(self.af.inva(), 1.0)
self.assertEqual(self.af.invb(), 0.0)
self.assertEqual(self.af.invc(), 0.0)
self.assertEqual(self.af.invd(), 1.0)
self.assertEqual(self.af.invp(), 0.0)
self.assertEqual(self.af.invq(), 0.0)
self.af.__init__()
self.af.add(147.515, -149.851, 166.412, -157.347)
self.af.add(-152.437, 150.203, -135.746, 140.264)
self.af.calc()
self.assertAlmostEqual(self.af.a(), 0.99960369)
self.assertAlmostEqual(self.af.b(), -0.00774819)
self.assertAlmostEqual(self.af.c(), 0.00774819)
self.assertAlmostEqual(self.af.d(), 0.99960369)
self.assertAlmostEqual(self.af.p(), 17.7944, 4)
self.assertAlmostEqual(self.af.q(), -8.6984, 4)
self.assertAlmostEqual(self.af.inva(), 1.00033637)
self.assertAlmostEqual(self.af.invb(), 0.00775387)
self.assertAlmostEqual(self.af.invc(), -0.00775387)
self.assertAlmostEqual(self.af.invd(), 1.00033637)
self.assertAlmostEqual(self.af.invp(), -17.7329, 4)
self.assertAlmostEqual(self.af.invq(), 8.8393, 4)
if __name__ == '__main__':
#### unittest:
unittest.main()
#### usage:
af = Affine()
print "3points AffineTrans"
af.add(0, 0, 0, 0)
af.add(1, 0, 0, 2)
af.add(0, 1, -2, 0)
af.calc()
af.printParams()
myxy = (2,4)
myXY = af.affTransTuple(myxy)
myxy2 = af.invTransTuple(myXY)
af.__init__()
print "2points HelmertTrans"
af.add(0, 0, 0, 0)
af.add(1, 0, 0, 2)
af.calc()
af.printParams()
af.__init__()
print "2points HelmertTrans 2"
af.add(147.515, -149.851, 166.412, -157.347)
af.add(-152.437, 150.203, -135.746, 140.264)
af.calc()
af.printParams()
af.__init__()
print "4points AffineTrans"
af.add(992, 161, -59382, 99769)
af.add(1153, 53952, -59505, 69904)
af.add(18811, 185, -49437, 99732)
af.add(18990, 53958, -49548, 69866)
af.calc()
af.printParams()
print "---"
print af.a(), af.b(), af.c(), af.d(), af.p(), af.q()
print af.inva(), af.invb(), af.invc(), af.invd(), af.invp(), af.invq()
print "---"
myxy = (992, 161)
myXY = af.affTransTuple(myxy)
myxy2 = af.invTransTuple(myXY)
print myxy2
af.__init__()
af.printParams()
print "---"
print af.a(), af.b(), af.c(), af.d(), af.p(), af.q()
print af.inva(), af.invb(), af.invc(), af.invd(), af.invp(), af.invq()
print "---"
Raw
gistfile1.txt
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 28 15:33:21 2015
@author: jyoshida-sci
!not-validated! !not-validated! !not-validated! !not-validated! !not-validated!
!not-validated! !not-validated! !not-validated! !not-validated! !not-validated!
"""
import itertools
decaymodes = []
#for p in itertools.permutations('PPPPPNNNNNNLL||', 15):
for p in itertools.permutations('PPPNNLL||', 3+2+2+2):
mystr = "".join(p)
if mystr[0] == '|':
continue
if mystr[7] == '|':
continue
if '||' in mystr:
continue
tmptuple = []
flag = 0
fragments = mystr.split('|')
for f in range(len(fragments)):
nucleus = map(str,fragments[f])
nucleus.sort()
tmptuple.append("".join(nucleus))
if nucleus.count('P') == 0:
flag += 1
if flag > 0:
continue
decaymodes.append(tmptuple)
for c in range(len(decaymodes)):
decaymodes[c].sort()
decaymodes.sort()
#result = list(set(decaymodes)) !!DEKINAI!!
possiblemodes = []
tmp = ['X','X','XXXX']
for c in range(len(decaymodes)):
if tmp != decaymodes[c]:
possiblemodes.append(decaymodes[c])
tmp = decaymodes[c]
for c in range(len(possiblemodes)):
print possiblemodes[c]
print len(possiblemodes)
Raw
kinematicfit_2body.py
# -*- coding: utf-8 -*-
"""
Created on Wed Jul 27 13:48:51 2016
@author: jyoshida-sci, MyintKyawSoe
"""
import numpy as np
if __name__ == '__main__':
a0 = []
a0r = []#radian
SM = []
f = open("momthetaphi.txt","r")#Input data is p1, p2, theta1, theta2, phi1, phi2, and their uncertainties.
for row in f:
print row
vals = row.split()
a0.append( float(vals[0]) )
a0r.append( np.deg2rad(float(vals[0])))
SM.append( float(vals[1]) )
f.close()
A0 = np.array([a0[0], a0[1], a0[2], a0[3], a0[4], a0[5]])
Vr0 = np.zeros((6,6))
for i in range(6):
Vr0[i,i] = SM[i]**2
D = np.array([
[
np.sin(a0r[1])*np.cos(a0r[2]),
a0[0]*np.cos(a0r[1])*np.cos(a0r[2]),
-a0[0]*np.sin(a0r[1])*np.sin(a0r[2]),
np.sin(a0r[4])*np.cos(a0r[5]),
a0[3]*np.cos(a0r[4])*np.cos(a0r[5]),
-a0[3]*np.sin(a0r[4])*np.sin(a0r[5])
],[
np.sin(a0r[1])*np.sin(a0r[2]),
a0[0]*np.cos(a0r[1])*np.sin(a0r[2]),
a0[0]*np.sin(a0r[1])*np.cos(a0r[2]),
np.sin(a0r[4])*np.sin(a0r[5]),
a0[3]*np.cos(a0r[4])*np.sin(a0r[5]),
a0[3]*np.sin(a0r[4])*np.cos(a0r[5])
],[
np.cos(a0r[1]),
-a0[0]*np.sin(a0r[1]),
0.0,
np.cos(a0r[4]),
-a0[3]*np.sin(a0r[4]),
0.0
]
])
d = np.array([
[a0[0]*np.sin(a0r[1])*np.cos(a0r[2])+a0[3]*np.sin(a0r[4])*np.cos(a0r[5])],
[a0[0]*np.sin(a0r[1])*np.sin(a0r[2])+a0[3]*np.sin(a0r[4])*np.sin(a0r[5])],
[a0[0]*np.cos(a0r[1])+a0[3]*np.cos(a0r[4])]
]);
Dt = D.transpose()
Vd0 = D.dot(Vr0.dot(Dt))
VD = np.linalg.inv(Vd0)
lmbda = VD.dot(d)
lmbda_t = lmbda.transpose()
a = A0 - Vr0.dot(Dt.dot(lmbda))
A = []
Ar = []#radian
for k in range(6):
A.append( a[0,k])
Ar.append( np.deg2rad(a[0,k]))
Vr = Vr0 - Vr0.dot(Dt.dot(VD.dot(D.dot(Vr0))))
chi2 = lmbda_t.dot(d)
dd = np.array([
[A[0]*np.sin(Ar[1])*np.cos(Ar[2])+A[3]*np.sin(Ar[4])*np.cos(Ar[5])],
[A[0]*np.sin(Ar[1])*np.sin(Ar[2])+A[3]*np.sin(Ar[4])*np.sin(Ar[5])],
[A[0]*np.cos(Ar[1])+A[3]*np.cos(Ar[4])]
])
print "track1 parameters"
print "P\t{0}\t+-{1}".format(a[0,0], np.sqrt(Vr[0,0]))
print "theta\t{0}\t+-{1}".format(a[0,1], np.sqrt(Vr[1,1]))
print "phi\t{0}\t+-{1}".format(a[0,2], np.sqrt(Vr[2,2]))
print ""
print "track2 parameters"
print "P\t{0}\t+-{1}".format(a[0,3],np.sqrt(Vr[3,3]))
print "theta\t{0}\t+-{1}".format( a[0,4], np.sqrt(Vr[4,4]))
print "phi\t{0}\t+-{1}".format(a[0,5], np.sqrt(Vr[5,5]))
print "";
print "the constraints equations";
print dd
print "chi2";
print chi2
'''
output of original CPP code
track2 parameters
P 100 +-0.707107
theta 90 +-0.707107
phi -1.7949e-009 +-0.707107
track3 parameters
P 100 +-0.707107
theta 90 +-0.707107
phi 180 +-0.707107
the constraints equations are
[0;
3.527139231042815e-007;
3.527139147759126e-007]
'''
Raw
kinematicfit_nagaraA.py
# -*- coding: utf-8 -*-
"""
Created on 20160806
@author: jyoshida-sci
"""
import numpy as np
from scipy import optimize
def chisq(m1, dispflag=False):
B_xi = 0.0
m_Nuclear = 11174.866 #Carbon
m_xi = 1321.71
err_m_xi = 0.07
initial_m1 = m1
m2 = 3727.380
m3 = 2808.922
mom1 = 170.961
mom2 = 82.661
mom3 = 166.049
theta1 = np.deg2rad(44.9)
theta2 = np.deg2rad(57.7)
theta3 = np.deg2rad(156.2)
phi1 = np.deg2rad(337.5)
phi2 = np.deg2rad(174.9)
phi3 = np.deg2rad(143.0)
err_mom1 = 3.506
err_mom2 = 5.855
err_mom3 = 0.666
err_theta1 = np.deg2rad(4.72)
err_theta2 = np.deg2rad(8.41)
err_theta3 = np.deg2rad(1.86)
err_phi1 = np.deg2rad(4.64)
err_phi2 = np.deg2rad(7.21)
err_phi3 = np.deg2rad(2.05)
err_theta1_phi1 = 1.74231E-5
err_theta2_phi2 = 7.87901E-4
err_theta3_phi3 = 5.28957E-6
B_xi=0 #0.13
#a-matrix: input value
a0 = np.array([
[mom1],
[theta1],
[phi1],
[mom2],
[theta2],
[phi2],
[mom3],
[theta3],
[phi3],
[m_xi]
])
#constraints derivation matrix
D11 = mom1/np.sqrt(m1**2+mom1**2)
D21 = np.sin(theta1)*np.cos(phi1)
D31 = np.sin(theta1)*np.sin(phi1)
D41 = np.cos(theta1)
D12 = 0.0
D22 = mom1*np.cos(theta1)*np.cos(phi1)
D32 = mom1*np.cos(theta1)*np.sin(phi1)
D42 = -mom1*np.sin(theta1)
D13 = 0.0
D23 = -mom1*np.sin(theta1)*np.sin(phi1)
D33 = mom1*np.sin(theta1)*np.cos(phi1)
D43 = 0.0
D14 = mom2/np.sqrt(m2**2+mom2**2)
D24 = np.sin(theta2)*np.cos(phi2)
D34 = np.sin(theta2)*np.sin(phi2)
D44 = np.cos(theta2)
D15 = 0.0
D25 = mom2*np.cos(theta2)*np.cos(phi2)
D35 = mom2*np.cos(theta2)*np.sin(phi2)
D45 = -mom2*np.sin(theta2)
D16 = 0.0
D26 = -mom2*np.sin(theta2)*np.sin(phi2)
D36 = mom2*np.sin(theta2)*np.cos(phi2)
D46 = 0.0
D17 = mom3/np.sqrt(m3**2+mom3**2)
D27 = np.sin(theta3)*np.cos(phi3)
D37 = np.sin(theta3)*np.sin(phi3)
D47 = np.cos(theta3)
D18 = 0.0
D28 = mom3*np.cos(theta3)*np.cos(phi3)
D38 = mom3*np.cos(theta3)*np.sin(phi3)
D48 = -mom3*np.sin(theta3)
D19 = 0.0
D29 = -mom3*np.sin(theta3)*np.sin(phi3)
D39 = mom3*np.sin(theta3)*np.cos(phi3)
D49 = 0.0
D110 = -1.0
D210 = 0.0
D310 = 0.0
D410 = 0.0
D = np.array([
[D11,D12,D13,D14,D15,D16,D17,D18,D19,D110],
[D21,D22,D23,D24,D25,D26,D27,D28,D29,D210],
[D31,D32,D33,D34,D35,D36,D37,D38,D39,D310],
[D41,D42,D43,D44,D45,D46,D47,D48,D49,D410]
])
#variance-covariance matrix
Va = np.zeros((10, 10))
Va[0,0] = err_mom1**2
Va[1,1] = err_theta1**2
Va[2,2] = err_phi1**2
Va[1,2] = err_theta1_phi1
Va[2,1] = err_theta1_phi1
Va[3,3] = err_mom2**2
Va[4,4] = err_theta2**2
Va[5,5] = err_phi2**2
Va[4,5] = err_theta2_phi2
Va[5,4] = err_theta2_phi2
Va[6,6] = err_mom3**2
Va[7,7] = err_theta3**2
Va[8,8] = err_phi3**2
Va[7,8] = err_theta3_phi3
Va[8,7] = err_theta3_phi3
Va[9,9] = err_m_xi**2
#constraints
H1 = np.sqrt(m1**2+mom1**2)+np.sqrt(m2**2+mom2**2)+np.sqrt(m3**2+mom3**2)-(m_Nuclear + m_xi - B_xi)
H2 = mom1*np.sin(theta1)*np.cos(phi1)+mom2*np.sin(theta2)*np.cos(phi2)+mom3*np.sin(theta3)*np.cos(phi3)
H3 = mom1*np.sin(theta1)*np.sin(phi1)+mom2*np.sin(theta2)*np.sin(phi2)+mom3*np.sin(theta3)*np.sin(phi3)
H4 = mom1*np.cos(theta1)+mom2*np.cos(theta2)+mom3*np.cos(theta3)
d = np.array([
[H1],
[H2],
[H3],
[H4]
])
#unknown parameter
E = np.array([[ m1/np.sqrt(m1**2 +mom1**2)],
[0.0],
[0.0],
[0.0]
])
#solve
Dt = D.transpose()#10x4matrix
Vd0 = D.dot(Va.dot(Dt))#4x4matrix
Vd = np.linalg.inv(Vd0)#4x4matrix
E_t = E.transpose()#
VE0 = E_t.dot(Vd.dot(E))
VE = np.linalg.inv(VE0)
lambda0 = Vd.dot(d)
z = -VE.dot( E_t.dot(lambda0))
lmbda = lambda0 + Vd.dot(E.dot(z))
lambda_t = lmbda.transpose()
Vd_i = np.linalg.inv(Vd)
chi_square = lambda_t.dot(Vd_i.dot(lmbda))
chi_square_value = chi_square[0,0]
#caluclated value
a = a0 - Va.dot(Dt.dot(lmbda))
Vz = VE
V_lambda = Vd - Vd.dot(E.dot( VE.dot( E_t.dot(Vd))))
V = Va - Va.dot(Dt.dot(V_lambda.dot(D.dot(Va))))
new_mom1 = a[0,0]
new_theta1 = a[1,0]
new_phi1 = a[2,0]
new_mom2 = a[3,0]
new_theta2 = a[4,0]
new_phi2 = a[5,0]
new_mom3 = a[6,0]
new_theta3 = a[7,0]
new_phi3 = a[8,0]
new_m_xi = a[9,0]
new_err_mom1 = np.sqrt(V[0,0])
new_err_theta1 = np.sqrt(V[1,1])
new_err_phi1 = np.sqrt(V[2,2])
new_err_mom2 = np.sqrt(V[3,3])
new_err_theta2 = np.sqrt(V[4,4])
new_err_phi2 = np.sqrt(V[5,5])
new_err_mom3 = np.sqrt(V[6,6])
new_err_theta3 = np.sqrt(V[7,7])
new_err_phi3 = np.sqrt(V[8,8])
new_err_m_xi = np.sqrt(V[9,9])
new_m1 = initial_m1 + z[0,0]
new_err_m1 = np.sqrt(Vz[0,0])
new_H1 = np.sqrt(new_m1**2+new_mom1**2)+np.sqrt(m2**2+new_mom2**2)+np.sqrt(m3**2+new_mom3**2)-(m_Nuclear+new_m_xi-B_xi)
new_H2 = new_mom1*np.sin(new_theta1)*np.cos(new_phi1)+new_mom2*np.sin(new_theta2)*np.cos(new_phi2)+new_mom3*np.sin(new_theta3)*np.cos(new_phi3)
new_H3 = new_mom1*np.sin(new_theta1)*np.sin(new_phi1)+new_mom2*np.sin(new_theta2)*np.sin(new_phi2)+new_mom3*np.sin(new_theta3)*np.sin(new_phi3)
new_H4 = new_mom1*np.cos(new_theta1)+new_mom2*np.cos(new_theta2)+new_mom3*np.cos(new_theta3)
delta_H1 = d[0,0]/np.sqrt(Vd_i[0,0])
delta_H2 = d[1,0]/np.sqrt(Vd_i[1,1])
delta_H3 = d[2,0]/np.sqrt(Vd_i[2,2])
delta_H4 = d[3,0]/np.sqrt(Vd_i[3,3])
new_err_BLL = np.sqrt(4*0.006**2 + new_err_m1**2)
if dispflag == True:
print "new_H1 {0:.4f} {1:.4f}".format(new_H1,delta_H1)
print "new_H2 {0:.4f} {1:.4f}".format(new_H2,delta_H2)
print "new_H3 {0:.4f} {1:.4f}".format(new_H3,delta_H3)
print "new_H4 {0:.4f} {1:.4f}".format(new_H4,delta_H4)
print ""
print "chi_square {0:.4f}".format(chi_square_value)
print "momentum1 = {0:.3f} +- {1:.3f}".format(new_mom1,new_err_mom1)
print "momentum2 = {0:.3f} +- {1:.3f}".format(new_mom2,new_err_mom2)
print "momentum3 = {0:.3f} +- {1:.3f}".format(new_mom3,new_err_mom3)
print "theta1 = {0:.3f} +- {1:.3f}".format(np.rad2deg(new_theta1),np.rad2deg(new_err_theta1))
print "phi1 = {0:.3f} +- {1:.3f}".format(np.rad2deg(new_phi1),np.rad2deg(new_err_phi1))
print "theta2 = {0:.3f} +- {1:.3f}".format(np.rad2deg(new_theta2),np.rad2deg(new_err_theta2))
print "phi2 = {0:.3f} +- {1:.3f}".format(np.rad2deg(new_phi2),np.rad2deg(new_err_phi2))
print "theta3 = {0:.3f} +- {1:.3f}".format(np.rad2deg(new_theta3),np.rad2deg(new_err_theta3))
print "phi3 = {0:.3f} +- {1:.3f}".format(np.rad2deg(new_phi3),np.rad2deg(new_err_phi3))
print "m_xi = {0:.3f} +- {1:.3f}".format(new_m_xi, new_err_m_xi)
print "m_D = {0:.3f}".format(m1)
print "BLL = {0:.3f}".format(1115.683 * 2.0 + 3727.380 - m1)
print "m_D(new) = {0:.4f} +- {1:.4f}".format(new_m1, new_err_m1)
print "BLL(new) = {0:.4f} +- {1:.4f}".format(1115.683 * 2 + 3727.380 - new_m1, new_err_BLL)
return chi_square_value
if __name__ == '__main__':
# x_min = optimize.brent(f)
m1_min = optimize.brent(chisq,brack=(5950.0, 5955.0), maxiter=100 )
print "m1_min={0}".format(m1_min)
chisq(m1_min, True)
# m1 = 5950.102826
# chisq(m1)
'''
Mishina's output
m1_min = 5950.102826
new_H1 0.0008 -10.3660
new_H2 -0.2159 -0.8452
new_H3 0.3284 0.0261
new_H4 -0.2957 0.8892
chi_square 1.2697
momentum1 = 170.873 ± 3.362
momentum2 = 80.295 ± 5.407
momentum3 = 166.064 ± 0.665
theta1 = 49.302 ± 2.331
phi1 = 338.882 ± 3.159
theta2 = 59.448 ± 2.829
phi2 = 173.929 ± 5.454
theta3 = 156.710 ± 1.739
phi3 = 142.807 ± 2.017
m_xi = 1321.710 ± 0.070
m_D = 5950.103
BLL = 8.643
m_D(new) = 5952.0532 ± 0.1763
BLL(new) = 6.6928 ± 0.1767
This output
m1_min = 5954.09830062
new_H1 0.0008 11.4780
new_H2 -0.2159 -0.8452
new_H3 0.3284 0.0261
new_H4 -0.2957 0.8892
chi_square 1.2697
momentum1 = 170.873 +- 3.362
momentum2 = 80.295 +- 5.407
momentum3 = 166.064 +- 0.665
theta1 = 49.302 +- 2.331
phi1 = 338.882 +- 3.159
theta2 = 59.448 +- 2.829
phi2 = 173.929 +- 5.454
theta3 = 156.710 +- 1.739
phi3 = 142.807 +- 2.017
m_xi = 1321.710 +- 0.070
m_D = 5954.098
BLL = 4.648
m_D(new) = 5952.0532 +- 0.1762
BLL(new) = 6.6928 +- 0.1766
'''
Raw
linear_gaussian_fitting.py
# -*- coding: utf-8 -*-
"""
Created on 20160803
@author: jyoshida-sci
linear fitting and gaussian fitting
"""
from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
def pol1(x, a, b):
return a*x + b
def gaus(x,k,sigma,mean):
normfactor = 1.0/(sigma*np.sqrt(2.0*np.pi))
expx2pow = np.exp(-(x-mean)**2/(2*sigma**2))
return k*normfactor*expx2pow
def pol1fitting():
xs = []
ys = []
for i in range(1000):
x = i*0.01 - 5.0
y = pol1(x, 2.0, 1.0) + np.random.normal(0.0, 1.0)
xs.append(x)
ys.append(y)
arrx = np.array(xs)
arry = np.array(ys)
param, cov = curve_fit(pol1, arrx, arry)
print param
print cov
fs = []
for i in range(1000):
x = i*0.01 - 5.0
y = pol1(x, param[0], param[1])
fs.append(y)
plt.plot(xs, ys,'bo')#blue circle
plt.plot(xs, fs,'r')#red line
plt.show()
def gausfitting():
xs = []
ys = []
for i in range(1000):
x = i*0.01 - 5.0
y = gaus(x, 2.0, 1.0, 0.0) + np.random.normal(0.0, 0.1)
xs.append(x)
ys.append(y)
arrx = np.array(xs)
arry = np.array(ys)
param, cov = curve_fit(gaus, arrx, arry)
print param
print cov
fs = []
for i in range(1000):
x = i*0.01 - 5.0
y = gaus(x, param[0], param[1], param[2])
fs.append(y)
plt.plot(xs, ys,'bo')#blue circle
plt.plot(xs, fs,'r')#red line
plt.show()
if __name__ == '__main__':
pol1fitting()
gausfitting()
'''
[ 1.9874254 1.01249703]
[[ 1.21347556e-04 6.06738409e-07]
[ 6.06738409e-07 1.01123164e-03]]
<linear graph>
[ 1.97578425 0.98679569 0.003557 ]
[[ 4.91074126e-04 1.63509746e-04 4.70664315e-11]
[ 1.63509746e-04 1.63328314e-04 -1.65566333e-11]
[ 4.70664315e-11 -1.65566333e-11 1.63328312e-04]]
<Gaussian graph>
'''
Raw
log.txt
2016/06/21 16:12:30.547 636021223505536758 actionA start
2016/06/21 16:12:31.554 636021223515547331 actionA end
2016/06/21 16:14:36.010 636021224760168519 actionA start
2016/06/21 16:14:37.017 636021224770179092 actionA end
2016/06/21 16:14:37.018 636021224770189092 actionA start
2016/06/21 16:14:38.018 636021224780189664 actionA end
2016/06/21 16:14:38.019 636021224780199665 actionA start
2016/06/21 16:14:39.021 636021224790210237 actionA end
2016/06/21 16:14:39.022 636021224790220238 actionA start
2016/06/21 16:14:40.023 636021224800230811 actionA end
2016/06/21 16:14:39.022 636021224790220238 actionB start
2016/06/21 16:14:40.024 636021224800240811 actionA start
2016/06/21 16:14:41.024 636021224810241383 actionA end
2016/06/21 16:14:41.025 636021224810251384 actionA start
2016/06/21 16:14:42.025 636021224820251956 actionA end
2016/06/21 16:14:42.026 636021224820261956 actionA start
2016/06/21 16:14:43.026 636021224830262528 actionA end
2016/06/21 16:14:43.027 636021224830272529 actionA start
2016/06/21 16:14:44.028 636021224840283101 actionA end
2016/06/21 16:14:44.029 636021224840293102 actionA start
2016/06/21 16:14:45.029 636021224850293674 actionA end
2016/06/21 16:14:45.030 636021224850303675 actionA start
2016/06/21 16:14:40.023 636021224800230811 actionB end
2016/06/21 16:14:46.030 636021224860304247 actionA end
2016/06/21 16:14:46.031 636021224860314247 actionA start
2016/06/21 16:14:47.032 636021224870324820 actionA end
2016/06/21 16:14:47.033 636021224870334820 actionA start
2016/06/21 16:14:48.033 636021224880335392 actionA end
2016/06/21 16:14:48.034 636021224880345393 actionA start
2016/06/21 16:14:49.035 636021224890355965 actionA end
2016/06/21 16:14:49.036 636021224890365966 actionA start
2016/06/21 16:14:50.036 636021224900366538 actionA end
2016/06/21 16:14:50.037 636021224900376539 actionA start
2016/06/21 16:14:51.038 636021224910387111 actionA end
2016/06/21 16:14:51.039 636021224910397112 actionA start
2016/06/21 16:14:52.040 636021224920407684 actionA end
2016/06/21 16:14:52.041 636021224920417685 actionA start
2016/06/21 16:14:53.042 636021224930428257 actionA end
2016/06/21 16:14:53.043 636021224930438258 actionA start
2016/06/21 16:14:54.043 636021224940438830 actionA end
2016/06/21 16:14:54.044 636021224940448831 actionA start
2016/06/21 16:14:55.044 636021224950449403 actionA end
2016/06/21 16:14:55.045 636021224950459403 actionA start
2016/06/21 16:14:56.046 636021224960469976 actionA end
2016/06/21 16:14:56.047 636021224960479976 actionA start
2016/06/21 16:14:57.049 636021224970490549 actionA end
2016/06/21 16:14:57.050 636021224970500549 actionA start
2016/06/21 16:14:58.050 636021224980501121 actionA end
2016/06/21 16:14:58.051 636021224980511122 actionA start
2016/06/21 16:14:59.051 636021224990511694 actionA end
2016/06/21 16:14:59.052 636021224990521695 actionA start
2016/06/21 16:15:00.052 636021225000522267 actionA end
2016/06/21 16:15:00.053 636021225000532267 actionA start
2016/06/21 16:15:01.054 636021225010542840 actionA end
2016/06/21 16:15:01.055 636021225010552840 actionA start
2016/06/21 16:15:02.055 636021225020553412 actionA end
2016/06/21 16:15:02.056 636021225020563413 actionA start
2016/06/21 16:15:03.057 636021225030573985 actionA end
2016/06/21 16:15:03.058 636021225030583986 actionA start
2016/06/21 16:15:04.059 636021225040594559 actionA end
2016/06/21 16:15:04.060 636021225040604559 actionA start
2016/06/21 16:15:05.060 636021225050605131 actionA end
2016/06/21 16:15:05.061 636021225050615132 actionA start
2016/06/21 16:15:06.062 636021225060625704 actionA end
2016/06/21 16:15:06.063 636021225060635705 actionA start
2016/06/21 16:15:07.064 636021225070646277 actionA end
2016/06/21 16:15:07.065 636021225070656278 actionA start
2016/06/21 16:15:08.066 636021225080666851 actionA end
2016/06/21 16:15:08.067 636021225080676851 actionA start
2016/06/21 16:15:09.068 636021225090687424 actionA end
2016/06/21 16:15:09.069 636021225090697424 actionA start
2016/06/21 16:15:10.070 636021225100707997 actionA end
2016/06/21 16:15:10.071 636021225100717997 actionA start
2016/06/21 16:15:11.072 636021225110728570 actionA end
2016/06/21 16:15:11.073 636021225110738571 actionA start
2016/06/21 16:15:12.074 636021225120749143 actionA end
2016/06/21 16:15:12.075 636021225120759144 actionA start
2016/06/21 16:15:13.076 636021225130769716 actionA end
2016/06/21 16:15:13.077 636021225130779717 actionA start
2016/06/21 16:15:14.079 636021225140790289 actionA end
2016/06/21 16:15:14.080 636021225140800290 actionA start
2016/06/21 16:15:15.081 636021225150810863 actionA end
2016/06/21 16:15:15.082 636021225150820863 actionA start
2016/06/21 16:15:16.083 636021225160831436 actionA end
2016/06/21 16:15:16.084 636021225160841436 actionA start
2016/06/21 16:15:17.085 636021225170852009 actionA end
2016/06/21 16:15:17.086 636021225170862009 actionA start
2016/06/21 16:15:18.087 636021225180872582 actionA end
2016/06/21 16:15:18.088 636021225180882583 actionA start
2016/06/21 16:15:19.089 636021225190893155 actionA end
2016/06/21 16:15:19.090 636021225190903156 actionA start
2016/06/21 16:15:20.091 636021225200913728 actionA end
2016/06/21 16:15:20.092 636021225200923729 actionA start
2016/06/21 16:15:21.093 636021225210934301 actionA end
2016/06/21 16:15:21.094 636021225210944302 actionA start
2016/06/21 16:15:22.095 636021225220954875 actionA end
2016/06/21 16:15:22.096 636021225220964875 actionA start
2016/06/21 16:15:23.097 636021225230975448 actionA end
2016/06/21 16:15:23.098 636021225230985448 actionA start
2016/06/21 16:15:24.099 636021225240996021 actionA end
2016/06/21 16:15:24.100 636021225241006021 actionA start
2016/06/21 16:15:25.101 636021225251016594 actionA end
2016/06/21 16:15:25.102 636021225251026595 actionA start
2016/06/21 16:15:26.103 636021225261037167 actionA end
2016/06/21 16:15:26.104 636021225261047168 actionA start
2016/06/21 16:15:27.105 636021225271057740 actionA end
2016/06/21 16:15:27.106 636021225271067741 actionA start
2016/06/21 16:15:28.107 636021225281078313 actionA end
2016/06/21 16:15:28.108 636021225281088314 actionA start
2016/06/21 16:15:29.109 636021225291098887 actionA end
2016/06/21 16:15:29.110 636021225291108887 actionA start
2016/06/21 16:15:30.111 636021225301119460 actionA end
2016/06/21 16:15:30.112 636021225301129460 actionA start
2016/06/21 16:15:31.114 636021225311140033 actionA end
2016/06/21 16:15:31.115 636021225311150033 actionA start
2016/06/21 16:15:32.116 636021225321160606 actionA end
2016/06/21 16:15:32.117 636021225321170607 actionA start
2016/06/21 16:15:33.118 636021225331181179 actionA end
2016/06/21 16:15:33.119 636021225331191180 actionA start
2016/06/21 16:15:34.121 636021225341211753 actionA end
2016/06/21 16:15:34.122 636021225341221753 actionA start
2016/06/21 16:15:35.123 636021225351232326 actionA end
2016/06/21 16:15:35.124 636021225351242327 actionA start
2016/06/21 16:15:36.125 636021225361252899 actionA end
2016/06/21 16:15:36.126 636021225361262900 actionA start
2016/06/21 16:15:37.151 636021225371513486 actionA end
2016/06/21 16:15:37.152 636021225371523487 actionA start
2016/06/21 16:15:38.153 636021225381534059 actionA end
2016/06/21 16:15:38.154 636021225381544060 actionA start
2016/06/21 16:15:39.155 636021225391554632 actionA end
2016/06/21 16:15:39.156 636021225391564633 actionA start
2016/06/21 16:15:40.157 636021225401575205 actionA end
2016/06/21 16:15:40.158 636021225401585206 actionA start
2016/06/21 16:15:41.159 636021225411595779 actionA end
2016/06/21 16:15:41.160 636021225411605779 actionA start
2016/06/21 16:15:42.161 636021225421616352 actionA end
2016/06/21 16:15:42.162 636021225421626352 actionA start
2016/06/21 16:15:43.163 636021225431636925 actionA end
2016/06/21 16:15:43.164 636021225431646925 actionA start
2016/06/21 16:15:44.165 636021225441657498 actionA end
2016/06/21 16:15:44.166 636021225441667499 actionA start
2016/06/21 16:15:45.167 636021225451678071 actionA end
2016/06/21 16:15:45.169 636021225451698072 actionA start
2016/06/21 16:15:46.170 636021225461708645 actionA end
2016/06/21 16:15:46.171 636021225461718645 actionA start
2016/06/21 16:15:47.172 636021225471729218 actionA end
2016/06/21 16:15:47.173 636021225471739219 actionA start
2016/06/21 16:15:48.174 636021225481749791 actionA end
2016/06/21 16:15:48.175 636021225481759792 actionA start
2016/06/21 16:15:49.177 636021225491770364 actionA end
2016/06/21 16:15:49.178 636021225491780365 actionA start
2016/06/21 16:15:50.179 636021225501790937 actionA end
2016/06/21 16:15:50.180 636021225501800938 actionA start
2016/06/21 16:15:51.181 636021225511811511 actionA end
2016/06/21 16:15:51.182 636021225511821511 actionA start
2016/06/21 16:15:52.183 636021225521832084 actionA end
2016/06/21 16:15:52.184 636021225521842084 actionA start
2016/06/21 16:15:53.185 636021225531852657 actionA end
2016/06/21 16:15:53.186 636021225531862657 actionA start
2016/06/21 16:15:54.187 636021225541873230 actionA end
2016/06/21 16:15:54.188 636021225541883231 actionA start
2016/06/21 16:15:55.189 636021225551893803 actionA end
2016/06/21 16:15:55.190 636021225551903804 actionA start
2016/06/21 16:15:56.191 636021225561914376 actionA end
2016/06/21 16:15:56.192 636021225561924377 actionA start
2016/06/21 16:15:57.193 636021225571934949 actionA end
2016/06/21 16:15:57.194 636021225571944950 actionA start
2016/06/21 16:15:58.195 636021225581955523 actionA end
2016/06/21 16:15:58.196 636021225581965523 actionA start
2016/06/21 16:15:59.197 636021225591976096 actionA end
2016/06/21 16:15:59.198 636021225591986096 actionA start
2016/06/21 16:16:00.199 636021225601996669 actionA end
2016/06/21 16:16:00.201 636021225602016670 actionA start
2016/06/21 16:16:01.203 636021225612037243 actionA end
2016/06/21 16:16:01.204 636021225612047244 actionA start
2016/06/21 16:16:02.205 636021225622057816 actionA end
2016/06/21 16:16:02.206 636021225622067817 actionA start
2016/06/21 16:16:03.207 636021225632078389 actionA end
2016/06/21 16:16:03.208 636021225632088390 actionA start
2016/06/21 16:16:04.209 636021225642098963 actionA end
2016/06/21 16:16:04.210 636021225642108963 actionA start
2016/06/21 16:16:05.211 636021225652119536 actionA end
2016/06/21 16:16:05.212 636021225652129536 actionA start
2016/06/21 16:16:06.214 636021225662140109 actionA end
2016/06/21 16:16:06.216 636021225662160110 actionA start
2016/06/21 16:16:07.217 636021225672170683 actionA end
2016/06/21 16:16:07.218 636021225672180683 actionA start
2016/06/21 16:16:08.219 636021225682191256 actionA end
2016/06/21 16:16:08.220 636021225682201256 actionA start
2016/06/21 16:16:09.221 636021225692211829 actionA end
2016/06/21 16:16:09.222 636021225692221829 actionA start
2016/06/21 16:16:10.223 636021225702232402 actionA end
2016/06/21 16:16:10.224 636021225702242403 actionA start
2016/06/21 16:16:11.225 636021225712252975 actionA end
2016/06/21 16:16:11.226 636021225712262976 actionA start
2016/06/21 16:16:12.227 636021225722273548 actionA end
2016/06/21 16:16:12.228 636021225722283549 actionA start
2016/06/21 16:16:13.229 636021225732294121 actionA end
2016/06/21 16:16:13.230 636021225732304122 actionA start
2016/06/21 16:16:14.231 636021225742314695 actionA end
2016/06/21 16:16:14.232 636021225742324695 actionA start
2016/06/21 16:16:15.233 636021225752335268 actionA end
2016/06/21 16:16:15.234 636021225752345268 actionA start
2016/06/21 16:16:16.235 636021225762355841 actionA end
2016/06/21 16:29:39.682 636021233796895390 actionB start
2016/06/21 16:29:39.691 636021233796915391 actionA start
2016/06/21 16:29:40.693 636021233806935965 actionA end
2016/06/21 16:29:40.694 636021233806945965 actionA start
2016/06/21 16:29:41.695 636021233816956538 actionA end
2016/06/21 16:29:41.696 636021233816966538 actionA start
2016/06/21 16:29:42.697 636021233826977111 actionA end
2016/06/21 16:29:42.698 636021233826987111 actionA start
2016/06/21 16:29:43.699 636021233836997684 actionA end
2016/06/21 16:29:43.700 636021233837007685 actionA start
2016/06/21 16:29:44.701 636021233847018257 actionA end
2016/06/21 16:29:44.702 636021233847028258 actionB end
2016/06/21 16:29:44.703 636021233847038258 actionB start
2016/06/21 16:29:44.704 636021233847048259 actionA start
2016/06/21 16:29:45.705 636021233857058831 actionA end
2016/06/21 16:29:45.706 636021233857068832 actionA start
2016/06/21 16:29:46.707 636021233867079405 actionA end
2016/06/21 16:29:46.708 636021233867089405 actionA start
2016/06/21 16:29:47.709 636021233877099978 actionA end
2016/06/21 16:29:47.710 636021233877109978 actionA start
2016/06/21 16:29:48.712 636021233887120551 actionA end
2016/06/21 16:29:48.713 636021233887130551 actionA start
2016/06/21 16:29:49.714 636021233897141124 actionA end
2016/06/21 16:29:49.715 636021233897151125 actionB end
2016/06/21 16:29:49.716 636021233897161125 actionB start
2016/06/21 16:29:49.717 636021233897171126 actionA start
2016/06/21 16:29:50.717 636021233907171698 actionA end
2016/06/21 16:29:50.718 636021233907181698 actionA start
2016/06/21 16:29:51.719 636021233917192271 actionA end
2016/06/21 16:29:51.720 636021233917202271 actionA start
2016/06/21 16:29:52.721 636021233927212844 actionA end
2016/06/21 16:29:52.722 636021233927222845 actionA start
2016/06/21 16:29:53.723 636021233937233417 actionA end
2016/06/21 16:29:53.724 636021233937243418 actionA start
2016/06/21 16:29:54.725 636021233947253990 actionA end
2016/06/21 16:29:54.726 636021233947263991 actionB end
2016/06/21 16:29:54.727 636021233947273991 actionB start
2016/06/21 16:29:54.728 636021233947283992 actionA start
2016/06/21 16:29:55.729 636021233957294565 actionA end
2016/06/21 16:29:55.730 636021233957304565 actionA start
2016/06/21 16:29:56.731 636021233967315138 actionA end
2016/06/21 16:29:56.732 636021233967325138 actionA start
2016/06/21 16:29:57.733 636021233977335711 actionA end
2016/06/21 16:29:57.734 636021233977345711 actionA start
2016/06/21 16:29:58.735 636021233987356284 actionA end
2016/06/21 16:29:58.736 636021233987366285 actionA start
2016/06/21 16:29:59.737 636021233997376857 actionA end
2016/06/21 16:29:59.738 636021233997386858 actionB end
2016/06/21 16:29:59.739 636021233997396858 actionB start
2016/06/21 16:29:59.740 636021233997406859 actionA start
2016/06/21 16:30:00.741 636021234007417431 actionA end
2016/06/21 16:30:00.742 636021234007427432 actionA start
2016/06/21 16:30:01.743 636021234017438005 actionA end
2016/06/21 16:30:01.744 636021234017448005 actionA start
2016/06/21 16:30:02.745 636021234027458578 actionA end
2016/06/21 16:30:02.746 636021234027468578 actionA start
2016/06/21 16:30:03.747 636021234037479151 actionA end
2016/06/21 16:30:03.748 636021234037489151 actionA start
2016/06/21 16:30:04.749 636021234047499724 actionA end
2016/06/21 16:30:04.750 636021234047509725 actionB end
2016/06/21 16:30:04.751 636021234047519725 actionB start
2016/06/21 16:30:04.752 636021234047529726 actionA start
2016/06/21 16:30:05.754 636021234057540298 actionA end
2016/06/21 16:30:05.755 636021234057550299 actionA start
2016/06/21 16:30:06.756 636021234067560871 actionA end
2016/06/21 16:30:06.757 636021234067570872 actionA start
2016/06/21 16:30:07.758 636021234077581445 actionA end
2016/06/21 16:30:07.759 636021234077591445 actionA start
2016/06/21 16:30:08.760 636021234087602018 actionA end
2016/06/21 16:30:08.761 636021234087612018 actionA start
2016/06/21 16:30:09.762 636021234097622591 actionA end
2016/06/21 16:30:09.763 636021234097632591 actionB end
2016/06/21 16:30:09.764 636021234097642592 actionB start
2016/06/21 16:30:09.765 636021234097652593 actionA start
2016/06/21 16:30:10.766 636021234107663165 actionA end
2016/06/21 16:30:10.767 636021234107673166 actionA start
2016/06/21 16:30:11.768 636021234117683738 actionA end
2016/06/21 16:30:11.769 636021234117693739 actionA start
2016/06/21 16:30:12.770 636021234127704311 actionA end
2016/06/21 16:30:12.771 636021234127714312 actionA start
2016/06/21 16:30:13.772 636021234137724885 actionA end
2016/06/21 16:30:13.773 636021234137734885 actionA start
2016/06/21 16:30:14.775 636021234147755458 actionA end
2016/06/21 16:30:14.776 636021234147765459 actionB end
2016/06/21 16:30:14.777 636021234147775459 actionB start
2016/06/21 16:30:14.778 636021234147785460 actionA start
2016/06/21 16:30:15.779 636021234157796033 actionA end
2016/06/21 16:30:15.780 636021234157806033 actionA start
2016/06/21 16:30:16.781 636021234167816606 actionA end
2016/06/21 16:30:16.782 636021234167826606 actionA start
2016/06/21 16:30:17.783 636021234177837179 actionA end
2016/06/21 16:30:17.784 636021234177847179 actionA start
2016/06/21 16:30:18.785 636021234187857752 actionA end
2016/06/21 16:30:18.786 636021234187867753 actionA start
2016/06/21 16:30:19.787 636021234197878325 actionA end
2016/06/21 16:30:19.788 636021234197888326 actionB end
2016/06/21 16:30:19.789 636021234197898326 actionB start
2016/06/21 16:30:19.790 636021234197908327 actionA start
2016/06/21 16:30:20.791 636021234207918899 actionA end
2016/06/21 16:30:20.792 636021234207928900 actionA start
2016/06/21 16:30:21.793 636021234217939473 actionA end
2016/06/21 16:30:21.794 636021234217949473 actionA start
2016/06/21 16:30:22.796 636021234227960046 actionA end
2016/06/21 16:30:22.797 636021234227970046 actionA start
2016/06/21 16:30:23.798 636021234237980619 actionA end
2016/06/21 16:30:23.799 636021234237990619 actionA start
2016/06/21 16:30:24.801 636021234248011193 actionA end
2016/06/21 16:30:24.802 636021234248021193 actionB end
2016/06/21 16:30:24.803 636021234248031194 actionB start
2016/06/21 16:30:24.804 636021234248041194 actionA start
2016/06/21 16:30:25.805 636021234258051767 actionA end
2016/06/21 16:30:25.806 636021234258061767 actionA start
2016/06/21 16:30:26.807 636021234268072340 actionA end
2016/06/21 16:30:26.808 636021234268082341 actionA start
2016/06/21 16:30:27.809 636021234278092913 actionA end
2016/06/21 16:30:27.810 636021234278102914 actionA start
2016/06/21 16:30:28.811 636021234288113486 actionA end
2016/06/21 16:30:28.812 636021234288123487 actionA start
2016/06/21 16:30:29.813 636021234298134060 actionA end
2016/06/21 16:30:29.814 636021234298144060 actionB end
2016/06/21 16:30:29.815 636021234298154061 actionB start
2016/06/21 16:30:29.816 636021234298164061 actionA start
2016/06/21 16:30:30.817 636021234308174634 actionA end
2016/06/21 16:30:30.818 636021234308184634 actionA start
2016/06/21 16:30:31.819 636021234318195207 actionA end
2016/06/21 16:30:31.820 636021234318205208 actionA start
2016/06/21 16:30:32.821 636021234328215780 actionA end
2016/06/21 16:30:32.822 636021234328225781 actionA start
2016/06/21 16:30:33.823 636021234338236353 actionA end
2016/06/21 16:30:33.824 636021234338246354 actionA start
2016/06/21 16:30:34.825 636021234348256926 actionA end
2016/06/21 16:30:34.826 636021234348266927 actionB end
2016/06/21 16:30:34.827 636021234348276928 actionB start
2016/06/21 16:30:34.828 636021234348286928 actionA start
2016/06/21 16:30:35.829 636021234358297501 actionA end
2016/06/21 16:30:35.888 636021234358887534 actionA start
2016/06/21 16:30:36.889 636021234368898107 actionA end
2016/06/21 16:30:36.890 636021234368908108 actionA start
2016/06/21 16:30:37.891 636021234378918680 actionA end
2016/06/21 16:30:37.892 636021234378928681 actionA start
2016/06/21 16:30:38.893 636021234388939253 actionA end
2016/06/21 16:30:38.894 636021234388949254 actionA start
2016/06/21 16:30:39.895 636021234398959826 actionA end
2016/06/21 16:30:39.896 636021234398969827 actionB end
2016/06/21 16:30:39.897 636021234398979828 actionB start
2016/06/21 16:30:39.898 636021234398989828 actionA start
2016/06/21 16:30:40.900 636021234409000401 actionA end
2016/06/21 16:30:40.901 636021234409010401 actionA start
2016/06/21 16:30:41.902 636021234419020974 actionA end
2016/06/21 16:30:41.903 636021234419030974 actionA start
2016/06/21 16:30:42.904 636021234429041547 actionA end
2016/06/21 16:30:42.905 636021234429051548 actionA start
2016/06/21 16:30:43.906 636021234439062120 actionA end
2016/06/21 16:30:43.907 636021234439072121 actionA start
2016/06/21 16:30:44.908 636021234449082693 actionA end
2016/06/21 16:30:44.909 636021234449092694 actionB end
2016/06/21 16:30:44.910 636021234449102694 actionB start
2016/06/21 16:30:44.911 636021234449112695 actionA start
2016/06/21 16:30:45.912 636021234459123268 actionA end
2016/06/21 16:30:45.914 636021234459143269 actionA start
2016/06/21 16:30:46.915 636021234469153841 actionA end
2016/06/21 16:30:46.916 636021234469163842 actionA start
2016/06/21 16:30:47.917 636021234479174414 actionA end
2016/06/21 16:30:47.918 636021234479184415 actionA start
2016/06/21 16:30:48.919 636021234489194988 actionA end
2016/06/21 16:30:48.920 636021234489204988 actionA start
2016/06/21 16:30:49.921 636021234499215561 actionA end
2016/06/21 16:30:49.922 636021234499225561 actionB end
2016/06/21 16:30:49.923 636021234499235562 actionB start
2016/06/21 16:30:49.924 636021234499245562 actionA start
2016/06/21 16:30:50.925 636021234509256135 actionA end
2016/06/21 16:30:50.926 636021234509266136 actionA start
2016/06/21 16:30:51.927 636021234519276708 actionA end
2016/06/21 16:30:51.928 636021234519286709 actionA start
2016/06/21 16:30:52.929 636021234529297281 actionA end
2016/06/21 16:30:52.930 636021234529307282 actionA start
2016/06/21 16:30:53.931 636021234539317854 actionA end
2016/06/21 16:30:53.932 636021234539327855 actionA start
2016/06/21 16:30:54.933 636021234549338428 actionA end
2016/06/21 16:30:54.934 636021234549348428 actionB end
2016/06/21 16:30:54.935 636021234549358429 actionB start
2016/06/21 16:30:54.936 636021234549368429 actionA start
2016/06/21 16:30:55.937 636021234559379002 actionA end
2016/06/21 16:30:55.938 636021234559389002 actionA start
2016/06/21 16:30:56.940 636021234569409576 actionA end
2016/06/21 16:30:56.941 636021234569419576 actionA start
2016/06/21 16:30:57.943 636021234579430149 actionA end
2016/06/21 16:30:57.944 636021234579440149 actionA start
2016/06/21 16:30:58.945 636021234589450722 actionA end
2016/06/21 16:30:58.946 636021234589460722 actionA start
2016/06/21 16:30:59.947 636021234599471295 actionA end
2016/06/21 16:30:59.948 636021234599481296 actionB end
2016/06/21 16:30:59.949 636021234599491296 actionB start
2016/06/21 16:30:59.950 636021234599501297 actionA start
2016/06/21 16:31:00.951 636021234609511869 actionA end
2016/06/21 16:31:00.952 636021234609521870 actionA start
2016/06/21 16:31:01.954 636021234619542443 actionA end
2016/06/21 16:31:01.955 636021234619552444 actionA start
2016/06/21 16:31:02.956 636021234629563016 actionA end
2016/06/21 16:31:02.957 636021234629573017 actionA start
2016/06/21 16:31:03.958 636021234639583589 actionA end
2016/06/21 16:31:03.959 636021234639593590 actionA start
2016/06/21 16:31:04.960 636021234649604162 actionA end
2016/06/21 16:31:04.961 636021234649614163 actionB end
2016/06/21 16:31:04.962 636021234649624164 actionB start
2016/06/21 16:31:04.963 636021234649634164 actionA start
2016/06/21 16:31:05.964 636021234659644737 actionA end
2016/06/21 16:31:05.965 636021234659654737 actionA start
2016/06/21 16:31:06.966 636021234669665310 actionA end
2016/06/21 16:31:06.967 636021234669675310 actionA start
2016/06/21 16:31:07.968 636021234679685883 actionA end
2016/06/21 16:31:07.969 636021234679695884 actionA start
2016/06/21 16:31:08.970 636021234689706456 actionA end
2016/06/21 16:31:08.971 636021234689716457 actionA start
2016/06/21 16:31:09.972 636021234699727029 actionA end
2016/06/21 16:31:09.973 636021234699737030 actionB end
2016/06/21 16:31:09.974 636021234699747030 actionB start
2016/06/21 16:31:09.975 636021234699757031 actionA start
2016/06/21 16:31:10.976 636021234709767604 actionA end
2016/06/21 16:31:10.977 636021234709777604 actionA start
2016/06/21 16:31:11.978 636021234719788177 actionA end
2016/06/21 16:31:11.979 636021234719798177 actionA start
2016/06/21 16:31:12.980 636021234729808750 actionA end
2016/06/21 16:31:12.981 636021234729818750 actionA start
2016/06/21 16:31:13.982 636021234739829323 actionA end
2016/06/21 16:31:13.983 636021234739839324 actionA start
2016/06/21 16:31:14.984 636021234749849896 actionA end
2016/06/21 16:31:14.985 636021234749859897 actionB end
2016/06/21 16:31:14.986 636021234749869897 actionB start
2016/06/21 16:31:14.987 636021234749879898 actionA start
2016/06/21 16:31:15.989 636021234759890470 actionA end
2016/06/21 16:31:15.990 636021234759900471 actionA start
2016/06/21 16:31:16.991 636021234769911044 actionA end
2016/06/21 16:31:16.992 636021234769921044 actionA start
2016/06/21 16:31:17.993 636021234779931617 actionA end
2016/06/21 16:31:17.994 636021234779941617 actionA start
2016/06/21 16:31:18.995 636021234789952190 actionA end
2016/06/21 16:31:18.996 636021234789962190 actionA start
2016/06/21 16:31:19.997 636021234799972763 actionA end
2016/06/21 16:31:19.998 636021234799982764 actionB end
Raw
logAnalyzer.py
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 21 2016
@author: jyoshida-sci
imputfile:
2016/06/21 16:12:30.547 636021223505536758 surface_recog start
2016/06/21 16:12:31.554 636021223515547331 surface_recog end
"""
from ROOT import ROOT, gROOT, gStyle, TCanvas, TH1D
import numpy as np
class logAnalyzer:
def __init__(self):
self.loglist = []
self.difflist = []# in sec units
def readLogFile(self, filename):
try:
for line in open(filename, 'r'):
l = line.split()
self.loglist.append(l)
except:
print "fileopen error: " + filename
def makeHistimg(self, actionname, imgfilename):
gROOT.Reset()
c1 = TCanvas( 'c1', 'Canvas', 800, 600 )
h = TH1D("h", actionname + ";time[sec]",50, 0, 5.5)
gStyle.SetOptStat(0)
start = 0
end = 0
self.difflist = []
for rec in self.loglist:
if "start" in rec[4] and actionname in rec[3]:
start = int(rec[2])
if "end" in rec[4] and actionname in rec[3]:
end = int(rec[2])
diff = (end-start)/10000.0/1000.0 #second
#print diff
self.difflist.append(diff)
h.Fill(diff)
h.Draw()
c1.Print(imgfilename)
def calcAverage(self):
return np.average(self.difflist)
def calcStdev(self):
return np.std(self.difflist)
if __name__ == '__main__':
la2 = logAnalyzer()
la2.readLogFile("aaaaa.txt")#not exist -> exception
la = logAnalyzer()
la.readLogFile("log.txt")
la.makeHistimg("tracking","hist.png")
la.makeHistimg("tracking","hist.pn")#invalid file will be produced if extention is invalid
print "{0:.2f}+-{1:.2f}".format(la.calcAverage(), la.calcStdev())
la.makeHistimg("surface_recog","hist2.png")
print "{0:.2f}+-{1:.2f}".format(la.calcAverage(), la.calcStdev())
Raw
logAnalyzer_app.py
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 21 2016
@author: jyoshida-sci
imputfile:
2016/06/21 16:12:30.547 636021223505536758 surface_recog start
2016/06/21 16:12:31.554 636021223515547331 surface_recog end
"""
import logAnalyzer
if __name__ == '__main__':
la = logAnalyzer.logAnalyzer()
la.readLogFile("log.txt")
la.makeHistimg("tracking","hist_trackingtime.png")
Raw
momthetaphi.txt
100 1.0 #particle1 |mom1|=100.0+-1.0MeV
90 1.0 #theta1 90.0+-1.0deg
0 1.0 #phi2 0.0+-1.0deg
100 1.0 #particle2 |mom1|=100.0+-1.0MeV
90 1.0 #theta2 90.0+-1.0deg
180 1.0 #phi2 180.0+-1.0deg
Raw
rollingmin_2darrayoperation.py
# -*- coding: utf-8 -*-
"""
Created on 20160802
@author: jyoshida-sci
"""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
if __name__ == '__main__':
#rolling minimum minimum of 3vals
#http://stackoverflow.com/questions/32436689/find-maxand-min-on-the-moving-interval-using-python
L = [5.5, 6.0, 6.0, 6.5, 6.0, 5.5, 5.5, 5.0, 4.5]
a = pd.DataFrame(L)
b = pd.rolling_min(a, 3)
print b
#part of matrix or image
myarray = np.random.randint(90,110,(5,6))#90~119random int numbers, 5x6 matrix
print "myarray"
print myarray
print ""
print "y=4 (5th row)"
print myarray[4,:]
print ""
print "x=2 (3rd column)"
print myarray[:,2]
# myarray.shape[0] = 5
# myarray.shape[1] = 6
for r in range(myarray.shape[0]):
verti = myarray[r,:]
horiz = range(myarray.shape[1])
plt.plot(horiz, verti, "o")
plt.ylim(0, 120)
plt.show()
for c in range(myarray.shape[1]):
verti = myarray[:,c]
horiz = range(myarray.shape[0])
plt.plot(horiz, verti, "o")
plt.ylim(0, 120)
plt.show()
Raw
SciPy_Constants.py
"""
Created on Wed Apr 15 00:07:54 2015
SciPy.constants
http://docs.scipy.org/doc/scipy/reference/constants.html#module-scipy.constants
@author: jyoshida-sci
"""
import numpy as np
import scipy.constants as const
#Physical constants
print(str.format('c: {0}', const.c))
print(str.format('the Planck constant: {0}', const.h))
print(str.format('standard acceleration of gravity: {0}', const.g))
print('\n')
#Constants database
names = np.array([
"electron mass energy equivalent in MeV",
"proton mass energy equivalent in MeV",
"neutron mass energy equivalent in MeV",
"deuteron mass energy equivalent in MeV",
"triton mass energy equivalent in MeV",
"alpha particle mass energy equivalent in MeV",
])
for i in range(len(names)):
a = const.physical_constants[names[i]]
print(str.format('{0}: ', names[i]))
print(str.format('{0:.5f} +- {1:.5f} [{2}]\n', a[0], a[2], a[1]))
"""
c: 299792458.0
the Planck constant: 6.62606957e-34
standard acceleration of gravity: 9.80665
electron mass energy equivalent in MeV:
0.51100 +- 0.00000 [MeV]
proton mass energy equivalent in MeV:
938.27205 +- 0.00002 [MeV]
neutron mass energy equivalent in MeV:
939.56538 +- 0.00002 [MeV]
deuteron mass energy equivalent in MeV:
1875.61286 +- 0.00004 [MeV]
triton mass energy equivalent in MeV:
2808.92101 +- 0.00006 [MeV]
alpha particle mass energy equivalent in MeV:
3727.37924 +- 0.00008 [MeV]
"""
Raw
solving_polyeq.py
#http://stackoverflow.com/questions/22742951/solve-an-equation-using-a-python-numerical-solver-in-numpy
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
# func0(x) = 2.0x^3 + 1.5x^2 + 1.0x^1 + 0.5
# func1(y) = 2.2y^3 + 1.7y^2 + 1.2y^1 + 0.7
# x+y=0.4
# when func0(x) = func1(y)?
func0 = lambda x : 2.0*x**3 + 1.5*x**2 + 1.0*x**1 + 0.5
func1 = lambda y : 2.2*y**3 + 1.7*y**2 + 1.2*y**1 + 0.7
func = lambda x : func0(x) - func1(-x+0.4)
# Plot it
x = np.linspace(-0.5, 1.5, 201)
plt.plot(x, func0(x))
plt.plot(x, func1(-x+0.4))
plt.plot(x, func(x))
plt.xlabel("range")
plt.ylabel("value")
plt.grid()
plt.show()
# Use the numerical solver to find the roots
r_initial_guess = 0.5
r_solution = fsolve(func, r_initial_guess)
print "The solution is tau = %f" % r_solution
print "at which the value of the expression is %f" % func(r_solution)
Raw
tbdkinematicfit.py
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 24 16:32:51 2015
non-validated non-validated non-validated non-validated non-validated non-validated
non-validated non-validated non-validated non-validated non-validated non-validated
non-validated non-validated non-validated non-validated non-validated non-validated
@author: jyoshida-sci
"""
import numpy as np
from numpy import sqrt,sin,cos
#data_in=fopen("data.txt","r")
mpi= 139.575#mass of pion
Mhy = 3922.565
m1 = 938.272
m2 = 2808.922
vx, vy = [], []
for l in open('001.txt').readlines():
data = l[:-1].split(' ')
vx += [float(data[0])]
vy += [float(data[1])]
idata = np.loadtxt('001.txt',delimiter=' ')
#a0 = np.array([33.012, 12.7, 314.8, 83.66, 47.9, 118.6, -34.6, 314.8])
#SM = np.array([1.385, 3.0, 15.0, 5.0, 10.0, 15.0, 2.0, 11.0])
#a0 = np.array([2.1073, 38.8518, 126.1315, 0.2061, 97.0867, 155.2786, 115.1204, 323.0121])
#SM = np.array([0.2, 1.0, 0.2, 0.3, 1.5, 0.3, 0.5, 0.1])
a0 = np.zeros(8)
SM = np.zeros(8)
a0[0] = vx[0]
a0[1] = np.deg2rad(vx[1])
a0[2] = np.deg2rad(vx[2])
a0[3] = vx[3]
a0[4] = np.deg2rad(vx[4])
a0[5] = np.deg2rad(vx[5])
a0[6] = np.deg2rad(vx[7])
a0[7] = np.deg2rad(vx[8])
SM[0] = vy[0]
SM[1] = np.deg2rad(vy[1])
SM[2] = np.deg2rad(vy[2])
SM[3] = vy[3]
SM[4] = np.deg2rad(vy[4])
SM[5] = np.deg2rad(vy[5])
SM[6] = np.deg2rad(vy[7])
SM[7] = np.deg2rad(vy[8])
#Pi momentum calculation from Energy conservation......................
Epi = Mhy - sqrt(a0[0]*a0[0]+m1*m1) - sqrt(a0[3]*a0[3]+m2*m2)
Ppi = sqrt(Epi*Epi-mpi*mpi)
# alpha0 matrix.............................
A0 = np.array([[a0[0]], [a0[1]], [a0[2]], [a0[3]], [a0[4]], [a0[5]], [a0[6]], [a0[7]]])
#V alpha matrix (covariance matrix)..........
Va = np.zeros((8,8))
for i in range(0, 8):
Va[i,i] = SM[i]*SM[i]
#D matrix (For Known parameters)..........
D = np.empty((4,8))
D[0,0]= sin(a0[1])*cos(a0[2])
D[0,1]= a0[0]*cos(a0[1])*cos(a0[2])
D[0,2]= -a0[0]*sin(a0[1])*sin(a0[2])
D[0,3]= sin(a0[4])*cos(a0[5])
D[0,4]= a0[3]*cos(a0[4])*cos(a0[5])
D[0,5]= -a0[3]*sin(a0[4])*sin(a0[5])
D[0,6]= Ppi*cos(a0[6])*cos(a0[7])
D[0,7]= -Ppi*sin(a0[4])*sin(a0[5])
D[1,0]= sin(a0[1])*sin(a0[2])
D[1,1]= a0[0]*cos(a0[1])*sin(a0[2])
D[1,2]= a0[0]*sin(a0[1])*cos(a0[2])
D[1,3]= sin(a0[4])*sin(a0[5])
D[1,4]= a0[3]*cos(a0[4])*sin(a0[5])
D[1,5]= a0[3]*sin(a0[4])*cos(a0[5])
D[1,6]= Ppi*cos(a0[6])*sin(a0[7])
D[1,7]= Ppi*sin(a0[4])*cos(a0[5])
D[2,0]= cos(a0[1])
D[2,1]= -a0[0]*sin(a0[1])
D[2,2]= 0
D[2,3]= cos(a0[4])
D[2,4]= -a0[3]*sin(a0[4])
D[2,5]= 0
D[2,6]= -Ppi*sin(a0[6])
D[2,7]= 0
D[3,0]= a0[0]/sqrt(m1*m1+a0[0]*a0[0])
D[3,1]= 0
D[3,2]= 0
D[3,3]= a0[3]/sqrt(m2*m2+a0[3]*a0[3])
D[3,4]= 0
D[3,5]= 0
D[3,6]= 0
D[3,7]= 0
# E matrix (for unknown paramenters)...............................
E = np.zeros((4,2))
E[0,0]= sin(a0[6])*cos(a0[7])
E[0,1]= 0
E[1,0]= sin(a0[6])*sin(a0[7])
E[1,1]= 0
E[2,0]= cos(a0[6])
E[2,1]= 0
E[3,0]= Ppi/sqrt(mpi*mpi+Ppi*Ppi)
E[3,1]= 1
#constraints equations....................................
d = np.zeros((4,1))
d[0,0]= a0[0]*sin(a0[1])*cos(a0[2])+a0[3]*sin(a0[4])*cos(a0[5])+Ppi*sin(a0[6])*cos(a0[7])
d[1,0]= a0[0]*sin(a0[1])*sin(a0[2])+a0[3]*sin(a0[4])*sin(a0[5])+Ppi*sin(a0[6])*sin(a0[7])
d[2,0]= a0[0]*cos(a0[1])+a0[3]*cos(a0[4])+Ppi*cos(a0[6])
d[3,0]= sqrt(a0[0]*a0[0]+m1*m1)+sqrt(a0[3]*a0[3]+m2*m2)+sqrt(Ppi*Ppi+mpi*mpi)-Mhy
#VD matrix
Dt = D.transpose()
Vd0 = D.dot(Va.dot(Dt))
Vd = np.linalg.inv(Vd0)
#VE matrix
E_T = E.transpose()
VE0 = E_T.dot(Vd.dot(E))
VE = np.linalg.inv(VE0)
#Largrium multipliers.............................
DdotA0 = D.dot(A0)
DdotA0minusd = np.add(DdotA0,d)
lambda0 = Vd.dot(DdotA0minusd)
# // Unknown parameters matrix.........................
z = -VE.dot(E_T.dot(lambda0))
lambd = lambda0 + Vd.dot(E.dot(z))
lambda_T = lambd.transpose()
Vd_i = np.linalg.inv(Vd)
# // chi2 ................................................
chi_2 = lambda_T.dot(Vd_i.dot(lambd))
Vadotlambda = Va.dot(Dt.dot(lambd))
A = A0 - Vadotlambda
#Vz = VE
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment