32from gblfit
import GblPoint, GblTrajectory
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
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)
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 -= cmp(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 -= cmp(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 getDirection(self)
Get (track) direction.
def getArcLength(self)
Get arc-length
def getCosIncidence(self)
Get cosine of incidence.
def getCurvilinearDirs(self)
Get curvilinear directions (U,V)
def getMeasPred(self)
Get measurement prediction.
def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos)
Constructor.
def getPos(self)
Get Position.
def generateHits(self, qbyp, genPar)
Generate hits on helix.
def getLayers(self)
Get layers.
def __init__(self, layers, bfac)
Constructor.
def isComposite(self)
Is composite?
__precision
precision (for reconstruction)
__measDirs
measurement directions
def intersectWithHelix(self, helix)
Intersect with helix.
def getSpacing(self)
Get spacing.
__ijkDirs
local alignment system (IJK = YZX)
def getRadiationLength(self)
Get radiation length.
def __init__(self, layer)
Constructor.
__uDir
measurement direction u
def getRigidBodyDerLocal(self, position, trackDir)
Get rigid body derivatives in local (alignment) frame.
def getPrecision(self)
Get precision.
def intersectWithHelix2(self, helix)
Intersect with helix (2nd sub layer)
__nDir
normal to measurement plane
__resolution
resolution (for simulation)
__vDir
measurement direction v
def getResolution(self)
Get resolution.
def getMeasSystemDirs(self)
Get directions of measurement system.
def getRigidBodyDerGlobal(self, position, trackDir)
Get rigid body derivatives in global frame.
__spacing
spacing (for composite layers)
__xRelCenter
XY circle parameter: X position of center / R.
__dir0
direction vector at point of closest approach (in XY)
__phi0
flight direction at point of closest approach (in XY)
def __init__(self, parameter)
Constructor.
__z0
Z position at distance of closest approach.
__dca
distance of closest approach in (XY)
def moveTo(self, newRefPoint)
Change reference point.
def getPrediction(self, refPos, uDir, vDir)
Get prediction.
__yRelCenter
XY circle parameter: Y position of center / R.
def getArcLengthXY(self, xPos, yPos)
Get (2D) arc length for given point.
def gblSimpleJacobian(ds, cosl, bfac)
Simple jacobian.
def gblMultipleScatteringError(qbyp, xbyx0)
Multiple scattering error.
def exampleSit()
Simulate and reconstruct helical tracks in silicon pixel and (1D or 2D) strip detectors.