Tutorial for DLD sectors alignment using photon peak#
Preparation#
Import necessary libraries#
[1]:
%load_ext autoreload
%autoreload 2
from pathlib import Path
import os
import numpy as np
from sed import SedProcessor
from sed.dataset import dataset
import matplotlib.pyplot as plt
%matplotlib widget
from lmfit.models import GaussianModel
Get data paths#
If it is your beamtime, you can read the raw data and write to the processed directory. For the public data, you can not write to the processed directory.
The paths are such that if you are on Maxwell, it uses those. Otherwise, data is downloaded in the current directory from Zenodo: https://zenodo.org/records/15011781
[2]:
beamtime_dir = "/asap3/flash/gpfs/pg2/2021/data/11010004" # on Maxwell
if os.path.exists(beamtime_dir) and os.access(beamtime_dir, os.R_OK):
path = beamtime_dir + "/raw/hdf/FL1USER3"
buffer_path = beamtime_dir + "/processed/tutorial/"
else:
# data_path can be defined and used to store the data in a specific location
dataset.get("Photon_peak") # Put in Path to a storage of at least 10 GByte free space.
path = dataset.dir
buffer_path = path + "/processed/"
INFO - Not downloading Photon_peak data as it already exists at "/home/runner/work/sed/sed/docs/tutorial/datasets/Photon_peak".
Set 'use_existing' to False if you want to download to a new location.
INFO - Using existing data path for "Photon_peak": "/home/runner/work/sed/sed/docs/tutorial/datasets/Photon_peak"
INFO - Photon_peak data is already present.
Config setup#
Here, we get the path to the config file and set up the relevant directories. This can also be done directly in the config file.
[3]:
# pick the default configuration file for hextof@FLASH
config_file = Path('../src/sed/config/flash_example_config.yaml')
assert config_file.exists()
[4]:
# here we setup a dictionary that will be used to override the path configuration
# a few setting changes are needed as well to work with older data
config_override = {
"core": {
"beamtime_id": 11010004,
"paths": {
"raw": path,
"processed": buffer_path
},
},
"dataframe": {
"ubid_offset": 0,
"channels": {
"timeStamp": {
"index_key": "/zraw/TIMINGINFO/TIME1.BUNCH_FIRST_INDEX.1/dGroup/index",
"dataset_key": "/zraw/TIMINGINFO/TIME1.BUNCH_FIRST_INDEX.1/dGroup/time",
},
"pulseId": {
"index_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/index",
"dataset_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/value",
},
"dldPosX": {
"index_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/index",
"dataset_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/value",
},
"dldPosY": {
"index_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/index",
"dataset_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/value",
},
"dldTimeSteps": {
"index_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/index",
"dataset_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/value",
},
"dldAux": {
"index_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/index",
"dataset_key": "/zraw/FLASH.FEL/HEXTOF.DAQ/DLD1/dGroup/value",
},
"bam": {
"index_key": "/zraw/FLASH.SDIAG/BAM.DAQ/4DBC3.HIGH_CHARGE_ARRIVAL_TIME/dGroup/index",
"dataset_key": "/zraw/FLASH.SDIAG/BAM.DAQ/4DBC3.HIGH_CHARGE_ARRIVAL_TIME/dGroup/value",
},
"delayStage": {
"index_key": "/zraw/FLASH.SYNC/LASER.LOCK.EXP/FLASH1.MOD1.PG.OSC/FMC0.MD22.1.ENCODER_POSITION.RD/dGroup/index",
"dataset_key": "/zraw/FLASH.SYNC/LASER.LOCK.EXP/FLASH1.MOD1.PG.OSC/FMC0.MD22.1.ENCODER_POSITION.RD/dGroup/value",
},
"opticalDiode": {
"format": "per_train",
"index_key": "/uncategorised/FLASH.LASER/FLACPUPGLASER1.PULSEENERGY/PG2_incoupl/PULSEENERGY.MEAN/index",
"dataset_key": "/uncategorised/FLASH.LASER/FLACPUPGLASER1.PULSEENERGY/PG2_incoupl/PULSEENERGY.MEAN/value",
},
},
},
}
Read data#
[5]:
run_number = 40887
sp_ph_peak = SedProcessor(runs=[run_number], config=config_override, system_config=config_file, verbose=True)
sp_ph_peak.add_jitter()
INFO - Folder config loaded from: [/home/runner/work/sed/sed/docs/tutorial/sed_config.yaml]
INFO - System config loaded from: [/home/runner/work/sed/sed/docs/src/sed/config/flash_example_config.yaml]
INFO - Default config loaded from: [/opt/hostedtoolcache/Python/3.10.17/x64/lib/python3.10/site-packages/sed/config/default.yaml]
INFO - Reading files: 0 new files of 17 total.
loading complete in 0.12 s
INFO - add_jitter: Added jitter to columns ['dldPosX', 'dldPosY', 'dldTimeSteps'].
Check which channels are included in the dataframe
[6]:
sp_ph_peak.dataframe[["dldTimeSteps", "dldSectorID"]].head()
[6]:
| dldTimeSteps | dldSectorID | |
|---|---|---|
| 0 | 2390.487828 | 3 |
| 1 | 2390.865981 | 0 |
| 2 | 4275.866290 | 7 |
| 3 | 4277.303092 | 6 |
| 4 | 4421.184809 | 7 |
Data w/o correction of quadrants in time#
First, we take a look at the photon peak and apply separation by single quadrants before any corrections. We plot the data in detector time (dldTimeSteps) as well as in detector position (dldPosX and dldPosY) coordinates, with additional separation by single sectors.
[7]:
axes = ['dldSectorID', 'dldTimeSteps','dldPosX','dldPosY']
ranges = [[0,8], [2360,2460], [435,885], [445,895]]
bins = [8,700,225,225]
res_ph_peak = sp_ph_peak.compute(bins=bins, axes=axes, ranges=ranges)
[8]:
res_ph_peak['dldPosX'].attrs['unit'] = 'pixel'
res_ph_peak['dldPosY'].attrs['unit'] = 'pixel'
fig,ax = plt.subplots(1,2,figsize=(6,2.25), layout='tight')
res_ph_peak.sum(('dldSectorID','dldPosX','dldPosY')).plot(ax=ax[0])
res_ph_peak.sel(dldTimeSteps=slice(2380,2400)).mean(('dldSectorID','dldTimeSteps')).plot(ax=ax[1], robust=True)
[8]:
<matplotlib.collections.QuadMesh at 0x7f80d8f85ba0>
Just photon peak itself without surrounding background
[9]:
ph_peak = res_ph_peak.sel(dldTimeSteps=slice(2380,2400)).sum(('dldPosX','dldPosY'))
plt.figure(figsize=(6,4))
ph_peak.sum('dldSectorID').plot()
[9]:
[<matplotlib.lines.Line2D at 0x7f80e03e1330>]
Let’s check the signal (photon peak) from every single sector
[10]:
plt.figure(figsize=(6,4))
ph_peak.plot()
[10]:
<matplotlib.collections.QuadMesh at 0x7f80e01983a0>
[11]:
plt.figure(figsize=(6,4))
for i, item in enumerate(ph_peak):
item.plot(label=f'S{i}')
plt.legend()
Position of the photon peak#
[12]:
Gauss_mod = GaussianModel()
x=ph_peak['dldTimeSteps']
y=ph_peak.sum('dldSectorID')
pars = Gauss_mod.make_params(amplitude=200, center=2390, sigma=1)
# pars = Gauss_mod.guess(y, x=x)
out = Gauss_mod.fit(y, pars, x=x)
print(out.fit_report())
plt.figure(figsize=(6,4))
plt.plot(x,y, 'rx')
plt.plot(x,out.best_fit, "b", label="FWHM = {:.3f}".format(out.values['fwhm']))
plt.title(f'Run {run_number}, full photon peak')
plt.legend(loc="best")
plt.xlabel("dldTimeSteps [step]")
[[Model]]
Model(gaussian)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 45
# data points = 141
# variables = 3
chi-square = 348712.352
reduced chi-square = 2526.90110
Akaike info crit = 1107.66723
Bayesian info crit = 1116.51351
R-squared = 0.93905197
[[Variables]]
amplitude: 1939.42461 +/- 45.8617460 (2.36%) (init = 200)
center: 2388.77841 +/- 0.02992332 (0.00%) (init = 2390)
sigma: 1.09575685 +/- 0.02992198 (2.73%) (init = 1)
fwhm: 2.58031015 +/- 0.07046088 (2.73%) == '2.3548200*sigma'
height: 706.104198 +/- 16.6986338 (2.36%) == '0.3989423*amplitude/max(1e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, sigma) = +0.5773
[12]:
Text(0.5, 0, 'dldTimeSteps [step]')
Width of the photon peak from every quadrant and they offset in ps in respect to 0 quadrant
[13]:
plt.figure(figsize=(6,4))
sector_delays = np.zeros(8)
for i, item in enumerate(ph_peak):
x=ph_peak['dldTimeSteps']
y=item
pars = Gauss_mod.make_params(amplitude=200, center=2390, sigma=1)
out = Gauss_mod.fit(y, pars, x=x)
Center = 2388.984276411258
Diff = "{:.3f}".format(Center - out.values['center'])
sector_delays[i] = (out.values['center'])
FWHM = "{:.3f}".format(out.values['fwhm'])
item.plot(label=f'S{i}={Diff}, FWHM = {FWHM}')
plt.title(f'Run {run_number}, individual sectors, not aligned')
plt.legend()
[14]:
sector_delays = sector_delays - np.mean(sector_delays)
sector_delays
[14]:
array([ 0.84143035, 0.63440337, -0.2350536 , -0.70526687, 0.31571148,
-0.48552155, -0.14427478, -0.2214284 ])
sector alignment#
as usual, first, we jitter, but here we also align in time the 8 sectors of the dld. This is done by finding the time of the maximum of the signal in each sector, and then shifting the signal in each sector by the difference between the maximum time and the time of the maximum in each sector.
[15]:
sp_ph_peak.align_dld_sectors(sector_delays=sector_delays)
INFO - Aligning 8s sectors of dataframe
INFO - Dask DataFrame Structure:
trainId pulseId electronId dldPosX dldPosY dldTimeSteps bam timeStamp delayStage opticalDiode monochromatorPhotonEnergy gmdBda sampleBias tofVoltage extractorVoltage extractorCurrent cryoTemperature sampleTemperature dldTimeBinSize dldSectorID
npartitions=17
uint32 int64 int64 float64 float64 float32 float32 float64 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 int8
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dask Name: assign, 16 graph layers
[16]:
sp_ph_peak.dataframe[["dldTimeSteps", "dldSectorID"]].head()
[16]:
| dldTimeSteps | dldSectorID | |
|---|---|---|
| 0 | 2391.074463 | 3 |
| 1 | 2390.388184 | 0 |
| 2 | 4276.004883 | 7 |
| 3 | 4276.660645 | 6 |
| 4 | 4421.423340 | 7 |
Width of the photon peak after sector alignment#
Now we can repeat the fit procedure for combined and sector-separated photon peaks to see the effect of sector alignment
[17]:
axes = ['dldSectorID', 'dldTimeSteps','dldPosX','dldPosY']
ranges = [[0,8], [2360,2460], [435,885], [445,895]]
bins = [8,700,225,225]
res_ph_peak_align = sp_ph_peak.compute(bins=bins, axes=axes, ranges=ranges)
ph_peak_align = res_ph_peak_align.sel(dldTimeSteps=slice(2380,2400)).sum(('dldPosX','dldPosY'))
fig,ax = plt.subplots(1,2,figsize=(6,3.25), layout='tight')
ph_peak_align.sum('dldSectorID').plot(ax=ax[0])
for i, item in enumerate(ph_peak_align):
item.plot(ax=ax[1], label=f'S{i}')
plt.legend()
[18]:
Gauss_mod = GaussianModel()
x=ph_peak_align['dldTimeSteps']
y=ph_peak_align.sum('dldSectorID')
pars = Gauss_mod.make_params(amplitude=200, center=2390, sigma=1)
out = Gauss_mod.fit(y, pars, x=x)
print(out.fit_report())
plt.figure(figsize=(6,4))
plt.plot(x,y, 'rx')
plt.plot(x,out.best_fit, "b", label="FWHM = {:.3f}".format(out.values['fwhm']))
plt.title(f'Run {run_number}, full photon peak, sectors aligned')
plt.legend(loc="best")
plt.xlabel("dldTimeSteps [step]")
plt.show()
[[Model]]
Model(gaussian)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 49
# data points = 141
# variables = 3
chi-square = 155431.599
reduced chi-square = 1126.31593
Akaike info crit = 993.733361
Bayesian info crit = 1002.57964
R-squared = 0.97271223
[[Variables]]
amplitude: 1950.51777 +/- 30.4190303 (1.56%) (init = 200)
center: 2388.97948 +/- 0.01947815 (0.00%) (init = 2390)
sigma: 1.08151832 +/- 0.01947726 (1.80%) (init = 1)
fwhm: 2.54678098 +/- 0.04586543 (1.80%) == '2.3548200*sigma'
height: 719.492242 +/- 11.2216435 (1.56%) == '0.3989423*amplitude/max(1e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, sigma) = +0.5773
As we can see from the result of the last fit, after sector alignment, we have improved the photon peak width by 0.058 steps.
The same check can be done for every single sector to see/check that all sectors were properly corrected in time by their difference from the 0 sector.
[19]:
plt.figure(figsize=(6,4))
for i, item in enumerate(ph_peak_align):
x=ph_peak_align['dldTimeSteps']
y=item
pars = Gauss_mod.make_params(amplitude=800.0, center=2390, sigma=1)
out = Gauss_mod.fit(y, pars, x=x)
Center = 2388.984276411258
Diff = "{:.3f}".format(Center - out.values['center'])
FWHM = "{:.3f}".format(out.values['fwhm'])
item.plot(label=f'S{i}={Diff}, FWHM = {FWHM} step')
plt.title(f'Run {run_number}, individual sectors, aligned')
plt.legend()
plt.show()
Now we can make an ns conversion and do a fitting procedure again to see, e.g. DLD time resolution.
[20]:
sp_ph_peak.append_tof_ns_axis()
INFO - Adding time-of-flight column in nanoseconds to dataframe.
INFO - Dask DataFrame Structure:
trainId pulseId electronId dldPosX dldPosY dldTimeSteps bam timeStamp delayStage opticalDiode monochromatorPhotonEnergy gmdBda sampleBias tofVoltage extractorVoltage extractorCurrent cryoTemperature sampleTemperature dldTimeBinSize dldSectorID dldTime
npartitions=17
uint32 int64 int64 float64 float64 float32 float32 float64 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 int8 float64
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dask Name: assign, 22 graph layers
[21]:
sp_ph_peak.dataframe[["dldTimeSteps", "dldTime", "dldSectorID"]].head()
[21]:
| dldTimeSteps | dldTime | dldSectorID | |
|---|---|---|---|
| 0 | 2390.302979 | 393.465517 | 3 |
| 1 | 2389.917969 | 393.402141 | 0 |
| 2 | 4276.472168 | 703.946046 | 7 |
| 3 | 4277.381836 | 704.095786 | 6 |
| 4 | 4421.681152 | 727.848760 | 7 |
[22]:
axes = ['dldSectorID', 'dldTime']
ranges = [[0,8], [390,397]]
bins = [8,350]
res_ph_peak_ns = sp_ph_peak.compute(bins=bins, axes=axes, ranges=ranges)
plt.figure()
res_ph_peak_ns.sum(('dldSectorID')).plot()
plt.show()
[23]:
Gauss_mod = GaussianModel()
x=res_ph_peak_ns['dldTime']
y=res_ph_peak_ns.sum('dldSectorID')
pars = Gauss_mod.make_params(amplitude=340.0, center=393.2, sigma=0.19)
# pars = Gauss_mod.guess(y, x=x)
out = Gauss_mod.fit(y, pars, x=x)
print(out.fit_report())
plt.figure()
plt.plot(x,y, 'rx')
plt.plot(x,out.best_fit, "b", label="FWHM = {:.3f} ns".format(out.values['fwhm']))
#plt.title(f'Run {runs}, full photon peak')
plt.legend(loc="best")
plt.xlabel("dldTime [ns]")
plt.title(f'Run {run_number}, all sectors aligned')
plt.show()
[[Model]]
Model(gaussian)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 25
# data points = 350
# variables = 3
chi-square = 142167.863
reduced chi-square = 409.705656
Akaike info crit = 2108.39072
Bayesian info crit = 2119.96451
R-squared = 0.97421364
[[Variables]]
amplitude: 273.025296 +/- 2.77706424 (1.02%) (init = 340)
center: 393.246191 +/- 0.00207909 (0.00%) (init = 393.2)
sigma: 0.17699854 +/- 0.00207899 (1.17%) (init = 0.19)
fwhm: 0.41679971 +/- 0.00489566 (1.17%) == '2.3548200*sigma'
height: 615.379867 +/- 6.25982720 (1.02%) == '0.3989423*amplitude/max(1e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
C(amplitude, sigma) = +0.5773
[ ]: