2Track fit with general broken lines.
32from .gblnum
import BorderedBandMatrix
33from .mille
import MilleRecord
55 def __init__(self, aMeasurement, minPrecision=0.):
72 if (aMeasurement[2].ndim == 2):
74 eigenVal, eigenVec = np.linalg.eigh(aMeasurement[2])
77 if (aMeasurement[0]
is None):
287 if (aScatterer[1].ndim == 2):
288 eigenVal, eigenVec = np.linalg.eigh(aScatterer[1])
289 scatTransformation = eigenVec.T
292 self.
__scatterer[1] = np.dot(scatTransformation, aScatterer[0])
321 scatPrec = np.zeros((4, 4))
322 scatPrec[2:, 2:] = np.linalg.inv(posCov)
324 eigenVal, eigenVec = np.linalg.eigh(scatPrec)
325 scatTransformation = eigenVec.T
326 return [scatTransformation, np.dot(scatTransformation, scatRes), eigenVal]
430 matJ = aJacobian[3:5, 3:5]
431 matS = aJacobian[3:5, 1:3]
432 vecD = aJacobian[3:5, 0:1]
435 matW = np.linalg.inv(matS)
436 return matW, np.dot(matW, matJ), np.dot(matW, vecD)
459 def __init__(self, aLabel=0, aType=0, aValue=0., aPrec=0., aMeas=0):
497 derLocal=None, labGlobal=None, derGlobal=None):
498 if (derLocal
is not None):
499 for i
in range(derLocal.shape[1]):
500 if (derLocal[iRow, i] != 0.):
504 for i
in range(len(labDer)):
505 if (labDer[i] != 0
and matDer[iRow, i] != 0.):
509 if (derGlobal
is not None):
510 for i
in range(derGlobal.shape[1]):
511 if (derGlobal[iRow, i] != 0.):
521 for i
in range(len(derExt)):
522 if (derExt[i] != 0.):
532 aMatrix = np.dot(aVector.T, aVector)
535 return self.
__parameters, aValue * aVector * aWeight, aMatrix * aWeight
562 if (scaledResidual < 4.6851):
563 aWeight = (1.0 - (scaledResidual / 4.6851) ** 2) ** 2
567 if (scaledResidual < 1.345):
570 aWeight = 1.345 / scaledResidual
572 aWeight = 1.0 / (1.0 + (scaledResidual / 2.3849) ** 2)
584 Chi2 = scaledResidual ** 2
586 if (scaledResidual < 4.6851):
587 Chi2 = 4.6851 ** 2 / 3. * (1. - (1. - (scaledResidual / 4.6851) ** 2) ** 3)
589 Chi2 = 4.6851 ** 2 / 3.
591 if (scaledResidual > 1.345):
592 Chi2 = 1.345 * (2.*scaledResidual - 1.345)
594 Chi2 = math.log(1. + (scaledResidual / 2.3849) ** 2) * 2.3849 ** 2
651 if (i < maxPar - maxBand):
653 return maxPar, maxBor
805 point.setLabel(label)
807 for m
in point.getMeasurements():
835 print(
"GblData blocks")
855 rec.addData(aData.toRecord())
857 rec.writeRecord(aFile)
865 rec.readRecord(aFile)
869 while (rec.moreData()):
870 aTag = rec.specialDataTag()
874 aData.fromRecord(rec.getData())
876 nPar, nBor = aData.analyzeData(mBand)
877 mPar = max(mPar, nPar)
878 mBor = max(mBor, nBor)
891 anIndex = abs(aLabel) - 1
906 nBorder = nCurv + nLocals
907 nParBRL = nBorder + 2 * nDim
908 nParLoc = nLocals + 5
909 aJacobian = np.zeros((nParLoc, nParBRL))
914 for i
in range(nLocals):
915 aJacobian[i + 5, i] = 1.0;
916 anIndex.append(i + 1);
922 anIndex.append(labDer[i]);
924 aJacobian[j, iCol] = matDer[j, i];
927 return anIndex, aJacobian
944 nOffset = aPoint.getOffset()
945 scatDim = aPoint.getScatDim()
946 anIndex = [0, 0, 0, 0, 0]
947 aJacobian = np.zeros((measDim, 5))
948 labOffset = measDim - 2
949 labSlope = measDim - 4
950 labCurv = measDim - 5
953 prevW, prevWJ, prevWd = aPoint.getDerivatives(0)
954 nextW, nextWJ, nextWd = aPoint.getDerivatives(1)
955 sumWJ = prevWJ + nextWJ
956 matN = np.linalg.inv(sumWJ)
960 prevNW = np.dot(matN, prevW)
961 nextNW = np.dot(matN, nextW)
962 prevNd = np.dot(matN, prevWd)
963 nextNd = np.dot(matN, nextWd)
964 iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1
966 aJacobian[labOffset:measDim, 0:1] = -prevNd - nextNd
967 anIndex[0] = nLocals + 1
968 aJacobian[labOffset:measDim, 1:3] = prevNW
969 aJacobian[labOffset:measDim, 3:5] = nextNW
970 for i
in range(nDim):
971 anIndex[1 + aDim[i]] = iOff + i
972 anIndex[3 + aDim[i]] = iOff + nDim + i
976 prevWPN = np.dot(nextWJ, prevNW)
977 nextWPN = np.dot(prevWJ, nextNW)
978 prevWNd = np.dot(nextWJ, prevNd)
979 nextWNd = np.dot(prevWJ, nextNd)
981 aJacobian[labSlope:labOffset, 0:1] = prevWNd - nextWNd
982 aJacobian[labSlope:labOffset, 1:3] = -prevWPN
983 aJacobian[labSlope:labOffset, 3:5] = nextWPN
992 iOff1 = nDim * nOffset + nCurv + nLocals + 1
993 index1 = 3 - 2 * nJacobian
994 iOff2 = iOff1 + nDim * (nJacobian * 2 - 1)
995 index2 = 1 + 2 * nJacobian
998 aJacobian[labOffset, index1] = 1.0
999 aJacobian[labOffset + 1, index1 + 1] = 1.0
1000 for i
in range(nDim):
1001 anIndex[index1 + aDim[i]] = iOff1 + i
1004 matW, matWJ, vecWd = aPoint.getDerivatives(nJacobian)
1005 sign = 2 * nJacobian - 1
1007 aJacobian[labSlope:labOffset, 0:1] = -sign * vecWd
1008 anIndex[0] = nLocals + 1
1009 aJacobian[labSlope:labOffset, index1:index1 + 2] = -sign * matWJ
1010 aJacobian[labSlope:labOffset, index2:index2 + 2] = sign * matW
1011 for i
in range(nDim):
1012 anIndex[index2 + aDim[i]] = iOff2 + i
1017 aJacobian[labCurv, labCurv] = 1.0
1019 return anIndex, aJacobian
1034 nOffset = aPoint.getOffset()
1035 anIndex = [0, 0, 0, 0, 0, 0, 0]
1036 aJacobian = np.zeros((2, 7))
1038 prevW, prevWJ, prevWd = aPoint.getDerivatives(0)
1039 nextW, nextWJ, nextWd = aPoint.getDerivatives(1)
1040 sumWJ = prevWJ + nextWJ
1041 sumWd = prevWd + nextWd
1042 iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1
1046 aJacobian[:, 0:1] = -sumWd
1047 anIndex[0] = nLocals + 1
1048 aJacobian[:, 1:3] = prevW
1049 aJacobian[:, 3:5] = -sumWJ
1050 aJacobian[:, 5:7] = nextW
1051 for i
in range(nDim):
1052 anIndex[1 + aDim[i]] = iOff + i
1053 anIndex[3 + aDim[i]] = iOff + nDim + i
1054 anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
1056 return anIndex, aJacobian
1073 nOffset = aPoint.getOffset()
1074 anIndex = [0, 0, 0, 0]
1075 aJacobian = np.zeros((4, 4))
1077 iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1
1080 aJacobian[2:,:2] = -np.eye(2)
1081 aJacobian[2:, 2:] = np.eye(2)
1083 for i
in range(nDim):
1084 anIndex[aDim[i]] = iOff + nDim + i
1085 anIndex[2 + aDim[i]] = iOff + nDim * 2 + i
1087 return anIndex, aJacobian
1104 nOffset = aPoint.getOffset()
1105 anIndex = [0, 0, 0, 0, 0, 0, 0, 0, 0]
1106 aJacobian = np.zeros((4, 9))
1108 prevW, prevWJ, prevWd = aPoint.getDerivatives(0)
1109 nextW, nextWJ, nextWd = aPoint.getDerivatives(1)
1110 iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1
1114 aJacobian[:2, 0:1] = -(prevWd + nextWd)
1115 anIndex[0] = nLocals + 1
1116 aJacobian[:2, 1:3] = prevW
1117 aJacobian[:2, 3:5] = -prevWJ
1118 aJacobian[:2, 5:7] = -nextWJ
1119 aJacobian[:2, 7:9] = nextW
1121 aJacobian[2:, 3:5] = -np.eye(2)
1122 aJacobian[2:, 5:7] = np.eye(2)
1124 for i
in range(nDim):
1125 anIndex[1 + aDim[i]] = iOff + i
1126 anIndex[3 + aDim[i]] = iOff + nDim + i
1127 anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
1128 anIndex[7 + aDim[i]] = iOff + nDim * 3 + i
1130 return anIndex, aJacobian
1133 aResidual, aMeasVar, aDownWeight, indLocal, derLocal = self.
__data[aData].
getResidual()
1134 aVec = np.array(derLocal)
1135 aMat = self.
__matrix.getBlockMatrix(indLocal)
1136 aFitVar = np.dot(aVec, np.dot(aMat, aVec.T))
1137 aFitVar *= aDownWeight
1138 aMeasError = math.sqrt(aMeasVar)
1140 aResError = math.sqrt(aMeasVar - aFitVar)
if aFitVar < aMeasVar
else 0.
1142 aResError = math.sqrt(aMeasVar + aFitVar)
1143 return aResidual, aMeasError, aResError, aDownWeight
1156 nParBRL = len(anIndex)
1157 aVec = np.empty(nParBRL)
1158 for i
in range(nParBRL):
1159 aVec[i] = self.
__vector[anIndex[i] - 1]
1160 aMat = self.
__matrix.getBlockMatrix(anIndex)
1161 locPar = np.dot(aJacobian, aVec)
1162 locCov = np.dot(aJacobian, np.dot(aMat, aJacobian.T))
1163 return locPar, locCov
1176 return 0,
None,
None,
None,
None
1178 aResiduals = np.empty(numData)
1179 aMeasErr = np.empty(numData)
1180 aResErr = np.empty(numData)
1181 aDownWeight = np.empty(numData)
1182 for i
in range(numData):
1184 return numData, aResiduals, aMeasErr, aResErr, aDownWeight
1197 return 0,
None,
None,
None,
None
1199 aResiduals = np.empty(numData)
1200 aMeasErr = np.empty(numData)
1201 aResErr = np.empty(numData)
1202 aDownWeight = np.empty(numData)
1203 for i
in range(numData):
1204 aResiduals[i], aMeasErr[i], aResErr[i], aDownWeight[i] = self.
__getResAndErr(firstData + i)
1205 return numData, aResiduals, aMeasErr, aResErr, aDownWeight
1220 def fit(self, optionList="", aLabel=0):
1223 def defineOffsets():
1231 scatDim = aPoint.getScatDim()
1233 aPoint.setOffset(nOffsets)
1239 aPoint.setOffset(-nOffsets)
1241 self.
__points[-1].setOffset(nOffsets)
1243 if (self.
__points[-1].getScatDim() == 4):
1251 def calcJacobians():
1252 scatJacobian = np.empty((5, 5))
1258 scatJacobian = self.
__points[iPoint].getP2pJacobian()
1260 scatJacobian = np.dot(self.
__points[iPoint].getP2pJacobian(), scatJacobian)
1262 self.
__points[iPoint].addPrevJacobian(scatJacobian)
1263 if (self.
__points[iPoint].getOffset() >= 0):
1264 self.
__points[lastPoint].addNextJacobian(scatJacobian)
1269 if (self.
__points[iPoint].getOffset() >= 0):
1270 scatJacobian = self.
__points[iPoint].getP2pJacobian()
1272 self.
__points[iPoint].addNextJacobian(scatJacobian)
1273 scatJacobian = np.dot(scatJacobian, self.
__points[iPoint].getP2pJacobian())
1281 if (aPoint.hasMeasurement()):
1282 nLabel = aPoint.getLabel()
1284 for m
in aPoint.getMeasurements():
1286 measDim = m.getMeasDim()
1287 measPrecision = m.getMeasMinPrec()
1288 localDer = m.getLocalDerivatives()
1289 globalLab = m.getGlobalLabels()
1290 globalDer = m.getGlobalDerivatives()
1291 matP, aMeas, aPrec = m.getMeasurement()
1292 nJacobian = 1
if aPoint.getOffset() < self.
__numOffsets - 1
else 0
1294 matPDer = matDer
if matP
is None else np.dot(matP, matDer)
1295 for i
in range(measDim):
1296 if (aPrec[i] > measPrecision):
1297 aData =
GblData(nLabel, 1, aMeas[i], aPrec[i], measIndex)
1298 aData.addDerivatives(i, labDer, matPDer, localDer, \
1299 globalLab, globalDer)
1300 self.
__data.append(aData)
1307 scatDim = aPoint.getScatDim()
1309 nLabel = aPoint.getLabel()
1312 if (aPoint.isLast()):
1314 matT, aMeas, aPrec = aPoint.getReducedScatterer()
1318 matT, aMeas, aPrec = aPoint.getScatterer()
1322 if (aPoint.isLast()):
1326 matT, aMeas, aPrec = aPoint.getScatterer()
1328 matTDer = matDer
if matT
is None else np.dot(matT, matDer)
1331 aData =
GblData(nLabel, 2, aMeas[i], aPrec[i])
1332 aData.addDerivatives(i, labDer, matTDer)
1333 self.
__data.append(aData)
1337 if (aPrec[i + 2] > 0.):
1338 aData =
GblData(nLabel, 2, aMeas[i + 2], aPrec[i + 2])
1339 aData.addDerivatives(i + 2, labDer, matTDer)
1340 self.
__data.append(aData)
1348 aMatrix = np.dot(eigenVec.T, aJacobian)
1349 for i
in range(len(eigenVec)):
1350 if (eigenVal[i] > 0.):
1351 externalDerivatives = []
1352 for j
in range(len(externalIndex)):
1353 externalDerivatives.append(aMatrix[i, j])
1355 aData.addExtDerivatives(externalIndex, externalDerivatives)
1356 self.
__data.append(aData)
1360 def buildLinearEquationSystem():
1364 for aData
in self.
__data:
1368 index, aVector, aMatrix = aData.getMatrices()
1369 for i
in range(len(index)):
1370 self.
__vector[ index[i] - 1 ] += aVector[0, i]
1371 self.
__matrix.addBlockMatrix(index, aMatrix)
1378 def downWeight(aMethod):
1380 for aData
in self.
__data:
1381 aLoss += (1. - aData.setDownWeighting(aMethod))
1386 for aData
in self.
__data:
1397 buildLinearEquationSystem()
1405 for o
in optionList:
1407 aMethod =
"THC".index(o.upper()) + 1
1408 lostWeight = downWeight(aMethod)
1409 buildLinearEquationSystem()
1417 for aData
in self.
__data:
1421 Chi2 += aData.getChi2()
1423 Chi2 /= [1.0, 0.8737, 0.9326, 0.8228 ][aMethod]
1424 return Chi2, Ndf, lostWeight
1426 except (ZeroDivisionError, np.linalg.linalg.LinAlgError):
Data (block) containing value, precision and derivatives for measurements, kinks and external seed.
__dwMethod
down weighting method; int
__globalDerivatives
derivatives (prediction vs global parameters); list(float)
def toRecord(self)
Get data components (for copying to MP binaty record)
__parameters
labels of fit parameters (with non zero derivative); list(int)
def setDownWeighting(self, aMethod)
Outlier down weighting with M-estimators.
def setPrediction(self, aVector)
Calculate prediction for data from fit.
def analyzeData(self, maxBand)
Analyze labels of fit parameters to determine number of parameters and border size with given maximal...
def getIndex(self)
Get index.
def addExtDerivatives(self, indexExt, derExt)
Add derivatives to data (block) from external seed.
__index
index for internal measurement
def getMatrices(self)
Calculate compressed matrix and right hand side from data.
def getType(self)
Get type.
__value
value (residual or kink); float
__downWeight
down weighting factor (M-estimators); float
def printData(self)
Print data.
def __init__(self, aLabel=0, aType=0, aValue=0., aPrec=0., aMeas=0)
Create new data.
__globalLabels
labels of global parameters; list(int)
__label
label of corresponding point; int
__derivatives
derivatives (prediction vs fit parameters); list(float)
def getLabel(self)
Get Label.
def getChi2(self)
Calculate Chi2 (contribution) from data.
def fromRecord(self, dataList)
Set data components (from MP binaty record)
__precision
precision (diagonal element of inverse covariance matrix); float
def addDerivatives(self, iRow, labDer, matDer, derLocal=None, labGlobal=None, derGlobal=None)
Add derivatives to data (block) from measurements and kinks.
__type
type of data (0: none, 1: internal measurement, 2: internal kink, 3: external seed,...
def getResidual(self)
Get data for residual (and errors).
__prediction
prediction (for value from fit); float
User supplied measurement at point on (initial) trajectory.
def getMeasDim(self)
Retrieve dimension of measurement of a point.
def addLocals(self, derivatives)
Add local derivatives.
def getMeasurement(self)
Retrieve measurement of a point.
__globalLabels
global labels; matrix(int)
def getMeasMinPrec(self)
Retrieve minimal precision to accept measurement.
__measTransformation
transformation (to eigen-vectors of precision matrix); matrix(float)
def printMeasurement(self)
Print Measurement.
__globalDerivatives
global derivatives; matrix(float)
__measMinPrec
minimal precision to accept measurement <
__measDim
dimension of measurement (2D, 4D or 5D); int
__localDerivatives
local derivatives; matrix(float)
def getNumLocals(self)
Get number of local derivatives.
def __init__(self, aMeasurement, minPrecision=0.)
Create new measurement.
def hasMeasurement(self)
Check point for a measurement.
def getLocalDerivatives(self)
Get local derivatives.
def addGlobals(self, labels, derivatives)
Add global derivatives.
__measurement
measurement at point: projection (dm/du), residuals (to initial trajectory), precision; list(matrix(f...
def getGlobalDerivatives(self)
Get global derivatives.
def getGlobalLabels(self)
Get global labels.
User supplied point on (initial) trajectory.
def hasMeasurement(self)
Check point for a measurement.
def addNextJacobian(self, aJacobian)
Add jacobian to next offset.
__measurements
measurements at point; list(GblMeasurement)
def printPoint(self)
Print point.
def getMeasurements(self)
Retrieve measurements of a point.
def setOffset(self, anOffset)
Define offset of a point and references to previous and next point with offset.
def addScatterer(self, aScatterer)
Add a (thin or thick) scatterer to a point.
def addPrevJacobian(self, aJacobian)
Add jacobian to previous offset.
def __init__(self, aJacobian)
Create new point.
def addLocals(self, derivatives)
Add local derivatives (to last measurement).
__label
label for referencing point (0,1,..,number(points)-1); int
def getLabel(self)
Retrieve label of a point.
__scatterer
scatterer at point: (transformation or None,) initial kinks, precision (inverse covariance matrix); l...
def addMeasurement(self, aMeasurement, minPrecision=0.)
Add a mesurement to a point.
def getScatDim(self)
Get scatterer dimension.
def setLabel(self, aLabel)
Define label of a point.
def getScatterer(self)
Retrieve scatterer of a point.
def addGlobals(self, labels, derivatives)
Add global derivatives (to last measurement).
def getDerivatives(self, index)
Get derivatives for locally linearized track model (backward or forward propagation).
__offset
>=0: offset number at point, <0: offset number at next point with offset; int
def isFirst(self)
Is first point?
__jacobians
jacobians for propagation to previous or next point with offsets; pair(matrix(float))
def getP2pJacobian(self)
Retrieve jacobian of a point.
__p2pJacobian
Point-to-point jacobian from previous point; matrix(float)
def isLast(self)
Is last point?
def getOffset(self)
Get offset of a point.
__scatDim
scatterer dimension (valid 2(thin) or 4(thick))
__type
type (-1: first, 0: inner, 1: last)
def getReducedScatterer(self)
Retrieve (reduced) scatterer of a point.
def setType(self, aType)
Define type of a point.
General Broken Lines Trajectory.
def __getFitToLocalJacobian(self, aPoint, measDim, nJacobian=1)
def milleIn(self, aFile)
Read (data blocks of) trajectory from MP (binary) file.
__numLocals
number of local parameters; int
__numCurvature
'curvature' is fit parameter (=1); int
def getNumPoints(self)
Get number of points on trajectory.
__numOffsets
number of (points with) offsets on trajectory; int
__dimensions
active components of offsets (both ([0,1]) or single ([0] or [1]); list(int)
__externalSeed
external seed (for local, fit parameters); matrix(float)
__measDataIndex
mapping points to data blocks from measurements; list(int)
def getResidual(self, iData)
Get residuals from data of trajectory.
def __getJacobian(self, aLabel)
Get jacobian for transformation from fit to track parameters at point.
def getMeasResults(self, aLabel)
Get residuals at point from measurement.
def fit(self, optionList="", aLabel=0)
Perform fit of trajectory.
def __getFitToKinkJacobian(self, aPoint)
Get jacobian for transformation from (trajectory) fit to kink parameters at point.
def getResults(self, aLabel)
Get results (corrections, covarinace matrix) at point in forward or backward direction.
def __getFitToStepJacobian(self, aPoint)
Get jacobian for transformation from (trajectory) fit to step parameters at point.
__scatDataIndex
mapping points to data blocks from scatterers; list(int)
__matrix
Calculate Jacobians to previous/next scatterer from point to point ones.
__numParameters
number fit parameters; int
__data
data (blocks) of trajectory; list(GblData)
__points
points on trajectory; list(GblPoint)
def addExternalSeed(self, aLabel, aSeed)
Add external seed to trajectory.
__skippedMeasLabel
label of point with measurements skipped in fit (for unbiased residuals)
def milleOut(self, aFile, doublePrec=False)
Write (data blocks of) trajectory to MP (binary) file (as float or double values).
__numPoints
number of points on trajectory; int
__externalPoint
label of point with external seed; int
def __init__(self, hasCurv=True, aDim=[0, 1])
Create new trajectory.
def printData(self)
Print data of trajectory.
def getScatResults(self, aLabel)
Get residuals at point from scatterer.
def getData(self)
Get data of trajectory.
def __getResAndErr(self, aData, used=True)
def printPoints(self)
Print points of trajectory.
def addPoint(self, point)
Add point to trajectory.
def __getFitToKinkAndStepJacobian(self, aPoint)
Get jacobian for transformation from (trajectory) fit to kink and step parameters at point.
(Symmetric) Bordered Band Matrix.
Millepede-II (binary) record.