diff --git a/autodeer/DEER_analysis.py b/autodeer/DEER_analysis.py index 81a141a..51bd9d7 100644 --- a/autodeer/DEER_analysis.py +++ b/autodeer/DEER_analysis.py @@ -14,9 +14,13 @@ import numbers import scipy.signal as sig from scipy.optimize import minimize,brute +import scipy.optimize as opt from autodeer.colors import primary_colors, ReIm_colors from autodeer.delay_optimise import calculate_optimal_tau import scipy.signal as signal +import dill + +opt._pickle = dill log = logging.getLogger('autoDEER.DEER') @@ -576,7 +580,7 @@ def DEERanalysis_plot_pub(results, ROI=None, fig=None, axs=None): axs[0].plot(results.t,background_func(results.t, results), '--',alpha=1,color='#42A399', lw=2) - axs[0].set_xlabel('Time / $\mu s$') + axs[0].set_xlabel(r'Time / $\mu s$') axs[0].xaxis.tick_top() axs[0].xaxis.set_label_position('top') @@ -815,10 +819,95 @@ def functional(f_axis,fieldsweep,A,B,filter=None,A_shift=0,B_shift=0): return -1*(Aspins * Bspins)/np.sqrt(Noise) +def shift_freq(pulse,shift): + # shift the frequency of a pulse by a given amount + if hasattr(pulse,'freq'): + pulse.freq.value += shift + else: + pulse.init_freq.value += shift + pulse.final_freq.value += shift + +def Functional_basic(x, pump_pulse, exc_pulse, ref_pulse, Fieldsweep, resonator, num_ref_pulses): + """ + Parameters + ---------- + x : list + The frequency shifts of the pump and excitation pulses. + pump_pulse : ad.Pulse + The pump pulse object + exc_pulse : ad.Pulse + The excitation pulse object + ref_pulse : ad.Pulse + The refocusing pulse object + Fieldsweep : ad.FieldSweepAnalysis + The FieldSweep analysis object + resonator : ad.ResonatorProfile + The resonator profile + num_ref_pulses : int + The number of refocusing pulses + + Returns + ------- + float + The functional value after optimisation + + """ + pump_freq = x[0] + obs_freq = x[1] + tmp_pump = pump_pulse.copy() + tmp_exc = exc_pulse.copy() + tmp_ref = ref_pulse.copy() + + shift_freq(tmp_pump, pump_freq) + shift_freq(tmp_exc, obs_freq) + shift_freq(tmp_ref, obs_freq) + + return calc_functional(Fieldsweep, tmp_pump, tmp_exc, tmp_ref, resonator=resonator, num_ref_pulses=num_ref_pulses) + +def Functional_spectrum_shift(x, pump_pulse, exc_pulse, ref_pulse, Fieldsweep, resonator, num_ref_pulses): + pump_freq = x[0] + obs_freq = x[1] + spectrum_shift = x[2] + + tmp_pump = pump_pulse.copy() + tmp_exc = exc_pulse.copy() + tmp_ref = ref_pulse.copy() + + shift_freq(tmp_pump, pump_freq) + shift_freq(tmp_exc, obs_freq) + shift_freq(tmp_ref, obs_freq) + + return calc_functional(Fieldsweep, tmp_pump, tmp_exc, tmp_ref, resonator=resonator, num_ref_pulses=num_ref_pulses, spectrum_shift=spectrum_shift) + + + +def optimise_pulses2(Fieldsweep, pump_pulse, exc_pulse, ref_pulse=None, resonator=None, num_ref_pulses=2, verbosity=0, method='basic', **kwargs): + + if method == 'basic': + res = minimize(lambda x: -Functional_basic(x, pump_pulse, exc_pulse, ref_pulse, Fieldsweep, resonator, num_ref_pulses), [0, 0], bounds=[(-0.15, 0.05), (-0.05, 0.1)]) + + elif method == 'spectrum_shift': + res = minimize(lambda x: -Functional_spectrum_shift(x, pump_pulse, exc_pulse, ref_pulse, Fieldsweep, resonator, num_ref_pulses), [0, 0, 0], bounds=[(-0.15, 0.05), (-0.05, 0.1), (-0.02, 0)]) + + if verbosity > 0: + print(f"Functional: {res.fun:.3f} \t Pump shift: {res.x[0]*1e3:.1f}MHz \t Exc/Ref shift: {res.x[1]*1e3:.1f}MHz") + if len(res.x) > 2: + print(f"Spectrum shift: {res.x[2]*1e3:.1f}MHz") + print(f"Number of function evaluations: {res.nfev}") + + new_pump = pump_pulse.copy() + new_exc = exc_pulse.copy() + new_ref = ref_pulse.copy() + shift_freq(new_pump, res.x[0]) + shift_freq(new_exc, res.x[1]) + shift_freq(new_ref, res.x[1]) + + return {'pump_pulse': new_pump, 'exc_pulse': new_exc, 'ref_pulse': new_ref}, res.fun + def optimise_pulses(Fieldsweep, pump_pulse, exc_pulse, ref_pulse=None, filter=None, verbosity=0, method='brute', nDEER=False, num_ref_pulses=2, full_output=False, resonator=None, **kwargs): - """Optimise the pulse positions to maximise the pump-exc overlap. + r"""Optimise the pulse positions to maximise the pump-exc overlap. Parameters ---------- @@ -945,9 +1034,147 @@ def optimise(filter): return new_pump_pulse, new_exc_pulse, new_ref_pulse, funct, grid, Jout else: return new_pump_pulse, new_exc_pulse, new_ref_pulse, +def shift_m11_01(x): + """ + Shifts the input array from the range [-1, 1] to the range [0, 1]. + + Parameters: + x (numpy.ndarray): Input array to be shifted. + + Returns: + numpy.ndarray: Shifted array in the range [0, 1]. + """ + return (-x + 1) / 2 +def build_profile(pulses,freqs=None,resonator=None): + + if freqs is None: + freqs = np.linspace(-0.3,0.3,100) + + excite_profile = np.ones_like(freqs) + + if isinstance(pulses, list): + for pulse in pulses: + tmp_excite_profile = pulse.exciteprofile(freqs=freqs,resonator=resonator)[:,2].real + if pulse.flipangle.value == np.pi: + tmp_excite_profile = shift_m11_01(tmp_excite_profile) + elif pulse.flipangle.value == np.pi/2: + tmp_excite_profile = -1*tmp_excite_profile+ 1 + excite_profile *= tmp_excite_profile + else: + excite_profile = pulses.exciteprofile(freqs=freqs,resonator=resonator)[:,2].real + if pulses.flipangle.value == np.pi: + excite_profile = shift_m11_01(excite_profile) + elif pulses.flipangle.value == np.pi/2: + excite_profile = -1*excite_profile+ 1 + + return excite_profile + +def calc_functional(Fieldsweep, pump_pulse, exc_pulse, ref_pulse,resonator=None, num_ref_pulses=2, spec_shift=0,**kwargs): + """ + + Parameters + ---------- + Fieldsweep : ad.FieldSweepAnalysis + The FieldSweep analysis object + pump_pulse : ad.Pulse + The pump pulse object + exc_pulse : ad.Pulse + The excitation pulse object + ref_pulse : ad.Pulse + The refocusing pulse object + respro : ad.ResonatorProfileAnalysis, optional + The resonator profile, by default None + num_ref_pulses : int, optional + The total number of refocusing pulses, by default 2 + Returns + ------- + float + The estimated modulation depth + """ + if 'respro' in kwargs: + resonator = kwargs['respro'] + print('Warning: respro is deprecated, use resonator instead') + + fieldsweep_fun = Fieldsweep.func_freq + f = np.linspace(-0.3,0.3,100) + fieldsweep_profile = fieldsweep_fun(f) + + fieldsweep_profile = resample_and_shift_vector(fieldsweep_profile, f, spec_shift) + fieldsweep_profile /= np.trapz(fieldsweep_profile,f) # normalise the fieldsweep profile + + + obs_sequence = [exc_pulse] + for i in range(num_ref_pulses): + obs_sequence.append(ref_pulse) + + pump_profile = build_profile(pump_pulse,f,resonator) + exc_profile = build_profile(obs_sequence,f,resonator) + dead_profile = pump_profile * exc_profile + pump_profile -= dead_profile + exc_profile -= dead_profile + pump_profile[pump_profile <0] = 0 + exc_profile[exc_profile <0] = 0 + + P_pump = np.trapz(pump_profile*fieldsweep_profile,f) / np.trapz(fieldsweep_profile,f) + P_obs = np.trapz(exc_profile*fieldsweep_profile,f) / np.trapz(fieldsweep_profile,f) + P_none = 1 - P_pump - P_obs + + return 2*P_pump*P_obs + +def calc_est_modulation_depth(Fieldsweep, pump_pulse, exc_pulse, ref_pulse,respro=None, num_ref_pulses=2): + """ + Calculate the estimated modulation depth from the EPR spectrum, the pulses and the resonator profile. + + .. math:: + \lambda = \frac{p_{\text{mod}}}{p_{\text{mod}}+p_{\text{un-mod}}} = \frac{2p_{\text{pump}}p_{\text{obs}}}{2p_{\text{pump}}p_{\text{obs}}+2p_{\text{obs}}p_{\text{none}} +p_{\text{obs}}^2 } -def plot_overlap(Fieldsweep, pump_pulse, exc_pulse, ref_pulse, filter=None, respro=None, num_ref_pulses=2, axs=None, fig=None): + + Parameters + ---------- + Fieldsweep : ad.FieldSweepAnalysis + The FieldSweep analysis object + pump_pulse : ad.Pulse + The pump pulse object + exc_pulse : ad.Pulse + The excitation pulse object + ref_pulse : ad.Pulse + The refocusing pulse object + respro : ad.ResonatorProfileAnalysis, optional + The resonator profile, by default None + num_ref_pulses : int, optional + The total number of refocusing pulses, by default 2 + + Returns + ------- + float + The estimated modulation depth + """ + resonator = respro + fieldsweep_fun = Fieldsweep.func_freq + f = np.linspace(-0.3,0.3,100) + fieldsweep_profile = fieldsweep_fun(f) + + obs_sequence = [exc_pulse] + for i in range(num_ref_pulses): + obs_sequence.append(ref_pulse) + + pump_profile = build_profile(pump_pulse,f,resonator) + exc_profile = build_profile(obs_sequence,f,resonator) + dead_profile = pump_profile * exc_profile + pump_profile -= dead_profile + exc_profile -= dead_profile + pump_profile[pump_profile <0] = 0 + exc_profile[exc_profile <0] = 0 + + P_pump = np.trapz(pump_profile*fieldsweep_profile,f) / np.trapz(fieldsweep_profile,f) + P_obs = np.trapz(exc_profile*fieldsweep_profile,f) / np.trapz(fieldsweep_profile,f) + P_none = 1 - P_pump - P_obs + + return (2*P_pump*P_obs)/(2*P_pump*P_obs + P_obs**2 + 2*P_obs*P_none) + + +def plot_overlap(Fieldsweep, pump_pulse, exc_pulse, ref_pulse, filter=None, resonator=None, num_ref_pulses=2, axs=None, fig=None,**kwargs): """Plots the pump and excitation profiles as well as the fieldsweep and filter profile. Parameters @@ -961,10 +1188,7 @@ def plot_overlap(Fieldsweep, pump_pulse, exc_pulse, ref_pulse, filter=None, resp The excitation pulse object ref_pulse : ad.Pulse, optional The refocusing pulse object, by default None - filter : str or number, optional - The filter profile if applicable, by default None. If it is a number a filter is generated with this cutoff frequency. - If the string 'Matched' is used a matched filter is used. - respro : ad.ResonatorProfileAnalysis, optional + resonator : ad.ResonatorProfileAnalysis, optional The resonator profile for fitting, by default None. The resonator profile must include the fit. num_ref_pulses : int, optional The total number of refocusing pulses, by default 2 @@ -975,47 +1199,43 @@ def plot_overlap(Fieldsweep, pump_pulse, exc_pulse, ref_pulse, filter=None, resp """ + if 'respro' in kwargs: + resonator = kwargs['respro'] + print('Warning: respro is deprecated, use resonator instead') - gyro = Fieldsweep.gyro + if filter is not None: + print('Warning: filter is deprecated, ignoring input') + fieldsweep_fun = Fieldsweep.func_freq - f = np.linspace(-0.4,0.4,100) - + f = np.linspace(-0.3,0.3,100) fieldsweep_profile = fieldsweep_fun(f) - - pump_Mz = normalise_01(-1*pump_pulse.exciteprofile(freqs=f)[:,2].real) - exc_Mz = normalise_01(-1*exc_pulse.exciteprofile(freqs=f)[:,2].real) - - if ref_pulse is not None: - for i in range(num_ref_pulses): - exc_Mz *= normalise_01(-1*ref_pulse.exciteprofile(freqs=f)[:,2].real) - + obs_sequence = [exc_pulse] + for i in range(num_ref_pulses): + obs_sequence.append(ref_pulse) - if filter == 'Matched': - # Matched filter - filter_profile = lambda f_new: np.interp(f_new,f,exc_Mz) - elif isinstance(filter, numbers.Number): - filter_profile = build__lowpass_butter_filter(filter) + pump_profile = build_profile(pump_pulse,f,resonator) + exc_profile = build_profile(obs_sequence,f,resonator) + dead_profile = pump_profile * exc_profile + pump_profile -= dead_profile + exc_profile -= dead_profile + pump_profile[pump_profile <0] = 0 + exc_profile[exc_profile <0] = 0 if axs is None and fig is None: fig, axs = plt.subplots(1,1,figsize=(5,5), layout='constrained') elif axs is None: - axs = fig.subplots(1,1,subplot_kw={},gridspec_kw={'hspace':0} - ) - - # Normalise the fieldsweep profile + axs = fig.subplots(1,1,subplot_kw={},gridspec_kw={'hspace':0}) + fieldsweep_profile /= np.abs(fieldsweep_profile).max() - - # Plot the profiles + axs.plot(f*1e3,fieldsweep_profile, label = 'Fieldsweep', c='k') - axs.fill_between(f*1e3,pump_Mz*fieldsweep_profile, label = 'Pump Profile', alpha=0.5,color='#D95B6F') - axs.fill_between(f*1e3,exc_Mz*fieldsweep_profile, label = 'Observer Profile',alpha=0.5,color='#42A399') - if filter is not None: - axs.plot(f*1e3,filter_profile(f),'--', label = 'Filter') + axs.fill_between(f*1e3,pump_profile*fieldsweep_profile, label = 'Pump Profile', alpha=0.5,color=primary_colors[0]) + axs.fill_between(f*1e3,exc_profile*fieldsweep_profile, label = 'Observer Profile',alpha=0.5,color=primary_colors[1]) - if respro is not None: - model_norm = respro.model / np.max(respro.model) - axs.plot((respro.model_x - Fieldsweep.LO)*1e3, model_norm,'--', label='Resonator Profile') + if resonator is not None: + model_norm = resonator.model / np.max(resonator.model) + axs.plot((resonator.model_x - Fieldsweep.LO)*1e3, model_norm,'--', label='Resonator Profile') fmin = f[~np.isclose(fieldsweep_profile,0)].min() fmax = f[~np.isclose(fieldsweep_profile,0)].max() @@ -1026,6 +1246,40 @@ def plot_overlap(Fieldsweep, pump_pulse, exc_pulse, ref_pulse, filter=None, resp return fig +def calc_dt_from_tau(deer_settings): + """ + Calculate the time step from the tau values in the DEER settings. + + Parameters + ---------- + deer_settings : dict + The DEER settings with keys 'tau1', 'tau2' and 'tau3' and 'ExpType'. + + Returns + ------- + deer_settings : dict + The DEER settings with the key 'dt' added. + """ + if deer_settings['ExpType'] == '4pDEER': + if deer_settings['tau2'] > 10: + deer_settings['dt'] = 12 + elif deer_settings['tau2'] > 20: + deer_settings['dt'] = 16 + else: + deer_settings['dt'] = 8 + + elif deer_settings['ExpType'] == '5pDEER': + if deer_settings['tau2']*2 > 10: + deer_settings['dt'] = 12 + elif deer_settings['tau2']*2 > 20: + deer_settings['dt'] = 16 + else: + deer_settings['dt'] = 8 + + return deer_settings + + + def calc_DEER_settings(relaxation_data,mode='auto', target_time=2, target_SNR=20, waveform_precision=2, corr_factor=1, rec_tau=None): @@ -1146,22 +1400,102 @@ def calc_DEER_settings(relaxation_data,mode='auto', target_time=2, deer_settings['tau1'] = rec_tau/2 log.info(f"Using recommended tau2: {rec_tau}us") - if deer_settings['ExpType'] == '4pDEER': - if deer_settings['tau2'] > 10: - deer_settings['dt'] = 12 - elif deer_settings['tau2'] > 20: - deer_settings['dt'] = 16 - else: - deer_settings['dt'] = 8 - elif deer_settings['ExpType'] == '5pDEER': - if deer_settings['tau2']*2 > 10: - deer_settings['dt'] = 12 - elif deer_settings['tau2']*2 > 20: - deer_settings['dt'] = 16 - else: - deer_settings['dt'] = 8 + calc_dt_from_tau(deer_settings) + return deer_settings + +def BackgroundAnalysis(dataset,model=dl.bg_hom3d,verbosity=0,**kwargs): + + Vexp:np.ndarray = dataset.data + + if np.iscomplexobj(Vexp): + Vexp,Vexp_im,_ = dl.correctphase(Vexp,full_output=True) + + Vexp /= Vexp.max() + + if 't' in dataset.coords: # If the dataset is a xarray and has sequence data + t = dataset['t'].data + if 'seq_name' in dataset.attrs: + if dataset.attrs['seq_name'] == "4pDEER": + exp_type = "4pDEER" + tau1 = dataset.attrs['tau1'] / 1e3 + tau2 = dataset.attrs['tau2'] / 1e3 + elif dataset.attrs['seq_name'] == "5pDEER": + exp_type = "5pDEER" + tau1 = dataset.attrs['tau1'] / 1e3 + tau2 = dataset.attrs['tau2'] / 1e3 + tau3 = dataset.attrs['tau3'] / 1e3 + elif dataset.attrs['seq_name'] == "3pDEER": + exp_type = "3pDEER" + tau1 = dataset.attrs['tau1'] / 1e3 + else: + if 'tau1' in dataset.attrs: + tau1 = dataset.attrs['tau1'] / 1e3 + elif 'tau1' in kwargs: + tau1 = kwargs.pop('tau1') + else: + raise ValueError("tau1 not found in dataset or kwargs") + if 'tau2' in dataset.attrs: + tau2 = dataset.attrs['tau2'] / 1e3 + elif 'tau2' in kwargs: + tau2 = kwargs.pop('tau2') + else: + raise ValueError("tau2 not found in dataset or kwargs") + if 'tau3' in dataset.attrs: + tau3 = dataset.attrs['tau3'] / 1e3 + exp_type = "5pDEER" + elif 'tau3' in kwargs: + tau3 = kwargs.pop('tau3') + exp_type = "5pDEER" + else: + exp_type = "4pDEER" + + else: + # Extract params from kwargs + if "tau1" in kwargs: + tau1 = kwargs.pop("tau1") + if "tau2" in kwargs: + tau2 = kwargs.pop("tau2") + if "tau3" in kwargs: + tau3 = kwargs.pop("tau3") + exp_type = kwargs.pop("exp_type") + # t = dataset.axes[0] + t = dataset['X'] + + + if hasattr(dataset,'attrs') and "pulse0_tp" in dataset.attrs: + # seach over all pulses to find the longest pulse using regex + pulselength = np.max([dataset.attrs[i] for i in dataset.attrs if re.match(r"pulse\d*_tp", i)]) + pulselength /= 1e3 # Convert to us + pulselength /= 2 # Account for the too long range permitted by deerlab + else: + if "pulselength" in kwargs: + pulselength = kwargs.pop("pulselength")/1e3 + else: + pulselength = 16/1e3 + + + if exp_type == "4pDEER": + pathways = [1] + experimentInfo = dl.ex_4pdeer(tau1=tau1,tau2=tau2,pathways=pathways,pulselength=pulselength) + elif exp_type == "5pDEER": + pathways = [1,5] + experimentInfo = dl.ex_fwd5pdeer(tau1=tau1,tau2=tau2,tau3=tau3,pathways=pathways,pulselength=pulselength) + elif exp_type == "3pDEER": + pathways = [1] + experimentInfo = dl.ex_3pdeer(tau=tau1,pathways=pathways,pulselength=pulselength) + + r= np.linspace(1,8,100) + Vmodel = dl.dipolarmodel(t, r, experiment=experimentInfo, Pmodel=model) + Vmodel.pathways = pathways + + + + + fit = dl.fit(Vmodel, Vexp, + verbose=verbosity, + **kwargs) diff --git a/autodeer/__init__.py b/autodeer/__init__.py index 42ee984..01d836b 100644 --- a/autodeer/__init__.py +++ b/autodeer/__init__.py @@ -1,5 +1,5 @@ from ._version import __version__ -from .DEER_analysis import DEERanalysis,DEERanalysis_plot, DEERanalysis_plot_pub,IdentifyROI,optimise_pulses,plot_overlap,normalise_01, calc_DEER_settings +from .DEER_analysis import DEERanalysis,DEERanalysis_plot, DEERanalysis_plot_pub,IdentifyROI,optimise_pulses,plot_overlap,normalise_01, calc_DEER_settings, calc_dt_from_tau from .delay_optimise import * from .sequences import * from .criteria import * diff --git a/autodeer/criteria.py b/autodeer/criteria.py index 4dcfeed..87ee74e 100644 --- a/autodeer/criteria.py +++ b/autodeer/criteria.py @@ -82,3 +82,43 @@ def test_func(data, verbosity=verbosity): return test super().__init__(name, test_func, description,**kwargs) + + +# class BackgroundCriteria(Criteria): + +# def __init__(self, SNR_target:float, model=None, verbosity=0, update_func=None, **kwargs)->None: +# """Criteria for terminating DEER background experiments. + +# Parameters +# ---------- +# SNR_target : float +# Target signal-to-noise ratio. + +# Returns +# ------- +# _type_ +# _description_ +# """ +# name = "BackgroundCriteria" +# description = "Criteria for terminating background experiments." + +# def test_func(data, verbosity=verbosity): + +# fit = BackgroundAnalysis(data, model=model, verbosity=verbosity,lin_maxiter=50,max_nfev=100) + +# if fit.SNR < SNR_target: +# test = False + +# if update_func is not None: +# update_func(fit) + +# test_msg = f"Test {self.name}: {test}\t - SNR:{fit.SNR}" +# log.debug(test_msg) +# if verbosity > 0: +# print(test_msg) + +# return test + + + +# super().__init__(name, test_func, description,**kwargs) \ No newline at end of file diff --git a/autodeer/gui/autoDEER_worker.py b/autodeer/gui/autoDEER_worker.py index 5d6e1a1..44f8c1e 100644 --- a/autodeer/gui/autoDEER_worker.py +++ b/autodeer/gui/autoDEER_worker.py @@ -61,7 +61,10 @@ class autoDEERWorker(QtCore.QRunnable): ''' - def __init__(self, interface, wait:QtCore.QWaitCondition, mutex:QtCore.QMutex,pulses:dict, results:dict, LO, gyro,AWG=True, user_inputs:dict = None, *args, **kwargs): + def __init__(self, interface, wait:QtCore.QWaitCondition, + mutex:QtCore.QMutex,pulses:dict, results:dict, LO, gyro, + AWG=True, user_inputs:dict = None, + operating_mode = 'auto',fixed_tau=None, *args, **kwargs): super(autoDEERWorker,self).__init__() # Store constructor arguments (re-used for processing) @@ -82,6 +85,8 @@ def __init__(self, interface, wait:QtCore.QWaitCondition, mutex:QtCore.QMutex,pu self.user_inputs = user_inputs self.stop_flag = False self.quick_deer_state = True + self.fixed_tau = fixed_tau + self.operating_mode = operating_mode self.max_tau = 3.5 print(f"Waveform Precision is: {get_waveform_precision()}") @@ -144,7 +149,7 @@ def run_fsweep(self): reptime = self.reptime p90, p180 = self.interface.tune_rectpulse(tp=self.tp, LO=LO, B=LO/gyro_N, reptime = reptime,shots=int(100*self.noise_mode)) shots = int(300*self.noise_mode) - shots = np.min([shots,100]) + shots = np.max([shots,100]) fsweep = FieldSweepSequence( B=LO/gyro_N, LO=LO,reptime=reptime,averages=50,shots=shots, Bwidth = 250, @@ -199,7 +204,7 @@ def run_CP_relax(self,dt=10,tmin=0.7,averages=30,autoStop=True,autoIFGain=True,* gyro = self.gyro reptime = self.reptime shots = int(100*self.noise_mode) - shots = np.min([shots,25]) + shots = np.max([shots,25]) relax = DEERSequence( B=LO/gyro, LO=LO,reptime=reptime,averages=averages,shots=shots, tau1=tmin/2, tau2=tmin/2, tau3=0.2, dt=15, @@ -235,7 +240,7 @@ def run_T2_relax(self,dt=10,tmin=0.5,averages=30,autoStop=True,autoIFGain=True,* gyro = self.gyro reptime = self.reptime shots = int(100*self.noise_mode) - shots = np.min([shots,25]) + shots = np.max([shots,25]) seq = T2RelaxationSequence( B=LO/gyro, LO=LO,reptime=reptime,averages=averages,shots=shots, @@ -299,6 +304,20 @@ def run_long_deer(self): end_signal = self.signals.longdeer_result.emit self.run_deer(total_crit,end_signal, dt=16,shot=50,averages=1e4) + def run_single_deer(self): + # Run a DEER experiment background measurement + self.signals.status.emit('Running DEER Background') + if 'autoStop' in self.deer_inputs and not self.deer_inputs['autoStop']: + SNR_crit = SNRCriteria(150,verbosity=2) + total_crit = [SNR_crit] + else: # autoStop is True + SNR_crit = SNRCriteria(150,verbosity=2) + total_crit = [SNR_crit, self.EndTimeCriteria] + end_signal = self.signals.longdeer_result.emit + self.run_deer(total_crit,end_signal, dt=16,shot=50,averages=1e4) + + + def run_deer(self,end_criteria,signal, dt=16,shot=50,averages=1000,): @@ -388,7 +407,7 @@ def run_reptime_opt(self): LO = self.LO p90, p180 = self.interface.tune_rectpulse(tp=self.tp, LO=LO, B=LO/self.gyro, reptime = reptime_guess,shots=int(100*self.noise_mode)) - n_shots = int(np.min([int(50*self.noise_mode),50])) + n_shots = int(np.max([int(50*self.noise_mode),50])) scan = ReptimeScan(B=LO/self.gyro, LO=LO,reptime=reptime_guess, reptime_max=20e3, averages=10, shots=n_shots, pi2_pulse=p90, pi_pulse=p180) self.interface.launch(scan,savename=f"{self.samplename}_reptimescan",) @@ -414,29 +433,27 @@ def tune_pulses(self): return 'skip' def _build_methods(self): - - seq = self.deer_inputs['ExpType'] - - if (self.deer_inputs['tau1'] == 0) and (self.deer_inputs['tau2'] == 0): - quick_deer = True - else: - quick_deer = False methods = [self.run_fsweep,self.run_reptime_opt,self.run_respro,self.run_fsweep, self.tune_pulses] - if (seq is None) or (seq == '5pDEER'): + if ( self.operating_mode is None) or ( self.operating_mode == '5pDEER'): methods.append(self.run_CP_relax) methods.append(self.run_T2_relax) - elif (seq == '4pDEER') or (seq == 'Ref2D'): + + elif ( self.operating_mode == '4pDEER') or ( self.operating_mode == 'Ref2D'): methods.append(self.run_CP_relax) methods.append(self.run_T2_relax) methods.append(self.run_2D_relax) + if self.operating_mode == 'Ref2D': + return methods - if seq == 'Ref2D': - return methods + elif self.operating_mode == 'single': + methods.append(self.run_T2_relax) + methods.append(self.run_single_deer) + return methods - if quick_deer: + if self.fixed_tau is None: methods.append(self.run_quick_deer) methods.append(self.run_long_deer) diff --git a/autodeer/gui/main.py b/autodeer/gui/main.py index e1c8639..7a6790d 100644 --- a/autodeer/gui/main.py +++ b/autodeer/gui/main.py @@ -222,7 +222,7 @@ def __init__(self): self.Min_tp=12 self.deer_settings = {'ESEEM':None, 'ExpType':'5pDEER'} - self.priorties = {'Auto': 150, 'MNR':300, 'Distance': 80} + self.priorties = {'Auto': 150, 'MNR':300, 'Distance': 80, 'Single':200} self.priotityComboBox.addItems(list(self.priorties.keys())) self.correction_factor=1 @@ -858,14 +858,22 @@ def refresh_relax(self, fitresult): def initialise_deer_settings(self): - aim_SNR = self.aim_MNR/(self.est_lambda*self.label_eff) - + # Assemble all relaxation data if (self.Exp_types.currentText() == '4pDEER'): exp = '4pDEER' else: exp = 'auto' + if self.userinput['priority'].lower() == 'single': + exp = '4pDEER' + aim_SNR = self.priorties[self.userinput['priority']] + self.aim_time = self.MaxTime.value() - ((time.time() - self.starttime) / (60*60)) # in hours + self.aim_MNR = aim_SNR + else: + aim_SNR = self.aim_MNR/(self.est_lambda*self.label_eff) + + relax_data = {} if 'relax' in self.current_results: relax_data['CP'] = self.current_results['relax'] @@ -910,9 +918,6 @@ def update_deer_settings(self): self.correction_factor = ad.calc_correction_factor(relax,self.current_results['quickdeer']) main_log.info(f"Correction factor {self.correction_factor:.3f}") - MNR_target = self.priorties[self.userinput['priority']] - SNR_target = MNR_target/(mod_depth) - main_log.info(f"SNR target {SNR_target:.2f}") # Assemble all relaxation data if (self.Exp_types.currentText() == '4pDEER'): @@ -920,6 +925,17 @@ def update_deer_settings(self): else: exp = 'auto' + if self.userinput['priority'].lower() == 'single': + SNR_target = self.priorties[self.userinput['priority']] + MNR_target = SNR_target + single_mode = True + exp = '4pDEER' + else: + MNR_target = self.priorties[self.userinput['priority']] + SNR_target = MNR_target/(mod_depth) + main_log.info(f"SNR target {SNR_target:.2f}") + + relax_data = {} if 'relax' in self.current_results: relax_data['CP'] = self.current_results['relax'] @@ -951,9 +967,10 @@ def update_tau_delays_figure(self, SNRs, MeasTimes, labels=None): fig = self.relax_canvas.figure axs = self.relax_ax[2] axs.cla() - # Only supports 5pDEER expand to 4pDEER - CP_analysis = self.current_results['relax'] - ad.plot_optimal_tau(CP_analysis,SNRs,MeasTimes,MaxMeasTime=36, labels=['5pDEER'],fig=fig,axs=axs,cmap=[epr.primary_colors[0]],corr_factor=self.correction_factor) + + if 'relax' in self.current_results: + CP_analysis = self.current_results['relax'] + ad.plot_optimal_tau(CP_analysis,SNRs,MeasTimes,MaxMeasTime=36, labels=['5pDEER'],fig=fig,axs=axs,cmap=[epr.primary_colors[0]],corr_factor=self.correction_factor) if 'relax2D' in self.current_results: ad.plot_optimal_tau(self.current_results['relax2D'],SNRs,MeasTimes,MaxMeasTime=36, labels=['4pDEER'],fig=fig,axs=axs,cmap=[epr.primary_colors[1]],corr_factor=self.correction_factor) @@ -1036,6 +1053,8 @@ def check_T2(self, fitresult): new_tmin = fitresult.axis[-1].values new_tmin += new_dt*1e-3 nAvgs = fitresult.dataset.attrs['nAvgs'] + + self.initialise_deer_settings() if self.worker is not None: self.worker.run_T2_relax(dt=new_dt,tmin=new_tmin,averages=nAvgs,autoStop=False) @@ -1109,7 +1128,7 @@ def update_longdeer(self, dataset=None): self.current_data['longdeer'] = dataset self.save_data(dataset,'DEER_final',folder='main') - self.Tab_widget.setCurrentIndex(4) + self.Tab_widget.setCurrentIndex(5) self.longDEER.current_data['quickdeer'] = dataset self.longDEER.update_inputs_from_dataset() self.longDEER.update_figure() @@ -1191,11 +1210,21 @@ def clear_all(self): def RunAutoDEER(self, advanced=False): self.clear_all() + self.operating_mode = 'auto' + self.fixed_tau = None + if self.spectromterInterface is None or self.connected is False: QMessageBox.about(self,'ERORR!', 'A interface needs to be connected first!') main_log.error('Could not run autoDEER. A interface needs to be connected first!') return None + + # Block the autoDEER buttons + self.FullyAutoButton.setEnabled(False) + self.AdvancedAutoButton.setEnabled(False) + self.resonatorComboBox.setEnabled(False) + self.fcDoubleSpinBox.setEnabled(False) + userinput = {} userinput['MaxTime'] = self.MaxTime.value() @@ -1211,13 +1240,40 @@ def RunAutoDEER(self, advanced=False): self.userinput = userinput - if advanced: - self.deer_settings['ExpType'] = self.Exp_types.currentText() - self.deer_settings['tau1'] = self.Tau1Value.value() - self.deer_settings['tau2'] = self.Tau2Value.value() - self.deer_settings['tau3'] = self.Tau3Value.value() + if self.priotityComboBox.currentText().lower() == 'single': + self.operating_mode = 'single' else: - self.deer_settings = {'ExpType':'5pDEER','tau1':0,'tau2':0,'tau3':0} + self.operating_mode = self.Exp_types.currentText() + + if advanced: + + if self.Tau1Value.value() > 0 or self.Tau2Value.value() > 0: + self.fixed_tau = {'tau1':self.Tau1Value.value(),'tau2':self.Tau2Value.value(),'tau3':self.Tau3Value.value()} + # Skip the relaxation data analysis and go straight to DEER + if self.Exp_types.currentText() == 'auto': + self.deer_settings['ExpType'] = '5pDEER' + elif self.Exp_types.currentText() == 'Ref2D': + print("ERROR: Ref2D not supported with fixed tau") + return None + else: + self.deer_settings['ExpType'] = self.Exp_types.currentText() + + self.deer_settings['tau1'] = self.Tau1Value.value() + self.deer_settings['tau2'] = self.Tau2Value.value() + self.deer_settings['tau3'] = self.Tau3Value.value() + self.deer_settings['criteria'] = self.priorties[self.userinput['priority']] + self.deer_settings['autoStop'] = self.Time_autoStop_checkbox.isChecked() + self.deer_settings = calc_dt_from_tau(self.deer_settings) + main_log.info(f"tau1 set to {self.deer_settings['tau1']:.2f} us") + main_log.info(f"tau2 set to {self.deer_settings['tau2']:.2f} us") + main_log.info(f"DEER Sequence set to {self.deer_settings['ExpType']}") + + else: + self.fixed_tau = None + + + + if not 'ESEEM' in self.deer_settings: self.deer_settings['ESEEM'] = None @@ -1226,11 +1282,6 @@ def RunAutoDEER(self, advanced=False): else: self.create_relax_figure(3) - # Block the autoDEER buttons - self.FullyAutoButton.setEnabled(False) - self.AdvancedAutoButton.setEnabled(False) - self.resonatorComboBox.setEnabled(False) - self.fcDoubleSpinBox.setEnabled(False) try: night_hours = self.config['autoDEER']['Night Hours'] @@ -1240,9 +1291,11 @@ def RunAutoDEER(self, advanced=False): self.waitCondition = QtCore.QWaitCondition() mutex = QtCore.QMutex() self.worker = autoDEERWorker( - self.spectromterInterface,wait=self.waitCondition,mutex=mutex, - pulses=self.pulses,results=self.current_results, AWG=self.AWG, LO=self.LO, gyro = self.gyro, - user_inputs=userinput, cores=self.cores,night_hours=night_hours) + self.spectromterInterface, wait=self.waitCondition, mutex=mutex, + pulses=self.pulses, results=self.current_results, AWG=self.AWG, LO=self.LO, gyro = self.gyro, + user_inputs=userinput, operating_mode = self.operating_mode, fixed_tau=self.fixed_tau, + cores=self.cores, night_hours=night_hours, + ) self.starttime = time.time() diff --git a/autodeer/gui/modetune.py b/autodeer/gui/modetune.py index 4878a07..60357f2 100644 --- a/autodeer/gui/modetune.py +++ b/autodeer/gui/modetune.py @@ -74,7 +74,7 @@ def start_mode_tune(self): seq.pulses[2].tp.value = 512 seq.evolution([freq]) - self.interface.launch(seq,savename="modetune",IFgain=1) + self.interface.launch(seq,savename="modetune") # Start thread self.worker = get_dataWorker(self.interface) @@ -121,34 +121,42 @@ def update_dip(self): x = dataset_pc.LO + dataset_pc.pulse0_freq y = dataset_pc.data - def gaussian(x, A, mu, sigma): - return A * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2)) - popt, pcov = curve_fit(gaussian, x, y, p0=[1, dataset_pc.LO, 0.1]) - fit_data = gaussian(x, *popt) - # check r^2 - residuals = y - fit_data - ss_res = np.sum(residuals ** 2) - ss_tot = np.sum((y - np.mean(y)) ** 2) - R2 = 1 - (ss_res / ss_tot) - - if R2 < 0.5: - #set QDoubleSpinBox color to red - self.ui.calc_freq.setStyleSheet("background-color: red") - elif R2 < 0.75: - #set QDoubleSpinBox color to yellow - self.ui.calc_freq.setStyleSheet("background-color: yellow") - else: - #set QDoubleSpinBox color to green - self.ui.calc_freq.setStyleSheet("background-color: green") - - self.ui.calc_freq.setValue(popt[1]) - - # update plot self.mode_ax.cla() - self.mode_ax.plot(x, gaussian(x, *popt), label='fit',color=primary_colors[0]) self.mode_ax.plot(x, y, '.', color='0.7', label='Data', ms=6) self.mode_ax.set_xlabel("Frequency (GHz)") self.mode_ax.set_ylabel("Intensity (a.u.)") + + def gaussian(x, A, mu, sigma): + return A * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2)) + try: + popt, pcov = curve_fit(gaussian, x, y, p0=[1, dataset_pc.LO, 0.1]) + except: + popt = None + + + if popt is not None: + fit_data = gaussian(x, *popt) + # check r^2 + residuals = y - fit_data + ss_res = np.sum(residuals ** 2) + ss_tot = np.sum((y - np.mean(y)) ** 2) + R2 = 1 - (ss_res / ss_tot) + + if R2 < 0.5: + #set QDoubleSpinBox color to red + self.ui.calc_freq.setStyleSheet("background-color: red") + elif R2 < 0.75: + #set QDoubleSpinBox color to yellow + self.ui.calc_freq.setStyleSheet("background-color: yellow") + else: + #set QDoubleSpinBox color to green + self.ui.calc_freq.setStyleSheet("background-color: green") + + self.ui.calc_freq.setValue(popt[1]) + + self.mode_ax.plot(x, gaussian(x, *popt), label='fit',color=primary_colors[0]) + + # update plot self.mode_ax.legend() self.modetTunecanvas.draw() pass