From d9af7044efff52792d9e4c82c35f26b278b33391 Mon Sep 17 00:00:00 2001 From: "xiangyu.xie" Date: Fri, 29 May 2026 16:06:11 +0200 Subject: [PATCH] Add hists and logics for 3-photon clusters --- UsefulFuncs.py | 96 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 95 insertions(+), 1 deletion(-) diff --git a/UsefulFuncs.py b/UsefulFuncs.py index f6e66bd..851df99 100644 --- a/UsefulFuncs.py +++ b/UsefulFuncs.py @@ -12,6 +12,7 @@ _aduToKev3DMap = None ### global caliFile _pedestalAduFrame = None ### global pedestalAduFrame _noiseEneFrame = None ### global noiseEneFrame nChunks = 16 +clusterSize3Photon = 9 clusterSize2Photon = 7 clusterSize1Photon = 3 @@ -111,26 +112,51 @@ def bookHistograms(energy, suffix='', energyBinWidth=0.1, isMC=False): ) histograms['h1_2PhotonClusterEnergy'].SetTitle('2-Photon Cluster Energy;Energy [keV];Counts') + histograms['h1_3PhotonClusterEnergy'] = TH1D( + f'h1_3PhotonClusterEnergy{suffix}', f'h1_3PhotonClusterEnergy{suffix}', nEnergyBins, -1, energy*6 + ) + histograms['h1_3PhotonClusterEnergy'].SetTitle('3-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_3PhotonClusterSize'] = TH1D( + f'h1_3PhotonClusterSize{suffix}', f'h1_3PhotonClusterSize{suffix}', 30, 0, 30 + ) + histograms['h1_3PhotonClusterSize'].SetTitle('3-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_3PhotonClusterDiameterX'] = TH1D( + f'h1_3PhotonClusterDiameterX{suffix}', f'h1_3PhotonClusterDiameterX{suffix}', 15, 0, 15 + ) + histograms['h1_3PhotonClusterDiameterX'].SetTitle('3-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_3PhotonClusterDiameterY'] = TH1D( + f'h1_3PhotonClusterDiameterY{suffix}', f'h1_3PhotonClusterDiameterY{suffix}', 15, 0, 15 + ) + histograms['h1_3PhotonClusterDiameterY'].SetTitle('3-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['h1_3PhotonClusterPixelEnergy'] = TH1D( + f'h1_3PhotonClusterPixelEnergy{suffix}', f'h1_3PhotonClusterPixelEnergy{suffix}', nEnergyBins, -1, energy*2.5 + ) + histograms['h1_3PhotonClusterPixelEnergy'].SetTitle('3-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 ) @@ -156,18 +182,28 @@ def bookHistograms(energy, suffix='', energyBinWidth=0.1, isMC=False): ) histograms['h2_2PhotonSumSample'].SetTitle('2-Photon Sum Sample;Sum Sample;Counts') + histograms['h2_3PhotonSumSample'] = TH2D( + f'h2_3PhotonSumSample{suffix}', f'h2_3PhotonSumSample{suffix}', clusterSize3Photon, 0, clusterSize3Photon, clusterSize3Photon, 0, clusterSize3Photon + ) + histograms['h2_3PhotonSumSample'].SetTitle('3-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_3PhotonSumFrame'] = TH2D( + f'h2_3PhotonSumFrame{suffix}', f'h2_3PhotonSumFrame{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1 + ) + histograms['h2_3PhotonSumFrame'].SetTitle('3-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 + f'h2_DiameterVsEnergy{suffix}', f'h2_DiameterVsEnergy{suffix}', nEnergyBins, -1, histMaxEnergy, 15, 0, 15 ) histograms['h2_DiameterVsEnergy'].SetTitle('Cluster Diameter vs Energy;Energy [keV];Cluster Diameter [pixel]') @@ -224,6 +260,19 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin ) cluster_2Photon_list = [] refpoint_2Photon_list = [] + + h5_3Photon_file = h5py.File(f'{_cfg["outputFolder"]}/3Photon_CS{clusterSize3Photon}_chunk{idxChunk}.h5', 'w') + dset_3Photon_clusters = h5_3Photon_file.create_dataset( + 'clusters', (0, clusterSize3Photon, clusterSize3Photon), maxshape=(None, clusterSize3Photon, clusterSize3Photon), dtype='f4', + chunks=True, compression='gzip' + ) + dset_3Photon_refs = h5_3Photon_file.create_dataset( + 'referencePoint', (0, 2), maxshape=(None, 2), dtype='i4', + chunks=True + ) + cluster_3Photon_list = [] + refpoint_3Photon_list = [] + BUFFER_THRESHOLD = 10000 def flush_h5_buffer(): if not cluster_1Photon_list: return # no new clusters to write @@ -236,6 +285,7 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin cluster_1Photon_list.clear() refpoint_1Photon_list.clear() + if not cluster_2Photon_list: return # no new clusters to write old_len = dset_2Photon_clusters.shape[0] new_len = old_len + len(cluster_2Photon_list) dset_2Photon_clusters.resize((new_len, clusterSize2Photon, clusterSize2Photon)) @@ -245,6 +295,16 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin cluster_2Photon_list.clear() refpoint_2Photon_list.clear() + if not cluster_3Photon_list: return # no new clusters to write + old_len = dset_3Photon_clusters.shape[0] + new_len = old_len + len(cluster_3Photon_list) + dset_3Photon_clusters.resize((new_len, clusterSize3Photon, clusterSize3Photon)) + dset_3Photon_refs.resize((new_len, 2)) + dset_3Photon_clusters[old_len:new_len] = cluster_3Photon_list + dset_3Photon_refs[old_len:new_len] = refpoint_3Photon_list + cluster_3Photon_list.clear() + refpoint_3Photon_list.clear() + for idxFrame in range(startFrame, endFrame): idxFile, idxFrame = divmod(idxFrame, nFramePerFile) try: @@ -357,6 +417,40 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin cluster_2Photon_list.append(cluster_2Photon) refpoint_2Photon_list.append([int(ref_x)+Roi[0], int(ref_y)+Roi[2]]) + if 3 * Energy - 3 * selectionRange < energy < 3 * Energy + 3 * selectionRange: + x_center = np.sum(xs * enes) / np.sum(enes) + y_center = np.sum(ys * enes) / np.sum(enes) + ref_x = int(x_center - clusterSize3Photon / 2) + 1 ### refered to the lower-left corner of the cluster + ref_y = int(y_center - clusterSize3Photon / 2) + 1 + + if int(ref_x) < 0 or int(ref_x)+clusterSize3Photon > signalEneFrame.shape[1] or int(ref_y) < 0 or int(ref_y)+clusterSize3Photon > signalEneFrame.shape[0]: + continue ### skip clusters too close to the border + for i in range(len(xs)): + _hists['h2_3PhotonSumFrame'].Fill(xs[i]+Roi[0], ys[i]+Roi[2], enes[i]) + _hists['h1_3PhotonClusterPixelEnergy'].Fill(enes[i]) + _hists['h1_3PhotonClusterSize'].Fill(clusterSize) + diameterX = max(xs) - min(xs) + 1 + diameterY = max(ys) - min(ys) + 1 + _hists['h1_3PhotonClusterDiameterX'].Fill(diameterX) + _hists['h1_3PhotonClusterDiameterY'].Fill(diameterY) + + ### filling the 3-photon cluster. Add background as well. + cluster_3Photon = np.zeros((clusterSize3Photon, clusterSize3Photon), 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 < clusterSize3Photon and 0 <= y_rel < clusterSize3Photon: + cluster_3Photon[y_rel, x_rel] = enes[i] + _hists['h2_3PhotonSumSample'].Fill(x_rel, y_rel, enes[i]) + cluster_3Photon += signalEneFrame[int(ref_y):int(ref_y)+clusterSize3Photon, int(ref_x):int(ref_x)+clusterSize3Photon] + _hists['h1_3PhotonClusterEnergy'].Fill(np.sum(cluster_3Photon)) + for i in range(clusterSize3Photon): + for j in range(clusterSize3Photon): + _hists['h2_3PhotonSumSample'].Fill(i, j, cluster_3Photon[i, j]) + if 'writeClusters' in _cfg and _cfg['writeClusters'] == True: + cluster_3Photon_list.append(cluster_3Photon) + refpoint_3Photon_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