gblfit.py 35.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
'''
Track fit with general broken lines.

Created on Jul 27, 2011

@author: kleinwrt
'''

import numpy as np
import math
from gblnum import BorderedBandMatrix
from mille import MilleRecord

class GblPoint(object):
  '''
  User supplied point on (initial) trajectory.
  
  Must have jacobians for propagation to previous or next point with offsets (first, last
  point, points with scatterer). May have:
  
    1. Measurement (1D or 2D)
    2. Scatterer (thin, 2D kinks)
    3. Additional local parameters (with derivatives). Fitted together with track parameters.
    4. Additional global parameters (with labels and derivatives). Not fitted, only passed
       on to (binary) file for fitting with Millepede-II. 
  '''
27
  def __init__(self, aJacobian):
28
29
    '''
    Create new point.
30
31
32
    
    @param aJacobian: jacobian from previous point
    @type matrix(float)
33
34
35
36
37
38
39
    '''
    self.__label = 0
    '''@ivar: label for referencing point (0,1,..,number(points)-1)
       @type: int'''
    self.__offset = 0
    '''@ivar: >=0: offset number at point, <0: offset number at next point with offset
       @type: int'''
40
41
42
43
    self.__p2pJacobian = aJacobian
    '''@ivar: Point-to-point jacobian from previous point
       @type: matrix(float)'''
    self.__jacobians = [ [], [] ]
44
    '''@ivar: jacobians for propagation to previous or next point with offsets
45
46
       @type: pair(matrix(float))'''
    self.__measurement = None
47
48
    '''@ivar: measurement at point: projection (dm/du), residuals (to initial trajectory), precision
       @type: list(matrix(float))'''
49
50
51
52
53
54
55
    self.__measDim = 0
    '''@ivar: dimension of measurement (2D, 4D or 5D)
       @type: int''' 
    self.__measTransformation = None    
    '''@ivar: transformation (to eigen-vectors of precision matrix)
       @type: matrix(float))''' 
    self.__scatterer = None
56
    '''@ivar: scatterer at point: initial kinks, precision (inverse covariance matrix)
57
58
       @type: matrix(float)'''
    self.__localDerivatives = None
59
    '''@ivar: local derivatives
60
61
       @type: matrix(float))'''    
    self.__globalLabels = None
62
    '''@ivar: global labels
63
64
       @type: matrix(int)'''  
    self.__globalDerivatives = None
65
    '''@ivar: global derivatives
66
       @type: matrix(float)'''  
67
68
69
70
71
72
73
74
75
76
      
# for extension to retrieval of residuals, pulls    
#    self.__dataMeas = [0, 0]
# for extension to retrieval of residuals, pulls    
#    self.__dataScat = [0, 0]

  def addMeasurement(self, aMeasurement):
    '''
    Add a mesurement to a point.
   
77
78
    @param aMeasurement: measurement (projection (or None), residuals, precision
                         (diagonal of or full matrix))
79
80
81
    @type aMeasurement: list(matrix(float))
    '''
    self.__measurement = aMeasurement
82
83
84
85
86
87
88
89
90
91
92
    self.__measDim = aMeasurement[1].shape[0]
    if (aMeasurement[2].ndim == 2): # full precision matrix, need to diagonalize
      eigenVal, eigenVec = np.linalg.eigh(aMeasurement[2])
      self.__measTransformation = eigenVec.T
#     transform measurement
      if (aMeasurement[0] is None):
        self.__measurement[0] = self.__measTransformation
      else:
        self.__measurement[0] = np.dot(self.__measTransformation, aMeasurement[0])
      self.__measurement[1] = np.dot(self.__measTransformation, aMeasurement[1])
      self.__measurement[2] = eigenVal
93
94
95
96
97
98
99

  def hasMeasurement(self):
    '''
    Check point for a measurement.
    
    @rtype: bool 
    '''
100
    return (self.__measurement is not None)
101
102
103
104
105
106
107
108
109

  def getMeasurement(self):
    '''
    Retrieve measurement of a point.
    
    @return: measurement (projection residuals, precision)
    @rtype: list(matrix(float))
    '''
    return self.__measurement
110
111
112
113
 
  def getMeasDim(self):
    '''
    Retrieve dimension of measurement of a point.
114
    
115
116
117
118
119
    @return: measurement dimension (2, 4 or 5)
    @rtype: int
    '''
    return self.__measDim
     
120
  def addScatterer(self, aScatterer):
121
122
    '''
    Add a (thin) scatterer to a point.
123
    
124
125
    Changes local track direction.
    
126
127
128
129
130
131
132
133
134
135
136
    @param aScatterer: scatterer (kinks, precision)
    @type aScatterer: list(matrix(float))
    '''
    self.__scatterer = aScatterer
  
  def hasScatterer(self):
    '''
    Check point for a scatterer.
    
    @rtype: bool     
    '''
137
    return (self.__scatterer is not None)
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152

  def getScatterer(self):
    '''
    Retrieve scatterer of a point.
    
    @return: scatterer (kinks, precision)
    @rtype: list(matrix(float))    
    '''
    return self.__scatterer
  
  def addLocals(self, derivatives):
    '''
    Add local derivatives.
    
    @param derivatives: local derivatives
153
    @type derivatives: matrix(float)    
154
    '''
155
156
157
158
159
    if (self.__measDim > 0):
      if (self.__measTransformation is None):
        self.__localDerivatives = derivatives
      else:
        self.__localDerivatives = np.dot(self.__measTransformation, derivatives)        
160
161
162
163
164
165

  def addGlobals(self, labels, derivatives):
    '''
    Add global derivatives.
    
    @param labels: global labels
166
    @type labels: matrix(int)
167
    @param derivatives: global derivatives
168
    @type derivatives: matrix(float)
169
    '''
170
171
172
173
174
175
    if (self.__measDim > 0):
      self.__globalLabels = labels
      if (self.__measTransformation is None):
        self.__globalDerivatives = derivatives
      else:
        self.__globalDerivatives = np.dot(self.__measTransformation, derivatives) 
176
177
178
179
180
181
182
183

  def getNumLocals(self):
    '''
    Get number of local derivatives.
    
    @return: Number of local derivatives at point
    @rtype: int
    '''
184
185
    if (self.__localDerivatives is not None):
      return self.__localDerivatives.shape[1]
186
187
188
189
190
191
192
193
194
195
    else:
      return 0
  
  def getLocalDerivatives(self):
    '''
    Get local derivatives.
    
    @return: local derivatives
    @rtype: matrix(float)
    '''
196
    return self.__localDerivatives
197
198
199
200
201
202
203
204

  def getGlobalLabels(self):
    '''
    Get global labels.
    
    @return: lglobal labels
    @rtype: matrix(int)   
    '''
205
    return self.__globalLabels
206
207
208
209
210
211
212
213
    
  def getGlobalDerivatives(self):
    '''
    Get global derivatives.
    
    @return: global derivatives
    @rtype: matrix(float)   
    '''
214
    return self.__globalDerivatives
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
    
  def setLabel(self, aLabel):
    '''
    Define label of a point.
    
    @param aLabel: label
    @type aLabel: int
    '''
    self.__label = aLabel

  def getLabel(self):
    '''
    Retrieve label of a point.
    
    @return: label
    @rtype: int
    '''
    return self.__label
 
234
  def setOffset(self, anOffset):
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
    '''
    Define offset of a point and references to previous and next point with offset.
    
    @param anOffset: offset number (at point (>=0) or at next point with offset (<0))
    @type anOffset: int
    '''
    self.__offset = anOffset

  def getOffset(self):
    '''
    Get offset of a point.
    
    @return: offset number (at point (>=0) or at next point with offset (<0))
    @rtype: int
    '''
    return self.__offset
  
#  def setDataMeas(self, aIndex, aData):
# for extension to retrieval of residuals, pulls    
#    self.__dataMeas[aIndex] = aData

#  def setDataScat(self, aIndex, aData):
# for extension to retrieval of residuals, pulls    
#    self.__dataScat[aIndex] = aData
259
260
 
  def getP2pJacobian(self):
261
    '''
262
    Retrieve jacobian of a point.
263
    
264
265
266
267
268
269
    @return: point-to-point jacobian (from previous point)
    @rtype: matrix(float)
    '''
    return self.__p2pJacobian         
 
  def addPrevJacobian(self, aJacobian):
270
    '''
271
    Add jacobian to previous offset.
272
    
273
274
275
276
277
278
    @param aJacobian: jacobian for propagation to previous point with offsets
    @type aJacobian: matrix(float)
    '''
    self.__jacobians[0] = np.linalg.inv(aJacobian)

  def addNextJacobian(self, aJacobian):
279
    '''
280
    Add jacobian to next offset.
281
    
282
283
    @param aJacobian: jacobian for propagation to next point with offsets
    @type aJacobian: matrix(float)
284
    '''
285
    self.__jacobians[1] = aJacobian
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
    
  def getDerivatives(self, index):
    '''
    Get derivatives for locally linearized track model (backward or forward propagation).
    
    @param index: 0 (previous) or 1 (next point with offsets)
    @type index: int
    @return: derivatives
    @rtype: list(matrix(float))
    '''
    aJacobian = self.__jacobians[index]
    matJ = aJacobian[3:5, 3:5] # J
    matS = aJacobian[3:5, 1:3] # S 
    vecD = aJacobian[3:5, 0:1] # d
    if (index < 1):
      matS = -matS
    matW = np.linalg.inv(matS) # W = +/- S^-1
