Skip to content

solver_interp Module

This module defines an interpolation based ODE solving-fitting routine. The ODE for the AT-TPC system is solved ahead of time for many different initial values. These solutions are then interpolated and fit to the data to determine a best set of parameters for a given trajectory. In this case the optimization used is L-BFGS-B.

distances(track, data, weights)

Calculate the average distance (error) of a track solution to the data

Loop over the data and approximate the error as the closest distance of a point in the track solution to the data. JIT-ed for speed

Parameters:

Name Type Description Default
track ndarray

The track solution

 required
data ndarray

The data to compare to

 required
weights ndarray

The weights due to data uncertainty (1/sigma)

 required

Returns:

Type Description
float

The average distance (error) weighted by uncertainty
Source code in src/spyral/solvers/solver_interp.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
@njit(fastmath=True, error_model="numpy", inline="always")
def distances(track: np.ndarray, data: np.ndarray, weights: np.ndarray) -> float:
    """Calculate the average distance (error) of a track solution to the data

    Loop over the data and approximate the error as the closest distance of a point in the
    track solution to the data. JIT-ed for speed

    Parameters
    ----------
    track: numpy.ndarray
        The track solution
    data: numpy.ndarray
        The data to compare to
    weights: numpy.ndarray
        The weights due to data uncertainty (1/sigma)

    Returns
    -------
    float
        The average distance (error) weighted by uncertainty
    """
    assert track.shape[1] == 3
    assert data.shape[1] == 3
    assert len(data) == len(weights)

    dists = np.zeros((len(data), len(track)))
    error = 0.0
    for i in prange(len(data)):
        for j in prange(len(track)):
            dists[i, j] = np.sqrt(
                (track[j, 0] - data[i, 0]) ** 2.0
                + (track[j, 1] - data[i, 1]) ** 2.0
                + (track[j, 2] - data[i, 2]) ** 2.0
            )
        error += np.min(dists[i]) * weights[i]
    return error / np.sum(weights)

fit_model_interp(cluster, guess, ejectile, interpolator, det_params, solver_params)

Used for jupyter notebooks examining the good-ness of the model

Parameters:

Name Type Description Default
cluster Cluster

the data to be fit

 required
guess Guess

the initial values of the parameters

 required
ejectile NucleusData

the data for the particle being tracked

 required
interpolator TrackInterpolator

the interpolation scheme to be used

 required
det_params DetectorParameters

Configuration parameters for detector characteristics

 required
solver_params SolverParameters

Configuration parameters for the solver

 required

Returns:

Type Description
Parameters | None

Returns the best fit Parameters upon success, or None upon failure
Source code in src/spyral/solvers/solver_interp.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def fit_model_interp(
    cluster: Cluster,
    guess: Guess,
    ejectile: NucleusData,
    interpolator: TrackInterpolator,
    det_params: DetectorParameters,
    solver_params: SolverParameters,
) -> Parameters | None:
    """Used for jupyter notebooks examining the good-ness of the model

    Parameters
    ----------
    cluster: Cluster
        the data to be fit
    guess: Guess
        the initial values of the parameters
    ejectile: spyral_utils.nuclear.NucleusData
        the data for the particle being tracked
    interpolator: TrackInterpolator
        the interpolation scheme to be used
    det_params: DetectorParameters
        Configuration parameters for detector characteristics
    solver_params: SolverParameters
        Configuration parameters for the solver

    Returns
    -------
    lmfit.Parameters | None
        Returns the best fit Parameters upon success, or None upon failure
    """
    traj_data = cluster.data[:, :3] * 0.001
    momentum = QBRHO_2_P * (guess.brho * float(ejectile.Z))
    kinetic_energy = np.sqrt(momentum**2.0 + ejectile.mass**2.0) - ejectile.mass
    if not interpolator.check_values_in_range(kinetic_energy, guess.polar):
        return None

    # Uncertainty due to TB resolution in meters
    z_error = (
        det_params.detector_length
        / float(det_params.window_time_bucket - det_params.micromegas_time_bucket)
        * 0.001
    ) * 0.5
    # uncertainty due to pad size, treat as bounding rectangle
    # Note that doesn't matter which side is short/long as we just use
    # total error (distance)
    x_error = cluster.data[:, 4] * BIG_PAD_HEIGHT * 0.5
    y_error = cluster.data[:, 4] * BIG_PAD_HEIGHT / np.sqrt(3.0)
    # total positional variance per point
    total_var = x_error**2.0 + y_error**2.0 + z_error**2.0
    weights = 1.0 / total_var

    fit_params = create_params(guess, ejectile, det_params, solver_params)

    result: MinimizerResult = minimize(
        objective_function,
        fit_params,
        args=(traj_data, weights, interpolator, ejectile),
        method="lbfgsb",
    )
    print(fit_report(result))

    return result.params  # type: ignore

