gbltst.py 12.6 KB
Newer Older
1
2
3
4
5
6
7
8
'''
Simple Test Program for General Broken Lines.

Created on Jul 27, 2011

@author: kleinwrt
'''

Claus Kleinwort's avatar
Claus Kleinwort committed
9
## \file
10
# (technical) examples
11
12
13
14
#
# \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
#
#  \copyright
Claus Kleinwort's avatar
Claus Kleinwort committed
15
#  Copyright (c) 2011 - 2021 Deutsches Elektronen-Synchroton,
16
17
18
19
20
21
22
23
24
25
26
27
28
#  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
#  This library is free software; you can redistribute it and/or modify
#  it under the terms of the GNU Library General Public License as
#  published by the Free Software Foundation; either version 2 of the
#  License, or (at your option) any later version. \n\n
#  This library is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU Library General Public License for more details. \n\n
#  You should have received a copy of the GNU Library General Public
#  License along with this program (see the file COPYING.LIB for more
#  details); if not, write to the Free Software Foundation, Inc.,
#  675 Mass Ave, Cambridge, MA 02139, USA.
Claus Kleinwort's avatar
Claus Kleinwort committed
29

30
31
32
33
34
import numpy as np
import math
import time
from gblfit import GblPoint, GblTrajectory
#
Claus Kleinwort's avatar
Claus Kleinwort committed
35

36

Claus Kleinwort's avatar
Claus Kleinwort committed
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
## 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).
#  Curvilinear system (T,U,V) or measurement system (I,J,K) as local
#  coordinate system and (q/p, slopes, offsets) as local track parameters.
#  
#  This example simulates and refits tracks in a system of planar detectors
#  with 2D measurements in a constant magnet field in Z direction using
#  the curvilinear system as local system. The true track parameters are
#  randomly smeared with respect to a (constant and straight) reference
#  trajectory with direction (lambda, phi) and are used (only) for the
#  on-the-fly simulation of the measurements and scatterers. The predictions
#  from the reference trajectory are therefore always zero and the residuals
#  needed (by addMeasurement) are equal to the measurements.
#
55
def example1():
Claus Kleinwort's avatar
Claus Kleinwort committed
56

57
#
Claus Kleinwort's avatar
Claus Kleinwort committed
58
59
  np.random.seed(47117)

60
61
  nTry = 1000  #: number of tries
  nLayer = 5  #: number of detector layers
Claus Kleinwort's avatar
Claus Kleinwort committed
62
  print " Gbltst $Id$ ", nTry, nLayer
63
64
65
66
67
68
  start = time.clock()
# track direction
  sinLambda = 0.3
  cosLambda = math.sqrt(1.0 - sinLambda ** 2)
  sinPhi = 0.
  cosPhi = math.sqrt(1.0 - sinPhi ** 2)
Claus Kleinwort's avatar
Claus Kleinwort committed
69
70
71
72
# Curvilinear system: track direction T, U = Z x T / |Z x T|, V = T x U
  tuvDir = np.array([[cosLambda * cosPhi, cosLambda * sinPhi, sinLambda], \
                     [-sinPhi, cosPhi, 0.], \
                     [-sinLambda * cosPhi, -sinLambda * sinPhi, cosLambda]])
73
# measurement resolution
74
  measErr = np.array([ 0.001, 0.001])  # 10 mu
75
76
  measPrec = 1.0 / measErr ** 2
# scattering error
77
  scatErr = 0.001  # 1 mrad
78
# RMS of CurviLinear track parameters (Q/P, slopes, offsets)
79
  clErr = np.array([0.001, -0.1, 0.2, -0.15, 0.25])
Claus Kleinwort's avatar
Claus Kleinwort committed
80
81
  # precision matrix for external seed (in local system)
  locSeed = None
82
  seedLabel = 0  # label of point with seed
83
84
  if seedLabel != 0:
    print " external seed at label ", seedLabel
85
#
86
87
  bfac = 0.2998  # Bz*c for Bz=1
  step = 1.5 / cosLambda  # constant steps in RPhi
88
89
90
91
92
93
94
95
#
  Chi2Sum = 0.
  NdfSum = 0
  LostSum = 0.
#
  binaryFile = open("milleBinaryISN.dat", "wb")
#
  for iTry in range(nTry):
96
# generate (CurviLinear) track parameters 
97
98
    clNorm = np.random.normal(0., 1., 5)  
    clPar = clErr * clNorm