303
    return matW, np.dot(matW, matJ), np.dot(matW, vecD)  # W, W*J, W*d
304
305
306
307
308
    
  def printPoint(self):
    '''
    Print point.
    '''
309
    print " point ", self.__label, self.__offset 
310
311
312
313
314

#------------------------------------------------------------------------------ 

class GblData(object):
  '''
315
316
  Data (block) containing value, precision and derivatives for measurements, kinks
  and external seed.
317
318
319
  
  Created from attributes of GblPoints, used to construct linear equation system for track fit.
  '''
320
  def __init__(self, aLabel=0, aValue=0., aPrec=0.):
321
322
    '''
    Create new data.
323
324
325
326
327
328
329
    
    @param aLabel: label of corresponding point
    @type aLabel: int
    @param aValue: value
    @type aValue: float
    @param aPrec: precision
    @type aPrec: float
330
    '''
331
    self.__label = aLabel
332
333
    '''@ivar: label of corresponding point
       @type: int'''
334
    self.__value = aValue
335
336
    '''@ivar: value (residual or kink)
       @type: float'''
337
    self.__precision = aPrec
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
    '''@ivar: precision (diagonal element of inverse covariance matrix)
       @type: float'''
    self.__downWeight = 1.
    '''@ivar: down weighting factor (M-estimators)
       @type: float'''
    self.__prediction = 0.
    '''@ivar: prediction (for value from fit)
       @type: float'''
    self.__parameters = []
    '''@ivar: labels of fit parameters (with non zero derivative)
       @type: list(int)'''
    self.__derivatives = []
    '''@ivar: derivatives (prediction vs fit parameters)
       @type: list(float)'''
    self.__globalLabels = []
    '''@ivar: labels of global parameters
       @type: list(int)'''
    self.__globalDerivatives = []
    '''@ivar: derivatives (prediction vs global parameters)
       @type: list(float)'''

#    self.__predictionVar = 0.
               
361
362
  def addDerivatives(self, iRow, labDer, matDer, \
                     derLocal=None, labGlobal=None, derGlobal=None):
363
    '''
364
    Add derivatives to data (block) from measurements and kinks. Generate lists of labels.
365
    
366
367
368
369
370
371
372
373
    @param iRow: row of measurement (vector)
    @type iRow: int
    @param labDer: labels of derivatives vs curvature and band parameters (offsets)
    @type labDer: vector(int)
    @param matDer: derivatives vs curvature and band parameters (offsets)
    @type matDer: matrix(float)
    @param derLocal: derivatives vs local parameters
    @type derLocal: list(float)
374
375
376
377
378
    @param labGlobal: labels of global parameters
    @type labGlobal: list(int)
    @param derGlobal: derivatives vs global parameters
    @type derGlobal: list(float)
    '''   
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
    if (derLocal is not None):
      for i in range(derLocal.shape[1]):           # local derivatives
        if (derLocal[iRow, i] != 0.):
          self.__derivatives.append(derLocal[iRow, i])
          self.__parameters.append(i + 1)
        
    for i in range(len(labDer)):           # curvature, offset derivatives
      if (labDer[i] != 0 and matDer[iRow, i] != 0.):
        self.__derivatives.append(matDer[iRow , i])
        self.__parameters.append(labDer[i])

    if (derGlobal is not None):  
      for i in range(derGlobal.shape[1]):          # global derivatives
        if (derGlobal[iRow, i] != 0.):
          self.__globalLabels.append(labGlobal[iRow, i])
          self.__globalDerivatives.append(derGlobal[iRow, i])

  def addExtDerivatives(self, indexExt, derExt):
    '''
    Add derivatives to data (block) from external seed. Generate lists of labels.
    
    @param indexExt: labels from exteranl seed
    @type indexExt: list(int)
    @param derExt: derivatives from exteranl seed
    @type derExt: list(float)
    '''       
    for i in range(len(derExt)):           # external derivatives
      if (derExt[i] != 0.):
        self.__derivatives.append(derExt[i])
        self.__parameters.append(indexExt[i])  
409
410
411
412
413
414
415
416
417
418
419
420

  def getMatrices(self):  
    '''
    Calculate compressed matrix and right hand side from data.
    
    @return: indices, compressed right hand side and matrix
    @rtype: list
    '''
    aVector = np.array([ self.__derivatives ])
    aMatrix = np.dot(aVector.T, aVector)
    aValue = self.__value
    aWeight = self.__precision * self.__downWeight