interpolate_trajectory(fit_params, interpolator, ejectile)

Use the interpolation scheme to generate a trajectory from the given fit parameters.

Parameters:

Name Type Description Default
fit_params Parameters

the set of lmfit Parameters

 required
interpolator TrackInterpolator

the interpolation scheme

 required
ejectile NucleusData

data for the particle being tracked

 required

Returns:

Type Description
ndarray | None

Returns a array of interpolated ODE trajectory data. Upon failure (typically an out of bounds for the interpolation scheme) returns None.
Source code in src/spyral/solvers/solver_interp.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def interpolate_trajectory(
    fit_params: Parameters, interpolator: TrackInterpolator, ejectile: NucleusData
) -> np.ndarray | None:
    """Use the interpolation scheme to generate a trajectory from the given fit parameters.

    Parameters
    ----------
    fit_params: lmfit.Parameters
        the set of lmfit Parameters
    interpolator: TrackInterpolator
        the interpolation scheme
    ejectile: NucleusData
        data for the particle being tracked

    Returns
    -------
    numpy.ndarray | None
        Returns a array of interpolated ODE trajectory data. Upon failure (typically an out of bounds for the interpolation scheme) returns None.
    """
    vertex_x = fit_params["vertex_x"].value
    vertex_y = fit_params["vertex_y"].value
    vertex_z = fit_params["vertex_z"].value
    momentum = QBRHO_2_P * (fit_params["brho"].value * float(ejectile.Z))
    kinetic_energy = np.sqrt(momentum**2.0 + ejectile.mass**2.0) - ejectile.mass
    polar = fit_params["polar"].value
    azimuthal = fit_params["azimuthal"].value

    return interpolator.get_trajectory(
        vertex_x,
        vertex_y,
        vertex_z,
        polar,
        azimuthal,
        kinetic_energy,
    )

objective_function(fit_params, x, weights, interpolator, ejectile)

Function to be minimized. Returns errors for data compared to estimated track.

Parameters:

Name Type Description Default
fit_params Parameters

the set of lmfit Parameters

 required
x ndarray

the data to be fit (x,y,z) coordinates in meters

 required
weights ndarray

The weights due to data uncertainty (1/sigma)

 required
interpolator TrackInterpolator

the interpolation scheme to be used

 required
ejectile NucleusData

the data for the particle being tracked

 required

Returns:

Type Description
float

the error between the estimate and the data
Source code in src/spyral/solvers/solver_interp.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
def objective_function(
    fit_params: Parameters,
    x: np.ndarray,
    weights: np.ndarray,
    interpolator: TrackInterpolator,
    ejectile: NucleusData,
) -> float:
    """Function to be minimized. Returns errors for data compared to estimated track.

    Parameters
    ----------
    fit_params: lmfit.Parameters
        the set of lmfit Parameters
    x: ndarray
        the data to be fit (x,y,z) coordinates in meters
    weights: numpy.ndarray
        The weights due to data uncertainty (1/sigma)
    interpolator: TrackInterpolator
        the interpolation scheme to be used
    ejectile: spyral_utils.nuclear.NucleusData
        the data for the particle being tracked

    Returns
    -------
    float
        the error between the estimate and the data
    """
    trajectory = interpolate_trajectory(fit_params, interpolator, ejectile)
    if trajectory is None:
        return 1.0e6
    return distances(trajectory, x, weights)

solve_physics_interp(cluster_index, orig_run, orig_event, cluster, guess, ejectile, interpolator, det_params, solver_params)

High level function to be called from the application.

Takes the Cluster and fits a trajectory to it using the initial Guess. It then writes the results to the dictionary.

Parameters:

Name Type Description Default
cluster_index int

