From f39127f7ea42465ae77dfd2dbdb4109e74915692 Mon Sep 17 00:00:00 2001 From: "xiangyu.xie" Date: Fri, 28 Nov 2025 08:31:44 +0100 Subject: [PATCH] Add 1photon sample outpu; add more hists --- UsefulFuncs.py | 143 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 124 insertions(+), 19 deletions(-) diff --git a/UsefulFuncs.py b/UsefulFuncs.py index 2bbdaf9..82b7f85 100644 --- a/UsefulFuncs.py +++ b/UsefulFuncs.py @@ -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))