#------------------------------
"""
:py:class:`UtilsJungfrau` contains utilities for Jungfrau detector correction
=============================================================================
Jungfrau gain range coding
bit: 15,14,...,0 Gain range, ind
0, 0 Normal, 0
0, 1 ForcedGain1, 1
1, 1 FixedGain2, 2
This software was developed for the SIT project.
If you use all or part of it, please give an appropriate acknowledgment.
Created on 2017-10-03 by Mikhail Dubrovin
"""
#------------------------------
import numpy as np
from time import time
from math import fabs
from Detector.GlobalUtils import print_ndarr, divide_protected
BW1 = 040000 # 16384 or 1<<14 (15-th bit starting from 1)
BW2 = 0100000 # 32768 or 2<<14 or 1<<15
BW3 = 0140000 # 49152 or 3<<14
MSK = 0x3fff # 16383 or (1<<14)-1 - 14-bit mask
#MSK = 037777 # 16383 or (1<<14)-1
#------------------------------
class Storage :
def __init__(self) :
self.offs = None
self.gfac = None
#------------------------------
store = Storage() # singleton
#------------------------------
#------------------------------
[docs]def calib_jungfrau(det, evt, src, cmpars=(7,3,100)) :
"""
Returns calibrated jungfrau data
- gets constants
- gets raw data
- evaluates (code - pedestal - offset)
- applys common mode correction if turned on
- apply gain factor
Parameters
- det (psana.Detector) - Detector object
- evt (psana.Event) - Event object
- src (psana.Source) - Source object
- cmpars (tuple) - common mode parameters
- cmpars[0] - algorithm # 7-for jungfrau
- cmpars[1] - control bit-word 1-in rows, 2-in columns
- cmpars[2] - maximal applied correction
"""
arr = det.raw(evt) # shape:(1, 512, 1024) dtype:uint16
#arr = np.array(det.raw(evt), dtype=np.float32)
peds = det.pedestals(evt) # - 4d pedestals shape:(3, 1, 512, 1024) dtype:float32
#mask = det.status_as_mask(evt, mode=0) # - 4d mask
gain = det.gain(evt) # - 4d gains
offs = det.offset(evt) # - 4d offset
cmp = det.common_mode(evt) if cmpars is None else cmpars
if gain is None : gain = np.ones_like(peds) # - 4d gains
if offs is None : offs = np.zeros_like(peds) # - 4d gains
gfac = store.gfac
if store.gfac is None :
store.gfac = gfac = divide_protected(np.ones_like(peds), gain)
#print_ndarr(cmp, 'XXX: common mode parameters ')
#print_ndarr(arr, 'XXX: calib_jungfrau arr ')
#print_ndarr(peds, 'XXX: calib_jungfrau peds')
#print_ndarr(gain, 'XXX: calib_jungfrau gain')
#print_ndarr(offs, 'XXX: calib_jungfrau offs')
# make bool arrays of gain ranges shaped as data
#abit15 = arr & BW1 # ~0.5ms
#abit16 = arr & BW2
#gr0 = np.logical_not(arr & BW3) #1.6ms
#gr1 = np.logical_and(abit15, np.logical_not(abit16))
#gr2 = np.logical_and(abit15, abit16)
#pro_den = np.select((den!=0,), (den,), default=1)
# Define bool arrays of ranges
# faster than bit operations
gr0 = arr < BW1 # 490 us
gr1 =(arr >= BW1) & (arr<BW2) # 714 us
gr2 = arr >= BW3 # 400 us
#print_ndarr(gr0, 'XXX: calib_jungfrau gr0')
#print_ndarr(gr1, 'XXX: calib_jungfrau gr1')
#print_ndarr(gr2, 'XXX: calib_jungfrau gr2')
# Subtract pedestals
arrf = np.array(arr & MSK, dtype=np.float32)
arrf[gr0] -= peds[0,gr0]
arrf[gr1] -= peds[1,gr1] #- arrf[gr1]
arrf[gr2] -= peds[2,gr2] #- arrf[gr2]
#a = np.select((gr0, gr1, gr2), (gain[0,:], gain[1,:], gain[2,:]), default=1) # 2msec
factor = np.select((gr0, gr1, gr2), (gfac[0,:], gfac[1,:], gfac[2,:]), default=1) # 2msec
offset = np.select((gr0, gr1, gr2), (offs[0,:], offs[1,:], offs[2,:]), default=0)
#print_ndarr(factor, 'XXX: calib_jungfrau factor')
#print_ndarr(offset, 'XXX: calib_jungfrau offset')
# Apply offset
arrf -= offset
t0_sec = time()
#if False :
if cmp is not None :
mode, cormax = int(cmp[1]), cmp[2]
#common_mode_2d(arrf, mask=gr0, cormax=cormax)
for s in range(arrf.shape[0]) :
if mode & 1 :
common_mode_rows(arrf[s,], mask=gr0[s,], cormax=cormax)
#common_mode_rows(arrf[s,:,:256], mask=gr0[s,:,:256], cormax=cormax)
#common_mode_rows(arrf[s,:,256:], mask=gr0[s,:,256:], cormax=cormax)
#common_mode_cols(arrf[s,:512,:], mask=gr0[s,:512,:], cormax=cormax)
if mode & 2 :
common_mode_cols(arrf[s,], mask=gr0[s,], cormax=cormax)
pass
#print '\nXXX: CM consumed time (sec) =', time()-t0_sec # 90-100msec total
# Apply gain
return arrf * factor
#------------------------------
[docs]def common_mode_rows(arr, mask=None, cormax=None) :
"""Defines and applys common mode correction to 2-d arr using the same shape mask in loop for rows.
"""
rows, cols = arr.shape
#print 'XXX.common_mode_rows.arr.shape', arr.shape
for r in range(rows) :
cmode = 0
cmode = np.median(arr[r,:][mask[r,:]]) if mask is not None else\
np.median(arr[r,:])
#print 'XXX.common_mode_2 cmode=%.3f, len(selected arr)=%d' % (cmode, len(arr[r,:][mask[r,:]]))
if cormax is None or fabs(cmode) < cormax :
if mask is not None : arr[r,:][mask[r,:]] -= cmode
else : arr[r,:] -= cmode
else :
return
#------------------------------
[docs]def common_mode_cols(arr, mask=None, cormax=None) :
"""Defines and applys common mode correction to 2-d arr using the same shape mask in loop for cols.
"""
rows, cols = arr.shape
#print 'XXX.common_mode_cols.arr.shape', arr.shape
for c in range(cols) :
cmode = 0
cmode = np.median(arr[:,c][mask[:,c]]) if mask is not None else\
np.median(arr[:,c])
#print 'XXX.common_mode_2 cmode=%.3f, len(selected arr)=%d' % (cmode, len(arr[:,c][mask[:,c]]))
if cormax is None or fabs(cmode) < cormax :
if mask is not None : arr[:,c][mask[:,c]] -= cmode
else : arr[:,c] -= cmode
else :
return
#------------------------------
[docs]def common_mode_2d(arr, mask=None, cormax=None) :
"""Defines and applys common mode correction to entire 2-d arr using the same shape mask.
"""
cmode = 0
cmode = np.median(arr[mask>0]) if mask is not None else\
np.median(arr)
#print 'XXX.common_mode_2 cmode=%.3f, len(selected arr)=%d' % (cmode, len(arr[mask>0]))
if cormax is None or fabs(cmode) < cormax :
if mask is not None : arr[mask>0] -= cmode
else : arr -= cmode
else :
return
#------------------------------
[docs]def common_mode_jungfrau(frame, cormax) :
"""
Parameters
- frame (np.array) - shape=(512, 1024)
"""
intmax = 100
rows = 512
cols = 1024
banks = 4
bsize = cols/banks
for r in range(rows):
col0 = 0
for b in range(banks):
try:
cmode = np.median(frame[r, col0:col0+bsize][frame[r, col0:col0+bsize]<cormax])
if not np.isnan(cmode):
## e.g. found no pixels below intmax
## print r, cmode, col0, b, bsize
if cmode<cormax-1 :
frame[r, col0:col0+bsize] -= cmode
except:
cmode = -666
print "cmode problem"
print frame[r, col0:col0 + bsize]
col0 += bsize
#------------------------------
#------------------------------
#------------------------------
#------------------------------
#------------------------------