diff --git a/Configs.py b/Configs.py new file mode 100644 index 0000000..f360c4e --- /dev/null +++ b/Configs.py @@ -0,0 +1,23 @@ +from glob import glob +import numpy as np + +Configs = {} + + +XrayFolder = f'/mnt/sls_det_storage/moench_data/2504_SoleilBeamtime/Moench040_g4_SimonStars_150V_50us' +Configs['SiemenStarUpperLeft'] = { + 'caliFile': np.load('/home/xie_x1/MLED/data-process/utils/BacksidePulsing_Calibration/Moench040/Moench040_AduToKevMapping_g4_50us_150V_250423.npy'), + 'pedestalFileName': f'{XrayFolder}/15keV_pedestal_d0_f0_0.raw', + 'signalFileNames': glob(f'{XrayFolder}/15keV_signal_d0_f*_0.raw'), + # 'etaxLUT': '/mnt/sls_det_storage/moench_data/MLXID/EtaInterpolation/Moench040_EtaxLUT_g4_50us_150V_250423.npy', + # 'etayLUT': '/mnt/sls_det_storage/moench_data/MLXID/EtaInterpolation/Moench040_EtayLUT_g4_50us_150V_250423.npy', + 'Roi': [140, 230, 120, 210], + 'NX': 400, 'NY': 400, + 'energy': 15, # keV + 'selectionRange': 2, + 'headerSize': 56, + # 'write3x3Cluster': True, + 'NFrame': 20_000_000, ### 20_000_000 in total + 'NThread': 16, + 'NChunks': 200, +} \ No newline at end of file diff --git a/DataAnalysis.py b/DataAnalysis.py index 97d7865..fa3c622 100644 --- a/DataAnalysis.py +++ b/DataAnalysis.py @@ -1,42 +1,8 @@ import UsefulFuncs -import imp -from glob import glob -from ROOT import * import numpy as np -from array import array -import matplotlib.pyplot as plt +from Configs import Configs -imp.reload(UsefulFuncs) -exptime = 50 -gain = 'g4' -moenchCode = 'Moench040' -date = '250423' -hv = 150 -Folder = f'/mnt/sls_det_storage/moench_data/xiangyu/{date}_BP_cali_{moenchCode}_{hv}V_{gain}_{exptime}us' -RootFileName = f'{moenchCode}_AduToKevMapping_{gain}_{exptime}us_{hv}V_{date}.root' -CaliLut3DFileName = f'/home/xie_x1/MLED/data-process/utils/BacksidePulsing_Calibration/Moench040/{moenchCode}_AduToKevMapping_{gain}_{exptime}us_{hv}V_{date}.npy' - -XrayFolder = f'/mnt/sls_det_storage/moench_data/2504_SoleilBeamtime/{moenchCode}_g4_SimonStars_150V_50us' -Roi = [140, 230, 120, 210] -Energy = 15 - -measurementConfig = { - 'caliFile': np.load(CaliLut3DFileName), - 'pedestalFileName': f'{XrayFolder}/{Energy}keV_pedestal_d0_f0_0.raw', - 'signalFileNames': glob(f'{XrayFolder}/{Energy}keV_signal_d0_f*_0.raw'), - 'etaxLUT': '/mnt/sls_det_storage/moench_data/MLXID/EtaInterpolation/Moench040_EtaxLUT_g4_50us_150V_250423.npy', - 'etayLUT': '/mnt/sls_det_storage/moench_data/MLXID/EtaInterpolation/Moench040_EtayLUT_g4_50us_150V_250423.npy', - 'Roi': Roi, - 'NX': 400, 'NY': 400, - 'energy': Energy, # keV - 'selectionRange': 2, - 'NFrame': 20_000_000, ### 20_000_000 in total - 'NThread': 16, - 'headerSize': 56, - 'clusterWidth': 3, ### cluster width, must be odd number - 'write3x3Cluster': True, - 'NChunks': 200, -} +measurementConfig = Configs['SiemenStarUpperLeft'] UsefulFuncs.init(measurementConfig) UsefulFuncs.getPedestalAndNoise_simplified() diff --git a/UsefulFuncs.py b/UsefulFuncs.py index a1cd2c5..1be273a 100644 --- a/UsefulFuncs.py +++ b/UsefulFuncs.py @@ -12,9 +12,6 @@ _cfg = {} ### global config _aduToKev3DMap = None ### global caliFile _pedestalAduFrame = None ### global pedestalAduFrame _noiseEneFrame = None ### global noiseEneFrame -h1_ClusterEnergy, h1_PixelEnergy, h1_CenterPixelEnergy, h2_Qc = None, None, None, None -sumFrame = np.zeros((400, 400), dtype=np.float32) ### global sumFrame for all frames -EtaxLUT, EtayLUT = None, None ### global LUT for charge weighted center nChunks = 16 def init(cfg): @@ -23,12 +20,6 @@ def init(cfg): _cfg = cfg _aduToKev3DMap = cfg['caliFile'] Roi = cfg['Roi'] - global sumFrame - sumFrame = sumFrame[Roi[2]:Roi[3], Roi[0]:Roi[1]] ### crop the sumFrame to the ROI - if 'etaxLUT' in cfg and 'etayLUT' in cfg: - global EtaxLUT, EtayLUT - EtaxLUT = np.load(cfg['etaxLUT']) - EtayLUT = np.load(cfg['etayLUT']) global nChunks nChunks = _cfg.get('NChunks', 16) @@ -58,37 +49,55 @@ def convertAduFrameToEnergyFrame(aduFrame): energyFrame *= _cfg['energyScalingFactor'] return energyFrame -def bookHistograms(energy, suffix = '', energyBinWidth = 0.1, isMC = False): +def bookHistograms(energy, suffix='', energyBinWidth=0.1, isMC=False): histMaxEnergy = energy * 5.5 - h1_ClusterEnergy = TH1D(f'h1_ClusterEnergy{suffix}', f'h1_ClusterEnergy{suffix}', int(histMaxEnergy//energyBinWidth), -1, histMaxEnergy) - h1_ClusterEnergy.SetDirectory(0) - h1_ClusterEnergy.SetTitle(f'Cluster Energy;Energy [keV];Counts') - h1_PixelEnergy = TH1D(f'h1_PixelEnergy{suffix}', f'h1_PixelEnergy{suffix}', int(histMaxEnergy//energyBinWidth), -1, histMaxEnergy) - h1_PixelEnergy.SetDirectory(0) - h1_PixelEnergy.SetTitle(f'Pixel Energy;Energy [keV];Counts') - h1_CenterPixelEnergy = TH1D(f'h1_CenterPixelEnergy{suffix}', f'h1_CenterPixelEnergy{suffix}', int(histMaxEnergy//energyBinWidth), -1, histMaxEnergy) - h1_CenterPixelEnergy.SetDirectory(0) - h1_CenterPixelEnergy.SetTitle(f'Center Pixel Energy;Energy [keV];Counts') - h2_Qc = TH2D(f'h2_Qc{suffix}', f'h2_Qc{suffix}', 201, 0, 1, 201, 0, 1) - h2_Qc.SetDirectory(0) - h2_Qc.SetTitle(f'Charge weighted center;#eta_x;#eta_y;Counts') - h2_counts = TH2D(f'h2_counts{suffix}', f'h2_counts{suffix}', _cfg['Roi'][1] - _cfg['Roi'][0], _cfg['Roi'][0], _cfg['Roi'][1], _cfg['Roi'][3] - _cfg['Roi'][2], _cfg['Roi'][2], _cfg['Roi'][3]) - h2_counts.SetDirectory(0) - h2_HitsSumFrame = TH2D(f'h2_HitsSumFrame{suffix}', f'h2_HitsSumFrame{suffix}', _cfg['Roi'][1] - _cfg['Roi'][0], _cfg['Roi'][0], _cfg['Roi'][1], _cfg['Roi'][3] - _cfg['Roi'][2], _cfg['Roi'][2], _cfg['Roi'][3]) - h2_HitsSumFrame.SetDirectory(0) + 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 - xBins = (_cfg['Roi'][1] - _cfg['Roi'][0]) * interpolationBins - yBins = (_cfg['Roi'][3] - _cfg['Roi'][2]) * interpolationBins - h2_interpolated = TH2D(f'h2_interpolated{suffix}', f'h2_interpolated{suffix}', xBins, _cfg['Roi'][0], _cfg['Roi'][1], yBins, _cfg['Roi'][2], _cfg['Roi'][3]) - h2_interpolated.SetDirectory(0) - h2_interpolatedSubpixel = TH2D(f'h2_interpolatedSubpixel{suffix}', f'h2_interpolatedSubpixel{suffix}', interpolationBins, 0, 1, interpolationBins, 0, 1) - h2_interpolatedSubpixel.SetDirectory(0) + + histograms = {} + histograms['h1_ClusterEnergy'] = TH1D( + f'h1_ClusterEnergy{suffix}', f'h1_ClusterEnergy{suffix}', nEnergyBins, -1, histMaxEnergy + ) + histograms['h1_PixelEnergy'] = TH1D( + f'h1_PixelEnergy{suffix}', f'h1_PixelEnergy{suffix}', nEnergyBins, -1, histMaxEnergy + ) + histograms['h1_CenterPixelEnergy'] = TH1D( + f'h1_CenterPixelEnergy{suffix}', f'h1_CenterPixelEnergy{suffix}', nEnergyBins, -1, histMaxEnergy + ) + histograms['h2_Qc'] = TH2D( + f'h2_Qc{suffix}', f'h2_Qc{suffix}', 201, 0, 1, 201, 0, 1 + ) + histograms['h2_counts'] = TH2D( + f'h2_counts{suffix}', f'h2_counts{suffix}', roi_width, roi_x0, roi_x1, roi_height, roi_y0, roi_y1 + ) + 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_interpolated'] = TH2D( + f'h2_interpolated{suffix}', f'h2_interpolated{suffix}', roi_width * interpolationBins, roi_x0, roi_x1, roi_height * interpolationBins, roi_y0, roi_y1 + ) + histograms['h2_interpolatedSubpixel'] = TH2D( + f'h2_interpolatedSubpixel{suffix}', f'h2_interpolatedSubpixel{suffix}', interpolationBins, 0, 1, interpolationBins, 0, 1 + ) + + for hist in histograms.values(): + hist.SetDirectory(0) + + histograms['h1_ClusterEnergy'].SetTitle('Cluster Energy;Energy [keV];Counts') + histograms['h1_PixelEnergy'].SetTitle('Pixel Energy;Energy [keV];Counts') + histograms['h1_CenterPixelEnergy'].SetTitle('Center Pixel Energy;Energy [keV];Counts') + histograms['h2_Qc'].SetTitle('Charge weighted center;#eta_x;#eta_y;Counts') if isMC: - h1_ChargeCollectionEfficiency = TH1D(f'h1_ChargeCollectionEfficiency{suffix}', f'h1_ChargeCollectionEfficiency{suffix}', 400, 0.8, 1.01) - return h1_ClusterEnergy.Clone(), h1_PixelEnergy.Clone(), h1_CenterPixelEnergy.Clone(), h2_Qc.Clone(), h1_ChargeCollectionEfficiency.Clone() + histograms['h1_ChargeCollectionEfficiency'] = TH1D(f'h1_ChargeCollectionEfficiency{suffix}', f'h1_ChargeCollectionEfficiency{suffix}', 400, 0.8, 1.01) + histograms['h1_ChargeCollectionEfficiency'].SetDirectory(0) + for key, hist in list(histograms.items()): + histograms[key] = hist.Clone() - return h1_ClusterEnergy, h1_PixelEnergy, h1_CenterPixelEnergy, h2_Qc, h2_counts, h2_HitsSumFrame, h2_interpolated, h2_interpolatedSubpixel + 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']) @@ -101,7 +110,7 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin headerSize = _cfg['headerSize'] Roi = _cfg['Roi'] nFramePerFile = os.path.getsize(f'{signalFileNames[0]}') // (NX * NY + 56) // 2 - _h1_PixelEnergy, _h1_ClusterEnergy, _h1_CenterPixelEnergy, _h2_Qc, _h2_counts, _h2_HitsSumFrame, _h2_interpolated, _h2_interpolatedSubpixel = bookHistograms(Energy, f'_{idxChunk}') + _hists = bookHistograms(Energy, f'_{idxChunk}') ### create a cluster finder nxROI = Roi[1] - Roi[0] @@ -155,7 +164,7 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin 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} from file {signalFileNames[idxFile]}: {e}, stop processing this chunk ({startFrame}-{endFrame})') - return _h1_ClusterEnergy.Clone(), _h1_PixelEnergy.Clone(), _h1_CenterPixelEnergy.Clone(), _h2_Qc.Clone(), _h2_counts.Clone(), _h2_HitsSumFrame.Clone(), _h2_interpolated.Clone(), _h2_interpolatedSubpixel.Clone() + return _hists signalAduFrame = signalAduFrame - _pedestalAduFrame signalEneFrame = convertAduFrameToEnergyFrame(signalAduFrame) @@ -164,41 +173,40 @@ def _processFrames(idxChunk): ### for both single and double photon events, usin clusters = CF.hits() for idxCluster in range(len(clusters)): energy = clusters['energy'][idxCluster] - _h1_ClusterEnergy.Fill(energy) + _hists['h1_ClusterEnergy'].Fill(energy) xs = clusters['cols'][idxCluster] ys = clusters['rows'][idxCluster] enes = clusters['enes'][idxCluster] if energy < 5: continue for i in range(len(xs)): - _h1_PixelEnergy.Fill(enes[i]) - _h2_HitsSumFrame.Fill(xs[i]+Roi[0], ys[i]+Roi[2], enes[i]) + _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: flush_h5_buffer() h5_file.close() - return _h1_ClusterEnergy, _h1_PixelEnergy, _h1_CenterPixelEnergy, _h2_Qc, _h2_counts, _h2_HitsSumFrame, _h2_interpolated, _h2_interpolatedSubpixel + return _hists def process(): with Pool(_cfg['NThread'], maxtasksperchild=1) as p: results = p.map(_processFrames, range(nChunks)) - global h1_ClusterEnergy, h1_PixelEnergy, h1_CenterPixelEnergy, h2_Qc, h2_counts, h2_HitsSumFrame, h2_interpolated, h2_interpolatedSubpixel - h1_ClusterEnergy, h1_PixelEnergy, h1_CenterPixelEnergy, h2_Qc, h2_counts, h2_HitsSumFrame, h2_interpolated, h2_interpolatedSubpixel = bookHistograms(_cfg['energy']) - for h in (h1_ClusterEnergy, h1_PixelEnergy, h1_CenterPixelEnergy, h2_Qc, h2_counts, h2_HitsSumFrame, h2_interpolated, h2_interpolatedSubpixel): + global hists + hists = bookHistograms(_cfg['energy']) + for h in hists.values(): h.SetDirectory(0) - global sumFrame - for _h1_ClusterEnergy, _h1_PixelEnergy, _h1_CenterPixelEnergy, _h2_Qc, _h2_counts, _h2_HitsSumFrame, _h2_interpolated, _h2_interpolatedSubpixel in results: - h1_ClusterEnergy.Add(_h1_ClusterEnergy) - h1_PixelEnergy.Add(_h1_PixelEnergy) - h1_CenterPixelEnergy.Add(_h1_CenterPixelEnergy) - h2_Qc.Add(_h2_Qc) - h2_interpolated.Add(_h2_interpolated) - h2_interpolatedSubpixel.Add(_h2_interpolatedSubpixel) - h2_counts.Add(_h2_counts) - h2_HitsSumFrame.Add(_h2_HitsSumFrame) - del _h1_ClusterEnergy, _h1_PixelEnergy, _h1_CenterPixelEnergy, _h2_Qc, _h2_counts, _h2_HitsSumFrame, _h2_interpolated, _h2_interpolatedSubpixel + for _hists in results: + hists['h1_ClusterEnergy'].Add(_hists['h1_ClusterEnergy']) + hists['h1_PixelEnergy'].Add(_hists['h1_PixelEnergy']) + hists['h1_CenterPixelEnergy'].Add(_hists['h1_CenterPixelEnergy']) + hists['h2_Qc'].Add(_hists['h2_Qc']) + hists['h2_interpolated'].Add(_hists['h2_interpolated']) + hists['h2_interpolatedSubpixel'].Add(_hists['h2_interpolatedSubpixel']) + hists['h2_counts'].Add(_hists['h2_counts']) + hists['h2_HitsSumFrame'].Add(_hists['h2_HitsSumFrame']) + del _hists def getPedestalAndNoise(): NX, NY = _cfg['NX'], _cfg['NY']