gbltst.py 6.11 KB
 Claus Kleinwort committed Nov 02, 2021 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 ``````''' Simple Test Program for General Broken Lines. Created on Jul 27, 2011 @author: kleinwrt ''' import numpy as np import math import time from gblfit import GblPoint, GblTrajectory # def example1(): ''' Create points on initial trajectory, create trajectory from points, fit and write trajectory to MP-II binary file, get track parameter corrections and covariance matrix at points. Equidistant measurement layers and thin scatterers, propagation with simple jacobian (quadratic in arc length differences). `````` Claus Kleinwort committed Nov 02, 2021 22 `````` Curvilinear system (U,V,T) as local coordinate system. `````` Claus Kleinwort committed Nov 02, 2021 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 `````` ''' def gblSimpleJacobian(ds, cosl, bfac): ''' Simple jacobian: quadratic in arc length difference. @param ds: arc length difference @type ds: float @param cosl: cos(lambda) @type cosl: float @param bfac: Bz*c @type bfac: float @return: jacobian to move by 'ds' on trajectory @rtype: matrix(float) ''' jac = np.eye(5) jac[1, 0] = -bfac * ds * cosl jac[3, 0] = -0.5 * bfac * ds * ds * cosl jac[3, 1] = ds jac[4, 2] = ds return jac # nTry = 10#00 #: number of tries nLayer = 5 #: number of detector layers `````` Claus Kleinwort committed Nov 02, 2021 46 `````` print " Gbltst \$Rev\$ ", nTry, nLayer `````` Claus Kleinwort committed Nov 02, 2021 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 `````` start = time.clock() # track direction sinLambda = 0.3 cosLambda = math.sqrt(1.0 - sinLambda ** 2) sinPhi = 0. cosPhi = math.sqrt(1.0 - sinPhi ** 2) # tDir = np.array([cosLambda * cosPhi, cosLambda * sinPhi, sinLambda]) # U = Z x T / |Z x T|, V = T x U uvDir = np.array([[-sinPhi, cosPhi, 0.], \ [-sinLambda * cosPhi, -sinLambda * sinPhi, cosLambda]]) # measurement resolution measErr = np.array([ 0.001, 0.001]) # 10 mu measPrec = 1.0 / measErr ** 2 # scattering error scatErr = np.array([ 0.001, 0.001]) # 1 mread scatPrec = 1.0 / scatErr ** 2 `````` Claus Kleinwort committed Nov 02, 2021 63 ``````# RMS of CurviLinear track parameters (Q/P, slopes, offsets) `````` Claus Kleinwort committed Nov 02, 2021 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 `````` clErr = np.array([0.001, -0.1, 0.2, -0.15, 0.25]) clSeed = np.eye(5) for i in range(5): clSeed[i, i] = 1.0 / clErr[i] ** 2 # bfac = 0.2998 # Bz*c for Bz=1 step = 1.5 / cosLambda # constant steps in RPhi # Chi2Sum = 0. NdfSum = 0 LostSum = 0. # binaryFile = open("milleBinaryISN.dat", "wb") # for iTry in range(nTry): `````` Claus Kleinwort committed Nov 02, 2021 79 ``````# generate (CurviLinear) track parameters `````` Claus Kleinwort committed Nov 02, 2021 80 81 82 83 84 `````` clNorm = np.random.normal(0., 1., 5) clPar = clErr * clNorm # arclength s = 0. sPoint = [] `````` Claus Kleinwort committed Nov 02, 2021 85 86 ``````# point-to-point jacobian (from previous point) jacPointToPoint = np.eye(5) `````` Claus Kleinwort committed Nov 02, 2021 87 88 ``````# additional (local or global) derivatives addDer = np.array([[1.0], [0.0]]) `````` Claus Kleinwort committed Nov 02, 2021 89 `````` labGlobal = np.array([[4711], [4711]]) `````` Claus Kleinwort committed Nov 02, 2021 90 91 92 93 94 95 96 97 ``````# create trajectory traj = GblTrajectory(bfac != 0.) for iLayer in range(nLayer): # measurement directions sinStereo = (0. if iLayer % 2 == 0 else 0.5) cosStereo = math.sqrt(1.0 - sinStereo ** 2) mDir = np.array([[sinStereo, cosStereo, 0.0], [0., 0, 1.]]) `````` Claus Kleinwort committed Nov 02, 2021 98 99 ``````# projection measurement to local (curvilinear uv) directions (duv/dm) proM2l = np.dot(uvDir, mDir.T) `````` Claus Kleinwort committed Nov 02, 2021 100 ``````# projection local (uv) to measurement directions (dm/duv) `````` Claus Kleinwort committed Nov 02, 2021 101 `````` proL2m = np.linalg.inv(proM2l) `````` Claus Kleinwort committed Nov 02, 2021 102 103 104 105 ``````# measurement - prediction in measurement system with error measNorm = np.random.normal(0., 1., 2) meas = np.dot(proL2m, clPar[3:5]) + measErr * measNorm # point with measurement `````` Claus Kleinwort committed Nov 02, 2021 106 `````` point = GblPoint(jacPointToPoint) `````` Claus Kleinwort committed Nov 02, 2021 107 108 `````` point.addMeasurement([proL2m, meas, measPrec]) # additional local parameters? `````` Claus Kleinwort committed Nov 02, 2021 109 ``````# point.addLocals(addDer) `````` Claus Kleinwort committed Nov 02, 2021 110 ``````# additional global parameters? `````` Claus Kleinwort committed Nov 02, 2021 111 `````` point.addGlobals(labGlobal, addDer) `````` Claus Kleinwort committed Nov 02, 2021 112 113 114 115 116 `````` addDer = -addDer # locDer flips sign every measurement # add point to trajectory iLabel = traj.addPoint(point) sPoint.append(s) # propagate to scatterer `````` Claus Kleinwort committed Nov 02, 2021 117 118 `````` jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac) clPar = np.dot(jacPointToPoint, clPar) `````` Claus Kleinwort committed Nov 02, 2021 119 120 121 122 `````` s += step if (iLayer < nLayer - 1): scat = np.array([0., 0.]) # point with scatterer `````` Claus Kleinwort committed Nov 02, 2021 123 `````` point = GblPoint(jacPointToPoint) `````` Claus Kleinwort committed Nov 02, 2021 124 125 126 127 128 129 130 `````` point.addScatterer([scat, scatPrec]) iLabel = traj.addPoint(point) sPoint.append(s) # scatter a little scatNorm = np.random.normal(0., 1., 2) clPar[1:3] = clPar[1:3] + scatErr * scatNorm # propagate to next measurement layer `````` Claus Kleinwort committed Nov 02, 2021 131 `````` clPar = np.dot(jacPointToPoint, clPar) `````` Claus Kleinwort committed Nov 02, 2021 132 133 134 `````` s += step # add external seed `````` Claus Kleinwort committed Nov 02, 2021 135 `````` traj.addExternalSeed(1, clSeed) `````` Claus Kleinwort committed Nov 02, 2021 136 137 138 139 ``````# dump trajectory # traj.dump() # fit trajectory `````` Claus Kleinwort committed Nov 02, 2021 140 `````` Chi2, Ndf, Lost = traj.fit() `````` Claus Kleinwort committed Nov 02, 2021 141 142 143 144 145 146 147 148 `````` print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost # write to MP binary file traj.milleOut(binaryFile) # sum up Chi2Sum += Chi2 NdfSum += Ndf LostSum += Lost # get corrections and covariance matrix at points `````` Claus Kleinwort committed Nov 02, 2021 149 150 151 152 153 154 155 156 157 158 `````` if (iTry == 0): for i in range(1, 1): locPar, locCov = traj.getResults(i) print " Point< ", i print " locPar ", locPar print " locCov ", locCov traj.getResults(-i) print " Point> ", i print " locPar ", locPar print " locCov ", locCov `````` Claus Kleinwort committed Nov 02, 2021 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 ``````# end = time.clock() print " Time [s] ", end - start print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum print " LostSum/nTry ", LostSum / nTry # def example2(): ''' Read trajectory from MP-II binary file and refit. ''' # binaryFile = open("milleBinaryISN.dat", "rb") nRec = 0 maxRec = 10 #: maximum number of records to read Chi2Sum = 0. NdfSum = 0 LostSum = 0. start = time.clock() try: while(nRec < maxRec): # create trajectory traj = GblTrajectory(0) # read from file `````` Claus Kleinwort committed Nov 02, 2021 183 `````` traj.milleIn(binaryFile) # get data blocks from file `````` Claus Kleinwort committed Nov 02, 2021 184 185 `````` nRec += 1 # fit trajectory `````` Claus Kleinwort committed Nov 02, 2021 186 `````` Chi2, Ndf, Lost = traj.fit() `````` Claus Kleinwort committed Nov 02, 2021 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 `````` print " Record, Chi2, Ndf, Lost", nRec, Chi2, Ndf, Lost # sum up Chi2Sum += Chi2 NdfSum += Ndf LostSum += Lost except EOFError: pass print " records read ", nRec end = time.clock() print " Time [s] ", end - start print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum print " LostSum/nTry ", LostSum / nRec # create points on initial trajectory, create trajectory from points, # fit and write trajectory to MP-II binary file # get track parameter corrections and covariance matrix at points example1() # read trajectory from MP-II binary file and refit example2()``````