A/E ctc correction

[fhagemann]

So, this is @verenaaur's and my approach to the charge trapping correction of the $(A/E) \text{classifier}$:

The idea is to determine a charge trapping factor $fct$ such that a linear correction of the $(A/E) \text{classifier}$
with respect to $Q\text{drift}/E$ results in corrected, $Q\text{drift}/E$-independent $(A/E) \text{classifier}$-values.

$$(A/E)_\text{classifier}^\text{corrected} = (A/E)^\text{raw}_\text{classifier} + fct * Q\text{drift} / E$$

The steps of the algorithm are as follows:

1. Create a mask to select all events that fill into the compton_bands provided to the function.
   
2. Create a histogram over all $(A/E)\text{classifier}$ in these compton_bands --> h_before
   
3. Apply different values of $fct$ and find the one resulting in the narrowest peak --> f_optimize_ctc
   
4. Apply the charge trapping correction for the optimal $fct$ and create a histogram --> h_norm
   
5. Fit h_norm and normalize the resulting peak again to μ = 0 and σ = 1 --> h_after
   
6. Return the value of $fct$ and a string to get to the charge-trapping-corrected, normalized $(A/E)\text{classifier}$ --> aoe_ctc_func.

Some things that would need some manual testing

There are some hard-coded values and some pseudo priors need to be improved for the fit to converge towards something reasonable. I will mark the hard coded values right after opening the PR.
These include histogram ranges, boundaries for $fct$ in the fit, etc.

How to use it

# determine ecal, aoe and qdrift/e
e_expression = :e_cusp
ecal_all = getproperty(hit_cal, :e_cusp_ctc_cal)
aoe_all  = ljl_propfunc(result_correction.func).(hit_cal)
qdrift_e_all = hit_cal.qdrift ./ getproperty(hit_cal, e_expression);

# for the processor
result, report = LegendSpecFits.ctc_aoe(aoe_all, ecal_all, qdrift_e_all, compton_bands,
    aoe_expression = result_correction.func, e_expression = e_expression,
    pseudo_prior = NamedTupleDist(B = LogUniform(0.01,10000000), B2 = LogUniform(0.01,1000000)),
    pseudo_prior_all = NamedTupleDist(B = LogUniform(5,1000000000), B2 = LogUniform(0.01,1000000000)));

[codecov]

Codecov Report

Attention: Patch coverage is 0% with 58 lines in your changes missing coverage. Please review.

Project coverage is 21.41%. Comparing base (fe962f7) to head (0a6dce2).

Files with missing lines	Patch %	Lines
src/aoe_ctc.jl	0.00%	58 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##              dev      #87      +/-   ##
==========================================
- Coverage   21.89%   21.41%   -0.49%     
==========================================
  Files          34       35       +1     
  Lines        2973     3031      +58     
==========================================
- Hits          651      649       -2     
- Misses       2322     2382      +60     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[fhagemann]

I had to adjust the pseudo priors on B and B2 in order to obtain a fit that looks like the distribution. The choice of prior on B and B2 seems to influence the fit result:

No prior on B2

B2 = LogUniform(0.01,100000)

B2 = LogUniform(0.01,1000000)

[fhagemann]

I put all default values from @verenaaur as default values of the arguments.

[fhagemann]

bin_width of the histograms need some engineering

[fhagemann]

hard-coded 0.001 to determine the maximum. Assuming a normalization to μ = 0 and σ = 1, this should be sufficiently fine.

[fhagemann]

The boundary values for fitting $fct$ might need some engineering

[fhagemann]

Here, the histogram edges are hard-coded to -50 and 30. I remember agreeing on something like -20 to 3, but the charge trapping correction introduces a bias towards higher values (that is afterwards normalized), so 3 might be a bit too little. Maybe let this depend on peakstats?

[fhagemann]

Again, -50 and 30 as histogram edges

[fhagemann]

Plotting code for now (very much inspired by the energy_ctc calibration in #67):

# This should become a plot recipe at some point
let aoe_final = ljl_propfunc(result.func).(hit_cal), _aoe = ljl_propfunc(result_correction.func).(hit_cal), _qdrift_e = hit_cal.qdrift ./ hit_cal.e_cusp
    sel = abs.(aoe_final) .< 100 #.&& mask
    plot(fit(Histogram, _aoe[sel], -9:0.1:9), fill = true, xlims = (-9,5), color = :darkgrey, subplot = 1, link = :x, framestyle = :semi, size = (1000,1000), margins = (0,:mm), layout = (2,1), grid = false, st = :stepbins, left_margin = (5,:mm), right_margin = (5,:mm), bottom_margin = (-4,:mm), label = "Before correction")
    plot!(fit(Histogram, aoe_final[sel], -9:0.1:9), fill = true, xlims = (-9,5), alpha = 0.5, color = :purple, subplot = 1, link = :x, framestyle = :semi, size = (1000,1000), margins = (0,:mm), layout = (2,1), grid = false, st = :stepbins, left_margin = (5,:mm), right_margin = (5,:mm), bottom_margin = (-4,:mm), label = "After correction", legend = :topleft, ylabel = "counts / 0.1")
    plot!(kde((_aoe[sel], (_qdrift_e)[sel])), subplot = 2, c = :binary, colorbar = :none, st = :line, fill = true, label = "After correction", yformatter = :plain, link = :x)
    plot!(kde((aoe_final[sel], (_qdrift_e)[sel])), subplot = 2, c = :plasma, link = :x, framestyle = :semi, colorbar = :none, st = :line, fill = false, label = "After correction", yformatter = :plain, xlims = (-9,5), ylims = (0,11), ylabel = "Eff. Drift time / Energy (a.u.)")
    plot!(xlabel = "A/E classifier", xtickfontsize = 12, xlabelfontsize = 14, ylabelfontsize = 14, ytickfontsize = 12, legendfontsize = 12, foreground_color_legend = :silver, background_color_legend = :white, fmt = :png)
end