GeneralBrokenLines V03-01-02
gblpy
gblsit.py
Go to the documentation of this file.
1'''
2Created on 28 Sep 2018
3
4@author: kleinwrt
5'''
6
7
27
28import numpy as np
29import math
30import random
31import time
32from gblfit import GblPoint, GblTrajectory
33
34
35
94
95 #
96 random.seed(47117)
97
98 # magnetic field
99 bfac = 0.003 # B*c for 1 T
100 #bfac = 0. # B*c for 0 T
101 # detector layers: name, position (x,y,z), thickness (X/X_0), (1 or 2) measurements (direction in YZ, resolution), spacing for composite layers
102 det = gblSiliconDet([ ['PIX1', (2.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ], # pixel
103 ['PIX2', (3.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ], # pixel
104 ['PIX3', (4.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ], # pixel
105 ['S2D4', (6.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)], 0.1 ], # strip 2D (composite), 5 deg stereo angle
106 ['S2D5', (8.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)], 0.1 ], # strip 2D (composite), -5 deg stereo angle
107 ['S2D6', (10.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)] ], # strip 2D (double sided), 5 deg stereo angle
108 ['S2D7', (12.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)] ], # strip 2D (double sided), -5 deg stereo angle
109 ['S1D8', (15.0, 0., 0.), 0.0033, [(0., 0.0040)] ] # strip 1D
110 ], bfac)
111
112 # Alignment with MillePede-II requires for alignables with a single 1D measurements to fix the unmeasured
113 # direction with a (linear equality) constraint (unless alignment equal to measurement system).
114 det.getMP2Constraints()
115
116 nTry = 1000 #: number of tries
117 qbyp = 0.2 # 5 GeV
118 binaryFile = open("milleBinary.dat", "wb")
119 #binaryFile = None
120 #
121 print " Gblsit $Id$ ", nTry
122 #
123 start = time.clock()
124 Chi2Sum = 0.
125 NdfSum = 0
126 LostSum = 0.
127 #
128 for iTry in range(nTry):
129 # helix parameter for track generation
130 dca = random.gauss(0., 0.1)
131 z0 = random.gauss(0., 0.1)
132 phi0 = 0.5 * random.uniform(-1., 1.)
133 dzds = 0.3 * random.uniform(-1., 1.)
134 curv = bfac * qbyp * math.sqrt(1. + dzds * dzds)
135 genPar = [curv, phi0, dca, dzds, z0]
136 #print " genPar ", iTry, genPar
137 # generate hits
138 genHits = det.generateHits(qbyp, genPar)
139 # seed (with true parameters)
140 seedPar = genPar
141 seed = gblSimpleHelix(seedPar)
142 sOld = 0.
143 cosLambda = 1. / math.sqrt(1. + seedPar[3] ** 2)
144 # construct GBL trajectory
145 traj = GblTrajectory(bfac != 0.)
146 # add GBL points
147 for l, layer in enumerate(det.getLayers()):
148 # prediction from seeding helix
149 pred = layer.intersectWithHelix(seed)
150 measPred = pred.getMeasPred()
151 sArc = pred.getArcLength()
152 # residuals
153 res = np.array([genHits[l][0] - measPred[0], genHits[l][1] - measPred[1]])
154 # measurement precision
155 measPrec = np.array(layer.getPrecision())
156 # Curvilinear system: track direction T, U = Z x T / |Z x T|, V = T x U
157 # as local system
158 curviDirs = pred.getCurvilinearDirs()
159 # projection matrix (local to measurement)
160 proL2m = np.linalg.inv(np.dot(curviDirs, np.linalg.inv(layer.getMeasSystemDirs())[:,:2]))
161 # propagation
162 jacPointToPoint = gblSimpleJacobian((sArc - sOld) / cosLambda, cosLambda, bfac)
163 sOld = sArc
164 # point with (independent) measurements (in measurement system)
165 point = GblPoint(jacPointToPoint)
166 # composite?
167 if layer.isComposite():
168 # 2nd prediction
169 pred = layer.intersectWithHelix2(seed)
170 measPred = pred.getMeasPred()
171 # 4D measurement
172 pro4D = np.zeros((4, 4)); pro4D[2:, 2:] = proL2m; pro4D[3,:2] = proL2m[1,:] * layer.getSpacing();
173 res4D = np.array([0., 0., res[0], genHits[l][1] - measPred[1]])
174 prec4D = np.array([0., 0., measPrec[0], measPrec[1]])
175 point.addMeasurement([pro4D, res4D, prec4D])
176 else:
177 # 2D measurement
178 point.addMeasurement([proL2m, res, measPrec])
179 # global parameters for rigid body alignment?
180 if binaryFile is not None:
181 # local (alignment) system, per layer
182 labels = [l * 10 + 1, l * 10 + 2, l * 10 + 3, l * 10 + 4, l * 10 + 5, l * 10 + 6]
183 labGlobal = np.array([labels, labels])
184 derGlobal = layer.getRigidBodyDerLocal(pred.getPos(), pred.getDirection())
185 # composite?
186 if layer.isComposite():
187 # 4D
188 labG4D = np.array([labels, labels, labels, labels])
189 derG4D = np.zeros((4, 6)); derG4D[2:] = derGlobal
190 point.addGlobals(labG4D, derG4D)
191 else:
192 # 2D
193 point.addGlobals(labGlobal, derGlobal)
194 # add scatterer to point
195 radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
196 scatErr = gblMultipleScatteringError(qbyp, radlen) # simple model
197 if scatErr > 0.:
198 scat = np.array([0., 0.])
199 scatP = np.array([1. / scatErr ** 2, 1. / scatErr ** 2])
200 # composite?
201 if layer.isComposite():
202 # two similar sub layers
203 scatP *= 0.5
204 point.addScatterer([scat, scatP])
205 # add point to trajectory
206 traj.addPoint(point)
207
208 # fit trajectory
209 Chi2, Ndf, Lost = traj.fit()
210 print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost
211 # sum up
212 Chi2Sum += Chi2
213 NdfSum += Ndf
214 LostSum += Lost
215 # write to binary file
216 if binaryFile is not None:
217 traj.milleOut(binaryFile)
218
219 end = time.clock()
220 print " Time [s] ", end - start
221 print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum
222 print " LostSum/nTry ", LostSum / nTry
223
224
225
235def gblSimpleJacobian(ds, cosl, bfac):
236 jac = np.eye(5)
237 jac[1, 0] = -bfac * ds * cosl
238 jac[3, 0] = -0.5 * bfac * ds * ds * cosl
239 jac[3, 1] = ds
240 jac[4, 2] = ds
241 return jac
242
243
244
252 return 0.015 * abs(qbyp) * math.sqrt(xbyx0)
253
254
255
256class gblSiliconLayer(object):
257
258
262 def __init__(self, layer):
263
264 self.__name = layer[0]
265
266 self.__center = np.array(layer[1])
267
268 self.__xbyx0 = layer[2]
269 # measurements (1D or 2D)
270 meas = layer[3]
271
272 self.__resolution = (meas[0][1], meas[1][1] if len(meas) > 1 else 0.)
273
274 self.__precision = (1. / meas[0][1] ** 2, 1. / meas[1][1] ** 2 if len(meas) > 1 else 0.)
275 # measurement angles
276 uPhi = meas[0][0] / 180. * math.pi
277 vPhi = meas[1][0] / 180. * math.pi if len(meas) > 1 else uPhi + 0.5 * math.pi
278
279 self.__uDir = np.array([0., math.cos(uPhi), math.sin(uPhi)])
280
281 self.__vDir = np.array([0., math.cos(vPhi), math.sin(vPhi)])
282
283 self.__nDir = np.array([1., 0., 0.])
284
285 self.__measDirs = np.array([self.__uDir, self.__vDir, self.__nDir])
286
287 self.__ijkDirs = np.array([[0., 1., 0.], [0., 0., 1.], [1., 0., 0.]])
288
289 self.__alignInMeasSys = np.array_equal(self.__measDirs, self.__ijkDirs)
290
291 self.__spacing = layer[4] if len(layer) > 4 else None
292
293
301 def getMP2Constraint(self, layer):
302 if self.__resolution[1] > 0. or self.__alignInMeasSys:
303 return
304 # transform vDir into alignment system
305 unMeasured = np.dot(self.__ijkDirs, self.__vDir)
306 print "Constraint 0. ! fix unmeasured direction in", self.__name
307 for i in range(3):
308 # 'zero' suppression
309 if abs(unMeasured[i]) > 1.0e-10:
310 print " ", layer * 10 + i + 1, unMeasured[i]
311
312
314 return self.__xbyx0
315
316
317 def getResolution(self):
318 return self.__resolution
319
320
321 def getPrecision(self):
322 return self.__precision
323
324
326 return self.__measDirs
327
328
333 def intersectWithHelix(self, helix):
334 return helix.getPrediction(self.__center, self.__uDir, self.__vDir)
335
336
341 def intersectWithHelix2(self, helix):
342 return helix.getPrediction(self.__center + self.__spacing * self.__nDir, self.__uDir, self.__vDir)
343
344
345 def isComposite(self):
346 return self.__spacing is not None
347
348
349 def getSpacing(self):
350 return self.__spacing
351
352
358 def getRigidBodyDerGlobal(self, position, trackDir):
359 # lever arms (for rotations)
360 dist = position
361 # dr/dm (residual vs measurement, 1-tdir*ndir^t/tdir*ndir)
362 drdm = np.eye(3) - np.outer(trackDir, self.__nDir) / np.dot(trackDir, self.__nDir)
363 # dm/dg (measurement vs 6 rigid body parameters)
364 dmdg = np.zeros((3, 6))
365 dmdg[0][0] = 1.; dmdg[0][4] = -dist[2]; dmdg[0][5] = dist[1]
366 dmdg[1][1] = 1.; dmdg[1][3] = dist[2]; dmdg[1][5] = -dist[0]
367 dmdg[2][2] = 1.; dmdg[2][3] = -dist[1]; dmdg[2][4] = dist[0]
368 # drl/drg (local vs global residuals)
369 drldrg = self.__measDirs
370 # drl/dg (local residuals vs rigid body parameters)
371 drldg = np.dot(drldrg, np.dot(drdm, dmdg))
372 return drldg
373
374
380 def getRigidBodyDerLocal(self, position, trackDir):
381 # track direction in local system
382 tLoc = np.dot(self.__ijkDirs, trackDir)
383 # local slopes
384 uSlope = tLoc[0] / tLoc[2]
385 vSlope = tLoc[1] / tLoc[2]
386 # (u,v) lever arms
387 uPos, vPos = np.dot(self.__ijkDirs, position - self.__center)[:2]
388 # wPos = 0 (in detector plane)
389 # drl/dg (local residuals vs rigid body parameters)
390 drldg = np.array([[1.0, 0.0, -uSlope, vPos * uSlope, -uPos * uSlope, vPos], \
391 [0.0, 1.0, -vSlope, vPos * vSlope, -uPos * vSlope, -uPos]])
392 # avoid numerics in case of unit transformation (below)
393 if self.__alignInMeasSys:
394 return drldg
395 # local (alignment) to measurement system
396 local2meas = np.dot(self.__measDirs, self.__ijkDirs.T)
397 return np.dot(local2meas[:2,:2], drldg)
398
399
400
401class gblSiliconDet(object):
402
403
408 def __init__(self, layers, bfac):
409
410 self.__layers = []
411
412 self.__bfac = bfac
413
414 for layer in layers:
415 self.__layers.append(gblSiliconLayer(layer))
416
417
418 def getLayers(self):
419 return self.__layers
420
421
423 print "! MillePede-II: constraints for alignables with SINGLE 1D measurements"
424 for l, layer in enumerate(self.__layers):
425 layer.getMP2Constraint(l)
426 print "! End of lines to be added to MillePede-II steering file"
427
428
434 def generateHits(self, qbyp, genPar):
435
436# list of hits
437 hits = []
438 localPar = genPar
439 # print " track ", helix
440 for layer in self.__layers:
441 # local constant (Bfield) helix
442 hlx = gblSimpleHelix(localPar)
443 # prediction from local helix
444 pred = layer.intersectWithHelix(hlx)
445 meas = [pred.getMeasPred()]
446 # scatter at intersection point
447 xPos, yPos = pred.getPos()[:2]
448 radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
449 errMs = gblMultipleScatteringError(qbyp, radlen) # simple model
450 cosLambda = 1. / math.sqrt(1. + localPar[3] ** 2)
451 # move to intersection point
452 newpar = hlx.moveTo((xPos, yPos))
453 newpar[1] += random.gauss(0., errMs / cosLambda) # phi
454 newpar[3] += random.gauss(0., errMs / cosLambda ** 2) # dzds
455 newhlx = gblSimpleHelix(newpar)
456 # move back
457 localPar = newhlx.moveTo((-xPos, -yPos))
458 # composite layer
459 if layer.isComposite():
460 # 2nd prediction from local helix
461 pred = layer.intersectWithHelix2(hlx)
462 meas.append(pred.getMeasPred())
463 # scatter at intersection point
464 xPos, yPos = pred.getPos()[:2]
465 cosLambda = 1. / math.sqrt(1. + localPar[3] ** 2)
466 # move to intersection point
467 newpar = hlx.moveTo((xPos, yPos))
468 newpar[1] += random.gauss(0., errMs / cosLambda) # phi
469 newpar[3] += random.gauss(0., errMs / cosLambda ** 2) # dzds
470 newhlx = gblSimpleHelix(newpar)
471 # move back
472 localPar = newhlx.moveTo((-xPos, -yPos))
473
474 # add (smeared) hit
475 sigma = layer.getResolution()
476 hits.append((random.gauss(meas[0][0], sigma[0]), random.gauss(meas[-1][1], sigma[1])))
477
478 return hits
479
480
481
485class gblSimpleHelix(object):
486
487
495 def __init__(self, parameter):
496
497 self.__rinv = parameter[0]
498
499 self.__phi0 = parameter[1]
500
501 self.__dir0 = (math.cos(self.__phi0), math.sin(self.__phi0))
502
503 self.__dca = parameter[2]
504
505 self.__dzds = parameter[3]
506
507 self.__z0 = parameter[4]
508
509 self.__xRelCenter = -(1. - self.__dca * self.__rinv) * self.__dir0[1]
510
511 self.__yRelCenter = (1. - self.__dca * self.__rinv) * self.__dir0[0]
512
513
522 def getPrediction(self, refPos, uDir, vDir):
523 # normal to (u,v) measurement plane
524 nDir = np.cross(uDir, vDir); nDir /= np.linalg.norm(nDir)
525 # ZS direction
526 cosLambda = 1. / math.sqrt(1. + self.__dzds * self.__dzds)
527 sinLambda = self.__dzds * cosLambda
528 # line (or helix)
529 if self.__rinv == 0.:
530 # track direction
531 tDir = np.array([cosLambda * self.__dir0[0], cosLambda * self.__dir0[1], sinLambda])
532 # distance (of point at dca to reference)
533 pca = np.array([ self.__dca * self.__dir0[1] , -self.__dca * self.__dir0[0], self.__z0])
534 dist = pca - refPos
535 # arc-length
536 sArc3D = -np.dot(dist, nDir) / np.dot(tDir, nDir); sArc2D = sArc3D * cosLambda
537 # distance (of point at sArc to reference)
538 pos = pca + sArc3D * tDir
539 dist = pos - refPos
540 else:
541 # initial guess of 2D arc-length
542 sArc2D = self.getArcLengthXY(refPos[0], refPos[1])
543 nIter = 0
544 while nIter < 10:
545 nIter += 1
546 # track direction
547 dPhi = sArc2D * self.__rinv
548 cosPhi = math.cos(self.__phi0 + dPhi); sinPhi = math.sin(self.__phi0 + dPhi)
549 tDir = np.array([cosLambda * cosPhi, cosLambda * sinPhi, sinLambda])
550 # distance (of point at sArc to reference)
551 pos = np.array([(self.__xRelCenter + sinPhi) / self.__rinv, (self.__yRelCenter - cosPhi) / self.__rinv, self.__z0 + self.__dzds * sArc2D])
552 dist = pos - refPos
553 # arc-length correction (linearizing helix at sNew)
554 sCorr3D = -np.dot(dist, nDir) / np.dot(tDir, nDir)
555 if abs(sCorr3D) > 0.0001:
556 sArc2D += sCorr3D * cosLambda
557 else:
558 break
559
560 # prediction in measurement directions
561 pred = [np.dot(dist, uDir), np.dot(dist, vDir)]
562 return gblHelixPrediction(sArc2D, pred, tDir, uDir, vDir, nDir, pos)
563
564
573 def getArcLengthXY(self, xPos, yPos):
574 # line
575 if self.__rinv == 0:
576 return self.__dir0[0] * xPos + self.__dir0[1] * yPos
577 # helix
578 dx = (xPos * self.__rinv - self.__xRelCenter)
579 dy = (yPos * self.__rinv - self.__yRelCenter)
580 dphi = math.atan2(dx, -dy) - self.__phi0
581 if (abs(dphi) > math.pi):
582 dphi -= cmp(dphi, 0.) * 2.0 * math.pi
583 return dphi / self.__rinv
584
585
591 def moveTo(self, newRefPoint):
592
593 rho = self.__rinv
594 phi = self.__phi0
595 dca = self.__dca
596 dzds = self.__dzds
597 z0 = self.__z0
598
599 u = 1. - rho * dca
600 dp = -newRefPoint[0] * self.__dir0[1] + newRefPoint[1] * self.__dir0[0] + dca
601 dl = newRefPoint[0] * self.__dir0[0] + newRefPoint[1] * self.__dir0[1]
602 sa = 2. * dp - rho * (dp * dp + dl * dl)
603 sb = rho * newRefPoint[0] + u * self.__dir0[1]
604 sc = -rho * newRefPoint[1] + u * self.__dir0[0]
605 sd = math.sqrt(1. - rho * sa)
606 # transformed parameters
607 if rho == 0.:
608 dca = dp
609 sArc = dl
610 newPar = [rho, phi, dca]
611 else:
612 phi = math.atan2(sb, sc)
613 dca = sa / (1. + sd)
614 dphi = phi - self.__phi0
615 if abs(dphi) > math.pi: dphi -= cmp(dphi, 0.) * 2.0 * math.pi
616 sArc = dphi / rho
617 newPar = [rho, phi, dca]
618 z0 += sArc * dzds
619 newPar += [dzds, z0]
620
621 return newPar
622
623
624
626class gblHelixPrediction(object):
627
628
638 def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos):
639
640 self.__sarc = sArc
641
642 self.__pred = pred
643
644 self.__tdir = tDir
645
646 self.__udir = uDir
647
648 self.__vdir = vDir
649
650 self.__ndir = nDir
651
652 self.__pos = pos
653 #
654
655
656 def getArcLength(self):
657 return self.__sarc
658
659
660 def getMeasPred(self):
661 return self.__pred
662
663
664 def getPos(self):
665 return self.__pos
666
667
668 def getDirection(self):
669 return self.__tdir
670
671
673 return np.dot(self.__tdir, self.__ndir)
674
675
680 cosTheta = self.__tdir[2]; sinTheta = math.sqrt(self.__tdir[0] ** 2 + self.__tdir[1] ** 2)
681 cosPhi = self.__tdir[0] / sinTheta; sinPhi = self.__tdir[1] / sinTheta
682 return np.array([[-sinPhi, cosPhi, 0.], [-cosPhi * cosTheta, -sinPhi * cosTheta, sinTheta]])
683
684
685if __name__ == '__main__':
686 exampleSit()
Prediction (from helix at measurement)
Definition: gblsit.py:626
__ndir
normal to (u,v)
Definition: gblsit.py:650
def getCosIncidence(self)
Get cosine of incidence.
Definition: gblsit.py:672
def getArcLength(self)
Get arc-length
Definition: gblsit.py:656
__tdir
track direction
Definition: gblsit.py:644
def getCurvilinearDirs(self)
Get curvilinear directions (U,V)
Definition: gblsit.py:679
def getDirection(self)
Get (track) direction.
Definition: gblsit.py:668
def getPos(self)
Get Position.
Definition: gblsit.py:664
def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos)
Constructor.
Definition: gblsit.py:638
def getMeasPred(self)
Get measurement prediction.
Definition: gblsit.py:660
Silicon detector.
Definition: gblsit.py:401
def generateHits(self, qbyp, genPar)
Generate hits on helix.
Definition: gblsit.py:434
def __init__(self, layers, bfac)
Constructor.
Definition: gblsit.py:408
def getMP2Constraints(self)
get MP2 constraints
Definition: gblsit.py:422
def getLayers(self)
Get layers.
Definition: gblsit.py:418
Silicon layer.
Definition: gblsit.py:256
__vDir
measurement direction v
Definition: gblsit.py:281
__xbyx0
radiation length
Definition: gblsit.py:268
def getSpacing(self)
Get spacing.
Definition: gblsit.py:349
__resolution
resolution (for simulation)
Definition: gblsit.py:272
__measDirs
measurement directions
Definition: gblsit.py:285
def getRigidBodyDerGlobal(self, position, trackDir)
Get rigid body derivatives in global frame.
Definition: gblsit.py:358
__spacing
spacing (for composite layers)
Definition: gblsit.py:291
def intersectWithHelix2(self, helix)
Intersect with helix (2nd sub layer)
Definition: gblsit.py:341
def getRigidBodyDerLocal(self, position, trackDir)
Get rigid body derivatives in local (alignment) frame.
Definition: gblsit.py:380
def isComposite(self)
Is composite?
Definition: gblsit.py:345
__uDir
measurement direction u
Definition: gblsit.py:279
def intersectWithHelix(self, helix)
Intersect with helix.
Definition: gblsit.py:333
__nDir
normal to measurement plane
Definition: gblsit.py:283
def getMP2Constraint(self, layer)
get MP2 constraint
Definition: gblsit.py:301
__alignInMeasSys
alignment == measurement system?
Definition: gblsit.py:289
def __init__(self, layer)
Constructor.
Definition: gblsit.py:262
__ijkDirs
local alignment system (IJK = YZX)
Definition: gblsit.py:287
def getResolution(self)
Get resolution.
Definition: gblsit.py:317
def getPrecision(self)
Get precision.
Definition: gblsit.py:321
def getRadiationLength(self)
Get radiation length.
Definition: gblsit.py:313
__precision
precision (for reconstruction)
Definition: gblsit.py:274
def getMeasSystemDirs(self)
Get directions of measurement system.
Definition: gblsit.py:325
Simple helix.
Definition: gblsit.py:485
__dca
distance of closest approach in (XY)
Definition: gblsit.py:503
def getPrediction(self, refPos, uDir, vDir)
Get prediction.
Definition: gblsit.py:522
__yRelCenter
XY circle parameter: Y position of center / R.
Definition: gblsit.py:511
def moveTo(self, newRefPoint)
Change reference point.
Definition: gblsit.py:591
def getArcLengthXY(self, xPos, yPos)
Get (2D) arc length for given point.
Definition: gblsit.py:573
__z0
Z position at distance of closest approach.
Definition: gblsit.py:507
__xRelCenter
XY circle parameter: X position of center / R.
Definition: gblsit.py:509
__rinv
curvature (in XY)
Definition: gblsit.py:497
def __init__(self, parameter)
Constructor.
Definition: gblsit.py:495
__dir0
direction vector at point of closest approach (in XY)
Definition: gblsit.py:501
__phi0
flight direction at point of closest approach (in XY)
Definition: gblsit.py:499
def gblMultipleScatteringError(qbyp, xbyx0)
Multiple scattering error.
Definition: gblsit.py:251
def exampleSit()
Simulate and reconstruct helical tracks in silicon pixel and (1D or 2D) strip detectors.
Definition: gblsit.py:93
def gblSimpleJacobian(ds, cosl, bfac)
Simple jacobian.
Definition: gblsit.py:235