421
    return self.__parameters, aValue * aVector * aWeight, aMatrix * aWeight
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
    
  def setPrediction(self, aVector):
    '''
    Calculate prediction for data from fit.
   
    @param aVector: values of fit parameters
    @type aVector: vector(float)
    '''
    self.__prediction = 0.
    for i in range(len(self.__parameters)):
      self.__prediction += self.__derivatives[i] * aVector[ self.__parameters[i] - 1 ]

#  def setPredictionVariance(self, aMatrix):
#    '''Calculate variance of prediction for data from fit.'''
# var(residual) = 1./precision**2 - var(prediction)
#    aBlockMatrix = aMatrix.getBlockMatrix(self.__parameters)
#    self.__predictionVar = np.dot(self.__derivatives.T, \
#                             np.dot(aBlockMatrix, self.__derivatives))
           
  def setDownWeighting(self, aMethod): 
    '''
    Outlier down weighting with M-estimators.
    
    @param aMethod: method (1=Tukey, 2=Huber, 3=Cauchy) 
    @type aMethod: int
    @return: weight (0..1)
    @rtype: float 
    '''
    scaledResidual = abs(self.__value - self.__prediction) * math.sqrt(self.__precision)   
    if (aMethod == 1):  # Tukey
      if (scaledResidual < 4.6851):
        aWeight = (1.0 - 0.045558 * scaledResidual ** 2) ** 2
      else:
        aWeight = 0.
    elif (aMethod == 2): # Huber
      if (scaledResidual < 1.345):
        aWeight = 1.
      else:
        aWeight = 1.345 / scaledResidual
    elif (aMethod == 3):  # Cauchy
      aWeight = 1.0 / (1.0 + (scaledResidual / 2.3849) ** 2)      
    self.__downWeight = aWeight
    return aWeight
        
  def getChi2(self):
    '''
    Calculate Chi2 (contribution) from data.
    
    @return: Chi2
    @rtype: float
    '''
    Chi2 = (self.__value - self.__prediction) ** 2 * self.__precision * self.__downWeight
    return Chi2
 
  def toRecord(self):
    '''
    Get data components (for copying to MP binaty record)
    
    @return: data components
    @rtype: list
    '''
483
484
    return self.__value, self.__precision, self.__parameters, self.__derivatives, \
            self.__globalLabels, self.__globalDerivatives
485
486
487
488
489
490
491
492

  def fromRecord(self, dataList):
    '''
    Set data components (from MP binaty record)
   
    @param dataList: data components
    @type dataList: list
    '''
493
494
    self.__value, self.__precision, self.__parameters, self.__derivatives, \
            self.__globalLabels, self.__globalDerivatives = dataList
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
            
  def analyzeData(self, maxBand):
    '''
    Analyze labels of fit parameters to determine number of parameters and
    border size with given maximal band width.
    
    @param maxBand: maximal band width
    @type maxBand: int   
    @return: number of parameters and border size (from this data)
    @rtype: pair(int) 
    '''
    maxPar = self.__parameters[-1]
    maxBor = 0
    for i in self.__parameters:
      if (i < maxPar - maxBand):
        maxBor = i
511
    return maxPar, maxBor
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
         
  def printData(self):
    '''
    Print data.
    '''
    print " meas. ", self.__label, self.__value, self.__precision
    print " param ", self.__parameters
    print " deriv ", self.__derivatives
 
#------------------------------------------------------------------------------ 
    
class GblTrajectory(object):
  '''
  General Broken Lines Trajectory.
  ================================
  
  For a track with an initial trajectory from a prefit of the
529
530
531
  (2D, 4D or 5D) measurements (internal seed) or an external 
  prediction(external seed) the description of multiple scattering
  is added by offsets in a local system. Along the initial
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
  trajectory points are defined with can describe a measurement
  or a (thin) scatterer or both. The refit provides corrections
  to the local track parameters (in the local system) and the 
  corresponding covariance matrix at any of those points.

  The broken lines trajectory is defined by (2D) offsets at the 
  first and last point and all points with a scatterer. The
  prediction for a measurement is obtained by interpolation of
  the enclosing offsets and for triplets of adjacent offsets
  kink angles are determined. This requires for all points the
  jacobians for propagation to the previous and next offset.

  Additional local or global parameters can be added and the
  trajectories can be written to special binary files for
  calibration and alignment with Millepede-II.
  (V. Blobel, NIM A, 566 (2006), pp. 5-13).

  The conventions for the coordinate systems follow:
  Derivation of Jacobians for the propagation of covariance
  matrices of track parameters in homogeneous magnetic fields
  A. Strandlie, W. Wittek, NIM A, 566 (2006) 687-698.
  
  Calling sequence:
  =================
    1. Create trajectory::
            traj = GblTrajectory()
    2. For all points on initial trajectory 
559
560
561
        - Create point (supply jacobian from previous point)::
            point = GblPoint(jacobian)
        - Optionally add measurement to point::    
562
            point.addMeasurement(..)
563
        - Optionally additional local or global parameters for measurement:: 
564
565
            point.addLocals(..)
            point.addGlobals(..)
566
567
        - Optionally add scatterer to point::    
            point.addScatterer(..)
568
569
        - Add point (ordered by arc length) to trajectory, get label of point::  
            label = traj.addPoint(point)
570
    3. Optionally add external seed::
571
            traj.addExternalSeed(..)
572
    4. Fit trajectory, bet Chi2, Ndf (and weight lost by M-estimators)::
573
            [..] = traj.fit()
574
    5. For any point on inital trajectory
575
576
        - Get corrections and covariance matrix for track parameters::
            [..] = traj.getResults(label) 
577
578
    6. Optionally write trajectory to MP binary file::
            traj.milleOut(..)
579
580
581
582
            
  Alternatively trajectories can by read from MP binary files and fitted. 
  As the points on the initial trajectory are not stored in this files results at
  points (corrections, covariance matrix) are not available.
583
584
585
586
587
588
589
590
591
  
  References:
  ===========  
    - V. Blobel, C. Kleinwort, F. Meier,
      Fast alignment of a complex tracking detector using advanced track models,
      Computer Phys. Communications (2011), doi:10.1016/j.cpc.2011.03.017
    - C. Kleinwort, General Broken Lines as advanced track fitting method,
      NIM A, 673 (2012), 107-110, doi:10.1016/j.nima.2012.01.024  
      
592
593
594
595
  ''' 
  def __init__(self, hasCurv=True, aDim=[0, 1]):
    '''
    Create new trajectory.
596
597
598
599
600
    
    @param hasCurv: flag for curvature
    @type hasCurv: bool
    @param aDim: active offset components (0:u_1, 1:u_2)
    @type aDim: list(int)
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
    '''  
    self.__numPoints = 0
    '''@ivar: number of points on trajectory
       @type: int'''
    self.__numOffsets = 0
    '''@ivar: number of (points with) offsets on trajectory
       @type: int'''
    self.__numCurvature = (1 if hasCurv else 0)
    '''@ivar: 'curvature' is fit parameter (=1)
       @type: int'''
    self.__numParameters = 0
    '''@ivar: number fit parameters
       @type: int'''
    self.__numLocals = 0
    '''@ivar: number of local parameters
       @type: int'''    
    self.__externalPoint = 0
    '''@ivar: label of point with external seed
       @type: int'''    
    self.__dimensions = aDim
621
    '''@ivar: active components of offsets (both ([0,1]) or single ([0] or [1])
622
623
624
625
626
627
628
       @type: list(int)'''
    self.__points = [] 
    '''@ivar: points on trajectory
       @type: list(GblPoint)'''
    self.__data = []
    '''@ivar: data (blocks) of trajectory
       @type: list(GblData)'''
629
    self.__externalSeed = None
630
    '''@ivar: external seed (for local, fit parameters)
631
       @type: matrix(float)'''
632
633
634
635
636
637
638

  def addPoint(self, point):
    '''
    Add point to trajectory. Points have to be ordered in arc length.
    
    @param point: point to be added
    @type point: GblPoint
639
    @return: label of added point (1..number(points))
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
    @rtype: int    
    '''
    self.__numPoints += 1
    label = self.__numPoints
    point.setLabel(label)
    self.__points.append(point)
    self.__numLocals = max(self.__numLocals, point.getNumLocals())
    return label
    
  def getNumPoints(self):
    '''
    Get number of points on trajectory.
    
    @return: number of points
    @rtype: int
    '''  
    return self.__numPoints
  
  def addExternalSeed(self, aLabel, aSeed):
    '''
    Add external seed to trajectory.
    
    @param aLabel: label of point with external seed
    @type aLabel: int 
664
    @param aSeed: seed (precision matrix of track parameters at point)
665
    @type aSeed: matrix(float) 
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
    '''
    self.__externalPoint = aLabel
    self.__externalSeed = aSeed
    
  def dump(self):
    '''
    Dump trajectory.
    '''
    for p in self.__points:
      p.printPoint() 
      
  def milleOut(self, aFile):
    '''
    Write (data blocks of) trajectory to MP (binary) file.
    
    @param aFile: MP file
    @type aFile: file
    '''
    rec = MilleRecord()
