Commit 8418fb31 authored by Claus Kleinwort's avatar Claus Kleinwort
Browse files

python: access to residuals of measurements/scatterers added

git-svn-id: http://svnsrv.desy.de/public/GeneralBrokenLines/trunk@116 281f6f2b-e318-4fd1-8bce-1a4ba7aab212
parent 8c4683a0
This diff is collapsed.
...@@ -76,7 +76,7 @@ class BorderedBandMatrix(object): ...@@ -76,7 +76,7 @@ class BorderedBandMatrix(object):
## size of border; int ## size of border; int
self.__numBorder = nBorder self.__numBorder = nBorder
## size of border; int ## size of border; int
self.__numBand = 0 # actual band width self.__numBand = 0 # actual band width
## size of band part of matrix; int ## size of band part of matrix; int
self.__numCol = nSizeBand self.__numCol = nSizeBand
## border part B; matrix(float) ## border part B; matrix(float)
...@@ -122,6 +122,7 @@ class BorderedBandMatrix(object): ...@@ -122,6 +122,7 @@ class BorderedBandMatrix(object):
# #
def getBlockMatrix(self, aIndex): def getBlockMatrix(self, aIndex):
nBorder = self.__numBorder nBorder = self.__numBorder
nBand = self.__numBand
nSize = len(aIndex) nSize = len(aIndex)
aMatrix = np.empty((nSize, nSize)) aMatrix = np.empty((nSize, nSize))
for i in range(nSize): for i in range(nSize):
...@@ -135,8 +136,11 @@ class BorderedBandMatrix(object): ...@@ -135,8 +136,11 @@ class BorderedBandMatrix(object):
elif (iMin < nBorder): elif (iMin < nBorder):
aMatrix[i, j] = self.__mixed[iMin, iMax - nBorder] aMatrix[i, j] = self.__mixed[iMin, iMax - nBorder]
else: else:
nBand = iIndex - jIndex iBand = iMax - iMin
aMatrix[i, j] = self.__band[nBand, jIndex - nBorder] if iBand <= nBand:
aMatrix[i, j] = self.__band[iBand, iMin - nBorder]
else:
aMatrix[i, j] = 0.
aMatrix[j, i] = aMatrix[i, j] aMatrix[j, i] = aMatrix[i, j]
return aMatrix return aMatrix
...@@ -173,7 +177,7 @@ class BorderedBandMatrix(object): ...@@ -173,7 +177,7 @@ class BorderedBandMatrix(object):
def decomposeBand(): def decomposeBand():
nRow = self.__numBand + 1 nRow = self.__numBand + 1
nCol = self.__numCol nCol = self.__numCol
auxVec = np.copy(self.__band[0]) * 16.0 # save diagonal elements auxVec = np.copy(self.__band[0]) * 16.0 # save diagonal elements
for i in range(nCol): for i in range(nCol):
if ((self.__band[0, i] + auxVec[i]) != self.__band[0, i]): if ((self.__band[0, i] + auxVec[i]) != self.__band[0, i]):
self.__band[0, i] = 1.0 / self.__band[0, i] self.__band[0, i] = 1.0 / self.__band[0, i]
...@@ -194,10 +198,10 @@ class BorderedBandMatrix(object): ...@@ -194,10 +198,10 @@ class BorderedBandMatrix(object):
nRow = self.__numBand + 1 nRow = self.__numBand + 1
nCol = self.__numCol nCol = self.__numCol
aSolution = np.copy(aRightHandSide) aSolution = np.copy(aRightHandSide)
for i in range(nCol): # forward substitution for i in range(nCol): # forward substitution
for j in range(min(nRow, nCol - i) - 1): for j in range(min(nRow, nCol - i) - 1):
aSolution[j + i + 1] -= self.__band[j + 1, i] * aSolution[i] aSolution[j + i + 1] -= self.__band[j + 1, i] * aSolution[i]
for i in range(nCol - 1, -1, -1): # backward substitution for i in range(nCol - 1, -1, -1): # backward substitution
rxw = self.__band[0, i] * aSolution[i] rxw = self.__band[0, i] * aSolution[i]
for j in range(min(nRow, nCol - i) - 1): for j in range(min(nRow, nCol - i) - 1):
rxw -= self.__band[j + 1, i] * aSolution[j + i + 1] rxw -= self.__band[j + 1, i] * aSolution[j + i + 1]
......
...@@ -56,8 +56,8 @@ def example1(): ...@@ -56,8 +56,8 @@ def example1():
# #
np.random.seed(47117) np.random.seed(47117)
nTry = 1000 #: number of tries nTry = 1000 #: number of tries
nLayer = 5 #: number of detector layers nLayer = 5 #: number of detector layers
print " Gbltst $Rev$ ", nTry, nLayer print " Gbltst $Rev$ ", nTry, nLayer
start = time.clock() start = time.clock()
# track direction # track direction
...@@ -70,20 +70,20 @@ def example1(): ...@@ -70,20 +70,20 @@ def example1():
[-sinPhi, cosPhi, 0.], \ [-sinPhi, cosPhi, 0.], \
[-sinLambda * cosPhi, -sinLambda * sinPhi, cosLambda]]) [-sinLambda * cosPhi, -sinLambda * sinPhi, cosLambda]])
# measurement resolution # measurement resolution
measErr = np.array([ 0.001, 0.001]) # 10 mu measErr = np.array([ 0.001, 0.001]) # 10 mu
measPrec = 1.0 / measErr ** 2 measPrec = 1.0 / measErr ** 2
# scattering error # scattering error
scatErr = 0.001 # 1 mread scatErr = 0.001 # 1 mread
# RMS of CurviLinear track parameters (Q/P, slopes, offsets) # RMS of CurviLinear track parameters (Q/P, slopes, offsets)
clErr = np.array([0.001, -0.1, 0.2, -0.15, 0.25]) clErr = np.array([0.001, -0.1, 0.2, -0.15, 0.25])
# precision matrix for external seed (in local system) # precision matrix for external seed (in local system)
locSeed = None locSeed = None
seedLabel = 0 # label of point with seed seedLabel = 0 # label of point with seed
if seedLabel != 0: if seedLabel != 0:
print " external seed at label ", seedLabel print " external seed at label ", seedLabel
# #
bfac = 0.2998 # Bz*c for Bz=1 bfac = 0.2998 # Bz*c for Bz=1
step = 1.5 / cosLambda # constant steps in RPhi step = 1.5 / cosLambda # constant steps in RPhi
# #
Chi2Sum = 0. Chi2Sum = 0.
NdfSum = 0 NdfSum = 0
...@@ -142,7 +142,7 @@ def example1(): ...@@ -142,7 +142,7 @@ def example1():
# point.addLocals(addDer) # point.addLocals(addDer)
# additional global parameters? # additional global parameters?
point.addGlobals(labGlobal, addDer) point.addGlobals(labGlobal, addDer)
addDer = -addDer # locDer flips sign every measurement addDer = -addDer # locDer flips sign every measurement
# add point to trajectory # add point to trajectory
iLabel = traj.addPoint(point) iLabel = traj.addPoint(point)
if iLabel == abs(seedLabel): if iLabel == abs(seedLabel):
...@@ -199,7 +199,15 @@ def example1(): ...@@ -199,7 +199,15 @@ def example1():
locPar, locCov = traj.getResults(i) locPar, locCov = traj.getResults(i)
print " Point> ", i print " Point> ", i
print " locPar ", locPar print " locPar ", locPar
#print " locCov ", locCov #print " locCov ", locCov
# check residuals
for i in range(traj.getNumPoints()):
numData, aResiduals, aMeasErr, aResErr, aDownWeight = traj.getMeasResults(i + 1)
for j in range(numData):
print " measRes " , i, j, aResiduals[j], aMeasErr[j], aResErr[j], aDownWeight[j]
numData, aResiduals, aMeasErr, aResErr, aDownWeight = traj.getScatResults(i + 1)
for j in range(numData):
print " scatRes " , i, j, aResiduals[j], aMeasErr[j], aResErr[j], aDownWeight[j]
# #
end = time.clock() end = time.clock()
print " Time [s] ", end - start print " Time [s] ", end - start
...@@ -212,7 +220,7 @@ def example2(): ...@@ -212,7 +220,7 @@ def example2():
# #
binaryFile = open("milleBinaryISN.dat", "rb") binaryFile = open("milleBinaryISN.dat", "rb")
nRec = 0 nRec = 0
maxRec = 10 #: maximum number of records to read maxRec = 10 #: maximum number of records to read
Chi2Sum = 0. Chi2Sum = 0.
NdfSum = 0 NdfSum = 0
LostSum = 0. LostSum = 0.
...@@ -223,7 +231,7 @@ def example2(): ...@@ -223,7 +231,7 @@ def example2():
# create trajectory # create trajectory
traj = GblTrajectory(0) traj = GblTrajectory(0)
# read from file # read from file
traj.milleIn(binaryFile) # get data blocks from file traj.milleIn(binaryFile) # get data blocks from file
nRec += 1 nRec += 1
# fit trajectory # fit trajectory
Chi2, Ndf, Lost = traj.fit() Chi2, Ndf, Lost = traj.fit()
...@@ -275,7 +283,7 @@ class gblMeasSystem(object): ...@@ -275,7 +283,7 @@ class gblMeasSystem(object):
## products (T,U,V) * (I,J,K) ## products (T,U,V) * (I,J,K)
self.__prod = np.dot(curviDir, measDir.T) self.__prod = np.dot(curviDir, measDir.T)
## Transformation from measurement to local system ## Transformation from local to measurement system
# #
# @return None (systems are identical) # @return None (systems are identical)
# #
...@@ -289,8 +297,8 @@ class gblMeasSystem(object): ...@@ -289,8 +297,8 @@ class gblMeasSystem(object):
def getTransLocalToCurvi(self): def getTransLocalToCurvi(self):
meas2crv = np.zeros((5, 5)) meas2crv = np.zeros((5, 5))
meas2crv[0, 0] = 1. meas2crv[0, 0] = 1.
meas2crv[1:3, 1:3] = self.__prod[1:3, 1:3] * self.__prod[0, 0] # (U,V)*(J,K) * T*I meas2crv[1:3, 1:3] = self.__prod[1:3, 1:3] * self.__prod[0, 0] # (U,V)*(J,K) * T*I
meas2crv[3:5, 3:5] = self.__prod[1:3, 1:3] # (U,V)*(J,K) meas2crv[3:5, 3:5] = self.__prod[1:3, 1:3] # (U,V)*(J,K)
return meas2crv return meas2crv
## Transformation from curvilinear to measurement system ## Transformation from curvilinear to measurement system
...@@ -298,15 +306,15 @@ class gblMeasSystem(object): ...@@ -298,15 +306,15 @@ class gblMeasSystem(object):
# @return Transformation for offsets; 2*2 matrix # @return Transformation for offsets; 2*2 matrix
# #
def getTransCurviToMeas(self): def getTransCurviToMeas(self):
return np.linalg.inv(self.__prod[1:3, 1:3]) #((U,V)*(J,K))^-1 return np.linalg.inv(self.__prod[1:3, 1:3]) #((U,V)*(J,K))^-1
## Scattering precision matrix in local system ## Scattering precision matrix in local system
# #
# @return Precision; 2*2 matrix # @return Precision; 2*2 matrix
# #
def getScatPrecision(self, scatErr): def getScatPrecision(self, scatErr):
c1 = self.__prod[0, 1] # T*J c1 = self.__prod[0, 1] # T*J
c2 = self.__prod[0, 2] # T*K c2 = self.__prod[0, 2] # T*K
fac = (1 - c1 * c1 - c2 * c2) / (scatErr * scatErr) fac = (1 - c1 * c1 - c2 * c2) / (scatErr * scatErr)
scatP = np.empty((2, 2)) scatP = np.empty((2, 2))
scatP[0, 0] = fac * (1 - c1 * c1) scatP[0, 0] = fac * (1 - c1 * c1)
...@@ -331,33 +339,33 @@ class gblCurviSystem(object): ...@@ -331,33 +339,33 @@ class gblCurviSystem(object):
## projection curvilinear to measurement system: ((U,V)*(J,K))^-1 ## projection curvilinear to measurement system: ((U,V)*(J,K))^-1
self.__c2m = np.linalg.inv(np.dot(curviDir[1:3, 1:3], measDir[1:3, 1:3].T)) self.__c2m = np.linalg.inv(np.dot(curviDir[1:3, 1:3], measDir[1:3, 1:3].T))
## Transformation from measurement to local system ## Transformation from local to measurement system
# #
# @return Transformation for offsets; 2*2 matrix # @return Transformation for offsets; 2*2 matrix
# #
def getTransLocalToMeas(self): def getTransLocalToMeas(self):
return self.__c2m #((U,V)*(J,K))^-1 return self.__c2m #((U,V)*(J,K))^-1
## Transformation of (q/p, slopes, offsets) from local to curvilinear system ## Transformation of (q/p, slopes, offsets) from local to curvilinear system
# #
# @return Transformation for track parameters; 5*5 matrix # @return Transformation for track parameters; 5*5 matrix
# #
def getTransLocalToCurvi(self): def getTransLocalToCurvi(self):
return np.eye(5) # unit matrix return np.eye(5) # unit matrix
## Transformation from curvilinear to measurement system ## Transformation from curvilinear to measurement system
# #
# @return Transformation for offsets; 2*2 matrix # @return Transformation for offsets; 2*2 matrix
# #
def getTransCurviToMeas(self): def getTransCurviToMeas(self):
return self.__c2m #((U,V)*(J,K))^-1 return self.__c2m #((U,V)*(J,K))^-1
## Scattering precision matrix in local system ## Scattering precision matrix in local system
# #
# @return Precision diagonal; vector(2) # @return Precision diagonal; vector(2)
# #
def getScatPrecision(self, scatErr): def getScatPrecision(self, scatErr):
return np.array([1., 1.]) / (scatErr * scatErr) # diagonal only return np.array([1., 1.]) / (scatErr * scatErr) # diagonal only
# create points on initial trajectory, create trajectory from points, # create points on initial trajectory, create trajectory from points,
# fit and write trajectory to MP-II binary file # fit and write trajectory to MP-II binary file
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment