MLXID/DataProcess
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add 1photon sample outpu; add more hists

Browse Source
This commit is contained in:
xie_x1
2025-11-28 08:31:44 +01:00
parent 09eae2231c
commit f39127f7ea
1 changed files with 124 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

143
UsefulFuncs.py
View File

@@ -69,15 +69,80 @@ def bookHistograms(energy, suffix='', energyBinWidth=0.1, isMC=False):
)
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_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_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_DoubleHitsSumFrame'] = TH2D(
f'h2_DoubleHitsSumFrame{suffix}', f'h2_DoubleHitsSumFrame{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1
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_DoubleHitsSumFrame'].SetTitle('Double Hits Sum Frame;X [pixel];Y [pixel];Energy Sum [keV]')
histograms['h2_1PhotonSumFrame'].SetTitle('1-Photon Hits Sum Frame;X [pixel];Y [pixel];Energy Sum [keV]')
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
@@ -115,7 +180,7 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin
CF.set_noiseMap(_noiseEneFrame)
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
h5_1Photon_file = h5py.File(f'{_cfg["outputFolder"]}/Clusters_1Photon_chunk{idxChunk}.h5', 'w')
h5_1Photon_file = h5py.File(f'{_cfg["outputFolder"]}/1Photon_CS3_chunk{idxChunk}.h5', 'w')
dset_1Photon_clusters = h5_1Photon_file.create_dataset(
'clusters', (0, 3, 3), maxshape=(None, 3, 3), dtype='f4',
chunks=True, compression='gzip'
@@ -128,7 +193,7 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin
cluster_1Photon_list = []
refpoint_1Photon_list = []
h5_2Photon_file = h5py.File(f'{_cfg["outputFolder"]}/Clusters_2Photon_CS{clusterSize2Photon}_chunk{idxChunk}.h5', 'w')
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'
@@ -160,8 +225,6 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin
refpoint_2Photon_list.clear()
for idxFrame in range(startFrame, endFrame):
if idxFrame % 10000 == 0:
print(f'Processing frame {idxFrame}...')
idxFile, idxFrame = divmod(idxFrame, nFramePerFile)
try:
if 'nFiber' in _cfg and _cfg['nFiber'] == 2:
@@ -192,22 +255,62 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin
ys = clusters['rows'][idxCluster][:clusterSize]
enes = clusters['enes'][idxCluster][:clusterSize]
if energy < 5:
continue
### single photon events
if Energy - selectionRange < energy < Energy + selectionRange:
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)
_hists['h1_1PhotonClusterEnergy'].Fill(energy)
### detect single photon events
maxX, maxY = clusters['col'][idxCluster], clusters['row'][idxCluster]
ref_x = maxX - 1 ### refered to the lower-left corner of the cluster
ref_y = maxY - 1
cluster_1Photon = np.zeros((3, 3), dtype=np.float32)
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 < 3 and 0 <= y_rel < 3:
cluster_1Photon[y_rel, x_rel] = enes[i]
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 - selectionRange < energy < 2 * Energy + selectionRange:
for i in range(len(xs)):
_hists['h2_DoubleHitsSumFrame'].Fill(xs[i]+Roi[0], ys[i]+Roi[2], enes[i])
_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)
ref_x = (min(xs) + max(xs)) // 2 - (clusterSize2Photon // 2) ### refered to the lower-left corner of the cluster
ref_y = (min(ys) + max(ys)) // 2 - (clusterSize2Photon // 2)
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])
if 'writeClusters' in _cfg and _cfg['writeClusters'] == True:
ref_x = (min(xs) + max(xs)) // 2 - (clusterSize2Photon // 2) ### refered to the lower-left corner of the cluster
ref_y = (min(ys) + max(ys)) // 2 - (clusterSize2Photon // 2)
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]
cluster_2Photon_list.append(cluster_2Photon)
refpoint_2Photon_list.append([int(ref_x)+Roi[0], int(ref_y)+Roi[2]])
@@ -227,6 +330,8 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin
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))