99
100
101
102
    # covariance matrix
    clCov = np.eye(5)
    for i in range(5):
      clCov[i, i] = clErr[i] ** 2
103
104
# arclength
    s = 0.
105
106
# point-to-point jacobian (from previous point)    
    jacPointToPoint = np.eye(5)
107
108
# additional (local or global) derivatives    
    addDer = np.array([[1.0], [0.0]])
109
    labGlobal = np.array([[4711], [4711]])
110
111
# create trajectory
    traj = GblTrajectory(bfac != 0.)
Claus Kleinwort's avatar
Claus Kleinwort committed
112
113
114
115

# at previous point: transformation from local to curvilinear system 
    oldL2c = np.eye(5)
   
116
    for iLayer in range(nLayer):
Claus Kleinwort's avatar
Claus Kleinwort committed
117
118
      #print " layer ", iLayer
# measurement directions (J,K) from stereo angle
119
      sinStereo = (0. if iLayer % 2 == 0 else 0.1) 
120
      cosStereo = math.sqrt(1.0 - sinStereo ** 2)    
121
# (orthogonal) measurement system: I, J ,K    
Claus Kleinwort's avatar
Claus Kleinwort committed
122
123
124
125
126
127
128
129
130
131
132
133
      ijkDir = np.array([[1., 0., 0.], \
                         [0., cosStereo, sinStereo], \
                         [0., -sinStereo, cosStereo]])
# local system: measurement or curvilinear 
      #local = gblMeasSystem(ijkDir, tuvDir)
      local = gblCurviSystem(ijkDir, tuvDir)
# projections
      proL2m = local.getTransLocalToMeas()
      proL2c = local.getTransLocalToCurvi()
      proC2l = np.linalg.inv(proL2c)
# projection curvilinear to measurement directions
      proC2m = local.getTransCurviToMeas()
134
135
# measurement - prediction in measurement system with error
      measNorm = np.random.normal(0., 1., 2)  
Claus Kleinwort's avatar
Claus Kleinwort committed
136
137
138
139
140
      meas = np.dot(proC2m, clPar[3:5]) + measErr * measNorm
# jacobian is calculated in curvilinear system and transformed
      jac = np.dot(proC2l, np.dot(jacPointToPoint, oldL2c))
# point with (independent) measurements (in measurement system)
      point = GblPoint(jac)
141
142
      point.addMeasurement([proL2m, meas, measPrec])
# additional local parameters?
143
#      point.addLocals(addDer)
144
# additional global parameters?
145
      point.addGlobals(labGlobal, addDer)
146
      addDer = -addDer  # locDer flips sign every measurement      
147
148
# add point to trajectory      
      iLabel = traj.addPoint(point)
149
150
      if iLabel == abs(seedLabel):
        clSeed = np.linalg.inv(clCov)
Claus Kleinwort's avatar
Claus Kleinwort committed
151
        locSeed = np.dot(proL2c, np.dot(clSeed, proL2c.T))
152
# propagate to scatterer
153
154
      jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac)
      clPar = np.dot(jacPointToPoint, clPar)
155
      clCov = np.dot(jacPointToPoint, np.dot(clCov, jacPointToPoint.T))
156
157
158
159
      s += step
      if (iLayer < nLayer - 1):
        scat = np.array([0., 0.])
# point with scatterer
Claus Kleinwort's avatar
Claus Kleinwort committed
160
161
        jac = np.dot(proC2l, np.dot(jacPointToPoint, proL2c))
        point = GblPoint(jac)
162
163
164
        if scatErr > 0:
          scatP = local.getScatPrecision(scatErr)
          point.addScatterer([scat, scatP])
165
        iLabel = traj.addPoint(point)
166
167
        if iLabel == abs(seedLabel):
          clSeed = np.linalg.inv(clCov)
Claus Kleinwort's avatar
Claus Kleinwort committed
168
          locSeed = np.dot(proL2c, np.dot(clSeed, proL2c.T))
169

170
171
172
173
# scatter a little    
        scatNorm = np.random.normal(0., 1., 2)  
        clPar[1:3] = clPar[1:3] + scatErr * scatNorm
# propagate to next measurement layer    
174
        clPar = np.dot(jacPointToPoint, clPar)
175
        clCov = np.dot(jacPointToPoint, np.dot(clCov, jacPointToPoint.T))