#   data: measurements and kinks        
    for aData in self.__data:       
      rec.addData(aData.toRecord())
                    
    rec.writeRecord(aFile)

  def milleIn(self, aFile):
    '''
    Read (data blocks of) trajectory from MP (binary) file.
    
    @param aFile: MP file
    @type aFile: file
    '''
    rec = MilleRecord()
    rec.readRecord(aFile)
    mPar = 0
    mBor = 0
    mBand = 3 * len(self.__dimensions) - 1  # max band width
    while (rec.moreData()):
      aTag = rec.specialDataTag()
      if (aTag < 0):
# get data
        aData = GblData()
        aData.fromRecord(rec.getData())
        self.__data.append(aData)
710
        nPar, nBor = aData.analyzeData(mBand)
711
712
713
714
715
716
717
718
719
720
        mPar = max(mPar, nPar)
        mBor = max(mBor, nBor)
        
    self.__numParameters = mPar
    self.__numLocals = mBor - self.__numCurvature
   
  def __getJacobian(self, aLabel):
    '''
    Get jacobian for transformation from fit to track parameters at point.
    
721
    @param aLabel: (signed) label of point
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
    @type aLabel: int
    @return: labels of required fit parameters and jacobian
    @rtype: list
    '''
    aDim = self.__dimensions
    nDim = len(aDim)
    anIndex = abs(aLabel) - 1
#   check consistency of (index, direction)    
    if (aLabel > 0):
      nJacobian = 1
      if (anIndex >= self.__numPoints - 1):
        anIndex = self.__numPoints - 1
        nJacobian = 0
    else:
      nJacobian = 0
      if (anIndex <= 0):
        anIndex = 0
        nJacobian = 1
# Jacobian broken lines (q/p,..,u_i,u_i+1..) to local (q/p,u',u) parameters   
    nCurv = self.__numCurvature
742
    nLocals = self.__numLocals   
743
744
745
746
747
    nBorder = nCurv + nLocals
    nParBRL = nBorder + 2 * nDim
    nParLoc = nLocals + 5
    aJacobian = np.zeros((nParLoc, nParBRL))
    aPoint = self.__points[anIndex]
748
749
750
    anIndex = []
    labDer, matDer = self.__getFitToLocalJacobian(aPoint, 5, nJacobian)
#   from local parameters
751
    for i in range(nLocals):
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
      aJacobian[i + 5, i] = 1.0;
      anIndex.append(i + 1);
  
#   from trajectory parameters
    iCol = nLocals;
    for i in range(5):
      if (labDer[i] > 0):
        anIndex.append(labDer[i]);
        for j in range(5):
          aJacobian[j, iCol] = matDer[j, i];
        iCol += 1

    return anIndex, aJacobian
   
  def __getFitToLocalJacobian(self, aPoint, measDim, nJacobian=1):
    '''
    Get (part of) jacobian for transformation from (trajectory) fit to track parameters at point.

    Jacobian broken lines (q/p,..,u_i,u_i+1..) to local (q/p,u',u) parameters.
    @param aPoint: point to use
    @type aPoint: GblPoint
    @param measDim: dimension of 'measurement' (2, 4 or 5)
    @type measDim: int
    @param nJacobian: direction (0: to previous offset, 1: to next offset)
    @type nJacobian: int
    @return: labels for fit parameters with non zero derivatives, corresponding transformation matrix
    @rtype: list(vector(int), matrix(float)) 
    '''
    aDim = self.__dimensions
    nDim = len(aDim)
    nCurv = self.__numCurvature
    nLocals = self.__numLocals      
    nOffset = aPoint.getOffset() 
    anIndex = [0, 0, 0, 0, 0]
    aJacobian = np.zeros((measDim, 5))
    labOffset = measDim - 2
    labSlope = measDim - 4
    labCurv = measDim - 5
    
791
    if (nOffset < 0): # need interpolation
792
793
      prevW, prevWJ, prevWd = aPoint.getDerivatives(0) # W-, W- * J-, W- * d-
      nextW, nextWJ, nextWd = aPoint.getDerivatives(1) # W+, W+ * J+, W+ * d+
794
795
      sumWJ = prevWJ + nextWJ
      matN = np.linalg.inv(sumWJ)      # N = (W- * J- + W+ * J+)^-1 
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
#     local offset
      if (labOffset >= 0):
#       derivatives for u_int      
        prevNW = np.dot(matN, prevW)     # N * W-
        nextNW = np.dot(matN, nextW)     # N * W+
        prevNd = np.dot(matN, prevWd)    # N * W- * d-
        nextNd = np.dot(matN, nextWd)    # N * W+ * d+
        iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1 # first offset ('i' in u_i)
        if (nCurv > 0):
          aJacobian[labOffset:measDim, 0:1] = -prevNd - nextNd # from curvature
          anIndex[0] = nLocals + 1
        aJacobian[labOffset:measDim, 1:3] = prevNW
        aJacobian[labOffset:measDim, 3:5] = nextNW
        for i in range(nDim):
          anIndex[1 + aDim[i]] = iOff + i
          anIndex[3 + aDim[i]] = iOff + nDim + i
