32from gblfit
import GblPoint, GblTrajectory
102 det =
gblSiliconDet([ [
'PIX1', (2.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],
103 [
'PIX2', (3.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],
104 [
'PIX3', (4.0, 0., 0.), 0.0033, [(0., 0.0010), (90., 0.0020)] ],
105 [
'S2D4', (6.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)], 0.1 ],
106 [
'S2D5', (8.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)], 0.1 ],
107 [
'S2D6', (10.0, 0., 0.), 0.0033, [(0., 0.0025), (5., 0.0025)] ],
108 [
'S2D7', (12.0, 0., 0.), 0.0033, [(0., 0.0025), (-5., 0.0025)] ],
109 [
'S1D8', (15.0, 0., 0.), 0.0033, [(0., 0.0040)] ]
114 det.getMP2Constraints()
118 binaryFile = open(
"milleBinary.dat",
"wb")
121 print " Gblsit $Id$ ", nTry
128 for iTry
in range(nTry):
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]
138 genHits = det.generateHits(qbyp, genPar)
143 cosLambda = 1. / math.sqrt(1. + seedPar[3] ** 2)
145 traj = GblTrajectory(bfac != 0.)
147 for l, layer
in enumerate(det.getLayers()):
149 pred = layer.intersectWithHelix(seed)
150 measPred = pred.getMeasPred()
151 sArc = pred.getArcLength()
153 res = np.array([genHits[l][0] - measPred[0], genHits[l][1] - measPred[1]])
155 measPrec = np.array(layer.getPrecision())
158 curviDirs = pred.getCurvilinearDirs()
160 proL2m = np.linalg.inv(np.dot(curviDirs, np.linalg.inv(layer.getMeasSystemDirs())[:,:2]))
165 point = GblPoint(jacPointToPoint)
167 if layer.isComposite():
169 pred = layer.intersectWithHelix2(seed)
170 measPred = pred.getMeasPred()
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])
178 point.addMeasurement([proL2m, res, measPrec])
180 if binaryFile
is not None:
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())
186 if layer.isComposite():
188 labG4D = np.array([labels, labels, labels, labels])
189 derG4D = np.zeros((4, 6)); derG4D[2:] = derGlobal
190 point.addGlobals(labG4D, derG4D)
193 point.addGlobals(labGlobal, derGlobal)
195 radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
198 scat = np.array([0., 0.])
199 scatP = np.array([1. / scatErr ** 2, 1. / scatErr ** 2])
201 if layer.isComposite():
204 point.addScatterer([scat, scatP])
209 Chi2, Ndf, Lost = traj.fit()
210 print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost
216 if binaryFile
is not None:
217 traj.milleOut(binaryFile)
220 print " Time [s] ", end - start
221 print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum
222 print " LostSum/nTry ", LostSum / nTry
237 jac[1, 0] = -bfac * ds * cosl
238 jac[3, 0] = -0.5 * bfac * ds * ds * cosl
252 return 0.015 * abs(qbyp) * math.sqrt(xbyx0)
274 self.
__precision = (1. / meas[0][1] ** 2, 1. / meas[1][1] ** 2
if len(meas) > 1
else 0.)
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
279 self.
__uDir = np.array([0., math.cos(uPhi), math.sin(uPhi)])
281 self.
__vDir = np.array([0., math.cos(vPhi), math.sin(vPhi)])
287 self.
__ijkDirs = np.array([[0., 1., 0.], [0., 0., 1.], [1., 0., 0.]])
306 print "Constraint 0. ! fix unmeasured direction in", self.
__name
309 if abs(unMeasured[i]) > 1.0e-10:
310 print " ", layer * 10 + i + 1, unMeasured[i]
362 drdm = np.eye(3) - np.outer(trackDir, self.
__nDir) / np.dot(trackDir, self.
__nDir)
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]
371 drldg = np.dot(drldrg, np.dot(drdm, dmdg))
384 uSlope = tLoc[0] / tLoc[2]
385 vSlope = tLoc[1] / tLoc[2]
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]])
397 return np.dot(local2meas[:2,:2], drldg)
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"
444 pred = layer.intersectWithHelix(hlx)
445 meas = [pred.getMeasPred()]
447 xPos, yPos = pred.getPos()[:2]
448 radlen = layer.getRadiationLength() / abs(pred.getCosIncidence())
450 cosLambda = 1. / math.sqrt(1. + localPar[3] ** 2)
452 newpar = hlx.moveTo((xPos, yPos))
453 newpar[1] += random.gauss(0., errMs / cosLambda)
454 newpar[3] += random.gauss(0., errMs / cosLambda ** 2)
457 localPar = newhlx.moveTo((-xPos, -yPos))
459 if layer.isComposite():
461 pred = layer.intersectWithHelix2(hlx)
462 meas.append(pred.getMeasPred())
464 xPos, yPos = pred.getPos()[:2]
465 cosLambda = 1. / math.sqrt(1. + localPar[3] ** 2)
467 newpar = hlx.moveTo((xPos, yPos))
468 newpar[1] += random.gauss(0., errMs / cosLambda)
469 newpar[3] += random.gauss(0., errMs / cosLambda ** 2)
472 localPar = newhlx.moveTo((-xPos, -yPos))
475 sigma = layer.getResolution()
476 hits.append((random.gauss(meas[0][0], sigma[0]), random.gauss(meas[-1][1], sigma[1])))
524 nDir = np.cross(uDir, vDir); nDir /= np.linalg.norm(nDir)
526 cosLambda = 1. / math.sqrt(1. + self.
__dzds * self.
__dzds)
527 sinLambda = self.
__dzds * cosLambda
531 tDir = np.array([cosLambda * self.
__dir0[0], cosLambda * self.
__dir0[1], sinLambda])
536 sArc3D = -np.dot(dist, nDir) / np.dot(tDir, nDir); sArc2D = sArc3D * cosLambda
538 pos = pca + sArc3D * tDir
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])
554 sCorr3D = -np.dot(dist, nDir) / np.dot(tDir, nDir)
555 if abs(sCorr3D) > 0.0001:
556 sArc2D += sCorr3D * cosLambda
561 pred = [np.dot(dist, uDir), np.dot(dist, vDir)]
580 dphi = math.atan2(dx, -dy) - self.
__phi0
581 if (abs(dphi) > math.pi):
582 dphi -= cmp(dphi, 0.) * 2.0 * math.pi
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)
610 newPar = [rho, phi, dca]
612 phi = math.atan2(sb, sc)
615 if abs(dphi) > math.pi: dphi -= cmp(dphi, 0.) * 2.0 * math.pi
617 newPar = [rho, phi, dca]
638 def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos):
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]])
685if __name__ ==
'__main__':
Prediction (from helix at measurement)
def getCosIncidence(self)
Get cosine of incidence.
def getArcLength(self)
Get arc-length
def getCurvilinearDirs(self)
Get curvilinear directions (U,V)
def getDirection(self)
Get (track) direction.
def getPos(self)
Get Position.
def __init__(self, sArc, pred, tDir, uDir, vDir, nDir, pos)
Constructor.
def getMeasPred(self)
Get measurement prediction.
def generateHits(self, qbyp, genPar)
Generate hits on helix.
def __init__(self, layers, bfac)
Constructor.
def getMP2Constraints(self)
get MP2 constraints
def getLayers(self)
Get layers.
__vDir
measurement direction v
def getSpacing(self)
Get spacing.
__resolution
resolution (for simulation)
__measDirs
measurement directions
def getRigidBodyDerGlobal(self, position, trackDir)
Get rigid body derivatives in global frame.
__spacing
spacing (for composite layers)
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?
__uDir
measurement direction u
def intersectWithHelix(self, helix)
Intersect with helix.
__nDir
normal to measurement plane
def getMP2Constraint(self, layer)
get MP2 constraint
__alignInMeasSys
alignment == measurement system?
def __init__(self, layer)
Constructor.
__ijkDirs
local alignment system (IJK = YZX)
def getResolution(self)
Get resolution.
def getPrecision(self)
Get precision.
def getRadiationLength(self)
Get radiation length.
__precision
precision (for reconstruction)
def getMeasSystemDirs(self)
Get directions of measurement system.
__dca
distance of closest approach in (XY)
def getPrediction(self, refPos, uDir, vDir)
Get prediction.
__yRelCenter
XY circle parameter: Y position of center / R.
def moveTo(self, newRefPoint)
Change reference point.
def getArcLengthXY(self, xPos, yPos)
Get (2D) arc length for given point.
__z0
Z position at distance of closest approach.
__xRelCenter
XY circle parameter: X position of center / R.
def __init__(self, parameter)
Constructor.
__dir0
direction vector at point of closest approach (in XY)
__phi0
flight direction at point of closest approach (in XY)
def gblMultipleScatteringError(qbyp, xbyx0)
Multiple scattering error.
def exampleSit()
Simulate and reconstruct helical tracks in silicon pixel and (1D or 2D) strip detectors.
def gblSimpleJacobian(ds, cosl, bfac)
Simple jacobian.