176
        s += step
Claus Kleinwort's avatar
Claus Kleinwort committed
177
      oldL2c = proL2c
178
 
179
# add external seed
Claus Kleinwort's avatar
Claus Kleinwort committed
180
181
    if locSeed is not None:    
      traj.addExternalSeed(seedLabel, locSeed)
182
183
  
# fit trajectory
184
    Chi2, Ndf, Lost = traj.fit()
185
    print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost
186
187
188
# dump trajectory
    #traj.printPoints()
    #traj.printData()
189
# write to MP binary file    
Claus Kleinwort's avatar
Claus Kleinwort committed
190
#    traj.milleOut(binaryFile)
191
192
193
194
195
# sum up    
    Chi2Sum += Chi2
    NdfSum += Ndf
    LostSum += Lost
# get corrections and covariance matrix at points 
196
    if (iTry == 0):
Claus Kleinwort's avatar
Claus Kleinwort committed
197
      for i in range(1, nLayer + 1):      
198
        locPar, locCov = traj.getResults(-i)
Claus Kleinwort's avatar
Claus Kleinwort committed
199
200
201
202
        print " >Point ", i
        print " locPar ", locPar
        #print " locCov ", locCov      
        locPar, locCov = traj.getResults(i)
203
204
        print " Point> ", i
        print " locPar ", locPar
205
206
207
208
209
210
211
212
213
        #print " locCov ", locCov
# check residuals        
      for i in range(traj.getNumPoints()):
        numData, aResiduals, aMeasErr, aResErr, aDownWeight = traj.getMeasResults(i + 1)
        for j in range(numData):
          print " measRes " , i, j, aResiduals[j], aMeasErr[j], aResErr[j], aDownWeight[j]   
        numData, aResiduals, aMeasErr, aResErr, aDownWeight = traj.getScatResults(i + 1)
        for j in range(numData):
          print " scatRes " , i, j, aResiduals[j], aMeasErr[j], aResErr[j], aDownWeight[j]   
214
215
216
217
218
#
  end = time.clock()
  print " Time [s] ", end - start
  print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum
  print " LostSum/nTry ", LostSum / nTry
Claus Kleinwort's avatar
Claus Kleinwort committed
219

220

Claus Kleinwort's avatar
Claus Kleinwort committed
221
## Read trajectory from MP-II binary file and refit.
222
223
224
225
226
#
def example2():
#  
  binaryFile = open("milleBinaryISN.dat", "rb")
  nRec = 0
227
  maxRec = 10  #: maximum number of records to read
228
229
230
231
232
233
234
235
236
237
  Chi2Sum = 0.
  NdfSum = 0
  LostSum = 0.
  start = time.clock()
  
  try:
    while(nRec < maxRec):
# create trajectory
      traj = GblTrajectory(0)
# read from file      
238
      traj.milleIn(binaryFile)  # get data blocks from file
239
240
      nRec += 1
# fit trajectory      
241
      Chi2, Ndf, Lost = traj.fit()
242
243
244
245
246
247
248
249
250
251
252
253
254
255
      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
256

Claus Kleinwort's avatar
Claus Kleinwort committed
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
 
## Simple jacobian.
#  
#  Simple jacobian for (q/p, slopes, offsets) in curvilinear system
#  quadratic in arc length difference.
#
#  @param ds arc length difference; float
#  @param cosl cos(lambda); float
#  @param bfac Bz*c; float
#  @return jacobian to move by 'ds' on trajectory matrix(float)
#
def gblSimpleJacobian(ds, cosl, bfac):
  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     

276

Claus Kleinwort's avatar
Claus Kleinwort committed
277
278
279
280
281
282
283
284
## Measurement system as local system
#
# In general the precision matrix for multiple scattering is not diagonal
#
class gblMeasSystem(object):
 
  ## Construct local system
  #
285
  #  @param  measDir (directions of orthogonal) measurement system
Claus Kleinwort's avatar
Claus Kleinwort committed
286
287
288
289
290
291
  #  @param  curviDir (directions of) curvilinear system
  #  
  def __init__(self, measDir, curviDir):
    ## products (T,U,V) * (I,J,K)
    self.__prod = np.dot(curviDir, measDir.T)

292
  ## Transformation from local to measurement system 
