90 det =
gblSiliconDet([ [
'PIX1', (2.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],
91 [
'PIX2', (3.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],
92 [
'PIX3', (4.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],
93 [
'S2D4', (6.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)], 0.1 ],
94 [
'S2D5', (8.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)], 0.1 ],
95 [
'S2D6', (10.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)] ],
96 [
'S2D7', (12.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)] ],
97 [
'S1D8', (15.0, 0., 0.), 0.0033, [(0., 0.0040)] ]
102 binaryFile = open(
"milleBinary.dat",
"wb")
105 print(
" Gblsit $Id$ ", nTry)
107 start = time.process_time()
112 for iTry
in range(nTry):
114 dca = random.gauss(0., 0.1)
115 z0 = random.gauss(0., 0.1)
116 phi0 = 0.5 * random.uniform(-1., 1.)
117 dzds = 0.3 * random.uniform(-1., 1.)
118 curv = bfac * qbyp * math.sqrt(1. + dzds * dzds)
119 genPar = [curv, phi0, dca, dzds, z0]
122 genHits = det.generateHits(qbyp, genPar)
127 cosLambda = 1. / math.sqrt(1. + seedPar[3] ** 2)
131 for l, layer
in enumerate(det.getLayers()):
133 pred = layer.intersectWithHelix(seed)
134 measPred = pred.getMeasPred()
135 sArc = pred.getArcLength()
137 res = np.array([genHits[l][0] - measPred[0], genHits[l][1] - measPred[1]])
139 measPrec = np.array(layer.getPrecision())
142 curviDirs = pred.getCurvilinearDirs()
144 proL2m = np.linalg.inv(np.dot(curviDirs, np.linalg.inv(layer.getMeasSystemDirs())[:,:2]))
151 if layer.isComposite():
153 pred = layer.intersectWithHelix2(seed)
154 measPred = pred.getMeasPred()
156 pro4D = np.zeros((4, 4)); pro4D[2:, 2:] = proL2m; pro4D[3,:2] = proL2m[1,:] * layer.getSpacing();
157 res4D = np.array([0., 0., res[0], genHits[l][1] - measPred[1]])
158 prec4D = np.array([0., 0., measPrec[0], measPrec[1]])
159 point.addMeasurement([pro4D, res4D, prec4D])
162 point.addMeasurement([proL2m, res, measPrec])
164 if binaryFile
is not None:
166 labels = [l * 10 + 1, l * 10 + 2, l * 10 + 3, l * 10 + 4, l * 10 + 5, l * 10 + 6]
167 labGlobal = np.array([labels, labels])
168 derGlobal = layer.getRigidBodyDerLocal(pred.getPos(), pred.getDirection())
170 if layer.isComposite():
172 labG4D = np.array([labels, labels, labels, labels])
173 derG4D = np.zeros((4, 6)); derG4D[2:] = derGlobal
174 point.addGlobals(labG4D, derG4D)
177 point.addGlobals(labGlobal, derGlobal)
179 radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
182 scat = np.array([0., 0.])
183 scatP = np.array([1. / scatErr ** 2, 1. / scatErr ** 2])
185 if layer.isComposite():
188 point.addScatterer([scat, scatP])
193 Chi2, Ndf, Lost = traj.fit()
194 print(
" Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost)
200 if binaryFile
is not None:
201 traj.milleOut(binaryFile)
203 end = time.process_time()
204 print(
" Time [s] ", end - start)
205 print(
" Chi2Sum/NdfSum ", Chi2Sum / NdfSum)
206 print(
" LostSum/nTry ", LostSum / nTry)
221 jac[1, 0] = -bfac * ds * cosl
222 jac[3, 0] = -0.5 * bfac * ds * ds * cosl
236 return 0.015 * abs(qbyp) * math.sqrt(xbyx0)
256 self.
__precision = (1. / meas[0][1] ** 2, 1. / meas[1][1] ** 2
if len(meas) > 1
else 0.)
258 uPhi = meas[0][0] / 180. * math.pi
259 vPhi = meas[1][0] / 180. * math.pi
if len(meas) > 1
else uPhi + 0.5 * math.pi
261 self.
__uDir = np.array([0., math.cos(uPhi), math.sin(uPhi)])
263 self.
__vDir = np.array([0., math.cos(vPhi), math.sin(vPhi)])
269 self.
__ijkDirs = np.array([[0., 1., 0.], [0., 0., 1.], [1., 0., 0.]])
323 drdm = np.eye(3) - np.outer(trackDir, self.
__nDir) / np.dot(trackDir, self.
__nDir)
325 dmdg = np.zeros((3, 6))
326 dmdg[0][0] = 1.; dmdg[0][4] = -dist[2]; dmdg[0][5] = dist[1]
327 dmdg[1][1] = 1.; dmdg[1][3] = dist[2]; dmdg[1][5] = -dist[0]
328 dmdg[2][2] = 1.; dmdg[2][3] = -dist[1]; dmdg[2][4] = dist[0]
332 drldg = np.dot(drldrg, np.dot(drdm, dmdg))
345 uSlope = tLoc[0] / tLoc[2]
346 vSlope = tLoc[1] / tLoc[2]
351 drldg = np.array([[1.0, 0.0, -uSlope, vPos * uSlope, -uPos * uSlope, vPos], \
352 [0.0, 1.0, -vSlope, vPos * vSlope, -uPos * vSlope, -uPos]])
393 pred = layer.intersectWithHelix(hlx)
394 meas = [pred.getMeasPred()]
396 xPos, yPos = pred.getPos()[:2]
397 radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
399 cosLambda = 1. / math.sqrt(1. + localPar[3] ** 2)
401 newpar = hlx.moveTo((xPos, yPos))
402 newpar[1] += random.gauss(0., errMs / cosLambda)
403 newpar[3] += random.gauss(0., errMs / cosLambda ** 2)
406 localPar = newhlx.moveTo((-xPos, -yPos))
408 if layer.isComposite():
410 pred = layer.intersectWithHelix2(hlx)
411 meas.append(pred.getMeasPred())
413 xPos, yPos = pred.getPos()[:2]
414 cosLambda = 1. / math.sqrt(1. + localPar[3] ** 2)
416 newpar = hlx.moveTo((xPos, yPos))
417 newpar[1] += random.gauss(0., errMs / cosLambda)
418 newpar[3] += random.gauss(0., errMs / cosLambda ** 2)
421 localPar = newhlx.moveTo((-xPos, -yPos))
424 sigma = layer.getResolution()
425 hits.append((random.gauss(meas[0][0], sigma[0]), random.gauss(meas[-1][1], sigma[1])))
473 nDir = np.cross(uDir, vDir); nDir /= np.linalg.norm(nDir)
475 cosLambda = 1. / math.sqrt(1. + self.
__dzds * self.
__dzds)
476 sinLambda = self.
__dzds * cosLambda
480 tDir = np.array([cosLambda * self.
__dir0[0], cosLambda * self.
__dir0[1], sinLambda])
485 sArc3D = -np.dot(dist, nDir) / np.dot(tDir, nDir); sArc2D = sArc3D * cosLambda
487 pos = pca + sArc3D * tDir
496 dPhi = sArc2D * self.
__rinv
497 cosPhi = math.cos(self.
__phi0 + dPhi); sinPhi = math.sin(self.
__phi0 + dPhi)
498 tDir = np.array([cosLambda * cosPhi, cosLambda * sinPhi, sinLambda])
503 sCorr3D = -np.dot(dist, nDir) / np.dot(tDir, nDir)
504 if abs(sCorr3D) > 0.0001:
505 sArc2D += sCorr3D * cosLambda
510 pred = [np.dot(dist, uDir), np.dot(dist, vDir)]
529 dphi = math.atan2(dx, -dy) - self.
__phi0
530 if (abs(dphi) > math.pi):
531 dphi -= (
lambda a, b:(a > b) - (a < b))(dphi, 0.) * 2.0 * math.pi
549 dp = -newRefPoint[0] * self.
__dir0[1] + newRefPoint[1] * self.
__dir0[0] + dca
550 dl = newRefPoint[0] * self.
__dir0[0] + newRefPoint[1] * self.
__dir0[1]
551 sa = 2. * dp - rho * (dp * dp + dl * dl)
552 sb = rho * newRefPoint[0] + u * self.
__dir0[1]
553 sc = -rho * newRefPoint[1] + u * self.
__dir0[0]
554 sd = math.sqrt(1. - rho * sa)
559 newPar = [rho, phi, dca]
561 phi = math.atan2(sb, sc)
564 if abs(dphi) > math.pi: dphi -= (
lambda a, b:(a > b) - (a < b))(dphi, 0.) * 2.0 * math.pi
566 newPar = [rho, phi, dca]
587 def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos):
629 cosTheta = self.
__tdir[2]; sinTheta = math.sqrt(self.
__tdir[0] ** 2 + self.
__tdir[1] ** 2)
630 cosPhi = self.
__tdir[0] / sinTheta; sinPhi = self.
__tdir[1] / sinTheta
631 return np.array([[-sinPhi, cosPhi, 0.], [-cosPhi * cosTheta, -sinPhi * cosTheta, sinTheta]])
634if __name__ ==
'__main__':
User supplied point on (initial) trajectory.
General Broken Lines Trajectory.
Prediction (from helix at measurement)
def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos)
Constructor.
def getArcLength(self)
Get arc-length
def getPos(self)
Get Position.
def getCurvilinearDirs(self)
Get curvilinear directions (U,V)
def getCosIncidence(self)
Get cosine of incidence.
def getMeasPred(self)
Get measurement prediction.
def getDirection(self)
Get (track) direction.
def generateHits(self, qbyp, genPar)
Generate hits on helix.
def __init__(self, layers, bfac)
Constructor.
def getLayers(self)
Get layers.
def getMeasSystemDirs(self)
Get directions of measurement system.
def getRigidBodyDerGlobal(self, position, trackDir)
Get rigid body derivatives in global frame.
__measDirs
measurement directions
def getRadiationLength(self)
Get radiation length.
__precision
precision (for reconstruction)
def __init__(self, layer)
Constructor.
def intersectWithHelix(self, helix)
Intersect with helix.
def getPrecision(self)
Get precision.
__nDir
normal to measurement plane
def getResolution(self)
Get resolution.
__spacing
spacing (for composite layers)
def getSpacing(self)
Get spacing.
__resolution
resolution (for simulation)
def intersectWithHelix2(self, helix)
Intersect with helix (2nd sub layer)
def getRigidBodyDerLocal(self, position, trackDir)
Get rigid body derivatives in local (alignment) frame.
def isComposite(self)
Is composite?
__vDir
measurement direction v
__ijkDirs
local alignment system (IJK = YZX)
__uDir
measurement direction u
__dca
distance of closest approach in (XY)
__phi0
flight direction at point of closest approach (in XY)
__dir0
direction vector at point of closest approach (in XY)
def getPrediction(self, refPos, uDir, vDir)
Get prediction.
def moveTo(self, newRefPoint)
Change reference point.
def getArcLengthXY(self, xPos, yPos)
Get (2D) arc length for given point.
__xRelCenter
XY circle parameter: X position of center / R.
__z0
Z position at distance of closest approach.
def __init__(self, parameter)
Constructor.
__yRelCenter
XY circle parameter: Y position of center / R.
def exampleSit()
Simulate and reconstruct helical tracks in silicon pixel and (1D or 2D) strip detectors.
def gblMultipleScatteringError(qbyp, xbyx0)
Multiple scattering error.
def gblSimpleJacobian(ds, cosl, bfac)
Simple jacobian.