#     local slope 
      if (labSlope >= 0):
#       derivatives for u'_int
        prevWPN = np.dot(nextWJ, prevNW) # W+ * J+ * N * W-
        nextWPN = np.dot(prevWJ, nextNW) # W- * J- * N * W+
        prevWNd = np.dot(nextWJ, prevNd) # W+ * J+ * N * W- * d-
        nextWNd = np.dot(prevWJ, nextNd) # W- * J- * N * W+ * d+
        if (nCurv > 0):
          aJacobian[labSlope:labOffset, 0:1] = prevWNd - nextWNd # from curvature
        aJacobian[labSlope:labOffset, 1:3] = -prevWPN
        aJacobian[labSlope:labOffset, 3:5] = nextWPN

    else: # at point               
#     anIndex must be sorted
#     forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
#     backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
      iOff1 = nDim * nOffset + nCurv + nLocals + 1 # first offset ('i' in u_i)
      index1 = 3 - 2 * nJacobian # index of first offset
      iOff2 = iOff1 + nDim * (nJacobian * 2 - 1) # second offset ('i' in u_i)
      index2 = 1 + 2 * nJacobian # index of second offset
#     local offset
      if (labOffset >= 0):
        aJacobian[labOffset, index1] = 1.0 # from 1st Offset
        aJacobian[labOffset + 1, index1 + 1] = 1.0
        for i in range(nDim):
          anIndex[index1 + aDim[i]] = iOff1 + i
#     local slope and curvature
      if (labSlope >= 0):
        matW, matWJ, vecWd = aPoint.getDerivatives(nJacobian) # W, W * J, W * d
841
        sign = 2 * nJacobian - 1
842
        if (nCurv > 0):
843
          aJacobian[labSlope:labOffset, 0:1] = -sign * vecWd # from curvature
844
          anIndex[0] = nLocals + 1
845
846
        aJacobian[labSlope:labOffset, index1:index1 + 2] = -sign * matWJ
        aJacobian[labSlope:labOffset, index2:index2 + 2] = sign * matW
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
        for i in range(nDim):
          anIndex[index2 + aDim[i]] = iOff2 + i  
          
#   local curvature
    if (labCurv >= 0):
      if (nCurv > 0):
        aJacobian[labCurv, labCurv] = 1.0  
                    
    return anIndex, aJacobian    
  
  def __getFitToKinkJacobian(self, aPoint):              
    ''' 
    Get jacobian for transformation from (trajectory) fit to kink parameters at point.

    Jacobian broken lines (q/p,..,u_i-1,u_i,u_i+1..) to kink (du') parameters.
    @param aPoint: point to use
    @type aPoint: GblPoint
    @return: labels for fit parameters with non zero derivatives, corresponding transformation matrix
    @rtype: list(vector(int), matrix(float)) 
    '''
    aDim = self.__dimensions
    nDim = len(aDim)
    nCurv = self.__numCurvature
    nLocals = self.__numLocals        
    nOffset = aPoint.getOffset() 
    anIndex = [0, 0, 0, 0, 0, 0, 0]
    aJacobian = np.zeros((2, 7))    

    prevW, prevWJ, prevWd = aPoint.getDerivatives(0) # W-, W- * J-, W- * d-
    nextW, nextWJ, nextWd = aPoint.getDerivatives(1) # W+, W+ * J+, W+ * d+
    sumWJ = prevWJ + nextWJ         # W- * J- + W+ * J+
    sumWd = prevWd + nextWd         # W+ * d+ + W- * d-
    iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1 # first offset ('i' in u_i)

#   local offset
    if (nCurv > 0):
      aJacobian[:, 0:1] = -sumWd # from curvature
      anIndex[0] = nLocals + 1      
    aJacobian[:, 1:3] = prevW # from 1st Offset
    aJacobian[:, 3:5] = -sumWJ # from 2nd Offset
    aJacobian[:, 5:7] = nextW # from 3rd Offset
    for i in range(nDim):
      anIndex[1 + aDim[i]] = iOff + i
      anIndex[3 + aDim[i]] = iOff + nDim + i
      anIndex[5 + aDim[i]] = iOff + nDim * 2 + i
        
    return anIndex, aJacobian    

