gbltst.py

Go to the documentation of this file.

'''
Simple Test Program for General Broken Lines.
 
Created on Jul 27, 2011
 
@author: kleinwrt
'''
 
 
 
import numpy as np
import math
import time
from gblfit import GblPoint, GblTrajectory
#
 
 
 
def example1():
 
#
  np.random.seed(47117)
 
  nTry = 1000  #: number of tries
  nLayer = 5  #: number of detector layers
  useThickScatterer = True
  print " Gbltst $Id$ ", nTry, nLayer, useThickScatterer
  start = time.clock()
# track direction
  sinLambda = 0.3
  cosLambda = math.sqrt(1.0 - sinLambda ** 2)
  sinPhi = 0.
  cosPhi = math.sqrt(1.0 - sinPhi ** 2)
# 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]])
# measurement resolution
  measErr = np.array([ 0.001, 0.001])  # 10 mu
  measPrec = 1.0 / measErr ** 2
# scattering error
  scatErr = 0.001  # 1 mrad
# scattering cov matrix for thick scatterers
  scatCov = np.zeros((4, 4))  
# RMS of CurviLinear track parameters (Q/P, slopes, offsets)
  clErr = np.array([0.001, -0.1, 0.2, -0.15, 0.25])
  # precision matrix for external seed (in local system)
  locSeed = None
  seedLabel = 0  # label of point with seed
  if seedLabel != 0:
    print " external seed at label ", seedLabel
#
  bfac = 0.2998  # Bz*c for Bz=1
  step = 1.5 / cosLambda  # constant steps in RPhi
#
  Chi2Sum = 0.
  NdfSum = 0
  LostSum = 0.
#
  binaryFile = open("milleBinaryISN.dat", "wb")
#
  for iTry in range(nTry):
# generate (CurviLinear) track parameters 
    clNorm = np.random.normal(0., 1., 5)  
    clPar = clErr * clNorm
    # covariance matrix
    clCov = np.eye(5)
    for i in range(5):
      clCov[i, i] = clErr[i] ** 2
# arclength
    s = 0.
# point-to-point jacobian (from previous point)    
    jacPointToPoint = np.eye(5)
# additional (local or global) derivatives    
    addDer = np.array([[1.0], [0.0]])
    labGlobal = np.array([[4711], [4711]])
# create trajectory
    traj = GblTrajectory(bfac != 0.)
 
# at previous point: transformation from local to curvilinear system 
    oldL2c = np.eye(5)
   
    for iLayer in range(nLayer):
      #print " layer ", iLayer
# measurement directions (J,K) from stereo angle
      sinStereo = (0. if iLayer % 2 == 0 else 0.1) 
      cosStereo = math.sqrt(1.0 - sinStereo ** 2)    
# (orthogonal) measurement system: I, J ,K    
      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()
# measurement - prediction in measurement system with error
      measNorm = np.random.normal(0., 1., 2)  
      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)
      point.addMeasurement([proL2m, meas, measPrec])
# additional local parameters?
#      point.addLocals(addDer)
# additional global parameters?
      point.addGlobals(labGlobal, addDer)
# add thick scatterer
      if useThickScatterer:
        if scatCov[0, 0] > 0:
          scat = np.zeros(4)
          scatPrec = np.linalg.inv(scatCov)
          #print " scatCov ", iLayer, scatCov
          #print " scatPrec ", iLayer, scatPrec
          point.addScatterer([scat, scatPrec])
        scatCov = np.zeros((4, 4))
# add 'noise', 10%
      measNorm = np.random.uniform(-5., 5., 2)  
      if np.random.uniform() < 0.:
        meas = np.dot(proC2m, clPar[3:5]) + measErr * measNorm 
        print " noise ", iLayer
        point.addMeasurement([proL2m, meas, measPrec])
        point.addGlobals(labGlobal, addDer)
# locDer flips sign every measurement      
      addDer = -addDer        
# add point to trajectory      
      iLabel = traj.addPoint(point)
      #print " meas. ", iLabel, s
      if iLabel == abs(seedLabel):
        clSeed = np.linalg.inv(clCov)
        locSeed = np.dot(proL2c, np.dot(clSeed, proL2c.T))
# propagate to scatterer
      jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac)
      clPar = np.dot(jacPointToPoint, clPar)
      clCov = np.dot(jacPointToPoint, np.dot(clCov, jacPointToPoint.T))
      s += step
# propagate scatCov
      scatCov = np.dot (jacPointToPoint[1:, 1:], np.dot(scatCov, jacPointToPoint[1:, 1:].T))      
      if (iLayer < nLayer - 1):
        scat = np.array([0., 0.])
# point with scatterer
        jac = np.dot(proC2l, np.dot(jacPointToPoint, proL2c))
        point = GblPoint(jac)
# add thin scatterer        
        if not useThickScatterer:
          if scatErr > 0:
            scatP = local.getScatPrecision(scatErr)
            point.addScatterer([scat, scatP])        
        iLabel = traj.addPoint(point)
        #print " scat1 ", iLabel, s
        if iLabel == abs(seedLabel):
          clSeed = np.linalg.inv(clCov)
          locSeed = np.dot(proL2c, np.dot(clSeed, proL2c.T))
 
