""" Tc_lib.py Critical temperature evaluation routines Lev Levitin, 24 May 2010 based on earlier code """ # Changelog # 24 May 2010 # fix error: load PLM only for a given time span, previously 5 minutes were added on either end # handle no data in superfluid/normal state correctly. # # 1 Jul 2010 # - bug fix in TcMeasTask.select_T_ranges() # - add an option of manually chosen output file in TcNMRTask # # 2 Jul 2010 # - bug fix in TInterp # # 5 Jul 2010 # add 'extended_output' option to Tc_lib.TcMeasTask.select_data, # change return format of Tc_lib.TcMeasTask.load_peak_data # # 11-14 Jul 2010 # various changes # - fixing misprints in the figure text # - fix a bug in spectra evaluation: the naive bulk peak used to be evaluated # using a spectrum with a slab Lorentzian fit subtracted, now the bare spectrum # is used instead. import aver, flib, deco, timelog, he3 import nmr, scopelog, timelog, nmrfit, massfft import numpy, pylab, os, matplotlib, matplotlib.dates, gc import datetime import plmfitlog skip_plm_points_at = [1.26499435e+09] class TInterp(object): """ time dependence of temperature interpolation: simple segment-linear interpolation """ def __init__(self, t, T, **vargs): """ t - time [sec since Unix Epoch] T - temperature [mK] """ self.t = t self.T = T def temperature(self, tstart, tend=None): """ Approximate temperature during periods tstart[k] <= t <= tend[k]. If 'tend' is not specified, it defaults to 'tstart'. Return a tuple of two arrays (T, dT), representing mean temperature and the uncertainty. """ if tend is None: tend = tstart Tmin = numpy.minimum(numpy.interp(tstart, self.t, self.T), numpy.interp(tend, self.t, self.T)) Tmax = numpy.maximum(numpy.interp(tstart, self.t, self.T), numpy.interp(tend, self.t, self.T)) T = 0.5 * (Tmax + Tmin) dT = 0.5 * abs(Tmax - Tmin) dT = numpy.maximum(dT, 0.001) return T, dT def ramp_rate(self, T): """ Return ramp rate while at temperature 'T'. This method calculates the mean ramp rate for the whole period being interpolated. """ return numpy.mean(numpy.diff(self.T) / numpy.diff(self.t)) @staticmethod def method_code(): "code of the temperature interpolation method for the log files"; return 0 def method_title(self): "the description of the interpolation method"; return "Linear Approximation" def need_aux_plot(self): "is there a need for an aux plot?"; return False def aux_plot_ylabel(self): "Y-axis label on the aux plot"; return '$T$ [mK]' def aux_plot_y(self, T): "function of temperature to represent in the aux plot"; return T def present(self, caption, tstart=None, tend=None, plm_sign=''): """ present the temperature interpolation in an A4 page. returns the figure object. 'title' - title for the page """ if len(plm_sign) > 0: plm_sign = ', ' + plm_sign fig = deco.portraitA4() if tstart is None or tend is None: tstart = tend = numpy.array([]) xmin = matplotlib.dates.epoch2num(tstart) xmax = matplotlib.dates.epoch2num(tend) ttt = numpy.concatenate([numpy.linspace(min(self.t), max(self.t), 200), self.t, tstart, tend]) ttt.sort() for sub, yfunc, ylabel, title in ( [(212, self.aux_plot_y, self.aux_plot_ylabel(), self.method_title()), (211, lambda y: y, 'PLM Temperature [mK]' + plm_sign, caption)] if self.need_aux_plot() else [(111, lambda y: y, 'PLM Temperature [mK]' + plm_sign, caption)]): pylab.subplot(sub) pylab.title(title) pylab.xlabel('Date and Time (GMT)') pylab.ylabel(ylabel) T, dT = self.temperature(tstart, tend) ymin = yfunc(T - dT) ymax = yfunc(T + dT) deco.plotdate(matplotlib.dates.epoch2num(ttt), yfunc(self.temperature(ttt)[0]), 'b-') if len(tstart) > 0: pylab.errorbar(0.5 * (xmin + xmax), 0.5 * (ymin + ymax), xerr=0.5 * (xmax - xmin), yerr=0.5 * (ymax - ymin), ls='none', color='b') deco.plotdate(matplotlib.dates.epoch2num(self.t), yfunc(self.T), 'r.') return fig class NaturalWarmupInterp(TInterp): """ natural warmup: 1/T_NS = 1/(T^2 - tau)^0.5 = At + B. tau = T^2 - T_NS^2 represents the thermal gradient between NS and SP, defaults to 0.07, can be specified as a 'tau' parameter of the constructor. """ def __init__(self, t, T, **vargs): if 'tau' in vargs: self.tau = vargs['tau'] else: self.tau = 0.07 self.t = t self.T = T self.t0 = min(t) self.Tfit = aver.LinearFit.fit(t - self.t0, 1 / (T**2 - self.tau)**0.5) def temperature(self, tstart, tend=None): if tend is None: tend = tstart Xs = self.Tfit.y_mean(tstart - self.t0) dXs = abs(self.Tfit.y_err(tstart - self.t0)) Xe = self.Tfit.y_mean(tend - self.t0) dXe = abs(self.Tfit.y_err(tend - self.t0)) Tmin = (1./numpy.maximum(Xs + dXs, Xe + dXe)**2 + self.tau)**0.5 Tmax = (1./numpy.minimum(Xs - dXs, Xe - dXe)**2 + self.tau)**0.5 T = 0.5 * (Tmax + Tmin) dT = 0.5 * abs(Tmax - Tmin) return T, dT def ramp_rate(self, T): return -self.Tfit.slope().mean / T * (T**2 - self.tau)**1.5 @staticmethod def method_code(): return 1 def method_title(self): return r"Natural Warmup: $1/T_{\mathrm{NS}} = 1/\sqrt{T^2 - %.3f\,\mathrm{mK}^2} = At + B$" % self.tau def need_aux_plot(self): return True def aux_plot_ylabel(self): return r'$1/T_{\mathrm{NS}}$ [1/mK]' def aux_plot_y(self, T): return 1 / (T**2 - self.tau)**0.5 class LinearHeatRampInterp(TInterp): """ time dependence of temperature interpolation: T^2 = At + B to be used on forced temperature ramps with linear ramp of the heat. """ def __init__(self, t, T, **vargs): self.t = t self.T = T self.t0 = min(t) self.Tfit = aver.LinearFit.fit(t - self.t0, T**2) def temperature(self, tstart, tend=None): if tend is None: tend = tstart T2s = self.Tfit.y_mean(tstart - self.t0) dT2s = abs(self.Tfit.y_err(tstart - self.t0)) T2e = self.Tfit.y_mean(tend - self.t0) dT2e = abs(self.Tfit.y_err(tend - self.t0)) Tmin = numpy.minimum(T2s - dT2s, T2e - dT2e)**0.5 Tmax = numpy.maximum(T2s + dT2s, T2e + dT2e)**0.5 T = 0.5 * (Tmax + Tmin) dT = 0.5 * abs(Tmax - Tmin) return T, dT def ramp_rate(self, T): return self.Tfit.slope().mean / (2*T) @staticmethod def method_code(): return 2 def method_title(self): return "Modelling a Linear S/P Power Ramp: $T^2 = At + B$" def need_aux_plot(self): return True def aux_plot_ylabel(self): return r'$T^2$ [mK$^2$]' def aux_plot_y(self, T): return T**2 class PolyHeatRampInterp(TInterp): """ time dependence of temperature interpolation: T^2 = A_0 + A_1*t + A_2*t^2 + ... + A_N*t^N to be used on forced temperature ramps with linear ramp of the heat. """ def __init__(self, t, T, **vargs): self.t = t self.T = T self.t0 = min(t) self.power = vargs['power'] if 'power' in vargs else 1 self.Tfit = numpy.polyfit(t - self.t0, T**2, self.power) def temperature(self, tstart, tend=None): if tend is None: tend = tstart T2s = numpy.polyval(self.Tfit, tstart - self.t0) dT2s = 0.000 * numpy.ones_like(tstart) T2e = numpy.polyval(self.Tfit, tend - self.t0) dT2e = 0.000 * numpy.ones_like(tend) Tmin = numpy.minimum(T2s - dT2s, T2e - dT2e)**0.5 Tmax = numpy.maximum(T2s + dT2s, T2e + dT2e)**0.5 T = 0.5 * (Tmax + Tmin) dT = 0.5 * abs(Tmax - Tmin) return T, dT def ramp_rate(self, T): TT = numpy.polyval(self.Tfit, self.t - self.t0) ii = TT.argsort() dt = numpy.interp(T, TT[ii], self.t[ii] - self.t0) return numpy.polyval(numpy.polyder(self.Tfit), dt) / (2*T) # def ramp_rate(self, T): # return self.Tfit.slope().mean / (2*T) @staticmethod def method_code(): return 2 def method_title(self): return "Modelling a Polynomial S/P Power Ramp: $T^2 = A_0 + A_1 t + A_2 t^2 + \dots + A_N t^N, N=%d$" % self.power def need_aux_plot(self): return True def aux_plot_ylabel(self): return r'$T^2$ [mK$^2$]' def aux_plot_y(self, T): return T**2 class TcNMRTask(object): """ Finding slab and bulk peaks around the Tc """ def __init__(self, demag, pressure, fft_settings, signals, indices, backgrounds, bg_indices, Naver,\ f_slab_min=1056.5e3, f_slab_max=1058.1e3, f_bulk_min=1058.1e3, f_bulk_max=1060.0e3, logfile_suffix='', smplrate=None, logdir=None, logfile=None, slab_tip_frange=30, bulk_tip_frange=100): self.demag = demag self.pressure = pressure self.fft_settings = fft_settings self.signals = signals self.indices = indices self.backgrounds = backgrounds self.bg_indices = bg_indices self.Naver = Naver self.f_slab_min = f_slab_min self.f_slab_max = f_slab_max self.f_bulk_min = f_bulk_min self.f_bulk_max = f_bulk_max self.logfile_suffix = logfile_suffix self.logfile = logfile self.logdir = logdir self.smplrate = smplrate self.slab_tip_frange = slab_tip_frange self.bulk_tip_frange = bulk_tip_frange def outdir(self): return ("/data/exp/run20/sf/Tc/demag%g_%.2fbar" % (self.demag, self.pressure)) if self.logdir is None else self.logdir def outfile(self): # v1 ## return os.path.join(self.outdir(), 'data', 'peaks_%s_%dav_%dus_%dms_%s.dat' % \ ## (self.signals.split('bar/')[1].split('/')[0], \ ## self.Naver, self.fft_settings.truncate * 1e6, \ ## self.fft_settings.T2filter * 1e3, self.logfile_suffix)) # v2 ## return os.path.join(self.outdir(), 'data', 'peaks_%s_%dav_%dus_%dms_%s.dat' % \ ## (self.signals.split('bar/')[-1].split('/')[0], \ ## self.Naver, self.fft_settings.truncate * 1e6, \ ## self.fft_settings.T2filter * 1e3, self.logfile_suffix)) # v3 if self.logfile is not None: return os.path.join(self.outdir(), self.logfile) basename = self.signals i = basename.lower().rfind('superfluid cell/') if i >= 0: basename = basename[i + len('superfluid cell/'):] if basename.count('bar/') > 0: basename = basename.split('bar/')[-1] basename = basename.split('/')[0] return os.path.join(self.outdir(), 'data', 'peaks_%s_%dav_%dus_%dms_%s.dat' % \ (basename, self.Naver, self.fft_settings.truncate * 1e6, \ self.fft_settings.T2filter * 1e3, self.logfile_suffix)) def analyse(self, append=False): """ if 'append' is True, then only the files not covered in the output file are analysed. """ gc.collect() outfile = self.outfile() if append: # check what files have been already analysed and skip those indices = numpy.array(self.indices, dtype=int) try: N = flib.loadascii(outfile, usecols=(0,), unpack=True) N = numpy.array(N, dtype=int) indices = numpy.array(tuple(set(indices).difference(set(N)))) except: pass else: indices = self.indices if len(indices) < 1: return if not os.path.exists(outfile) or not append: if not os.path.exists(self.outdir()): os.mkdir(self.outdir()) if not os.path.exists(os.path.join(self.outdir(), 'data')): os.mkdir(os.path.join(self.outdir(), 'data')) out = open(outfile, 'w') try: out.write('# Slab and Bulk peaks\n') out.write('# Signals: %s (%s).\n' % (self.signals, flib.a2str(numpy.array(self.indices, dtype=int), maxlen=4)))# Signal directory and data range out.write('# Backgrounds: %s (%s).\n' % (self.backgrounds, flib.a2str(numpy.array(self.bg_indices, dtype=int), maxlen=4)))# Background directory and data range out.write('# Slab peak window: %.4f MHz to %.4f MHz.\n' % (self.f_slab_min/1e6, self.f_slab_max/1e6))# Slab peak window out.write('# Bulk peak window: %.4f MHz to %.4f MHz.\n' % (self.f_bulk_min/1e6, self.f_bulk_max/1e6))# Bulk peak window out.write('# Rolling average of %d signals per point.\n' % (self.Naver))# Number of averages used out.write(self.fft_settings.report_settings() + '\n') out.write('#\n') out.write('# [0] First Signal.\n') out.write('# [1] Last Signal.\n') out.write('# [2] Start time of [0] (seconds since 1st Jan 1970).\n') out.write('# [3] End time of [1] (seconds since 1st Jan 1970).\n') out.write('# [4] Average PLM temperature, XMIT4, Gain 10 (mK).\n') out.write('# [5] Half PLM temperature span (mK).\n') out.write('# [6] Naive Slab peak frequency (Hz).\n') out.write('# [7] Naive Slab peak amplitude.\n') out.write('# [8] Naive Bulk peak frequency (Hz).\n') out.write('# [9] Naive Bulk peak amplitude.\n') out.write('# Double Lorentzian fit results (if only the slab peak is fitted, bulk peak values are all NaN)\n') out.write('# [10-13] Slab peak f0,A,T2,phi\n') out.write('# [14-17] Bulk peak f0,A,T2,phi\n') out.write('# [18] 1 if slab peak lorentzian fit converged, 2 if double lorenzian fit converged, 0 if fit did not converge\n') out.write('# [19] Mean pressure (mbar)\n') out.write('# Single Lorentzian fit to slab peak tip (+-%d Hz around naive slab peak)\n' % self.slab_tip_frange) out.write('# [20-23] f0,A,T2,phi\n') out.write('# Single Lorentzian fit to bulk peak tip (+-%d Hz around naive bulk peak)\n' % self.bulk_tip_frange) out.write('# [24-27] f0,A,T2,phi\n') out.write('# [28,29] min and max pressure (mbar)\n') out.write('#\n') finally: out.close() # get start and final time for each line in the log start = numpy.nan * numpy.zeros(len(indices)) end = numpy.nan * numpy.zeros(len(indices)) for i in range(len(indices)): info=scopelog.getinfo(self.signals % indices[i]) if info is not None: start[i] = info.starttime end[i] = info.finaltime # load PLM and p31k logs for a suitable time span to cover all signals. PLM data is for XMIT4, Gain 10 plm = timelog.PLMLog.loadrun(20, start = min(start[numpy.isfinite(start)]) - 1800, end = max(end[numpy.isfinite(start)]) + 1800).at(xmit=4, gain=10) p31k = timelog.ParoLog.loadHighP(20, start = min(start[numpy.isfinite(start)]) - 1800, end = max(end[numpy.isfinite(start)]) + 1800) p31k = p31k.remove_spikes(10) # load backgrounds b = massfft.getfft(self.fft_settings, self.backgrounds, self.bg_indices, smplrate=self.smplrate) for n in indices: try: # load signals first = n last = n + self.Naver - 1 s = massfft.getfft(self.fft_settings, self.signals, pylab.frange(first, last), smplrate=self.smplrate) # Subtract backgrounds and compute FFT w = 1e6*(s-b) w = w.frange(self.f_slab_min, self.f_bulk_max) wSlab = w.frange(self.f_slab_min, self.f_slab_max) wBulk = w.frange(self.f_bulk_min, self.f_bulk_max) try: slab_fit = nmrfit.Lorentzian.fit(wSlab) wBulkNoSlab = wBulk - slab_fit.model(wBulk) if abs(wBulkNoSlab.peak()) > 0.08: bulk_fit = nmrfit.Lorentzian.fit(wBulkNoSlab) f = nmrfit.DoubleLorentzian.fit(w, slab_fit, bulk_fit) good_fit = 2 else: f = nmrfit.DoubleLorentzian(slab_fit, nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan)) good_fit = 1 except nmrfit.FitError, e: f = nmrfit.DoubleLorentzian(nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan), nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan)) good_fit = 0 try: if self.slab_tip_frange > 0: slab_tip_fit = nmrfit.Lorentzian.fit(w.frange(wSlab.peakf() - self.slab_tip_frange, wSlab.peakf() + self.slab_tip_frange)) else: slab_tip_fit = nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan) except nmrfit.FitError, e: slab_tip_fit = nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan) try: if self.bulk_tip_frange > 0: bulk_tip_fit = nmrfit.Lorentzian.fit(w.frange(wBulk.peakf() - self.bulk_tip_frange, wBulk.peakf() + self.bulk_tip_frange)) else: bulk_tip_fit = nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan) except nmrfit.FitError, e: bulk_tip_fit = nmrfit.Lorentzian(numpy.nan,numpy.nan,numpy.nan,numpy.nan) # Calculate PLM reading at the start and end of the signal(s) and then compute the mean reading and error bar. T_start = plm.interpolate(s.info.starttime) T_end = plm.interpolate(s.info.finaltime) T = 0.5*(T_start + T_end) dT = 0.5*abs(T_start - T_end) #Calculate mean pressure between start and end of signals P = p31k.period(s.info.starttime, s.info.finaltime).P() # Write data to the output file # [0] First Signal. # [1] Last Signal. # [2] Start time of [0] (seconds since 1st Jan 1970). # [3] End time of [1] (seconds since 1st Jan 1970). # [4] Average PLM temperature, XMIT4, Gain 10 (mK). # [5] Half PLM temperature span (mK). # [6] Naive Slab peak frequency (Hz). # [7] Naive Slab peak amplitude. # [8] Naive Bulk peak frequency (Hz). # [9] Naive Bulk peak amplitude. # Double Lorentzian fit results (if only the slab peak is fitted, bulk peak values are all NaN # [10-13] Slab peak f0,A,T2,phi # [11-17] Bulk peak f0,A,T2,phi # [18] 1 if slab peak lorentzian fit converged, 2 if double lorenzian fit converged, 0 if fit did not converge # [19] Mean pressure (mbar) # Single Lorentzian fit to slab peak tip (+-%d Hz around naive slab peak) # [20-23] f0,A,T2,phi # Single Lorentzian fit to bulk peak tip (+-%d Hz around naive bulk peak) # [24-27] f0,A,T2,phi # [28,29] min and max pressure (mbar) out = open(outfile, 'a') try: out.write('\t'.join([ '%d' % first, '%d' % last, '%.1f' % s.info.starttime, '%.1f' % s.info.finaltime, '%.4f' % T, '%.4f' % dT, '%.1f' % wSlab.peakf(), '%.4e' % abs(wSlab.peak()), '%.1f' % wBulk.peakf(), '%.4e' % abs(wBulk.peak()), '%.1f\t%.4e\t%.4e\t%.2f' % (f.a.f0, f.a.A, f.a.T2, f.a.phase_deg()), '%.1f\t%.4e\t%.4e\t%.2f' % (f.b.f0, f.b.A, f.b.T2, f.b.phase_deg()), '%d' % good_fit, '%.1f' % numpy.mean(P), '%.1f\t%.4e\t%.4e\t%.2f' % (slab_tip_fit.f0, slab_tip_fit.A, slab_tip_fit.T2, slab_tip_fit.phase_deg()), '%.1f\t%.4e\t%.4e\t%.2f' % (bulk_tip_fit.f0, bulk_tip_fit.A, bulk_tip_fit.T2, bulk_tip_fit.phase_deg()), '%.1f\t%.1f' % (numpy.min(P), numpy.max(P)), ]) + '\n') finally: out.close() except Exception, e: print 'signals %d-%d: error: %s' % (first, last, e) gc.collect() def plm_filter(plmlog): """ choose to use xmit 4 gain 10 or 20 depending on which is better represented in the data. correct xmit 4 gain 20 readings """ ii4_10 = (plmlog.xmit() == 4) & (plmlog.gain() == 10) n4_10 = len(ii4_10.nonzero()[0]) ii4_20 = (plmlog.xmit() == 4) & (plmlog.gain() == 20) n4_20 = len(ii4_20.nonzero()[0]) ii8_10 = (plmlog.xmit() == 8) & (plmlog.gain() == 10) n8_10 = len(ii8_10.nonzero()[0]) T = plmlog.T() if n8_10 > 2 * (n4_10 + n4_20): T[ii8_10] *= 1.004 subset = ii8_10 elif n4_10 > 0.6*n4_20: subset = ii4_10 else: T[ii4_20] *= (1 - 0.29 * numpy.exp(-1657/plmlog.M0())[ii4_20]) subset = ii4_20 plmlog.set_T(T) return plmlog.subset(subset) class TcMeasTask(object): """ The parameters (task) for Tc evaluation after peaks in NMR spectra have been located (see 'TcNMRTask' class) """ def __init__(self, demag, name, pressure, filename, start, end, Tc_star, Ts_sf_min, Ts_sf_max, Ts_norm_min, Ts_norm_max, Tb_sf_min, Tb_sf_max, Tb_norm_min, Tb_norm_max, Tfit_method=TInterp, Tfit_vargs={}, plm_filter=plm_filter, freq_shift_relaxation=None, use_Tc_slab=True, use_Tc_bulk=True, use_slab_slope=True, use_bulk_slope=True, bulk_fit_class=0): self.demag = demag self.name = name self.pressure = pressure self.filename = filename self.start = flib.tm(start) self.end = flib.tm(end) self.Tc_star = Tc_star self.Ts_sf_min = Ts_sf_min self.Ts_sf_max = Ts_sf_max self.Ts_norm_min = Ts_norm_min self.Ts_norm_max = Ts_norm_max self.Tb_sf_min = Tb_sf_min self.Tb_sf_max = Tb_sf_max self.Tb_norm_min = Tb_norm_min self.Tb_norm_max = Tb_norm_max self.Tfit_method = Tfit_method self.Tfit_vargs = Tfit_vargs self.plm_filter = plm_filter self.freq_shift_relaxation = freq_shift_relaxation self.use_Tc_slab = use_Tc_slab self.use_Tc_bulk = use_Tc_bulk self.use_slab_slope = use_slab_slope self.use_bulk_slope = use_bulk_slope self.bulk_fit_class = bulk_fit_class def analyse(self, slab_sf_T_range, bulk_sf_T_range, use_Lorentzian_fits=False, save_results=True, plot_thermometry=True, plot_fshifts=True): Tfit, tstart, tend, T, dT, f_slab, good_slab, f_bulk, good_bulk, plm_xmit, plm_gain, minP, maxP = self.select_data2(use_Lorentzian_fits) analysis_type = 'lorentzian' if use_Lorentzian_fits else 'naive' plm_sign = '' if numpy.isfinite(plm_xmit): plm_sign += 'xmit %d' % plm_xmit else: plm_sign += 'multiple xmits' plm_sign += ', ' if numpy.isfinite(plm_gain): plm_sign += 'gain %d' % plm_gain else: plm_sign += 'multiple gains' if plot_thermometry: Tpage = Tfit.present('ULT Run 20 Demag %s. SF Transition at %.2f bar: %s' % (self.demag, self.pressure, self.name), tstart, tend, plm_sign) else: Tpage = None if save_results: dirnames = ['/data/exp/run20/sf/Tc/demag%s_%.2fbar/' % (self.demag, self.pressure), '/data/exp/run20/sf/Tc/demag%s_%.2fbar/f2shifts' % (self.demag, self.pressure), '/data/exp/run20/sf/Tc/demag%s_%.2fbar/results' % (self.demag, self.pressure)] if plot_fshifts: dirnames.append('/data/exp/run20/sf/Tc/demag%s_%.2fbar/plots' % (self.demag, self.pressure)) for dirname in dirnames: if not os.path.exists(dirname): os.mkdir(dirname) dirname = '/data/exp/run20/sf/Tc/demag%s_%.2fbar/f2shifts' % (self.demag, self.pressure) if plot_fshifts: SFpage = deco.portraitA4() pylab.subplot(211) pylab.title('ULT Run 20 Demag %s. SF Transition at %.2f bar: %s' % (self.demag, self.pressure, self.name)) pylab.xlabel('PLM Temperature [mK], ' + plm_sign) pylab.ylabel(r'$f_{\mathrm{bulk}}^2 - f_0^2$ [$10^{10}\,\mathrm{Hz}^2$] (%s)' % analysis_type) xx = numpy.array([0, 0.05]) pylab.plot(self.Tc_star.mean * (1. - xx), pylab.polyval([0.135, 0.3], self.pressure) * xx, 'r-') ii = good_bulk & (f_bulk >= numpy.mean(f_bulk[(T > self.Tb_norm_min) & (T < self.Tb_norm_max)]) - 10) Tc_bulk, f0_bulk, slope_bulk, chi2_bulk, min_T_sf_bulk, max_T_sf_bulk = self.find_Tc(T[ii], dT[ii], f_bulk[ii], \ self.Tb_sf_min, self.Tb_sf_max, self.Tb_norm_min, self.Tb_norm_max, \ 'bulk', self.Tc_star, self.Ts_sf_min, None, plot_fshifts=False) try: Tb_sf_min = max(self.Tb_sf_min, Tc_bulk.mean - bulk_sf_T_range * Tc_bulk.mean) except: Tb_sf_min = self.Tb_sf_min bulk_fit_classes = {1: "the A phase", 2: "the A phase, unclear", 0: "the B phase", -1: "unclear"} bulk_fit_class = bulk_fit_classes[self.bulk_fit_class if self.bulk_fit_class in bulk_fit_classes.keys() else -1] Tc_bulk, f0_bulk, slope_bulk, chi2_bulk, min_T_sf_bulk, max_T_sf_bulk = self.find_Tc(T[ii], dT[ii], f_bulk[ii], \ Tb_sf_min, self.Tb_sf_max, self.Tb_norm_min, self.Tb_norm_max, \ 'bulk', self.Tc_star, self.Ts_sf_min, \ os.path.join(dirname, 'bulk_%s_%s.dat' % (self.name.lower().replace(' ', '_'), analysis_type)) if save_results else None, plot_fshifts=plot_fshifts, extra_summary="The slope is attributed to: %s" % bulk_fit_class) if plot_fshifts: pylab.subplot(212) pylab.xlabel('PLM Temperature [mK], ' + plm_sign) pylab.ylabel(r'$f_{\mathrm{slab}}^2 - f_0^2$ [$10^{10}\,\mathrm{Hz}^2$] (%s)' % analysis_type) ii = good_slab & (f_slab <= numpy.mean(f_slab[(T > self.Ts_norm_min) & (T < self.Ts_norm_max)]) + 10) Tc_slab, f0_slab, slope_slab, chi2_slab, min_T_sf_slab, max_T_sf_slab = self.find_Tc(T[ii], dT[ii], f_slab[ii], self.Ts_sf_min, self.Ts_sf_max, self.Ts_norm_min, self.Ts_norm_max, 'slab', self.Tc_star, self.Ts_sf_min, None, plot_fshifts=False) try: Ts_sf_min = max(self.Ts_sf_min, Tc_slab.mean - slab_sf_T_range * Tc_slab.mean) except: Ts_sf_min = self.Ts_sf_min Tc_slab, f0_slab, slope_slab, chi2_slab, min_T_sf_slab, max_T_sf_slab = self.find_Tc(T[ii], dT[ii], f_slab[ii], Ts_sf_min, self.Ts_sf_max, self.Ts_norm_min, self.Ts_norm_max, 'slab', self.Tc_star, self.Ts_sf_min, \ os.path.join(dirname, 'slab_%s_%s.dat' % (self.name.lower().replace(' ', '_'), analysis_type)) if save_results else None, plot_fshifts=plot_fshifts) if plot_fshifts: deco.eqlims(pylab.gcf(), x=True) pylab.draw() Tramp_rate_star = Tfit.ramp_rate(self.Tc_star.mean) * 1e3*60*60 Tramp_rate_bulk = Tfit.ramp_rate(Tc_bulk.mean) * 1e3*60*60 Tramp_rate_slab = Tfit.ramp_rate(Tc_slab.mean) * 1e3*60*60 if plot_fshifts and plot_thermometry: pylab.figure(Tpage.number) xlim = pylab.xlim() ylim = pylab.ylim() text = r'\begin{flushleft}' if numpy.isfinite(Tramp_rate_star): text += r'$\dot T = %.0f\,$uK/hour @ $T_c^{*}$,\\' % Tramp_rate_star if numpy.isfinite(Tramp_rate_bulk): text += r'$\dot T = %.0f\,$uK/hour @ $T_c^{\mathrm{bulk}}$,\\' % Tramp_rate_bulk if numpy.isfinite(Tramp_rate_slab): text += r'$\dot T = %.0f\,$uK/hour @ $T_c^{\mathrm{slab}}$' % Tramp_rate_slab text += r'\end{flushleft}' pylab.text(0.5 * (xlim[0] + xlim[1]), ylim[0] + 0.97 * (ylim[1] - ylim[0]), text, va='top', ha='center') if save_results: dirname = '/data/exp/run20/sf/Tc/demag%s_%.2fbar/plots' % (self.demag, self.pressure) if plot_fshifts: pages = [SFpage] if plot_thermometry: pages.append(Tpage) deco.save_pdf_book(os.path.join(dirname, 'Tc_%s_%s.pdf' % (self.name.lower().replace(' ', '_'), analysis_type)), pages) if not self.use_Tc_bulk: Tc_bulk.mean = Tc_bulk.error = numpy.nan if not self.use_Tc_slab: Tc_slab.mean = Tc_slab.error = numpy.nan if not self.use_bulk_slope: slope_bulk.mean = slope_bulk.error = numpy.nan if not self.use_slab_slope: slope_slab.mean = slope_slab.error = numpy.nan dirname = '/data/exp/run20/sf/Tc/demag%s_%.2fbar/results' % (self.demag, self.pressure) flib.saveascii(os.path.join(dirname, 'Tc_%s_%s.dat' % (self.name.lower().replace(' ', '_'), analysis_type)), (self.pressure, self.demag, self.Tc_star.mean, self.Tc_star.error, Tc_bulk.mean, Tc_bulk.error, f0_bulk.mean, f0_bulk.error, slope_bulk.mean, slope_bulk.error, chi2_bulk, Tc_slab.mean, Tc_slab.error, f0_slab.mean, f0_slab.error, slope_slab.mean, slope_slab.error, chi2_slab, Tramp_rate_star, Tramp_rate_bulk, Tramp_rate_slab, min_T_sf_bulk, max_T_sf_bulk, min_T_sf_slab, max_T_sf_slab, self.bulk_fit_class, plm_xmit, plm_gain, minP, maxP), header="""# Extracted Tc Parameters: # [0] pressure [bar] # [1] demag # Tc, f0 and slope are pairs: mean and error # [2,3] Tc_star (manually picked bulk Tc) [mK] # [4,5] Tc_bulk [mK] # [6,7] f0_bulk [Hz] # [8,9] df_bulk^2 / d(1 - T/Tc_bulk) [10^10 Hz^2], linear fit # [10] chi2_bulk - chi^2 factor for the linear fit # [11,12] Tc_slab [mK] # [13,14] f0_slab [Hz] # [15,16] df_slab^2 / d(1 - T/Tc_bulk) [10^10 Hz^2], linear fit # [17] chi2_slab - chi^2 factor for the linear fit # [18,19,20] dT/dt [uK/hour] at Tc_star, Tc_bulk and Tc_slab # [21,22] bulk fit SF temperature range # [23,24] slab fit SF temperature range # [25] bulk fit class: -1: unclear, 0: B phase, 1 A phase, 2: A phase, unclear # [26,27] PLM xmit and gain (NaN if multipe settings are used) # [28,29] min and max P31k [mbar]""") def find_Tc(self, T, dT, f, T_sf_min, T_sf_max, T_norm_min, T_norm_max, name, Tc_star, T_plot_min, outfile=None, plot_fshifts=True, extra_summary=""): """ Find a Tc based on f(T+-dT) temperature dependency of the frequency. Use: T_sf_min < T < T_sf_max range to evaluate the slope """ dT = abs(dT) ii_sf = (T_sf_min <= T - dT) & (T_sf_max >= T + dT) ii_norm = (T_norm_min <= T - dT) & (T_norm_max >= T + dT) ii_plot = (T_plot_min <= T - dT) & (T_norm_max >= T + dT) ii_discard = ii_plot & ~ii_sf & ~ii_norm if not numpy.any(ii_norm): if plot_fshifts: pylab.plot([], [], 'w-', label=r"No data in the normal state, can't proceed") pylab.legend(loc='center') return aver.Mean(numpy.NaN, numpy.NaN), aver.Mean(numpy.NaN, numpy.NaN), aver.Mean(numpy.NaN, numpy.NaN), numpy.NaN, numpy.NaN, numpy.NaN # evaluate Larmour frequency f0 = aver.average(f[ii_norm]).be() f0_2_err = 2 * f0.std * f0.mean / 1e10 df2 = (f**2 - f0.mean**2) / 1e10 if plot_fshifts: # plot f^2 - f0^2 vs. T if numpy.any(ii_sf): pylab.errorbar(T[ii_sf], df2[ii_sf], xerr=dT[ii_sf], color='b', ls='None', zorder=-50) if numpy.any(ii_norm): pylab.errorbar(T[ii_norm], df2[ii_norm], xerr=dT[ii_norm], color='g', ls='None', zorder=-50) if numpy.any(ii_discard): pylab.errorbar(T[ii_discard], df2[ii_discard], xerr=dT[ii_discard], color='#777777', ls='None', zorder=-50) if not numpy.any(ii_sf): if plot_fshifts: pylab.plot([], [], 'w-', label=r"No data in the superfluid state, can't proceed") pylab.legend(loc='center') y = numpy.diff(pylab.ylim())[0] * 0.08 if numpy.isfinite(Tc_star.mean): pylab.fill([Tc_star.mean - Tc_star.error, Tc_star.mean + Tc_star.error, Tc_star.mean + Tc_star.error, Tc_star.mean - Tc_star.error],\ [-y,-y,y,y], fc='y', ec='y', alpha=0.4, zorder=-60) return aver.Mean(numpy.NaN, numpy.NaN), f0, aver.Mean(numpy.NaN, numpy.NaN), numpy.NaN, numpy.NaN, numpy.NaN # fit linearly T vs df2 sf_fit = aver.LinearFit.fit(df2[ii_sf], T[ii_sf], dT[ii_sf]) Tc_n = sf_fit.y(-f0_2_err) Tc_p = sf_fit.y(f0_2_err) TcTc = [Tc_n.mean - Tc_n.error, Tc_n.mean + Tc_n.error, Tc_p.mean + Tc_p.error, Tc_p.mean - Tc_p.error] Tc = aver.average(min(TcTc), max(TcTc)) slope = -Tc / sf_fit.slope() if plot_fshifts: f2_min = min(df2[ii_sf | ii_norm]) f2_max = max(df2[ii_sf | ii_norm]) ff2 = numpy.linspace(f2_min, f2_max, 100) TT = sf_fit.y(ff2) contour_x = numpy.concatenate([[t.mean - t.std for t in TT], [t.mean + t.std for t in TT][-1::-1]]) contour_y = numpy.concatenate([ff2, ff2[-1::-1]]) pylab.fill(contour_x, contour_y, fc='b', ec='b', alpha=0.3, zorder=-60) f0_2e = [-f0_2_err, -f0_2_err, f0_2_err, f0_2_err] pylab.fill([T_sf_max, T_norm_max, T_norm_max, T_sf_max], f0_2e, fc='g', ec='g', alpha=0.3, zorder=-60) # pylab.fill([Tc.mean - Tc.error, Tc.mean + Tc.error, Tc.mean + Tc.error, Tc.mean - Tc.error], f0_2e, fc='r', ec='r') if numpy.isfinite(Tc_star.mean): y = numpy.diff(pylab.ylim())[0] * 0.08 pylab.fill([Tc_star.mean - Tc_star.error, Tc_star.mean + Tc_star.error, Tc_star.mean + Tc_star.error, Tc_star.mean - Tc_star.error],\ [-y,-y,y,y], fc='y', ec='y', alpha=0.4, zorder=-60) if numpy.isfinite(Tc.mean): y = numpy.diff(pylab.ylim())[0] * 0.05 pylab.fill([Tc.mean - Tc.error, Tc.mean + Tc.error, Tc.mean + Tc.error, Tc.mean - Tc.error], \ [-y,-y,y,y], fc='m', ec='m', alpha=0.4, zorder=-60) # pylab.axvspan(Tc_star.mean - Tc_star.error, Tc_star.mean + Tc_star.error, ec='none', fc='yellow', alpha=0.5) # pylab.axvspan(Tc.mean - Tc.error, Tc.mean + Tc.error, ec='none', fc='red', alpha=0.5) summary = r"""$f^{\mathrm{%s}}_0 = %.1f \pm %.1f\,$Hz $T_c^{\mathrm{%s}} = %.4f \pm %.4f\,$mK $\partial (\Delta f^{\mathrm{%s}})^2/\partial (1-T/T_c^{\mathrm{%s}}) = %.3f \pm %.3f \times 10^{10}\,$Hz$^2$ $\tilde \chi^2$ = %.3f""" % (name, f0.mean, f0.std, name, Tc.mean, Tc.error, name, name, slope.mean, slope.error, sf_fit.chi2) if numpy.isfinite(Tc_star.mean): Tc_sup = Tc / Tc_star summary += r""" \rule{0pt}{2pt} $T_c^{\mathrm{%s}} / T_c^{*} = %.4f \pm %.4f$""" % (name, Tc_sup.mean, Tc_sup.error) Tc_star_Gr = Tc_star / he3.Tc(self.pressure) summary += r""" Use $T_c^{*} = %.4f \pm %.4f$ $T_c^{*} / T_c^{\mathrm{Greywall}}(%.3f\,\mathrm{bar}) = %.4f \pm %.4f$""" % \ (Tc_star.mean, Tc_star.error, self.pressure, Tc_star_Gr.mean, Tc_star_Gr.error) if len(extra_summary) > 0: summary += r""" \rule{0pt}{2pt} %s""" % extra_summary xx = pylab.xlim() yy = pylab.ylim() if slope.mean < 0: pylab.text(xx[1] - 0.02*(xx[1] - xx[0]), yy[0] + 0.02*(yy[1] - yy[0]), r'\noindent ' + summary.replace('\n', r'\\'), ha='right', va='bottom') else: pylab.text(xx[1] - 0.02*(xx[1] - xx[0]), yy[1] - 0.02*(yy[1] - yy[0]), r'\noindent ' + summary.replace('\n', r'\\'), ha='right', va='top') if outfile is not None: mode = -numpy.ones(len(T)) mode[ii_sf] = 1.0 mode[ii_norm] = 0.0 flib.saveascii(outfile, (T, dT, df2, mode), header="""# NMR Frequency shift near Tc: # [0] T [mK] # [1] T uncertainty [mK] # [2] f^2 - fL^2 [10^10 Hz^2] # [3] mode: 1 - superfluid, 0 - normal, -1 - transitional or superfluid but not used for fits""") return Tc, f0, slope, sf_fit.chi2, min(T[ii_sf] - dT[ii_sf]), max(T[ii_sf] + dT[ii_sf]) @staticmethod def load_peak_data(filename, start=None, end=None): if isinstance(filename, str): if os.path.isdir(filename): dirname = filename filenames = [os.path.join(dirname, f) for f in os.listdir(dirname) if os.path.isfile(os.path.join(dirname, f))] if len(filenames) < 1: dirname = os.path.join(filename, 'data') filenames = [os.path.join(dirname, f) for f in os.listdir(dirname) if os.path.isfile(os.path.join(dirname, f))] else: filenames = [filename] else: filenames = filename # Write data to the output file # [0] First Signal. # [1] Last Signal. # [2] Start time of [0] (seconds since 1st Jan 1970). # [3] End time of [1] (seconds since 1st Jan 1970). # [4] Average PLM temperature, XMIT4, Gain 10 (mK). # [5] Half PLM temperature span (mK). # [6] Naive Slab peak frequency (Hz). # [7] Naive Slab peak amplitude. # [8] Naive Bulk peak frequency (Hz). # [9] Naive Bulk peak amplitude. # Double Lorentzian fit results (if only the slab peak is fitted, bulk peak values are all NaN # [10-13] Slab peak f0,A,T2,phi # [11-17] Bulk peak f0,A,T2,phi # [18] 1 if slab peak lorentzian fit converged, 2 if double lorenzian fit converged, 0 if fit did not converge # [19] Mean pressure (mbar) # Single Lorentzian fit to slab peak tip (+-30 Hz around naive slab peak) # [20-23] f0,A,T2,phi # Single Lorentzian fit to bulk peak tip (+-100 Hz around naive bulk peak) # [24-27] f0,A,T2,phi try: data = numpy.concatenate([ # flib.loadascii(f, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,14,15,16,18,19), unpack=False, keep2D=True) flib.loadascii(f, usecols=(0,1,2,3,4,5,6,7,8,9,20,21,22,24,25,26,18,19), unpack=False, keep2D=True) for f in filenames]).transpose() except: data = numpy.concatenate([ flib.loadascii(f, usecols=(0,1,2,3,4,5,6,7,8,9,10,11,12,14,15,16,18,19), unpack=False, keep2D=True) for f in filenames]).transpose() Nfirst,Nlast,tstart,tend,T,dT,fSn,aSn,fBn,aBn,fS,aS,T2S,fB,aB,T2B,good,P = data data = data[:,numpy.argsort(tstart)] Nfirst,Nlast,tstart,tend,T,dT,fSn,aSn,fBn,aBn,fS,aS,T2S,fB,aB,T2B,good,P = data if start is None: start = min(tstart) if end is None: end = max(tend) data = data[:,(tstart >= flib.tm(start)) & (tend <= flib.tm(end))] gc.collect() return data @staticmethod def select_data(Tfit_method, Tfit_vargs, filename, start, end, use_Lorentzian_fits, plm_filter=plm_filter, freq_shift_relaxation=None, extended_output=False): Nfirst,Nlast,tstart,tend,T,dT,fSn,aSn,fBn,aBn,fS,aS,T2S,fB,aB,T2B,good,P = TcMeasTask.load_peak_data(filename, start, end) if len(tstart) <1: start = end = 0 ## else: ## if flib.tm(start) < min(tstart): ## start = min(tstart) ## if flib.tm(end) > max(tend): ## end = min(tend) plm = plm_filter(timelog.PLMLog.loadrun(20, start, end))#min(tstart) - 300, max(tend) + 300)) mask = numpy.ones(len(plm), dtype=bool) for t in skip_plm_points_at: mask &= abs(plm.t() - t) > 30 plm = plm.subset(mask) plm_xmit = list(set(plm.xmit())) plm_xmit = plm_xmit[0] if len(plm_xmit) == 1 else numpy.nan plm_gain = list(set(plm.gain())) plm_gain = plm_gain[0] if len(plm_gain) == 1 else numpy.nan Tfit = Tfit_method(plm.t(), plm.T(), **Tfit_vargs) T, dT = Tfit.temperature(tstart, tend) if use_Lorentzian_fits: # use Lorentzian fit frequencies f_slab = fS a_slab = aS T2_slab = T2S good_slab = (good > 0) & (fS > 1050e3) & (fS < 1060e3) & (T2S > 0.5e-3) & (T2S < 15e-3) f_bulk = fB a_bulk = aB T2_bulk = T2B good_bulk = (good == 2) & (fB > 1050e3) & (fB < 1060e3) & (T2B > 0.5e-3) & (T2B < 15e-3) else: # use naive peaks f_slab = fSn a_slab = aSn T2_slab = numpy.nan*numpy.ones(aSn.shape) good_slab = (fSn > 1050e3) & (fSn < 1060e3) f_bulk = fBn a_bulk = aBn T2_bulk = numpy.nan*numpy.ones(aBn.shape) good_bulk = (fBn > 1050e3) & (fBn < 1060e3) if freq_shift_relaxation is not None: f_shift = freq_shift_relaxation(0.5 * (tstart + tend)) f_slab -= f_shift f_bulk -= f_shift if len(P) < 1: # this is done to make 'min' and 'max' happy P = numpy.array([numpy.nan]) if extended_output: return Tfit, tstart, tend, T, dT, f_slab, a_slab, T2_slab, good_slab, f_bulk, a_bulk, T2_bulk, good_bulk, plm_xmit, plm_gain, min(P), max(P) else: return Tfit, tstart, tend, T, dT, f_slab, good_slab, f_bulk, good_bulk, plm_xmit, plm_gain, min(P), max(P) def select_data2(self, use_Lorentzian_fits): return TcMeasTask.select_data(self.Tfit_method, self.Tfit_vargs, self.filename, self.start, self.end, use_Lorentzian_fits, self.plm_filter, self.freq_shift_relaxation) @staticmethod def select_ramp_periods(filename, lead_time=10*60, shortest_steady=3*60): Nfirst,Nlast,tstart,tend,T,dT,fSn,aSn,fBn,aBn,fS,aS,T2S,fB,aB,T2B,good,P = TcMeasTask.load_peak_data(filename) p = timelog.PIDLog.loadTemperaturePID(20, min(tstart) - 300, max(tstart) + 300) if len(p) < 1 or max(p.power_total()) == 0.0: periods = [flib.Period(min(tstart), max(tstart), direction=0.0)] else: periods = p.ramp_periods(shortest_steady) for p in periods: p.start += lead_time periods = [p for p in periods if p.end - p.start > 0 and any((tstart >= p.start) & (tend <= p.end))] for p in periods: # print "'%s', '%s', '%s'" % ('warmup' if p.direction >= 0 else 'cooldown', flib.gmt2str(p.start), flib.gmt2str(p.end)) print "flib.Period('%s', '%s', direction=%d)," % (flib.gmt2str(p.start), flib.gmt2str(p.end), p.direction) return periods @staticmethod def select_T_ranges(demag, pressure, filename, periods=None, lead_time=None, shortest_steady=3*60, \ Tfit_methods={-1.0: LinearHeatRampInterp, 1.0: LinearHeatRampInterp, 0.0: NaturalWarmupInterp}, Tfit_vargs={}, use_Lorentzian_fits=False, plm_filter=plm_filter, freq_shift_relaxation=None): if periods is None: vargs = {} if lead_time is not None: vargs['lead_time'] = lead_time if shortest_steady is not None: vargs['shortest_steady'] = shortest_steady periods = TcMeasTask.select_ramp_periods(filename, **vargs) warmup = 1 cooldown = 1 Tpage = SFpage = None tasks = [] for p in periods: if p.direction >= 0: name = 'warmup %d' % warmup warmup += 1 else: name = 'cooldown %d' % cooldown cooldown += 1 try: Tfit, tstart, tend, T, dT, f_slab, good_slab, f_bulk, good_bulk, plm_xmit, plm_gain, minP, maxP = TcMeasTask.select_data(\ Tfit_methods[p.direction], Tfit_vargs, filename, p.start, p.end, use_Lorentzian_fits, plm_filter, freq_shift_relaxation) plm_sign = '' if numpy.isfinite(plm_xmit): plm_sign += 'xmit %d' % plm_xmit else: plm_sign += 'multiple xmits' plm_sign += ', ' if numpy.isfinite(plm_gain): plm_sign += 'gain %d' % plm_gain else: plm_sign += 'multiple gains' while True: print 'Picking temperature ranges for "%s". You can move on to the next measurement at any time by pressing "r" and [Enter]' % name Tpage = Tfit.present('%s: Temperature approximation' % name, tstart, tend, plm_sign) SFpage = deco.portraitA4() pylab.subplot(211) pylab.xlabel('PLM Temperature [mK], ' + plm_sign) pylab.ylabel(r'$f_{\mathrm{bulk}}$ [Hz] - 1057.9 kHz') pylab.errorbar(T[good_bulk], f_bulk[good_bulk] - 1057.9e3, xerr=dT[good_bulk], color='b', ls='None') b_xlim = pylab.xlim() b_ylim = pylab.ylim() pylab.subplot(212) pylab.xlabel('PLM Temperature [mK], ' + plm_sign) pylab.ylabel(r'$f_{\mathrm{slab}}$ [Hz] - 1057.9 kHz') pylab.errorbar(T[good_slab], f_slab[good_slab] - 1057.9e3, xerr=dT[good_slab], color='g', ls='None') s_xlim = pylab.xlim() s_ylim = pylab.ylim() pylab.subplot(211) pylab.title('%s: Please choose Tc* on the top plot' % name) SFpage.canvas.toolbar.zoom(True) if raw_input('please select Tc* on the top plot and press [Enter]').lower() == 'r': break pylab.subplot(211) Tc_star = aver.average(pylab.xlim()) print 'T* = %.4f +- %.4f' % (Tc_star.mean, Tc_star.error) pylab.subplot(211) pylab.axvspan(Tc_star.mean - Tc_star.error, Tc_star.mean + Tc_star.error, alpha=0.5, zorder=-99, fc='y', ec='y') pylab.xlim(b_xlim) pylab.ylim(b_ylim) pylab.subplot(212) pylab.axvspan(Tc_star.mean - Tc_star.error, Tc_star.mean + Tc_star.error, alpha=0.5, zorder=-99, fc='y', ec='y') pylab.xlim(s_xlim) pylab.ylim(s_ylim) pylab.subplot(211) pylab.title('%s: Please choose superfluid T ranges on both plots' % name) if raw_input('please choose superfluid T ranges on both plots and press [Enter]').lower() == 'r': break pylab.subplot(212) Ts_sf_min, Ts_sf_max = pylab.xlim() pylab.subplot(211) Tb_sf_min, Tb_sf_max = pylab.xlim() print 'Slab superfluid at %.4f-%.4f mK, Bulk superfluid at %.4f-%.4f mK' % (Ts_sf_min, Ts_sf_max, Tb_sf_min, Tb_sf_max) pylab.subplot(211) pylab.axvspan(Tb_sf_min, Tb_sf_max, alpha=0.5, zorder=-100, fc='b', ec='b') pylab.xlim(b_xlim) pylab.ylim(b_ylim) pylab.subplot(212) pylab.axvspan(Ts_sf_min, Ts_sf_max, alpha=0.5, zorder=-100, fc='b', ec='b') pylab.xlim(s_xlim) pylab.ylim(s_ylim) pylab.subplot(211) pylab.title('%s: Please choose normal T ranges on both plots' % name) if raw_input('please choose normal T ranges on both plots and press [Enter]').lower() == 'r': break pylab.subplot(212) Ts_norm_min, Ts_norm_max = pylab.xlim() pylab.subplot(211) Tb_norm_min, Tb_norm_max = pylab.xlim() print 'Slab normal at %.4f-%.4f mK, Bulk normal at %.4f-%.4f mK' % (Ts_norm_min, Ts_norm_max, Tb_norm_min, Tb_norm_max) pylab.subplot(211) pylab.axvspan(Tb_norm_min, Tb_norm_max, alpha=0.5, zorder=-100, fc='g', ec='g') pylab.xlim(b_xlim) pylab.ylim(b_ylim) pylab.subplot(212) pylab.axvspan(Ts_norm_min, Ts_norm_max, alpha=0.5, zorder=-100, fc='g', ec='g') pylab.xlim(s_xlim) pylab.ylim(s_ylim) pylab.subplot(211) pylab.title('%s: Do you like this?' % name) print 'Measurement name: "%s"' % name print '' print 'Select one of the following option and press [Enter]' print ' y - accept these parameters and continue' print ' r - reject the current measurement and move on' print ' q - stop picking the parameters' print ' any other key - try selecting the parameters again' choice = raw_input('what would you like? ').lower() print 'your choice is ' + choice if choice == 'q': return if choice == 'y': t = TcMeasTask(demag, name, pressure, filename, p.start, p.end, Tc_star, Ts_sf_min, Ts_sf_max, Ts_norm_min, Ts_norm_max, Tb_sf_min, Tb_sf_max, Tb_norm_min, Tb_norm_max, Tfit_method=Tfit_methods[p.direction], Tfit_vargs=Tfit_vargs) tasks.append(t) break else: print 'try again' except Exception, e: print 'Processing "%s": %s' % (name, e) finally: if Tpage is not None: pylab.close(Tpage.number) Tpage = None if SFpage is not None: pylab.close(SFpage.number) SFpage = None print '[' for t in tasks: print "Tc_lib.TcMeasTask(demag=%g, name='%s', pressure=%g," % (t.demag, t.name, t.pressure) print " filename='%s'," % (t.filename) print " start='%s', end='%s', Tc_star=aver.Mean(%.4f, %.4f)," % (flib.gmt2str(t.start), flib.gmt2str(t.end), t.Tc_star.mean, t.Tc_star.error) print " Ts_sf_min=%.4f, Ts_sf_max=%.4f, Ts_norm_min=%.4f, Ts_norm_max=%.4f," % (t.Ts_sf_min, t.Ts_sf_max, t.Ts_norm_min, t.Ts_norm_max) print " Tb_sf_min=%.4f, Tb_sf_max=%.4f, Tb_norm_min=%.4f, Tb_norm_max=%.4f," % (t.Tb_sf_min, t.Tb_sf_max, t.Tb_norm_min, t.Tb_norm_max) print " Tfit_method=Tc_lib.%s, Tfit_vargs=%s)," % (t.Tfit_method.__name__, t.Tfit_vargs) print ']' return tasks # global functions for batch processing of sets of measurements def plot_normal_fshifts(measurements, use_Lorentzian_fits=False, f_offset=1057.9e3, fit_line=True): deco.portraitA4() pylab.subplot(211) pylab.xlabel('Date Time (UTC)') pylab.ylabel('Bulk Peak Frequency [Hz], -%.1f kHz' % (f_offset * 1e-3)) pylab.subplot(212) pylab.xlabel('Date Time (UTC)') pylab.ylabel('Slab Peak Frequency [Hz], -%.1f kHz' % (f_offset * 1e-3)) is_int = pylab.isinteractive() pylab.interactive(False) try: tt_slab = [] ff_slab = [] tt_bulk = [] ff_bulk = [] for m in measurements: Tfit, tstart, tend, T, dT, f_slab, good_slab, f_bulk, good_bulk, plm_xmit, plm_gain, minP, maxP = TcMeasTask.select_data(\ m.Tfit_method, m.Tfit_vargs, m.filename, m.start, m.end, use_Lorentzian_fits, m.plm_filter, m.freq_shift_relaxation) ii = good_bulk & (T - dT >= m.Tb_norm_min) f = f_bulk - f_offset f[~ii] = numpy.nan pylab.subplot(211) deco.plotdate(matplotlib.dates.epoch2num(0.5 * (tstart + tend)), f) tt_bulk = numpy.concatenate([tt_bulk, (0.5 * (tstart + tend))[ii]]) ff_bulk = numpy.concatenate([ff_bulk, f[ii]]) ii = good_slab & (T - dT >= m.Ts_norm_min) f = f_slab - f_offset f[~ii] = numpy.nan tt_slab = numpy.concatenate([tt_slab, (0.5 * (tstart + tend))[ii]]) ff_slab = numpy.concatenate([ff_slab, f[ii]]) pylab.subplot(212) deco.plotdate(matplotlib.dates.epoch2num(0.5 * (tstart + tend)), f) pylab.draw() pylab.interactive(is_int) if fit_line: t0 = min([min(tt_slab), min(tt_bulk)]) fit_bulk = aver.LinearFit.fit(tt_bulk - t0, ff_bulk) print 'f_bulk(t) = %.2f %+4g * (t - %.2f)' % (fit_bulk.y(0).mean, fit_bulk.slope().mean, t0) pylab.subplot(211) deco.plotdate(matplotlib.dates.epoch2num(tt_bulk), fit_bulk.y_mean(tt_bulk - t0), 'k--') fit_slab = aver.LinearFit.fit(tt_slab - t0, ff_slab) print 'f_slab(t) = %.2f %+4g * (t - %.2f)' % (fit_slab.y(0).mean, fit_slab.slope().mean, t0) pylab.subplot(212) deco.plotdate(matplotlib.dates.epoch2num(tt_slab), fit_slab.y_mean(tt_slab - t0), 'k--') finally: pylab.interactive(is_int) return tt_bulk, ff_bulk, tt_slab, ff_slab def analyse_and_plot(measurements, slab_sf_T_range=0.03, bulk_sf_T_range=0.10, save_results=True, use_Lorentzian_fits=False, \ plot_thermometry=False, close_windows=False, freq_shift_relaxation=None): for m in measurements: if freq_shift_relaxation is not None: m.freq_shift_relaxation = freq_shift_relaxation if m.Ts_sf_min > m.Ts_sf_max * slab_sf_T_range: m.Ts_sf_min = m.Ts_sf_max * slab_sf_T_range m.analyse(slab_sf_T_range=slab_sf_T_range, bulk_sf_T_range=bulk_sf_T_range, \ save_results=save_results, use_Lorentzian_fits=use_Lorentzian_fits, plot_thermometry=plot_thermometry) if close_windows: pylab.close('all')