Histogram Analysis Functions

(*
Return measurement metadata in the format:
{time, chan0singles, chan1singles, PinTapdBm, PoutTapdBm, Voltage, currentmA}
*)

GetFileExtension[str_] :=
StringCases[str,RegularExpression[".*\\.(\\w+)"]-> {"$1"}][[1,1]];

LoadMetadata[filename_] := Module[{f = filename},
ext = GetFileExtension[f];
If[ext =="zip",
((Import[zip,f][[start;;-1, 1;;histStartIdx-1]])^\[Transpose]),
((Import[f][[start;;-1, 1;;histStartIdx-1]])^\[Transpose])
]
];

(*Load array of histograms*)
LoadHistograms[filename_] := Module[{f = filename},
ext = GetFileExtension[f];
If[ext =="zip",
((Import[zip,f][[start;;-1, histStartIdx;;-2]])^\[Transpose]),
((Import[f][[start;;-1, histStartIdx;;-2]])^\[Transpose])
]
];

(*Add all the histograms together*)
SumHistograms[histograms_] := Module[{hists = histograms},
nHists = Dimensions[hists][[2]];
Sum[hists[[All, i]], {i, 1, nHists}]
]

(*Get full measurement histogram*)
LoadTotalHistogram[filename_] := Module[{fn = filename},
hist = LoadHistograms[fn];
SumHistograms[hist]
]



Gauss[x_, A_,μ_, σ_] := A Exp[-((((x-μ))^2)/(2(σ^2)))];
(*Fit a Gaussian to histogram given inital guess parameters*)
FitGauss [histogram_,times_, Ag_, μg_, σg_] :=Module[{hist = histogram,timeVector = times,  Ai= Ag, μi = μg, σi = σg},
nmf = NonlinearModelFit[(({timeVector, hist})^\[Transpose])[[1;;-2]],
 Gauss[x, A, μ, σ],
{
{A, Ai},
{μ, μi},
{σ,σi}},
x];

(*Extract Coincidences by 3σ integrating around peak*)
{Af, μf, σf}={A, μ, σ} /.nmf["BestFitParameters"] ;

(*Return fit parameters*)
{Af, μf, σf}
]

(*Sum the counts around a given peak*)
nIntσ=3;
IntegratePeak[histogram_,times_,μ_, σ_] := Module[{hist = histogram, timeVector = times, μf = μ, σf = σ,totalCounts,idxLow,idxHigh},

idxLow = Ordering[Abs[timeVector - (μf-nIntσ*σf)], 1][[1]];
idxHigh = Ordering[Abs[timeVector - (μf+nIntσ*σf)], 1][[1]];

totalCounts = Total[hist[[idxLow;;idxHigh]]];

totalCounts
]

(*Identify the peaks in the histogram*)
FindCoincidencePeaks[histogram_, times_, downsampleFactor_:20,smoothing_:"Gaussian"] := Module[{hist = histogram,data, timeVector = times,dsf=downsampleFactor, thresh,diffs,xs,is},
(*  Hacked by JWS 05/2019 *)

(* Smooth data *)
If[smoothing=="Gaussian",
data={times,GaussianFilter[hist,dsf],hist}//Transpose;];

If[smoothing=="MovingAverage",
data={times[[Ceiling[dsf/2];;-Floor[dsf/2]-1]],MovingAverage[hist,dsf],hist[[Ceiling[dsf/2];;-Floor[dsf/2]-1]]}//Transpose;];

(* Set threshold as one σ above the mean *)
thresh=Mean@data[[All,2]]+StandardDeviation@data[[All,2]];
(* Select only data above threshold *)
data=Select[data,(#[[2]]>thresh)&];


diffs=Differences@data[[All,1]];(* Find discontinuities in thresholded data *)
is=Pick[Range[Length@data-1],#>Mean@diffs&/@diffs];(* Isolate indices of discontinuity *)
is=Join[{1},is,{Length@data-1}];(* Keep the first and last points *)

data=Table[data[[is[[i]]+1;;is[[i+1]]]],{i,Length@is-1}];(* Break data into sublists containing each peak *)

(* Return the point with the largest original hist value for each sublist *)
(* 3rd element is the original hist value *)
data=Table[MaximalBy[data[[i,All,{1,3}]],#[[2]]&,1][[1]],{i,Length@data}];

Return[data];
];


(*Get Gaussian fit parameters for peaks in a histogram*)
FitAllPeaks[histogram_, times_, σ_:0.2]:= Module[{hist = histogram, timeVector = times, σg = σ},
{ts, cs} = FindCoincidencePeaks[hist, timeVector]//Transpose;

Table[
(*Initial guess parameters*)
{Ag, μg} = {cs[[i]],ts[[i]]};

(*Fit Gaussians to our peaks *)
{As,μs,σs} =FitGauss[hist, timeVector, Ag, μg, σg],

{i, 1, Length@cs}
]
]

(*Integrate the peaks given the fit parameters*)
IntegrateAllPeaks[histogram_, times_, fitParameters_] := Table[
IntegratePeak[histogram, times, fitParameters[[i, 2]], fitParameters[[i, 3]]],
{i, 1, Length@fitParameters}
]


(*Given a histogram, fit and integrate all the peaks*)
IntegrateAndFitAllPeaks[histogram_, times_, σ_:0.2] := Module[{hist = histogram, timeVector = times, σg = σ},
(*Identify the peaks*)
{ts, cs} = FindCoincidencePeaks[hist, timeVector]//Transpose;

Table[
(*Initial guess parameters*)
{Ag, μg} = {cs[[i]],ts[[i]]};
(*Fit Gaussians to our peaks *)
{As,μs,σs} =FitGauss[hist, timeVector, Ag, μg, σg];

(*Return the integrated peak*)
IntegratePeak[hist, timeVector, μs, σs],
{i, 1, Length@ts}]
]

(*Find the peak index of the coincidence peak*)
CoincidencePeakIndex[histogram_, times_, fitParameters_] := Module[{hist = histogram, timeVector = times, fits = fitParameters},
allCounts = IntegrateAllPeaks[hist, timeVector, fits];
coincidencePeakIdx = Position[allCounts, Max@allCounts][[1,1]]
]

(*Return the integrated coincidence peak and the list of integrated accidentals*)
GetCoincidencesAccidentals[histogram_, times_, fitParameters_, coincidencePeakIdx_] := Module[{hist = histogram, timeVector = times, fits = fitParameters},
allCounts = IntegrateAllPeaks[hist, timeVector, fits];
coincidences = allCounts[[coincidencePeakIdx]];
accidentals = Drop[allCounts, {coincidencePeakIdx}];

{coincidences, Mean@accidentals}
]

(*Calculate net coincidences from histogram*)
NetCoincidences[histogram_, times_, fitParameters_, coincidencePeakIdx_]:= Module[{hist = histogram, timeVector = times, fits = fitParameters, coinPeakIdx = coincidencePeakIdx},
{gross, accs} =  GetCoincidencesAccidentals[hist, timeVector, fits, coinPeakIdx];
gross - accs
]