Read up on morphing of different functions:

• in HEP: https://indico.cern.ch/event/507948/contributions/2028505/attachments/1262169/1866169/atlas-hcomb-morphwshop-intro-v1.pdf
• https://mathematica.stackexchange.com/questions/208990/morphing-between-two-functions
• https://mathematica.stackexchange.com/questions/209039/convert-symbolic-to-numeric-code-speed-up-morphing

Aside from morphing, the theory of optimal transport seems to be directly related to such problems:

• https://de.wikipedia.org/wiki/Optimaler_Transport (funny, this can be described by topology using GR lingo; there's no English article on this)
• https://en.wikipedia.org/wiki/Transportation_theory_(mathematics) see in particular:

This seems to imply that given some functions \(f(x)\), \(g(x)\) we are looking for the transport function \(T\) which maps \(T: f(x) \rightarrow g(x)\) in the language of transportation theory.

See the linked article about the Wasserstein metric:

• https://en.wikipedia.org/wiki/Wasserstein_metric in particular the section about its connection to the optimal transport problem

It describes a distance metric between two probability distributions. In that sense the distance between two distributions to be transported is directly related to the Wasserstein distance.

One of the major constraints of transportation theory is that the transportation has to preserve the integral of the transported function. Technically, this is not the case for our application to the CDL data, due to different amount of data available for each target. However, of course we normalize the CDL data and assume that the given data actually is a PDF. At that point each distribution is normalized to 1 and thus each morphed function has to be normalized to 1 as well. This is a decent check for a morph result. If a morphing technique does not satisfy this property we need to renormalize the result.

On the other hand, considering the slides by Verkerke (indico.cern.ch link above), the morphing between two functions can also be interpreted as a simple interpolation problem.

In that sense there are multiple approaches to compute an intermediate step of the CDL distributions.

1. visualize the CDL data as a 2D heatmap:
   • value of the logL variable is x
   • energy of each fluorescence line is y
2. linear interpolation at a specific energy \(E\) based on the two neighboring CDL distributions (interpolation thus only along y axis)
3. spline interpolation over all energy ranges
4. KDE along the energy axis (only along y) Extend KDE to 2D?
5. bicubic interpolation (-> problematic, because our energies/variables are not spread on a rectilinear grid (energy is not spread evenly)
6. other distance based interpolations, i.e. KD-tree? Simply perform an interpolation based on all neighboring points in a certain distance?

Of these 2 and 4 seem to be the easiest implementations. KD-tree would also be easy, provided I finish the implementation finally.

We will investigate different ideas in ./../CastData/ExternCode/TimepixAnalysis/Tools/cdlMorphing/cdlMorphing.nim.

  let E = ... # given as argument to function
  let lineEnergies = getXrayFluorescenceLines()
  let refLowT = df.filter(f{string -> bool: `Dset` == refLow})["Hist", float]
  let refHighT = df.filter(f{string -> bool: `Dset` == refHigh})["Hist", float]
  result = zeros[float](refLowT.size.int)
  let deltaE = abs(lineEnergies[idx] - lineEnergies[idx+offset])
  # walk over each bin and compute linear interpolation between
  for i in 0 ..< refLowT.size:
    result[i] = refLowT[i] * (1 - (abs(lineEnergies[idx] - E)) / deltaE) +
      refHighT[i] * (1 - (abs(lineEnergies[idx+offset] - E)) / deltaE)

Using a KDE is problematic, because our data is already pre-binned of course. This leads to a very sparse phase space, which either makes the local prediction around a known distribution good, but fails miserably in between them (small bandwidth) or gives decent predictions in between, but pretty bad reconstruction of the known distributions (larger bandwidth).

There is also a strong conflict in bandwidth selection, due to the non-linear steps in energy between the different CDL distributions. This leads to a too large / too small bandwidth at either end of the energy range.

Fig. 24, 25, 26 show the default bandwidth (Silverman's rule of thumb). In comparison Fig. 27, 28, 29 show the same plot using a much smaller custom bandwidth of 0.3 keV. The agreement is much better, but the actual prediction between the different distributions becomes much worse. Compare fig. 30 (default bandwidth) to fig. 31. The latter has regions of almost no counts, which is obviously wrong.

Note that also the fig. 30 is problematic. An effect of a bad KDE input is visible, namely that the bandwidth vs. number of datapoints is such that the center region (in energy) has higher values than the edges, due to the fact that predictions near the boundaries see no signal (this boundary effect could be corrected for by assuming a suitable boundary conditions, e.g. just extending the first/last distributions in the respective ranges. It is not clear however in what spacing such a distribution should be placed etc.

\clearpage

\clearpage

Another idea is to use a spline interpolation. This has the advantage that the existing distributions will be correctly predicted (as for the linear interpolation), but possibly yields better results between distributions (or in the case of predicting a known distributions.

Fig. 32, 33, 34 show the prediction using a spline. Same as for the linear interpolation each morphed distribution was computed by excluding that distribution from the spline definition and then predicting the energy of the respective fluorescence line.

The result looks somewhat better in certain areas than the linear interpolation, but has unphysical artifacts in other areas (negative values) while also deviating quite a bit. For that reason it seems like simpler is better in case of CDL morphing (at least if it's done bin-wise).

\clearpage

For the time being we will use the linear interpolation method and see where this leads us. Should definitely be a big improvement over the current interval based option.

For the results of applying linear interpolation based morphing to the likelihood analysis see section 28.6.5.

Thoughts on the implememntition in likelihood.nim or CDL morphing.

1. Add interpolation code from cdlMorphing.nim in private/cdl_cuts.nim.
2. Add a field to config.nim that describes the morphing technique to be used.
3. Add an enum for the possible morphing techniques, MorphingKind with fields mkNone, mkLinear
4. In calcCutValueTab we currently return a Table[string, float] mapping target/filter combinations to cut values. This needs to be modified such that we have something that hides away input -> output and yields what we need. Define an CutValueInterpolator type, which is returned instead. It will be a variant object with case kind: MorphingKind. This object will allow access to cut values based on:
   • string: a target/filter combination.
     • mkNone: access internal Table[string, float] as done currently
     • mkLinear: raise exception, since does not make sense
   • float: an energy in keV:
     • mkNone: convert energy to a target/filter combination and access internal Table
     • mkLinear: access the closest energy distribution and return its cut value
5. in filterClustersByLogL replace cutTab name and access by energy of cluster instead of converted to target/filter dataset

With these steps we should have a working interpolation routine. The code used in the cdlMorphing.nim test script needs to be added of course to provide the linearly interpolated logic (see step 0).