442 lines
23 KiB
Python
442 lines
23 KiB
Python
|
|
import os
|
|
import numpy as np
|
|
from ROOT import TH1D, TH2D, TCanvas, TF1, TFile
|
|
from multiprocessing import Pool
|
|
from array import array
|
|
import h5py
|
|
from multiprocessing import Pool
|
|
import aare
|
|
|
|
_cfg = {} ### global config
|
|
_aduToKev3DMap = None ### global caliFile
|
|
_pedestalAduFrame = None ### global pedestalAduFrame
|
|
_noiseEneFrame = None ### global noiseEneFrame
|
|
nChunks = 16
|
|
clusterSize2Photon = 7
|
|
clusterSize1Photon = 3
|
|
|
|
def init(cfg):
|
|
global _cfg
|
|
global _aduToKev3DMap
|
|
_cfg = cfg
|
|
NX, NY = cfg['NX'], cfg['NY']
|
|
St_x, St_y = cfg.get('StCorner', [0, 0]) ### with default corner at (0, 0)
|
|
_aduToKev3DMap = np.load(cfg['caliFileName'])[:, St_y:St_y+NY, St_x:St_x+NX]
|
|
Roi = cfg['Roi']
|
|
global nChunks
|
|
nChunks = _cfg.get('NChunks', 16)
|
|
|
|
def convertAduFrameToEnergyFrame(aduFrame):
|
|
if _aduToKev3DMap is None:
|
|
return aduFrame
|
|
NX, NY = _cfg['NX'], _cfg['NY']
|
|
if _aduToKev3DMap.shape[0] == 1640:
|
|
idxFrame = (np.abs(aduFrame//10)).astype(np.int32)
|
|
idxFrame[idxFrame > _aduToKev3DMap.shape[0]-2] = _aduToKev3DMap.shape[0]-2
|
|
NY, NX = _aduToKev3DMap.shape[1:3]
|
|
_energyFrame0 = _aduToKev3DMap[idxFrame, np.arange(NY).reshape(-1, 1), np.arange(NX)]
|
|
_energyFrame1 = _aduToKev3DMap[idxFrame+1, np.arange(NY).reshape(-1, 1), np.arange(NX)]
|
|
energyFrame = _energyFrame0 + (_energyFrame1 - _energyFrame0)/10 * (np.abs(aduFrame) - idxFrame*10)
|
|
energyFrame *= np.sign(aduFrame)
|
|
|
|
else:
|
|
aduPerBin = 10
|
|
nExtraBins = 1000 // aduPerBin
|
|
idxFrame = (aduFrame//aduPerBin).astype(np.int32) + nExtraBins
|
|
idxFrame[idxFrame > _aduToKev3DMap.shape[0]-2] = _aduToKev3DMap.shape[0]-2
|
|
NY, NX = _aduToKev3DMap.shape[1:3]
|
|
_energyFrame0 = _aduToKev3DMap[idxFrame, np.arange(NY).reshape(-1, 1), np.arange(NX)]
|
|
_energyFrame1 = _aduToKev3DMap[idxFrame+1, np.arange(NY).reshape(-1, 1), np.arange(NX)]
|
|
energyFrame = _energyFrame0 + (_energyFrame1 - _energyFrame0)/aduPerBin * (aduFrame - (idxFrame-nExtraBins)*aduPerBin)
|
|
if 'energyScalingFactor' in _cfg:
|
|
energyFrame *= _cfg['energyScalingFactor']
|
|
return energyFrame
|
|
|
|
def bookHistograms(energy, suffix='', energyBinWidth=0.1, isMC=False):
|
|
histMaxEnergy = energy * 5.5
|
|
nEnergyBins = int(histMaxEnergy // energyBinWidth)
|
|
roi_x0, roi_x1, roi_y0, roi_y1 = _cfg['Roi'][0], _cfg['Roi'][1], _cfg['Roi'][2], _cfg['Roi'][3]
|
|
roi_width = roi_x1 - roi_x0
|
|
roi_height = roi_y1 - roi_y0
|
|
interpolationBins = 10
|
|
|
|
histograms = {}
|
|
histograms['h1_ClusterEnergy'] = TH1D(
|
|
f'h1_ClusterEnergy{suffix}', f'h1_ClusterEnergy{suffix}', nEnergyBins, -1, histMaxEnergy
|
|
)
|
|
histograms['h1_ClusterEnergy'].SetTitle('Cluster Energy;Energy [keV];Counts')
|
|
|
|
histograms['h1_PixelEnergy'] = TH1D(
|
|
f'h1_PixelEnergy{suffix}', f'h1_PixelEnergy{suffix}', nEnergyBins, -1, histMaxEnergy
|
|
)
|
|
histograms['h1_PixelEnergy'].SetTitle('Pixel Energy;Energy [keV];Counts')
|
|
|
|
histograms['h1_1PhotonClusterSize'] = TH1D(
|
|
f'h1_1PhotonClusterSize{suffix}', f'h1_1PhotonClusterSize{suffix}', 30, 0, 30
|
|
)
|
|
histograms['h1_1PhotonClusterSize'].SetTitle('1-Photon Cluster Size;Cluster Size [pixel];Counts')
|
|
|
|
histograms['h1_1PhotonClusterDiameterX'] = TH1D(
|
|
f'h1_1PhotonClusterDiameterX{suffix}', f'h1_1PhotonClusterDiameterX{suffix}', 10, 0, 10
|
|
)
|
|
histograms['h1_1PhotonClusterDiameterX'].SetTitle('1-Photon Cluster Diameter X;Cluster Diameter X [pixel];Counts')
|
|
|
|
histograms['h1_1PhotonClusterDiameterY'] = TH1D(
|
|
f'h1_1PhotonClusterDiameterY{suffix}', f'h1_1PhotonClusterDiameterY{suffix}', 10, 0, 10
|
|
)
|
|
histograms['h1_1PhotonClusterDiameterY'].SetTitle('1-Photon Cluster Diameter Y;Cluster Diameter Y [pixel];Counts')
|
|
|
|
histograms['h1_1PhotonClusterPixelEnergy'] = TH1D(
|
|
f'h1_1PhotonClusterPixelEnergy{suffix}', f'h1_1PhotonClusterPixelEnergy{suffix}', nEnergyBins, -1, energy*1.5
|
|
)
|
|
histograms['h1_1PhotonClusterPixelEnergy'].SetTitle('1-Photon Cluster Pixel Energy;Energy [keV];Counts')
|
|
|
|
histograms['h1_1PhotonClusterEnergy'] = TH1D(
|
|
f'h1_1PhotonClusterEnergy{suffix}', f'h1_1PhotonClusterEnergy{suffix}', nEnergyBins, -1, energy*1.5
|
|
)
|
|
histograms['h1_1PhotonClusterEnergy'].SetTitle('1-Photon Cluster Energy;Energy [keV];Counts')
|
|
|
|
histograms['h1_1PhotonEta_X'] = TH1D(
|
|
f'h1_1PhotonEta_X{suffix}', f'h1_1PhotonEta_X{suffix}', 100, -0.5, 0.5
|
|
)
|
|
histograms['h1_1PhotonEta_X'].SetTitle('1-Photon Eta X;Eta X;Counts')
|
|
|
|
histograms['h1_1PhotonEta_Y'] = TH1D(
|
|
f'h1_1PhotonEta_Y{suffix}', f'h1_1PhotonEta_Y{suffix}', 100, -0.5, 0.5
|
|
)
|
|
histograms['h1_1PhotonEta_Y'].SetTitle('1-Photon Eta Y;Eta Y;Counts')
|
|
|
|
histograms['h1_2PhotonClusterEnergy'] = TH1D(
|
|
f'h1_2PhotonClusterEnergy{suffix}', f'h1_2PhotonClusterEnergy{suffix}', nEnergyBins, -1, energy*4
|
|
)
|
|
histograms['h1_2PhotonClusterEnergy'].SetTitle('2-Photon Cluster Energy;Energy [keV];Counts')
|
|
|
|
histograms['h1_2PhotonClusterSize'] = TH1D(
|
|
f'h1_2PhotonClusterSize{suffix}', f'h1_2PhotonClusterSize{suffix}', 30, 0, 30
|
|
)
|
|
histograms['h1_2PhotonClusterSize'].SetTitle('2-Photon Cluster Size;Cluster Size [pixel];Counts')
|
|
|
|
histograms['h1_2PhotonClusterDiameterX'] = TH1D(
|
|
f'h1_2PhotonClusterDiameterX{suffix}', f'h1_2PhotonClusterDiameterX{suffix}', 10, 0, 10
|
|
)
|
|
histograms['h1_2PhotonClusterDiameterX'].SetTitle('2-Photon Cluster Diameter X;Cluster Diameter X [pixel];Counts')
|
|
|
|
histograms['h1_2PhotonClusterDiameterY'] = TH1D(
|
|
f'h1_2PhotonClusterDiameterY{suffix}', f'h1_2PhotonClusterDiameterY{suffix}', 10, 0, 10
|
|
)
|
|
histograms['h1_2PhotonClusterDiameterY'].SetTitle('2-Photon Cluster Diameter Y;Cluster Diameter Y [pixel];Counts')
|
|
|
|
histograms['h1_2PhotonClusterPixelEnergy'] = TH1D(
|
|
f'h1_2PhotonClusterPixelEnergy{suffix}', f'h1_2PhotonClusterPixelEnergy{suffix}', nEnergyBins, -1, energy*2.5
|
|
)
|
|
histograms['h1_2PhotonClusterPixelEnergy'].SetTitle('2-Photon Cluster Pixel Energy;Energy [keV];Counts')
|
|
|
|
histograms['h2_SumFrame'] = TH2D(
|
|
f'h2_SumFrame{suffix}', f'h2_SumFrame{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1
|
|
)
|
|
histograms['h2_SumFrame'].SetTitle('Sum Frame;X [pixel];Y [pixel];Energy Sum [keV]')
|
|
|
|
histograms['h2_HitsSumFrame'] = TH2D(
|
|
f'h2_HitsSumFrame{suffix}', f'h2_HitsSumFrame{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1
|
|
)
|
|
histograms['h2_HitsSumFrame'].SetTitle('Hits Sum Frame;X [pixel];Y [pixel];Energy Sum [keV]')
|
|
|
|
histograms['h2_1PhotonSumFrame'] = TH2D(
|
|
f'h2_1PhotonSumFrame{suffix}', f'h2_1PhotonSumFrame{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1
|
|
)
|
|
histograms['h2_1PhotonSumFrame'].SetTitle('1-Photon Hits Sum Frame;X [pixel];Y [pixel];Energy Sum [keV]')
|
|
|
|
histograms['h2_1PhotonSumSample'] = TH2D(
|
|
f'h2_1PhotonSumSample{suffix}', f'h2_1PhotonSumSample{suffix}', clusterSize1Photon, 0, clusterSize1Photon, clusterSize1Photon, 0, clusterSize1Photon
|
|
)
|
|
histograms['h2_1PhotonSumSample'].SetTitle('1-Photon Sum Sample;Sum Sample;Counts')
|
|
|
|
histograms['h2_2PhotonSumSample'] = TH2D(
|
|
f'h2_2PhotonSumSample{suffix}', f'h2_2PhotonSumSample{suffix}', clusterSize2Photon, 0, clusterSize2Photon, clusterSize2Photon, 0, clusterSize2Photon
|
|
)
|
|
histograms['h2_2PhotonSumSample'].SetTitle('2-Photon Sum Sample;Sum Sample;Counts')
|
|
|
|
histograms['h2_2PhotonSumFrame'] = TH2D(
|
|
f'h2_2PhotonSumFrame{suffix}', f'h2_2PhotonSumFrame{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1
|
|
)
|
|
histograms['h2_2PhotonSumFrame'].SetTitle('2-Photon Hits Sum Frame;X [pixel];Y [pixel];Energy Sum [keV]')
|
|
|
|
histograms['h2_ClusterSizeVsEnergy'] = TH2D(
|
|
f'h2_ClusterSizeVsEnergy{suffix}', f'h2_ClusterSizeVsEnergy{suffix}', nEnergyBins, -1, histMaxEnergy, 30, 0, 30
|
|
)
|
|
histograms['h2_ClusterSizeVsEnergy'].SetTitle('Cluster Size vs Energy;Energy [keV];Cluster Size [pixel]')
|
|
|
|
histograms['h2_DiameterVsEnergy'] = TH2D(
|
|
f'h2_DiameterVsEnergy{suffix}', f'h2_DiameterVsEnergy{suffix}', nEnergyBins, -1, histMaxEnergy, 10, 0, 10
|
|
)
|
|
histograms['h2_DiameterVsEnergy'].SetTitle('Cluster Diameter vs Energy;Energy [keV];Cluster Diameter [pixel]')
|
|
|
|
for hist in histograms.values():
|
|
hist.SetDirectory(0)
|
|
|
|
return histograms
|
|
|
|
def _processFrames(idxChunk): ### for both single and double photon events, using
|
|
startFrame, endFrame = _cfg['NFrame'] * idxChunk // nChunks, min(_cfg['NFrame'] * (idxChunk + 1) // nChunks, _cfg['NFrame'])
|
|
print(f'Processing frames from {startFrame} to {endFrame}...')
|
|
Energy = _cfg['energy']
|
|
selectionRange = _cfg['selectionRange']
|
|
signalFileNames = _cfg['signalFileNames']
|
|
NX = _cfg['NX']
|
|
NY = _cfg['NY']
|
|
headerSize = _cfg['headerSize']
|
|
Roi = _cfg['Roi']
|
|
nFramePerFile = os.path.getsize(f'{signalFileNames[0]}') // (NX * NY + 56) // 2
|
|
_hists = bookHistograms(Energy, f'_{idxChunk}')
|
|
sumFrame = np.zeros((NY, NX), dtype=np.float64)
|
|
|
|
### create a cluster finder
|
|
nxROI = Roi[1] - Roi[0]
|
|
nyROI = Roi[3] - Roi[2]
|
|
CF = aare.VarClusterFinder((nyROI, nxROI), 5) ### arg1: frame size, arg2: threshold
|
|
CF.set_peripheralThresholdFactor(3)
|
|
CF.set_noiseMap(_noiseEneFrame)
|
|
|
|
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
|
|
h5_1Photon_file = h5py.File(f'{_cfg["outputFolder"]}/1Photon_CS{clusterSize1Photon}_chunk{idxChunk}.h5', 'w')
|
|
dset_1Photon_clusters = h5_1Photon_file.create_dataset(
|
|
'clusters', (0, clusterSize1Photon, clusterSize1Photon), maxshape=(None, clusterSize1Photon, clusterSize1Photon), dtype='f4',
|
|
chunks=True, compression='gzip'
|
|
)
|
|
dset_1photon_refs = h5_1Photon_file.create_dataset(
|
|
'referencePoint', (0, 2), maxshape=(None, 2), dtype='i4',
|
|
chunks=True
|
|
)
|
|
|
|
cluster_1Photon_list = []
|
|
refpoint_1Photon_list = []
|
|
|
|
h5_2Photon_file = h5py.File(f'{_cfg["outputFolder"]}/2Photon_CS{clusterSize2Photon}_chunk{idxChunk}.h5', 'w')
|
|
dset_2Photon_clusters = h5_2Photon_file.create_dataset(
|
|
'clusters', (0, clusterSize2Photon, clusterSize2Photon), maxshape=(None, clusterSize2Photon, clusterSize2Photon), dtype='f4',
|
|
chunks=True, compression='gzip'
|
|
)
|
|
dset_2Photon_refs = h5_2Photon_file.create_dataset(
|
|
'referencePoint', (0, 2), maxshape=(None, 2), dtype='i4',
|
|
chunks=True
|
|
)
|
|
cluster_2Photon_list = []
|
|
refpoint_2Photon_list = []
|
|
BUFFER_THRESHOLD = 10000
|
|
def flush_h5_buffer():
|
|
if not cluster_1Photon_list: return # no new clusters to write
|
|
old_len = dset_1Photon_clusters.shape[0]
|
|
new_len = old_len + len(cluster_1Photon_list)
|
|
dset_1Photon_clusters.resize((new_len, clusterSize1Photon, clusterSize1Photon))
|
|
dset_1photon_refs.resize((new_len, 2))
|
|
dset_1Photon_clusters[old_len:new_len] = cluster_1Photon_list
|
|
dset_1photon_refs[old_len:new_len] = refpoint_1Photon_list
|
|
cluster_1Photon_list.clear()
|
|
refpoint_1Photon_list.clear()
|
|
|
|
old_len = dset_2Photon_clusters.shape[0]
|
|
new_len = old_len + len(cluster_2Photon_list)
|
|
dset_2Photon_clusters.resize((new_len, clusterSize2Photon, clusterSize2Photon))
|
|
dset_2Photon_refs.resize((new_len, 2))
|
|
dset_2Photon_clusters[old_len:new_len] = cluster_2Photon_list
|
|
dset_2Photon_refs[old_len:new_len] = refpoint_2Photon_list
|
|
cluster_2Photon_list.clear()
|
|
refpoint_2Photon_list.clear()
|
|
|
|
for idxFrame in range(startFrame, endFrame):
|
|
idxFile, idxFrame = divmod(idxFrame, nFramePerFile)
|
|
try:
|
|
if 'nFiber' in _cfg and _cfg['nFiber'] == 2:
|
|
offset = (headerSize + idxFrame * (NX * NY//2 + headerSize)) * 2
|
|
fileName1 = signalFileNames[idxFile]
|
|
fileName2 = fileName1.replace('d0', 'd1')
|
|
_signalAduFrame1 = np.fromfile(f'{fileName1}', dtype=np.uint16, offset=offset, count=NX * NY //2).astype(np.int32).reshape(NY//2, NX)
|
|
_signalAduFrame2 = np.fromfile(f'{fileName2}', dtype=np.uint16, offset=offset, count=NX * NY //2).astype(np.int32).reshape(NY//2, NX)
|
|
signalAduFrame = np.concatenate([_signalAduFrame1, _signalAduFrame2], axis=0)
|
|
del _signalAduFrame1, _signalAduFrame2
|
|
else:
|
|
offset = (headerSize + idxFrame * (NX * NY + headerSize)) * 2
|
|
signalAduFrame = np.fromfile(f'{signalFileNames[idxFile]}', dtype=np.uint16, offset=offset, count=NX * NY).astype(np.int32).reshape(NX, NY)
|
|
except Exception as e:
|
|
print(f'Error reading frame {idxFrame+startFrame} from file {signalFileNames[idxFile]}: {e}, stop processing this chunk ({startFrame}-{endFrame})')
|
|
return _hists
|
|
|
|
signalAduFrame = signalAduFrame - _pedestalAduFrame
|
|
signalEneFrame = convertAduFrameToEnergyFrame(signalAduFrame)
|
|
sumFrame += signalEneFrame
|
|
|
|
signalEneFrame = signalEneFrame[Roi[2]:Roi[3], Roi[0]:Roi[1]]
|
|
_signalEneFrame = np.copy(signalEneFrame) ### make a copy of the original energy frame for cluster energy correction later, since the cluster finder will modify the input energy frame
|
|
CF.find_clusters_X(signalEneFrame)
|
|
clusters = CF.hits()
|
|
for idxCluster in range(len(clusters)):
|
|
clusterSize = clusters['size'][idxCluster]
|
|
energy = clusters['energy'][idxCluster]
|
|
_hists['h1_ClusterEnergy'].Fill(energy)
|
|
xs = clusters['cols'][idxCluster][:clusterSize]
|
|
ys = clusters['rows'][idxCluster][:clusterSize]
|
|
enes = clusters['enes'][idxCluster][:clusterSize]
|
|
|
|
### single photon events
|
|
if Energy - selectionRange < energy < Energy + selectionRange:
|
|
ref_x = clusters['col'][idxCluster] - clusterSize1Photon//2 ### refered to the lower-left corner of the cluster
|
|
ref_y = clusters['row'][idxCluster] - clusterSize1Photon//2
|
|
if int(ref_x) < 0 or int(ref_x)+clusterSize1Photon > signalEneFrame.shape[1] or int(ref_y) < 0 or int(ref_y)+clusterSize1Photon > signalEneFrame.shape[0]:
|
|
continue ### skip clusters too close to the border
|
|
for i in range(len(xs)):
|
|
_hists['h2_1PhotonSumFrame'].Fill(xs[i]+Roi[0], ys[i]+Roi[2], enes[i])
|
|
_hists['h1_1PhotonClusterSize'].Fill(clusterSize)
|
|
diameterX = max(xs) - min(xs) + 1
|
|
diameterY = max(ys) - min(ys) + 1
|
|
_hists['h1_1PhotonClusterDiameterX'].Fill(diameterX)
|
|
_hists['h1_1PhotonClusterDiameterY'].Fill(diameterY)
|
|
|
|
cluster_1Photon = np.zeros((clusterSize1Photon, clusterSize1Photon), dtype=np.float32)
|
|
### var cluster is less accurate than the 3x3 fixed cluster, since the zero pixels are incontinuous, not friendly for ML.
|
|
# for i in range(len(xs)):
|
|
# x_rel = xs[i] - int(ref_x)
|
|
# y_rel = ys[i] - int(ref_y)
|
|
# # _hists['h1_1PhotonClusterPixelEnergy'].Fill(enes[i])
|
|
# if 0 <= x_rel < clusterSize1Photon and 0 <= y_rel < clusterSize1Photon:
|
|
# cluster_1Photon[y_rel, x_rel] = enes[i]
|
|
|
|
# use the sum of 3x3 cluster as the cluster energy for single photon events, which is more accurate than the original cluster energy given by the cluster finder, since the cluster finder may miss some pixels with low energy far from the center pixel
|
|
cluster_1Photon = _signalEneFrame[int(ref_y):int(ref_y)+clusterSize1Photon, int(ref_x):int(ref_x)+clusterSize1Photon]
|
|
_hists['h1_1PhotonClusterEnergy'].Fill(np.sum(cluster_1Photon))
|
|
for i in range(clusterSize1Photon):
|
|
for j in range(clusterSize1Photon):
|
|
_hists['h2_1PhotonSumSample'].Fill(i, j, cluster_1Photon[i, j])
|
|
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
|
|
cluster_1Photon_list.append(cluster_1Photon)
|
|
refpoint_1Photon_list.append([int(ref_x)+Roi[0], int(ref_y)+Roi[2]])
|
|
|
|
sumEne = np.sum(cluster_1Photon)
|
|
### calculate eta
|
|
projectedX = np.sum(cluster_1Photon, axis=0)
|
|
etaX = (projectedX[2] - projectedX[0]) / sumEne
|
|
projectedY = np.sum(cluster_1Photon, axis=1)
|
|
etaY = (projectedY[2] - projectedY[0]) / sumEne
|
|
_hists['h1_1PhotonEta_X'].Fill(etaX)
|
|
_hists['h1_1PhotonEta_Y'].Fill(etaY)
|
|
|
|
### double photon events
|
|
if 2 * Energy - 2 * selectionRange < energy < 2 * Energy + 2 * selectionRange:
|
|
x_center = np.sum(xs * enes) / np.sum(enes)
|
|
y_center = np.sum(ys * enes) / np.sum(enes)
|
|
ref_x = int(x_center - clusterSize2Photon / 2) + 1 ### refered to the lower-left corner of the cluster
|
|
ref_y = int(y_center - clusterSize2Photon / 2) + 1
|
|
|
|
if int(ref_x) < 0 or int(ref_x)+clusterSize2Photon > signalEneFrame.shape[1] or int(ref_y) < 0 or int(ref_y)+clusterSize2Photon > signalEneFrame.shape[0]:
|
|
continue ### skip clusters too close to the border
|
|
for i in range(len(xs)):
|
|
_hists['h2_2PhotonSumFrame'].Fill(xs[i]+Roi[0], ys[i]+Roi[2], enes[i])
|
|
_hists['h1_2PhotonClusterPixelEnergy'].Fill(enes[i])
|
|
_hists['h1_2PhotonClusterSize'].Fill(clusterSize)
|
|
diameterX = max(xs) - min(xs) + 1
|
|
diameterY = max(ys) - min(ys) + 1
|
|
_hists['h1_2PhotonClusterDiameterX'].Fill(diameterX)
|
|
_hists['h1_2PhotonClusterDiameterY'].Fill(diameterY)
|
|
|
|
cluster_2Photon = np.zeros((clusterSize2Photon, clusterSize2Photon), dtype=np.float32)
|
|
for i in range(len(xs)):
|
|
x_rel = xs[i] - int(ref_x)
|
|
y_rel = ys[i] - int(ref_y)
|
|
if 0 <= x_rel < clusterSize2Photon and 0 <= y_rel < clusterSize2Photon:
|
|
cluster_2Photon[y_rel, x_rel] = enes[i]
|
|
_hists['h2_2PhotonSumSample'].Fill(x_rel, y_rel, enes[i])
|
|
# cluster_2Photon = signalEneFrame[int(ref_y):int(ref_y)+clusterSize2Photon, int(ref_x):int(ref_x)+clusterSize2Photon]
|
|
_hists['h1_2PhotonClusterEnergy'].Fill(np.sum(cluster_2Photon))
|
|
for i in range(clusterSize2Photon):
|
|
for j in range(clusterSize2Photon):
|
|
_hists['h2_2PhotonSumSample'].Fill(i, j, cluster_2Photon[i, j])
|
|
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
|
|
cluster_2Photon_list.append(cluster_2Photon)
|
|
refpoint_2Photon_list.append([int(ref_x)+Roi[0], int(ref_y)+Roi[2]])
|
|
|
|
_hists['h2_ClusterSizeVsEnergy'].Fill(energy, clusterSize)
|
|
_diameterX = max(xs) - min(xs) + 1
|
|
_diameterY = max(ys) - min(ys) + 1
|
|
_diameter = max(_diameterX, _diameterY)
|
|
_hists['h2_DiameterVsEnergy'].Fill(energy, _diameter)
|
|
for i in range(len(xs)):
|
|
_hists['h1_PixelEnergy'].Fill(enes[i])
|
|
_hists['h2_HitsSumFrame'].Fill(xs[i]+Roi[0], ys[i]+Roi[2], enes[i])
|
|
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
|
|
if len(cluster_1Photon_list) > BUFFER_THRESHOLD or len(cluster_2Photon_list) > BUFFER_THRESHOLD:
|
|
flush_h5_buffer()
|
|
|
|
for _x in range(Roi[0], Roi[1]):
|
|
for _y in range(Roi[2], Roi[3]):
|
|
_hists['h2_SumFrame'].Fill(_x, _y, sumFrame[_y, _x])
|
|
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
|
|
flush_h5_buffer()
|
|
h5_1Photon_file.close()
|
|
h5_2Photon_file.close()
|
|
return _hists
|
|
|
|
def process():
|
|
import os
|
|
os.makedirs(_cfg['outputFolder'], exist_ok=True)
|
|
with Pool(_cfg['NThread'], maxtasksperchild=1) as p:
|
|
results = p.map(_processFrames, range(nChunks))
|
|
|
|
global hists
|
|
hists = bookHistograms(_cfg['energy'])
|
|
for h in hists.values():
|
|
h.SetDirectory(0)
|
|
|
|
for _hists in results:
|
|
for key in hists.keys():
|
|
hists[key].Add(_hists[key])
|
|
del _hists
|
|
|
|
outputTFileName = f'/home/xie_x1/MLXID/DataProcess/Results/{_cfg["runName"]}.root'
|
|
outputTFile = TFile(outputTFileName, 'RECREATE')
|
|
for hist in hists.values():
|
|
hist.Write()
|
|
outputTFile.Close()
|
|
|
|
def getPedestalAndNoise():
|
|
NX, NY = _cfg['NX'], _cfg['NY']
|
|
pedestalFileName = _cfg['pedestalFileName']
|
|
if 'nFiber' in _cfg and _cfg['nFiber'] == 2:
|
|
pedestalFileName2 = pedestalFileName.replace('d0', 'd1')
|
|
_pedestalAduFrames1 = np.fromfile(f'{pedestalFileName}', dtype=np.uint16).astype(np.int32)
|
|
_pedestalAduFrames1 = _pedestalAduFrames1.reshape(-1, NX * NY //2 + 56)
|
|
_pedestalAduFrames1 = _pedestalAduFrames1[:, 56:].reshape(-1, NY//2, NX)
|
|
_pedestalAduFrames2 = np.fromfile(f'{pedestalFileName2}', dtype=np.uint16).astype(np.int32)
|
|
_pedestalAduFrames2 = _pedestalAduFrames2.reshape(-1, NX * NY //2 + 56)
|
|
_pedestalAduFrames2 = _pedestalAduFrames2[:, 56:].reshape(-1, NY//2, NX)
|
|
pedestalAduFrames = np.concatenate([_pedestalAduFrames1, _pedestalAduFrames2], axis=1)
|
|
del _pedestalAduFrames1, _pedestalAduFrames2
|
|
else:
|
|
pedestalAduFrames = np.fromfile(f'{pedestalFileName}', dtype=np.uint16).astype(np.int32)
|
|
pedestalAduFrames = pedestalAduFrames.reshape(-1, NX * NY + 56)
|
|
pedestalAduFrames = pedestalAduFrames[:, 56:].reshape(-1, NX, NY)
|
|
|
|
pedestalAduFrames = pedestalAduFrames[pedestalAduFrames.shape[0]//10:] # skip the first 10% frames
|
|
global _pedestalAduFrame, _noiseEneFrame
|
|
_pedestalAduFrame = np.mean(pedestalAduFrames, axis=0)
|
|
noiseAduFrame = np.std(pedestalAduFrames, axis=0, ddof=1)
|
|
with Pool(16) as p:
|
|
pedestalEneFrames = p.map(convertAduFrameToEnergyFrame, pedestalAduFrames)
|
|
_noiseEneFrame = np.std(pedestalEneFrames, axis=0, ddof=1)
|
|
print(f'Average noise = {np.mean(_noiseEneFrame):.3f} keV; {np.mean(noiseAduFrame):.3f} ADU')
|
|
del pedestalAduFrames, pedestalEneFrames
|
|
|
|
def getPedestalAndNoise_simplified():
|
|
NX, NY = _cfg['NX'], _cfg['NY']
|
|
pedestalFileName = _cfg['pedestalFileName']
|
|
pedestalAduFrames = np.fromfile(f'{pedestalFileName}', dtype=np.uint16).astype(np.int32)
|
|
pedestalAduFrames = pedestalAduFrames.reshape(-1, NX * NY + 56)
|
|
pedestalAduFrames = pedestalAduFrames[:, 56:].reshape(-1, NX, NY)
|
|
pedestalAduFrames = pedestalAduFrames[pedestalAduFrames.shape[0]//10:] # skip the first 10% frames
|
|
global _pedestalAduFrame, _noiseEneFrame
|
|
_pedestalAduFrame = np.mean(pedestalAduFrames, axis=0)
|
|
noiseAduFrame = np.std(pedestalAduFrames, axis=0, ddof=1)
|
|
_noiseEneFrame = convertAduFrameToEnergyFrame(noiseAduFrame)
|
|
|
|
print(f'Average noise = {np.mean(_noiseEneFrame):.3f} keV; {np.mean(noiseAduFrame):.3f} ADU')
|
|
del pedestalAduFrames
|