diff --git a/scripts/jungfrau_create_pedestals.py b/scripts/jungfrau_create_pedestals.py new file mode 100644 index 0000000..1dcb9e4 --- /dev/null +++ b/scripts/jungfrau_create_pedestals.py @@ -0,0 +1,225 @@ +import argparse +import sys +import os +import numpy as np +import h5py +import logging + +ch = logging.StreamHandler() +ch.setFormatter(logging.Formatter('[%(levelname)s] %(message)s')) + +log = logging.getLogger("create_pedestals") +log.addHandler(ch) + + +def h5_printname(name): + print(" {}".format(name)) + + +def forcedGainValue(i, n0, n1, n2, n3): + if i <= n0 - 1: + return 0 + if i <= (n0 + n1) - 1: + return 1 + if i <= (n0 + n1 + n2) - 1: + return 3 + if i <= (n0 + n1 + n2 + n3) - 1: + return 4 + return 2 + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--verbosity", default=None, help="log verbosity level INFO/DEBUG/WARN/ERROR/CRITICAL") + parser.add_argument("--filename", default="pedestal.h5", help="pedestal file") + parser.add_argument("--X_test_pixel", type=int, default=0, help="x position of the test pixel") + parser.add_argument("--Y_test_pixel", type=int, default=0, help="y position of the test pixel") + parser.add_argument("--nFramesPede", type=int, default=1000, help="number of pedestal frames to average pedestal value") + parser.add_argument("--frames_G0", type=int, default=0, help="force to treat pedestal run as first frames_G0 taken in gain0, then frames_G1 in gain1, and frames_G2 in gain2 and HG0") + parser.add_argument("--frames_G1", type=int, default=0, help="force to treat pedestal run as first frames_G0 taken in gain0, then frames_G1 in gain1, and frames_G2 in gain2 and HG0") + parser.add_argument("--frames_G2", type=int, default=0, help="force to treat pedestal run as first frames_G0 taken in gain0, then frames_G1 in gain1, and frames_G2 in gain2 and HG0") + parser.add_argument("--frames_HG0", type=int, default=0, help="force to treat pedestal run as first frames_G0 taken in gain0, then frames_G1 in gain1, and frames_G2 in gain2 and HG0") + parser.add_argument("--number_frames", type=int, default=1000000, help="analyze only first number_frames frames") + parser.add_argument("--frames_average", type=int, default=1000, help="for pedestal in each gain average over last frames_average frames, reducing weight of previous") + parser.add_argument("--directory", default="./", help="Output directory where to store pixelmask and gain file") + parser.add_argument("--gain_check", type=int, default=1, help="check that gain setting in each of the module corresponds to the general gain switch, (0 - dont check)") + parser.add_argument("--add_pixel_mask", default=None, help="add additional masked pixels from external, specified file") + parser.add_argument("--number_bad_modules", type=int, default=0, help="Number of bad modules in detector") + args = parser.parse_args() + + if not (os.path.isfile(args.filename) and os.access(args.filename, os.R_OK)): + print("Pedestal file {} not found, exit".format(args.filename)) + exit() + + if args.verbosity: + log.setLevel(getattr(logging, args.verbosity.upper(), None)) + + overwriteGain = False + if (args.frames_G0 + args.frames_G1 + args.frames_G2) > 0: + log.info("Treat this run as taken with {} frames in gain0, then {} frames in gain1 and {} frames in gain2".format(args.frames_G0, args.frames_G1, args.frames_G2)) + overwriteGain = True + + f = h5py.File(args.filename, "r") + + #detector_name = (f.get("general/detector_name").value).decode('UTF-8') + detector_name = (f.get("general/detector_name")[()]).decode('UTF-8') + #n_bad_modules = f.get("general/n_bad_modules").value + n_bad_modules = args.number_bad_modules + + data_location = "data/" + detector_name + "/data" + daq_recs_location = "data/" + detector_name + "/daq_rec" + is_good_frame_location = "data/" + detector_name + "/is_good_frame" + + + numberOfFrames = len(f[data_location]) + (sh_y, sh_x) = f[data_location][0].shape + nModules = (sh_x * sh_y) // (1024 * 512) + if (nModules * 1024 * 512) != (sh_x * sh_y): + log.error(" {} : Something very strange in the data, Jungfrau consists of (1024x512) modules, while data has {}x{}".format(detector_name, sh_x, sh_y)) + exit() + + (tX, tY) = (args.X_test_pixel, args.Y_test_pixel) + if tX < 0 or tX > (sh_x - 1): + tX = 0 + if tY < 0 or tY > (sh_y - 1): + tY = 0 + + log.debug(" {} : test pixel is ( x y ): {}x{}".format(detector_name, tX, tY)) + log.info(" {} : In pedestal file {} there are {} frames".format(detector_name, args.filename, numberOfFrames + 1)) +# log.debug("Following groups are available:") +# if args.verbosity >= 3: +# f.visit(h5_printname) + log.debug(" {} : data has the following shape: {}, type: {}, {} modules ({} bad modules)".format(detector_name, f[data_location][0].shape, f[data_location][0].dtype, nModules, n_bad_modules)) + + pixelMask = np.zeros((sh_y, sh_x), dtype=np.int) + + adcValuesN = np.zeros((5, sh_y, sh_x)) + adcValuesNN = np.zeros((5, sh_y, sh_x)) + + + averagePedestalFrames = args.frames_average + + nMgain = [0] * 5 + + gainCheck = -1 + highG0Check = 0 + printFalseGain = False + nGoodFrames = 0 + nGoodFramesGain = 0 + + analyzeFrames = min(numberOfFrames, args.number_frames) + + for n in range(analyzeFrames): + + if not f[is_good_frame_location][n]: + continue + + nGoodFrames += 1 + + daq_rec = (f[daq_recs_location][n])[0] + + image = f[data_location][n][:] + frameData = (np.bitwise_and(image, 0b0011111111111111)) + gainData = np.bitwise_and(image, 0b1100000000000000) >> 14 + trueGain = forcedGainValue(n, args.framesG0, args.framesG1, args.framesG2, args.framesHG0) if overwriteGain else ( (daq_rec & 0b11000000000000) >> 12 ) + highG0 = (daq_rec & 0b1) + + gainGoodAllModules = True + if args.gain_check > 0: + daq_recs = f[daq_recs_location][n] + for i in range(len(daq_recs)): + if trueGain != ((daq_recs[i] & 0b11000000000000) >> 12) or highG0 != (daq_recs[i] & 0b1): + gainGoodAllModules = False + + if highG0 == 1 and trueGain != 0: + gainGoodAllModules = False + log.info(" {} : Jungfrau is in the high G0 mode ({}), but gain settings is strange: {}".format( detector_name, highG0, trueGain)) + + nFramesGain = np.sum(gainData==(trueGain)) + if nFramesGain < (nModules - 0.5 - n_bad_modules) * (1024 * 512): # make sure that most are the modules are in correct gain + gainGoodAllModules = False + log.debug(" {} : Too many bad pixels, skip the frame {}, true gain: {}(highG0: {}) ({}); gain0 : {}; gain1 : {}; gain2 : {}; undefined gain : {}".format( detector_name, n, trueGain, highG0, nFramesGain, np.sum(gainData==0), np.sum(gainData==1), np.sum(gainData==3), np.sum(gainData==2))) + + if not gainGoodAllModules: + log.debug(" {} : In Frame Number {} : mismatch in modules and general settings, Gain: {} vs {}; HighG0: {} vs {} (or too many bad pixels)".format( detector_name, n, trueGain, ((daq_recs & 0b11000000000000) >> 12), highG0, (daq_recs & 0b1))) + continue + nGoodFramesGain += 1 + + if gainData[tY][tX] != trueGain: + if not printFalseGain: + log.info(" {} : Gain wrong for channel ({}x{}) should be {}, but {}. Frame {}. {} {}".format( detector_name, tX, tY, trueGain, gainData[tY][tX], n, trueGain, daq_rec)) + printFalseGain = True + else: + if gainCheck != -1 and printFalseGain: + log.info(" {} : Gain was wrong for channel ({}x{}) in previous frames, but now correct : {}. Frame {}.".format( detector_name, tX, tY, gainData[tY, tX], n)) + printFalseGain = False + + if gainData[tY][tX] != gainCheck or highG0Check != highG0: + log.info(" {} : Gain changed for ({}x{}) channel {} -> {} (highG0 setting: {} -> {}), frame number {}, match: {}".format( detector_name, tX, tY, gainCheck, gainData[tY][tX], highG0Check, highG0, n, gainData[tY][tX] == trueGain)) + gainCheck = gainData[tY][tX] + highG0Check = highG0 + + if gainGoodAllModules: + + pixelMask[gainData != trueGain] |= (1 << (trueGain+4*highG0)) + + trueGain += 4 * highG0 + + + nMgain[trueGain] += 1 + + if nMgain[trueGain] > averagePedestalFrames: + adcValuesN[trueGain] -= adcValuesN[trueGain] / averagePedestalFrames + adcValuesNN[trueGain] -= adcValuesNN[trueGain] / averagePedestalFrames + + adcValuesN[trueGain] += frameData + adcValuesNN[trueGain] += np.float_power(frameData, 2) + + + log.info(" {} : {} frames analyzed, {} good frames, {} frames without settings mismatch. Gain frames distribution (0,1,2,3,HG0) : ({})".format( detector_name, analyzeFrames, nGoodFrames, nGoodFramesGain, nMgain)) + + if args.add_pixel_mask != None: + if (os.path.isfile(args.add_pixel_mask) and os.access(args.add_pixel_mask, os.R_OK)): + additional_pixel_mask_file = h5py.File(args.add_pixel_mask, "r") + additional_pixel_mask = np.array(additional_pixel_mask_file["pixel_mask"]) + log.info("Will add additional masked pixels from file %s , number %d " % (args.add_pixel_mask, np.sum(additional_pixel_mask == 1))) + if additional_pixel_mask.shape == pixelMask.shape: + pixelMask[additional_pixel_mask == 1] |= (1 << 5) + else: + log.error(" shape of additional pixel mask ({}) doesn't match current ({})".format( additional_pixel_mask.shape, pixelMask.shape)) + else: + log.error(" Specified addition file with pixel mask not found or not reachable {}".format( args.add_pixel_mask)) + + fileNameIn = os.path.splitext(os.path.basename(args.filename))[0] + full_fileNameOut = args.directory + "/" + fileNameIn + ".res.h5" + log.info(" {} : Output file with pedestal corrections in: {}".format( detector_name, full_fileNameOut)) + outFile = h5py.File(full_fileNameOut, "w") + + gains = [None] * 4 + gainsRMS = [None] * 4 + + for gain in range(5): + numberFramesAverage = max(1, min(averagePedestalFrames, nMgain[gain])) + mean = adcValuesN[gain] / float(numberFramesAverage) + mean2 = adcValuesNN[gain] / float(numberFramesAverage) + variance = mean2 - np.float_power(mean, 2) + stdDeviation = np.sqrt(variance) + log.debug(" {} : gain {} values results (pixel ({},{}) : {} {}".format( detector_name, gain, tY, tX, mean[tY][tX], stdDeviation[tY][tX])) + if gain != 2: + g = gain if gain < 3 else (gain-1) + gains[g] = mean + gainsRMS[g] = stdDeviation + + pixelMask[np.isclose(stdDeviation,0)] |= (1 << (6 + g)) + + dset = outFile.create_dataset('pixel_mask', data=pixelMask) + dset = outFile.create_dataset('gains', data=gains) + dset = outFile.create_dataset('gainsRMS', data=gainsRMS) + + outFile.close() + + log.info(" {} : Number of good pixels: {} from {} in total ({} bad pixels)".format( detector_name, np.sum(pixelMask == 0), sh_x * sh_y, (sh_x * sh_y - np.sum(pixelMask == 0)))) + + +if __name__ == "__main__": + main() diff --git a/scripts/retrieve_detector_data.sh b/scripts/retrieve_detector_data.sh index 1583a40..f1c7822 100755 --- a/scripts/retrieve_detector_data.sh +++ b/scripts/retrieve_detector_data.sh @@ -55,6 +55,10 @@ case ${DETECTOR} in NM=9 DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF02.json ;; +'JF04T01V01') + NM=1 + DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF04.json + ;; 'JF06T32V02') NM=32 DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF06.json @@ -67,6 +71,14 @@ case ${DETECTOR} in NM=32 DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF07.json ;; +'JF09T01V01') + NM=1 + DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF09.json + ;; +'JF10T01V01') + NM=1 + DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF10.json + ;; 'JF13T01V01') NM=1 DET_CONFIG_FILE=/gpfs/photonics/swissfel/buffer/config/stream-JF13.json @@ -141,6 +153,51 @@ echo $((date3-date2)) | awk '{print int($1/60)":"int($1%60)}' if [ ${JF_CONVERSION} == 0 ] then echo "File is written in raw format, no compression" + + dir_name=`dirname ${OUTFILE_RAW}` + base_name=`basename ${dir_name}` + + if [ ${base_name} == "JF_pedestals" ] + then + echo "Pedestal run will make conversion" + + export PATH=/home/dbe/miniconda3/bin:$PATH + + source /home/dbe/miniconda3/etc/profile.d/conda.sh + + conda deactivate + conda activate sf-daq + + if [ ${DETECTOR} == "JF07T32V01" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --add_pixel_mask /sf/bernina/config/jungfrau/pixel_mask/JF07T32V01/pixel_mask_13_full.h5 + elif [ ${DETECTOR} == "JF03T01V02" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --add_pixel_mask /sf/bernina/config/jungfrau/pixel_mask/JF03T01V02/pixel_mask_half_chip.h5 + elif [ ${DETECTOR} == "JF02T09V02" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --number_bad_modules=1 + elif [ ${DETECTOR} == "JF06T08V02" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --add_pixel_mask /sf/alvra/config/jungfrau/pixel_mask/JF06T08V01/mask_2lines_module3.h5 +# elif [ ${DETECTOR} == "JF06T32V02" ] +# then +# time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --add_pixel_mask /sf/alvra/config/jungfrau/pixel_mask/JF06T32V02/mask_noise_in_28.h5 + elif [ ${DETECTOR} == "JF13T01V01" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --add_pixel_mask /sf/bernina/config/jungfrau/pixel_mask/JF13T01V01/pixel_mask_bad_rb_22.09.2020.h5 + elif [ ${DETECTOR} == "JF11T04V01" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --number_bad_modules=2 + elif [ ${DETECTOR} == "JF11T04V01" ] + then + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG --number_bad_modules=1 + else + time python /home/dbe/git/sf_daq_buffer/scripts/jungfrau_create_pedestals.py --filename ${OUTFILE_RAW} --directory ${dir_name} --verbosity DEBUG + fi + + fi + else echo "Will call compression/convertion ${OUTFILE_RAW} --> ${OUTFILE}"