Claus Kleinwort's avatar
Claus Kleinwort committed
293
294
295
296
297
298
299
300
301
302
303
304
305
  #
  #  @return None (systems are identical)
  #   
  def getTransLocalToMeas(self):
    return None

  ## Transformation of (q/p, slopes, offsets) from local to curvilinear system
  #
  #  @return Transformation for track parameters; 5*5 matrix
  #   
  def getTransLocalToCurvi(self):
    meas2crv = np.zeros((5, 5))
    meas2crv[0, 0] = 1.
306
307
    meas2crv[1:3, 1:3] = self.__prod[1:3, 1:3] * self.__prod[0, 0]  # (U,V)*(J,K) * T*I
    meas2crv[3:5, 3:5] = self.__prod[1:3, 1:3]  # (U,V)*(J,K)
Claus Kleinwort's avatar
Claus Kleinwort committed
308
309
310
311
312
313
314
    return meas2crv
 
  ## Transformation from curvilinear to measurement system 
  #
  #  @return Transformation for offsets; 2*2 matrix
  #   
  def getTransCurviToMeas(self):
315
    return np.linalg.inv(self.__prod[1:3, 1:3])  #((U,V)*(J,K))^-1
Claus Kleinwort's avatar
Claus Kleinwort committed
316
317
318
319
320
321
 
  ## Scattering precision matrix in local system  
  #
  #  @return Precision; 2*2 matrix
  #   
  def getScatPrecision(self, scatErr):
322
323
    c1 = self.__prod[0, 1]  # T*J
    c2 = self.__prod[0, 2]  # T*K
Claus Kleinwort's avatar
Claus Kleinwort committed
324
325
326
327
328
329
330
    fac = (1 - c1 * c1 - c2 * c2) / (scatErr * scatErr)
    scatP = np.empty((2, 2))
    scatP[0, 0] = fac * (1 - c1 * c1)
    scatP[0, 1] = fac * (-c1 * c2)
    scatP[1, 0] = fac * (-c1 * c2)
    scatP[1, 1] = fac * (1 - c2 * c2)
    return scatP       
331

Claus Kleinwort's avatar
Claus Kleinwort committed
332
333
334
335
336
337
338
339
340
341
 
## Curvilinear system as local system
#
# The precision matrix for multiple scattering (for slopes)
# is diagonal (multiple of unit matrix).
#
class gblCurviSystem(object):

  ## Construct local system
  #
342
  #  @param  measDir (directions of orthogonal) measurement system
Claus Kleinwort's avatar
Claus Kleinwort committed
343
344
345
346
347
348
  #  @param  curviDir (directions of) curvilinear system
  #    
  def __init__(self, measDir, curviDir):
    ## projection curvilinear to measurement system: ((U,V)*(J,K))^-1
    self.__c2m = np.linalg.inv(np.dot(curviDir[1:3, 1:3], measDir[1:3, 1:3].T))

349
  ## Transformation from local to measurement system 
Claus Kleinwort's avatar
Claus Kleinwort committed
350
351
352
353
  #   
  #  @return Transformation for offsets; 2*2 matrix
  #   
  def getTransLocalToMeas(self):
354
    return self.__c2m  #((U,V)*(J,K))^-1
Claus Kleinwort's avatar
Claus Kleinwort committed
355
356
357
358
359
360

  ## Transformation of (q/p, slopes, offsets) from local to curvilinear system
  #
  #  @return Transformation for track parameters; 5*5 matrix
  # 
  def getTransLocalToCurvi(self):
361
    return np.eye(5)  # unit matrix 
Claus Kleinwort's avatar
Claus Kleinwort committed
362
363
364
365
366
367
 
  ## Transformation from curvilinear to measurement system 
  #   
  #  @return Transformation for offsets; 2*2 matrix
  #   
  def getTransCurviToMeas(self):
368
    return self.__c2m  #((U,V)*(J,K))^-1
Claus Kleinwort's avatar
Claus Kleinwort committed
369
370
371
372
373
374
 
  ## Scattering precision matrix in local system  
  #
  #  @return Precision diagonal; vector(2)
  #
  def getScatPrecision(self, scatErr):
375
    return np.array([1., 1.]) / (scatErr * scatErr)  # diagonal only     
376

Claus Kleinwort's avatar
Claus Kleinwort committed
377
        
378
379
380
381
382
# 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
Claus Kleinwort's avatar
Claus Kleinwort committed
383
#example2()