Index of the cluster in the hdf5 scheme. Used only for debugging

 required
orig_run int

The original run number

 required
orig_event int

The original event number

 required
cluster Cluster

the data to be fit

 required
guess Guess

the initial values of the parameters

 required
ejectile NucleusData

the data for the particle being tracked

 required
interpolator TrackInterpolator

the interpolation scheme to be used

 required
det_params DetectorParameters

Configuration parameters for detector characteristics

 required
solver_params SolverParameters

Configuration parameters for the solver

 required

Returns:

Type Description
SolverResult | None

The best-fit result of the solver, or None if it failed
Source code in src/spyral/solvers/solver_interp.py 
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
def solve_physics_interp(
    cluster_index: int,
    orig_run: int,
    orig_event: int,
    cluster: Cluster,
    guess: Guess,
    ejectile: NucleusData,
    interpolator: TrackInterpolator,
    det_params: DetectorParameters,
    solver_params: SolverParameters,
) -> SolverResult | None:
    """High level function to be called from the application.

    Takes the Cluster and fits a trajectory to it using the initial Guess. It then writes the results to the dictionary.

    Parameters
    ----------
    cluster_index: int
        Index of the cluster in the hdf5 scheme. Used only for debugging
    orig_run: int
        The original run number
    orig_event: int
        The original event number
    cluster: Cluster
        the data to be fit
    guess: Guess
        the initial values of the parameters
    ejectile: spyral_utils.nuclear.NucleusData
        the data for the particle being tracked
    interpolator: TrackInterpolator
        the interpolation scheme to be used
    det_params: DetectorParameters
        Configuration parameters for detector characteristics
    solver_params: SolverParameters
        Configuration parameters for the solver

    Returns
    -------
    SolverResult | None
        The best-fit result of the solver, or None if it failed
    """
    traj_data = cluster.data[:, :3] * 0.001
    momentum = QBRHO_2_P * (guess.brho * float(ejectile.Z))
    kinetic_energy = np.sqrt(momentum**2.0 + ejectile.mass**2.0) - ejectile.mass
    if not interpolator.check_values_in_range(kinetic_energy, guess.polar):
        return

    # Uncertainty due to TB resolution in meters
    z_error = (
        det_params.detector_length
        / float(det_params.window_time_bucket - det_params.micromegas_time_bucket)
        * 0.001
    ) * 0.5
    # uncertainty due to pad size, treat as bounding rectangle
    # Note that doesn't matter which side is short/long as we just use
    # total error (distance)
    x_error = cluster.data[:, 4] * BIG_PAD_HEIGHT * 0.5
    y_error = cluster.data[:, 4] * BIG_PAD_HEIGHT / np.sqrt(3.0)
    # total positional variance per point
    total_var = x_error**2.0 + y_error**2.0 + z_error**2.0
    weights = 1.0 / total_var

    fit_params = create_params(guess, ejectile, det_params, solver_params)

    best_fit: MinimizerResult = minimize(
        objective_function,
        fit_params,
        args=(traj_data, weights, interpolator, ejectile),
        method="lbfgsb",
    )

    p = best_fit.params["brho"].value * QBRHO_2_P * ejectile.Z  # type: ignore
    ke = np.sqrt(p**2.0 + ejectile.mass**2.0) - ejectile.mass

    return SolverResult(
        event=cluster.event,
        cluster_index=cluster_index,
        cluster_label=cluster.label,
        orig_run=orig_run,
        orig_event=orig_event,
        vertex_x=best_fit.params["vertex_x"].value,  # type: ignore
        sigma_vx=1.0e6,
        vertex_y=best_fit.params["vertex_y"].value,  # type: ignore
        sigma_vy=1.0e6,
        vertex_z=best_fit.params["vertex_z"].value,  # type: ignore
        sigma_vz=1.0e6,
        brho=best_fit.params["brho"].value,  # type: ignore
        sigma_brho=1.0e6,
        ke=ke,
        sigma_ke=1.0e6,
        polar=best_fit.params["polar"].value,  # type: ignore
        sigma_polar=1.0e6,
        azimuthal=best_fit.params["azimuthal"].value,  # type: ignore
        sigma_azimuthal=1.0e6,
        redchisq=best_fit.redchi,  # type: ignore
    )