From 81d4ed71220849eb18860e1f9649b1e8085a8663 Mon Sep 17 00:00:00 2001 From: Assmann Greta Marie Date: Thu, 27 Mar 2025 14:13:32 +0100 Subject: [PATCH] add new geom file --- src/10.cell | 12 ++++ src/10_jf.geom | 38 +++++++++++++ src/clara.py | 151 ++++++++++++++++++++++++++++++++++++++----------- src/results.py | 70 +++++++++++++++++------ 4 files changed, 221 insertions(+), 50 deletions(-) create mode 100644 src/10.cell create mode 100644 src/10_jf.geom diff --git a/src/10.cell b/src/10.cell new file mode 100644 index 0000000..9c72dbc --- /dev/null +++ b/src/10.cell @@ -0,0 +1,12 @@ +CrystFEL unit cell file version 1.0 + +lattice_type = cubic +centering = P + + +a = 1.0 A +b = 2.0 A +c = 3.0 A +al = 90.0 deg +be = 90.0 deg +ga = 90.0 deg diff --git a/src/10_jf.geom b/src/10_jf.geom new file mode 100644 index 0000000..9e3d886 --- /dev/null +++ b/src/10_jf.geom @@ -0,0 +1,38 @@ +; PSI JF8M + + +; Camera length (in m) and photon energy (eV) +clen = 0.15150000154972076 +photon_energy = 12000.0 +flag_lessthan = -32768 + +adu_per_eV = 0.00008285 +res = 13333.3 ; 75.00000356230885 micron pixel size + +rigid_group_0 = 0 +rigid_group_collection_0 = 0 + +; These lines describe the data layout for the JF native multi-event files +dim0 = % +dim1 = ss +dim2 = fs +data = /entry/data/data + + +;mask_file =/sf/cristallina/data/p21981/raw/run0010-6084triggers10wait/data/acq_datmaster.h5 +;mask = /entry/instrument/detector/pixel_mask +;mask_good = 0x0 +;mask_bad = 0xFFFFFFFF + +; corner_{x,y} set the position of the corner of the detector (in pixels) +; relative to the beam + +0/min_fs = 0 +0/min_ss = 0 +0/max_fs =3212 +0/max_ss =3332 +0/corner_x = -1613.0 +0/corner_y = -1666.0 +0/fs = x +0/ss = y + diff --git a/src/clara.py b/src/clara.py index e1731b8..125b2a7 100644 --- a/src/clara.py +++ b/src/clara.py @@ -11,14 +11,19 @@ import subprocess as sub import sys import time import zmq +import re from pathlib import Path from loguru import logger import receive_msg + + +#sys.path.insert(0, os.path.expanduser('/sf/cristallina/applications/mx/jungfraujoch_openapi_python/sf_datafiles')) +#from sfdata import SFDataFiles #define log file place: -LOG_FILENAME = time.strftime("/sf/cristallina/data/p21981/res/clara_%Y%m.log") +LOG_FILENAME = time.strftime("/sf/cristallina/data/p22217/res/clara_%Y%m.log") logger.add(LOG_FILENAME, level="INFO", rotation="100MB") @@ -113,7 +118,7 @@ class CollectedH5: except Exception as e: logger.info("Could not cd into processing directory: {}".format(e)) - return None + return out_path def convert_spg_num(self, sg: int): @@ -304,8 +309,8 @@ class CollectedH5: f2.write("flag_lessthan = " + str(self.message["underload"]) + "\n") f2.write("\n") #f2.write("adu_per_eV = 0.00008065\n") - f2.write("adu_per_eV = 0.00008285\n") - # f2.write("adu_per_photon = 1\n") + #f2.write("adu_per_eV = 1.0\n") + f2.write("adu_per_photon = 1\n") f2.write("res = 13333.3 ; " + str(self.message["pixel_size_m"] *1000000) + " micron pixel size\n") f2.write("\n") f2.write("rigid_group_0 = 0 \n") @@ -353,6 +358,36 @@ class CollectedH5: f2.close() return None + def copy_and_edit_geom_file_from_master(self, out_path): + output_file_path = os.path.join(out_path, 'JF17T16V01.geom') + geom_file_location = str(self.message["user_data"]["geometry_file_location"]) + pgroup = str(self.message["experiment_group"]) + res = "res" + geom_file_path = pa / pgroup / res / geom_file_location + #geom_file_location = '/sf/cristallina/data/p21981/work/crystfel_test/JF17T16V01_jfj.geom' + with open(geom_file_path, 'r') as file: + lines=file.read() + detector_distance = str(float(self.message["detector_distance_m"])-0.0002) + photon_energy = str(self.message["incident_energy_eV"]) + pattern_energy =r'photon_energy\s=\s(\d+)\seV\n' + energy = f'photon_energy = {photon_energy} eV\n' + lines_new = re.sub(pattern_energy,energy, lines) + pattern_clen =r'clen\s=\s(\d+.\d+)\n' + clen = f'clen = {detector_distance}\n' + lines_final = re.sub(pattern_clen, clen, lines_new) + #clen = str(self.message["detector_distance_m"]) + #photon_energy = str(self.message["incident_energy_eV"]) + #lines[0] = f'photon_energy = {photon_energy} eV\n' + #lines[1] = f'clen = {clen}\n' + + #lines[0] = 'photon_energy = 11000 eV\n' + #lines[1] = 'clen = 0.116\n' + + with open(output_file_path, 'w') as file: + file.writelines(lines_final) + + return output_file_path + def create_list_file(self): """ Function to generate a list file with the path of the input H5 file @@ -377,12 +412,14 @@ class CollectedH5: bsdata_name = Path(filen[:-15]+'.BSDATA.h5') bsdata_path = slash / bsdata_name - - print(bsdata_path) + + logger.info("Waiting for bs_data to exist") + while not bsdata_path.exists(): time.sleep(1) - print('survived the while loop') try: + time.sleep(10) + #bsdata = SFDataFiles(bsdata_path) bsdata = h5py.File(bsdata_path, "r") #r"/sf/cristallina/data/p21630/raw/run0065-lov_movestop_normal_1/data/acq0001.BSDATA.h5" except Exception as e: print(f"didn't open bsdata due to error {e}") #_logger.error(f"Cannot open {data_file} (due to {e})") @@ -396,6 +433,7 @@ class CollectedH5: try: x = h5py.File(jf_path, "r") + #data = SFDataFiles(bsdata_path, jf_path) except Exception as e: print(f"didn't open JF17T16V01.h5 due to error {e}") #_logger.error(f"Cannot open {data_file} (due to {e})") return @@ -406,29 +444,71 @@ class CollectedH5: index_dark = [] index_light = [] blanks = [] - - for i in range(n_pulse_id): + #print(evt_set) - p = pulseids_JF[i] - q = pulseids_BS[i] - if p != q: - event_i = np.where(pulseids_JF == q)[0] - event_i = i - else: - event_i=i - - events=evt_set[event_i] + _inters, ind_JF, ind_BS = np.intersect1d(pulseids_JF, pulseids_BS, return_indices=True) + #evt_set = evt_set[ind_BS] + #print(len(evt_set)) + + #print(min(ind_BS),max(ind_BS),len(ind_BS)) + for idx_JF, idx_BS in zip(ind_JF, ind_BS): + events = evt_set[idx_BS] + trigger_event = int(self.message["user_data"]["trigger_event"]) - if events[trigger_event] and events[200]: - index_light.append(i) + if trigger_event == 63: + + if events[trigger_event] and events[200]: + index_dark.append(idx_JF) - elif events[200]: - index_dark.append(i) + elif events[200]: + index_light.append(idx_JF) - else: - blanks.append(i) + else: + blanks.append(idx_JF) + + else: + + if events[trigger_event] and events[200]: + index_light.append(idx_JF) + + elif events[200]: + index_dark.append(idx_JF) + + else: + blanks.append(idx_JF) + + + #for i in range(n_pulse_id): + + # p = pulseids_JF[i] + # q = pulseids_BS[i] + # if p != q: + # #event_i = np.where(pulseids_JF == q)[0] for normal detector + # event_i = np.where(pulseids_BS == p)[0][0]#[::2] #[0][0] + # #event_i = i + # else: + # event_i=i# + + # events=evt_set[event_i] + #print(events) + #print(len(events)) + #print(event_i) + #print(p,q) + #print(i) + + # trigger_event = int(self.message["user_data"]["trigger_event"]) + + # if self.message["user_data"]["trigger_flag"]: + # if events[trigger_event] and events[200]: + # index_light.append(i) + + # elif events[200]: + # index_dark.append(i) + + # else: + # blanks.append(i) bsdata.close() x.close() @@ -470,7 +550,7 @@ class CollectedH5: return None - def create_slurm_script(self,trigger): + def create_slurm_script(self,trigger, geometry_file_path): """ Creates the input SLURM file with the following info: SLURM parameters ( CPUS , nodes, etc) @@ -499,6 +579,7 @@ class CollectedH5: #f.write("#SBATCH --nodelist=sf-cn-21 \n") # uncomment if on RA # f.write("#SBATCH --partition=hour \n") + #f.write("#SBATCH --account=p22217 \n") f.write("#SBATCH --cpus-per-task=32 \n") # f.write("#SBATCH --output=" + LOG_FILENAME + "\n") # f.write("#SBATCH --open-mode=append \n") @@ -508,6 +589,7 @@ class CollectedH5: f.write("module purge \n") f.write("module use MX unstable \n") f.write("module load crystfel/0.10.2-rhel8 \n") + #f.write("module load crystfel/0.11.1 \n") # f.write( # "module load crystfel/0.10.1-2 xgandalf/2018.01 HDF5_bitshuffle/2018.05 HDF5_LZ4/2018.05 gcc/4.8.5 hdf5_serial/1.10.3 \n" # ) @@ -539,8 +621,11 @@ class CollectedH5: + ".list -o " + data_file_name + ".stream -g " - + str(self.message["run_number"]) - + "_jf.geom" + + f"{geometry_file_path}" + #+ " --mille" + #+ " --mille-dir=mille-dir" + #+ str(self.message["run_number"]) + #+ "_jf.geom" + " -j `nproc` --min-res=75 " ) if self.message["user_data"]["retry"]: @@ -563,9 +648,9 @@ class CollectedH5: #+ "python /sls/MX/applications/git/vdp/src/results.py " + data_file_name + ".stream " - # + data_file_name - # + ".log " - ) + + trigger + + ) f.close() @@ -586,6 +671,7 @@ class CollectedH5: try: slurm_out = sub.run(["sbatch", "run_SLURM_" + trigger ], capture_output=True) + time.sleep(1) txt = slurm_out.stdout.decode().split() # grep the slurm number logger.info(f"submitted batch job number: {txt[-1]}") @@ -706,11 +792,12 @@ if __name__ == "__main__": mess_dec = json.loads(message.decode("utf-8")) logger.info(f"Received message is: {mess_dec}") mess_inp = CollectedH5(mess_dec) - mess_inp.mk_cd_output_dir_bl() + out_path = mess_inp.mk_cd_output_dir_bl() logger.info("output folder created") mess_inp.create_cell_file() logger.info("cell file created") mess_inp.create_geom_from_master() + geometry_file_path = mess_inp.copy_and_edit_geom_file_from_master(out_path) logger.info("geom file created") mess_inp.create_list_file() logger.info("list files created") @@ -720,7 +807,7 @@ if __name__ == "__main__": trigger_list.append("on") print(trigger_list) for i in trigger_list: - mess_inp.create_slurm_script(i) + mess_inp.create_slurm_script(i, geometry_file_path) mess_inp.submit_job_to_slurm(i) diff --git a/src/results.py b/src/results.py index a16b0ad..1c4c589 100644 --- a/src/results.py +++ b/src/results.py @@ -14,7 +14,7 @@ import MxdbVdpTools # import matplotlib.pyplot as plt #LOG_FILENAME = time.strftime("/sf/cristallina/applications/mx/clara_tools/log/clara_%Y%m.log") -LOG_FILENAME = time.strftime("/sf/cristallina/data/p21981/res/clara_%Y%m.log") +LOG_FILENAME = time.strftime("/sf/cristallina/data/p22226/res/clara_%Y%m.log") logger.add(LOG_FILENAME, level="INFO", rotation="100MB") @@ -152,8 +152,8 @@ def stream_to_dictionary(streamfile): ln += 1 # loop over the next 5 lines to get the diffraction resolution - line = loop_over_next_n_lines(text_file, 5).split() - ln += 5 + line = loop_over_next_n_lines(text_file, 5).split() # 5 for crystfel 10 # 7 for crystfel 11 + ln += 7 # multiple lattices not done yet if line[0] == "predict_refine/det_shift": @@ -161,9 +161,10 @@ def stream_to_dictionary(streamfile): tmp_crystal["det_shift_y"] = line[6] line = loop_over_next_n_lines(text_file, 1).split() ln += 1 - - tmp_crystal["diffraction_resolution_limit [nm^-1]"] = float(line[2]) - tmp_crystal["diffraction_resolution_limit [A]"] = float(line[5]) + + tmp_crystal["diffraction_resolution_limit_nm_inverse"] = float(line[2]) + #print(line) + tmp_crystal["diffraction_resolution_limit_angstrom"] = float(line[5]) #[5] for crystfel 0.11.1 [5] fpr 0.10.2-rehl8 # get the number of predicted reflections num_reflections = int(text_file.readline().split()[-1]) @@ -243,7 +244,10 @@ def get_data_from_streamfiles(): # parse stream into dict #TODO get on or off from streamfile naming and put it in output message parsed_stream = stream_to_dictionary(sys.argv[1]) + old_message["trigger_status"]= sys.argv[2] old_message["numberOfImages"] = len(parsed_stream) + + old_message["user_tag"] = old_message["user_data"]["file_prefix"].split('-')[1] # get number of indexable frames- not accounted for multilattice # print(parsed_stream[0].keys()) @@ -255,10 +259,13 @@ def get_data_from_streamfiles(): spots_per_frame = np.zeros((len(parsed_stream))) indexable_lattices = 0 for i in range(0, len(parsed_stream)): - # get spots per indexable frames: - spots_per_frame[i] = parsed_stream[i]["num_peaks"] - # get total number of indexable lattices - indexable_lattices = indexable_lattices + parsed_stream[i]["multiple_lattices"] + try: + # get spots per indexable frames: + spots_per_frame[i] = parsed_stream[i]["num_peaks"] + # get total number of indexable lattices + indexable_lattices = indexable_lattices + parsed_stream[i]["multiple_lattices"] + except: + logger.warning(f"Error: missing image from stream file: {i}") # put into dict for results, convert to list to be json serializable old_message["numberOfSpotsPerImage"] = (spots_per_frame.astype(int)).tolist() @@ -267,24 +274,32 @@ def get_data_from_streamfiles(): # get number of indexed lattices per pattern lattices_per_indx_frame = {} for i in indexable_frames: - lattices_per_indx_frame[int(i)] = parsed_stream[i]["multiple_lattices"] + lattices_per_indx_frame[int(i)] = parsed_stream[i]["multiple_lattices"] + old_message["numberOfLatticesPerImage"] = lattices_per_indx_frame # mean beam center shift X and Y list_x = [] list_y = [] + resolution_limit_list = [] + reflections_list = [] # define np.array uc_array = np.zeros((indexable_lattices, 6)) b = 0 for i in indexable_frames: for x in range(0, len(parsed_stream[i]["crystals"])): # det shift in x and y - list_x.append((parsed_stream[i]["crystals"][x]["det_shift_x"])) - list_y.append((parsed_stream[i]["crystals"][x]["det_shift_y"])) - # unit cell constants - uc_array[b] = np.asarray((parsed_stream[i]["crystals"][x]["Cell parameters"])) - b = b + 1 - + try: + list_x.append((parsed_stream[i]["crystals"][x]["det_shift_x"])) + list_y.append((parsed_stream[i]["crystals"][x]["det_shift_y"])) + # unit cell constants + uc_array[b] = np.asarray((parsed_stream[i]["crystals"][x]["Cell parameters"])) + resolution_limit_list.append((parsed_stream[i]["crystals"][x]["diffraction_resolution_limit_angstrom"])) + reflections_list.append((parsed_stream[i]["crystals"][x]["num_predicted_reflections"])) + b = b + 1 + except: + logger.warning(f"Error: missing image from stream file: {i}") + if len(indexable_frames) == 0: old_message["beamShiftMeanX_pxl"] = 0 old_message["beamShiftMeanY_pxl"] = 0 @@ -292,6 +307,10 @@ def get_data_from_streamfiles(): old_message["beamShiftStdY_pxl"] = 0 old_message["unitCellIndexingMean"] = 0 old_message["unitCellIndexingStd"] = 0 + old_message["resolutionLimitMean"] = 0 + old_message["resolutionLimitStd"] = 0 + old_message["numberReflectionsMean"] = 0 + old_message["numberReflectionsStd"] = 0 return old_message # ------ DET SHIFT MEAN and STD------- @@ -304,6 +323,7 @@ def get_data_from_streamfiles(): std_x = np.around(np.std(np.asarray(list_x).astype(float)), 4) mean_y = np.around(np.mean(np.asarray(list_y).astype(float)), 4) std_y = np.around(np.std(np.asarray(list_y).astype(float)), 4) + # convert to pixel unit # 0.075 mm = 1 pixel =75 um # print(mean_x, mean_x * (1/0.075), std_x, std_x * (1/0.075), "x") @@ -320,7 +340,21 @@ def get_data_from_streamfiles(): #convert to list to be json serializable old_message["unitCellIndexingMean"] = (np.around(mean_uc, 3)).tolist() old_message["unitCellIndexingStd"] = (np.around(std_uc, 3)).tolist() - + + #calculate diffraction resolution limit + mean_res = np.around(np.mean(np.asarray(resolution_limit_list).astype(float)), 4) + std_res = np.around(np.std(np.asarray(resolution_limit_list).astype(float)), 4) + old_message["resolutionLimitMean"] = mean_res + old_message["resolutionLimitStd"] = std_res + + #number of predicted reflections + mean_refs = np.around(np.mean(np.asarray(reflections_list).astype(float)), 4) + std_refs = np.around(np.std(np.asarray(reflections_list).astype(float)), 4) + old_message["numberReflectionsMean"] = mean_refs + old_message["numberReflectionsStd"] = std_refs + + #number of spots/peaks per image + # print(old_message) return old_message