# scatter a little    
        scatNorm = np.random.normal(0., 1., 2)  
        clPar[1:3] = clPar[1:3] + scatErr * scatNorm
        scatCov[0, 0] += scatErr ** 2
        scatCov[1, 1] += scatErr ** 2
# propagate to 2nd scatterer
      jacPointToPoint = gblSimpleJacobian(step, cosLambda, bfac)
      clPar = np.dot(jacPointToPoint, clPar)
      clCov = np.dot(jacPointToPoint, np.dot(clCov, jacPointToPoint.T))
      s += step
# propagate scatCov
      scatCov = np.dot (jacPointToPoint[1:, 1:], np.dot(scatCov, jacPointToPoint[1:, 1:].T))        
      if (iLayer < nLayer - 1):
        scat = np.array([0., 0.])
# point with scatterer
        jac = np.dot(proC2l, np.dot(jacPointToPoint, proL2c))
        point = GblPoint(jac)
# add thin scatterer        
        if not useThickScatterer:
          if scatErr > 0:
            scatP = local.getScatPrecision(scatErr)
            point.addScatterer([scat, scatP])
        iLabel = traj.addPoint(point)
        #print " scat2 ", iLabel, s
        if iLabel == abs(seedLabel):
          clSeed = np.linalg.inv(clCov)
          locSeed = np.dot(proL2c, np.dot(clSeed, proL2c.T))
 
# scatter a little    
        scatNorm = np.random.normal(0., 1., 2)  
        clPar[1:3] = clPar[1:3] + scatErr * scatNorm  
        scatCov[0, 0] += scatErr ** 2
        scatCov[1, 1] += scatErr ** 2      
# propagate to next measurement layer    
        clPar = np.dot(jacPointToPoint, clPar)
        clCov = np.dot(jacPointToPoint, np.dot(clCov, jacPointToPoint.T))
        s += step
        # propagate scatCov
        scatCov = np.dot (jacPointToPoint[1:, 1:], np.dot(scatCov, jacPointToPoint[1:, 1:].T))  
      oldL2c = proL2c
 
# add external seed
    if locSeed is not None: 
      traj.addExternalSeed(seedLabel, locSeed)
  
# fit trajectory
    Chi2, Ndf, Lost = traj.fit()
    #print " Record, Chi2, Ndf, Lost", iTry, Chi2, Ndf, Lost
# dump trajectory
    #traj.printPoints()
    #traj.printData()
# write to MP binary file    
    traj.milleOut(binaryFile)
# sum up    
    Chi2Sum += Chi2
    NdfSum += Ndf
    LostSum += Lost
# get corrections and covariance matrix at points 
    if (iTry == 0):
      for i in range(1, nLayer + 1): 
        locPar, locCov = traj.getResults(-i)
        print " >Point ", i
        print " locPar ", locPar
        #print " locCov ", locCov      
        locPar, locCov = traj.getResults(i)
        print " Point> ", i
        print " locPar ", locPar
        #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]   
#
  end = time.clock()
  print " Time [s] ", end - start
  print " Chi2Sum/NdfSum ", Chi2Sum / NdfSum
  print " LostSum/nTry ", LostSum / nTry
 
 
 
def example2():
#  
  binaryFile = open("milleBinaryISN.dat", "rb")
  nRec = 0
  maxRec = 10  #: maximum number of records to read
  Chi2Sum = 0.
  NdfSum = 0
  LostSum = 0.
  start = time.clock()
  
  try:
    while(nRec < maxRec):
# create trajectory
      traj = GblTrajectory(0)
# read from file      
      traj.milleIn(binaryFile)  # get data blocks from file
      nRec += 1
# fit trajectory      
      Chi2, Ndf, Lost = traj.fit()
      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
 
 
 
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     
 
 
 
class gblMeasSystem(object):
 
  
  def __init__(self, measDir, curviDir):
    
    self.__prod = np.dot(curviDir, measDir.T)
 
  
  def getTransLocalToMeas(self):
    return None
 
  
  def getTransLocalToCurvi(self):
    meas2crv = np.zeros((5, 5))
    meas2crv[0, 0] = 1.
    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)
    return meas2crv
 
  
  def getTransCurviToMeas(self):
    return np.linalg.inv(self.__prod[1:3, 1:3])  #((U,V)*(J,K))^-1
 
  
  def getScatPrecision(self, scatErr):
    c1 = self.__prod[0, 1]  # T*J
    c2 = self.__prod[0, 2]  # T*K
    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       
 
 
 
class gblCurviSystem(object):
 
  
  def __init__(self, measDir, curviDir):
    
    self.__c2m = np.linalg.inv(np.dot(curviDir[1:3, 1:3], measDir[1:3, 1:3].T))
 
  
  def getTransLocalToMeas(self):
    return self.__c2m  #((U,V)*(J,K))^-1
 
  
  def getTransLocalToCurvi(self):
    return np.eye(5)  # unit matrix 
 
  
  def getTransCurviToMeas(self):
    return self.__c2m  #((U,V)*(J,K))^-1
 
  
  def getScatPrecision(self, scatErr):
    return np.array([1., 1.]) / (scatErr * scatErr)  # diagonal only     
 
        
# 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
example2()