895
896
897
898
899
900
901
902
903
  def getResults(self, aLabel):
    '''
    Get results (corrections, covarinace matrix) at point in forward or backward direction.
    
    @param aLabel: signed label of point (<0 backward, >0 forward)
    @type aLabel: int
    @return: correction vector, covarinace matrix for track parameters
    @rtype: list
    '''
904
    anIndex, aJacobian = self.__getJacobian(aLabel)
905
906
907
    nParBRL = len(anIndex)
    aVec = np.empty(nParBRL)
    for i in range(nParBRL):
908
      aVec[i] = self.__vector[anIndex[i] - 1]         # compressed vector
909
910
911
    aMat = self.__matrix.getBlockMatrix(anIndex)  # compressed matrix
    locPar = np.dot(aJacobian, aVec) 
    locCov = np.dot(aJacobian, np.dot(aMat, aJacobian.T))
912
    return  locPar, locCov 
913
914
915
916
917
918
919
920
921
922
  
  def fit(self, optionList=""):
    '''
    Perform fit of trajectory.
    
    @param optionList: M-estimators to be used (one iteration per character)
    @type optionList: string
    @return: Chi2, Ndf, loss of weight from fit ([0., -1, 0.] if fit failed)
    @rtype: list
    '''
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
 
    def defineOffsets():
      '''
      Define offsets from list of points.
      '''
# set labels for previous/next offsets
#     first point is offset    
      self.__points[0].setOffset(0)
      nOffsets = 1
#     intermediate scatterers are offsets    
      for aPoint in self.__points[1:-1]:
        if (aPoint.hasScatterer()):
          aPoint.setOffset(nOffsets)
          nOffsets += 1
        else:  
          aPoint.setOffset(-nOffsets)
#     last point is offset    
      self.__points[-1].setOffset(nOffsets)
      self.__numOffsets = nOffsets + 1
      self.__numParameters = self.__numOffsets * len(self.__dimensions) \
                           + self.__numCurvature + self.__numLocals
 
    def calcJacobians():
      '''
      Calculate Jacobians to previous/next scatterer from point to point ones.
      '''
      scatJacobian = np.empty((5, 5)) 
#     forward propagation (all)
      lastPoint = 0;
      numStep = 0;
      for iPoint in range(1, self.__numPoints):
        if (numStep == 0):
          scatJacobian = self.__points[iPoint].getP2pJacobian()
        else: 
          scatJacobian = np.dot(self.__points[iPoint].getP2pJacobian(), scatJacobian)
        numStep += 1
        self.__points[iPoint].addPrevJacobian(scatJacobian) # aPoint -> previous scatterer
        if (self.__points[iPoint].getOffset() >= 0):
          self.__points[lastPoint].addNextJacobian(scatJacobian) # lastPoint -> next scatterer
          numStep = 0;
          lastPoint = iPoint
#     backward propagation (without scatterers)
965
      for iPoint in range(self.__numPoints - 1, 0, -1):
966
967
        if (self.__points[iPoint].getOffset() >= 0):
          scatJacobian = self.__points[iPoint].getP2pJacobian()
968
          continue # skip offsets
969
        self.__points[iPoint].addNextJacobian(scatJacobian) # iPoint -> next scatterer
970
        scatJacobian = np.dot(scatJacobian, self.__points[iPoint].getP2pJacobian())
971
      
972
973
974
975
976
977
978
979
980
    def prepare():
      '''
      Prepare fit; generate data from points.
      '''
      aDim = self.__dimensions
# measurements
      for aPoint in self.__points:
        if (aPoint.hasMeasurement()):
          nLabel = aPoint.getLabel()
981
          measDim = aPoint.getMeasDim()
982
983
984
          localDer = aPoint.getLocalDerivatives()
          globalLab = aPoint.getGlobalLabels()
          globalDer = aPoint.getGlobalDerivatives()
985
986
987
988
989
990
991
992
993
          matP, aMeas, aPrec = aPoint.getMeasurement()
          labDer, matDer = self.__getFitToLocalJacobian(aPoint, measDim)
          matPDer = np.dot(matP, matDer)
          for i in range(measDim):
            if (aPrec[i] > 0.):
              aData = GblData(nLabel, aMeas[i], aPrec[i])
              aData.addDerivatives(i, labDer, matPDer, localDer, \
                                   globalLab, globalDer)
              self.__data.append(aData)
994
995
996
997
998
#                aPoint.setDataMeas(i, len(self.__data)) 
# pseudo measurements from kinks
      for aPoint in self.__points[1:-1]:
        if (aPoint.hasScatterer()):
          nLabel = aPoint.getLabel()        
999
1000
          aMeas, aPrec = aPoint.getScatterer()
          labDer, matDer = self.__getFitToKinkJacobian(aPoint)
For faster browsing, not all history is shown. View entire blame