% Hbb_ReadCal reads in calibration files collected with the Hyper-bb and
% computes specific calibration factors used in Hyper-bb data processing
% 
% Usage:
%   cal = Hbb_ReadCal  
%   OR
%   cal = Hbb_ReadCal(PMT_Cal_FilePath, Plaque_Cal_FilePath, Dark_Cal_FilePath, PlotResults)
%
% Inputs:
%   Hbb_ReadCal with no arguments will open up a file browser that will
%   prompt for selection of the PMT cal file, plaque cal file, and dark file. 
%
%   Otherwise, the file paths to the cal files are specified as input
%   parameters to the function. Plotting of the results can be turned on/off
%   by passing 1 or 0 as the fourth parameter.
%
% Output:
%   Hbb_ReadCal outputs a structure containing the following data. The cal 
%   structure is automatically saved in a .mat file named with the cal
%   timestamp.
%       pmtGamma:           Exponent of power function fit between pmt gain vs signal
%       pmtRefGain:         Reference gain for power function fit
%    	darkCalWavelength:  Wavelength values from dark cal file
%       darkCalPmtGain:     PMT gain values from dark cal file
%    	darkCalScat1:       Dark offset values in scat channel 1 (wl x gain)
%    	darkCalScat2:       Dark offset values in scat channel 2 (wl x gain)
%    	darkCalScat3:       Dark offset values in scat channel 3 (wl x gain)
%       gain12:             Median gain ratio between raw and low gain
%       gain23:             Median gain ratio between low and high gain
%       H:                  Distance between centers of transmit and receive beams (mm)
%       rho:                Used in calculation of muFactors
%       muWavelengths:      Wavelengths (nm)
%       muFactors:          Instrument response factors
%       muLedTemp:          LED temperature during mu calibration
%       timestamp:          Timestamp of calibration
%
% Sequoia Scientific, Inc. - 10/23/2020

function [cal, figuresHandles] =  Hbb_ReadCal(varargin)
%% Parse function parameters

cal = [];
plotResults = 1;

if isempty(varargin) % no inputs, open GUI browser to select files
    
    [file,path] = uigetfile('*.txt','Select PMT Cal File');
    PMT_Cal_Filename = fullfile(path,file);
    
    [file,path] = uigetfile(fullfile(path,'*.txt'),'Select Plaque Cal File');
    Plaque_Cal_Filename = fullfile(path,file);
    
    [file,path] = uigetfile(fullfile(path,'*.txt'),'Select Dark Cal File');
    Dark_Cal_Filename = fullfile(path,file);
    
elseif length(varargin) == 4 % filenames provided as function parameters
    
    PMT_Cal_Filename = varargin{1};
    Plaque_Cal_Filename = varargin{2};
    Dark_Cal_Filename = varargin{3};
    plotResults = varargin{4};
else
    
    error('Hbb_ReadCal function requires no parameters or four parameters: PMT Cal Filename, Plaque Cal Filename, Dark Cal Filename, Plotting (0|1)');
end

%% PMT CAL - Power-law scaling of PMT output vs. PmtGain from fixed plaque, stepping PmtGain

gDat = [];
try
    gDat = Hbb_ReadFile(PMT_Cal_Filename, 'IncludeRawData',true);
end

if isempty(gDat)
    warning('Unable to read PMT Cal File');
else
    
    % Use wavelength closest to 470nm null
    [~, ix] = min(abs(gDat.wl-470));
    nullWl = gDat.wl(ix);
    
    % remove other wavelengths, low signals, saturated signals
    gDat(gDat.wl~=nullWl,:) = [];
    gDat(gDat.NetSig1 < 5, :) = [];
    gDat(gDat.SigOn1 > 3800, :) = [];
    gDat(isnan(gDat.SigOn1),:) = [];
    
    % scale and fit power function
    pmtRefGain = 2000;
    scaleFact = 1./(gDat.Scat1./nanmedian(gDat.Scat1(gDat.PmtGain==pmtRefGain)));
    
    [scaleSigRef, pmtGamma] = FitPmtGainCurve(gDat.PmtGain, scaleFact, pmtRefGain);
    scaleHat = scaleSigRef.*(gDat.PmtGain./pmtRefGain).^pmtGamma;
    
    % save data in cal struct
    cal.pmtGamma = pmtGamma;
    cal.pmtRefGain = pmtRefGain;
    
    if plotResults
        
        figuresHandles(1) = figure;
        subplot(121)
        hold on
        plot(gDat.PmtGain,gDat.Scat1, 'k.', 'markersize',8)
        set(gca, 'YScale','log')
        xlabel('PMT Gain Control Setting'), ylabel(['PMT Signal at ' num2str(nullWl) 'nm'])
        box on
        grid on
        
        subplot(122), hold on,
        plot(gDat.PmtGain,scaleFact, 'k.','markersize',8)
        plot(gDat.PmtGain, scaleHat,'b-')
        set(gca, 'YScale','log')
        xlabel('PMT Gain Control Setting'), ylabel(['PMT Gain Factor, G_{PMT} at ' num2str(nullWl) 'nm'])
        legend('Measured Signal',['pmtGain^{' num2str(pmtGamma) '}'],'fontsize',10)
        box on
        grid on
    end
    
end

%% Instrument dark calibration

darkFile = [];
try
    darkFile = Hbb_ReadFile(Dark_Cal_Filename, 'IncludeRawData',true);
end
if isempty(darkFile)
    warning('Unable to read Dark Cal File');
else
    
    gains = sort(unique(darkFile.PmtGain));
    colors = jet(length(gains));
    
    if plotResults
        figuresHandles(2) = figure;
        for k = 1:numel(gains)
            subplot(311)
            plot(darkFile.wl(darkFile.PmtGain==gains(k),:),darkFile.Scat1(darkFile.PmtGain==gains(k),:),'.-', 'Color',colors(k,:))
            hold on
            title('Dark Calibration Scat 1')
            xlabel('Wavelength (nm)')
            ylabel('Scat 1')
            subplot(312)
            plot(darkFile.wl(darkFile.PmtGain==gains(k),:),darkFile.Scat2(darkFile.PmtGain==gains(k),:),'.-', 'Color',colors(k,:))
            hold on
            title('Dark Calibration Scat 2')
            xlabel('Wavelength (nm)')
            ylabel('Scat 2')
            subplot(313)
            plot(darkFile.wl(darkFile.PmtGain==gains(k),:),darkFile.Scat3(darkFile.PmtGain==gains(k),:),'.-', 'Color',colors(k,:))
            hold on
            title('Dark Calibration Scat 3')
            xlabel('Wavelength (nm)')
            ylabel('Scat 3')
        end
    end
    
    % save dark data for each gain level as matrix (wl x pmtGain)
    cal.darkCalWavelength = sort(unique(darkFile.wl));
    cal.darkCalPmtGain = gains;
    cal.darkCalScat1 = reshape(darkFile.Scat1,[length(cal.darkCalWavelength),length(gains)]);
    cal.darkCalScat2 = reshape(darkFile.Scat2,[length(cal.darkCalWavelength),length(gains)]);
    cal.darkCalScat3 = reshape(darkFile.Scat3,[length(cal.darkCalWavelength),length(gains)]);
    
    
end


%% Analog front end gains, Gfe12 and Gfe23 from the z-resolved response measurements

muDat = [];
try
    muDat = Hbb_ReadFile(Plaque_Cal_Filename, 'IncludeRawData',true);
end

if isempty(muDat)
    warning('Unable to read Plaque Cal File');
else
    if isempty(cal) || ~isfield(cal,'darkCalScat1')
        warning('Instrument response vs distance could not be completed, missing dark calibration data');
    else
        
        % raw to low gain ratio
        msk12 = muDat.wl>500 & muDat.wl<650 & muDat.NetSig2>100 & muDat.NetSig1>100;
         
        % find dark offsets
        darkOffset_scat1 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat1,muDat.PmtGain(msk12),muDat.wl(msk12), 'linear');
        darkOffset_scat2 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat2,muDat.PmtGain(msk12),muDat.wl(msk12), 'linear');

        gain12 = nanmedian( (muDat.Scat2(msk12) - darkOffset_scat2) ./ (muDat.Scat1(msk12) - darkOffset_scat1) );
        
        % low to high gain ratio
        msk23 = muDat.wl>500 & muDat.wl<650 & muDat.NetSig3>100 & muDat.NetSig2>100;
        
        % find dark offsets
        darkOffset_scat2 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat2,muDat.PmtGain(msk23),muDat.wl(msk23), 'linear');
        darkOffset_scat3 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat3,muDat.PmtGain(msk23),muDat.wl(msk23), 'linear');

        gain23 = nanmedian( (muDat.Scat3(msk23) - darkOffset_scat3) ./ (muDat.Scat2(msk23) - darkOffset_scat2) );
        
        
        % save data in cal struct
        cal.gain12 = gain12;
        cal.gain23 = gain23;
        
        if plotResults
            
            figuresHandles(3) = figure;
            subplot(1,2,1), hold on,
            plot(muDat.wl(msk12),muDat.Scat2(msk12)./muDat.Scat1(msk12), 'bd')
            plot([min(muDat.wl) max(muDat.wl)], [gain12 gain12], 'k-','LineWidth',2)
            xlabel('Wavelength [nm]'), ylabel('Front-End Gain, G_{12}')
            box on
            grid on
            
            subplot(1,2,2), hold on,
            plot(muDat.wl(msk23),muDat.Scat3(msk23)./muDat.Scat2(msk23), 'ro')
            plot([min(muDat.wl) max(muDat.wl)], [gain23 gain23], 'k-','LineWidth',2)
            xlabel('Wavelength [nm]'), ylabel('Front-End Gain, G_{23}')
            box on
            grid on
            
        end
        
    end
    
    
end

%% Instrument response (mu) calibration

muDat = [];
try
    muDat = Hbb_ReadFile(Plaque_Cal_Filename, 'IncludeRawData',true);
end
if isempty(muDat)
    warning('Unable to read Plaque Cal File');
else
    
    if isempty(cal) || ~isfield(cal,'pmtGamma') || ~isfield(cal,'gain12') || ~isfield(cal,'darkCalScat1')
        warning('Instrument response vs distance could not be completed, missing PMT or gain ratio calibration data');
    else
        
        muDat(muDat.zDistance<6,:) = [];
        
        % Subtract Dark Offset
        darkOffset_scat1 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat1,muDat.PmtGain,muDat.wl, 'linear');
        darkOffset_scat2 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat2,muDat.PmtGain,muDat.wl, 'linear');
        darkOffset_scat3 = interp2(cal.darkCalPmtGain,cal.darkCalWavelength,cal.darkCalScat3,muDat.PmtGain,muDat.wl, 'linear');

        scat1_darkRemoved = muDat.Scat1 - darkOffset_scat1;
        scat2_darkRemoved = muDat.Scat2 - darkOffset_scat2;
        scat3_darkRemoved = muDat.Scat3 - darkOffset_scat3;
        
        % Apply PMT and front end gain factors
        scalePmt = (muDat.PmtGain./cal.pmtRefGain).^cal.pmtGamma;
        muDat = addvars(muDat,scat1_darkRemoved.*cal.gain12.*cal.gain23.*scalePmt, 'After','Scat3', 'NewVariableNames', 'ScatCor1');
        muDat = addvars(muDat,scat2_darkRemoved.*cal.gain23.*scalePmt,         'After','ScatCor1', 'NewVariableNames', 'ScatCor2');
        muDat = addvars(muDat,scat3_darkRemoved.*scalePmt,               'After','ScatCor2', 'NewVariableNames', 'ScatCor3');
        
        zDistance = unique(muDat.zDistance);
        muWavelengths = sort(unique(muDat.wl))';
        
        [gridZ, gridWL] = meshgrid(zDistance,muWavelengths);
        gridPmtSig = griddata(muDat.wl,muDat.zDistance,muDat.ScatCor2,gridWL,gridZ, 'nearest');
        
        H = 31.4; % distance between centers of transmit and receive beams in mm (based on CAD model)
        rho = 1.1;
        
        zDelta = (zDistance(2)-zDistance(1))*1e-3;
        muFactors = NaN(size(muWavelengths));
        muFactors2 = NaN(size(muWavelengths));
        muFactors_temp = NaN(size(muWavelengths));
        colors = jet(length(muWavelengths));
        
        if plotResults
            figuresHandles(4) = figure;
            subplot(2,2,1)
            hold on
        end
        
        for kw = 1:length(muWavelengths)
            pmtSig = muDat.ScatCor2(muDat.wl==muWavelengths(kw));
            if plotResults
                plot(zDistance,pmtSig, '.-', 'Color',colors(kw,:))
            end
            muFactors(kw) = rho/(pi*sum((pmtSig./cos(atan(H./(2.*zDistance))))*zDelta));
            muFactors2(kw) = rho/(pi*sum((pmtSig./(cos(atan(H./(2.*zDistance)))).^2)*zDelta));
            muFactors_temp(kw) = mean(muDat.LedTemp(muDat.wl==muWavelengths(kw)));
        end
        
        % save data in cal struct
        cal.H = H;
        cal.rho = rho;
        cal.muWavelengths = muWavelengths;
        cal.muFactors = muFactors;
        cal.muLedTemp = muFactors_temp;
        cal.timestamp = muDat.Timestamp(1);
        
        if plotResults
            xlabel('Distance, z [mm]'), ylabel('Instrument Response [au]')
            box on
            grid on
            
            subplot(2,2,2), hold on
            plot(muWavelengths, muFactors, 'bo-')
            %plot(muWavelengths, muFactors2, 'rs-')
            xlabel('Wavelength [nm]'), ylabel('Calibration Factor, \mu(\lambda)')
            box on
            grid on
            
            subplot(2,2,[3 4]), mesh(gridWL, gridZ, gridPmtSig)
            xlabel('Wavelength [nm]'), ylabel('Distance, z [mm]')
            zlabel('Instrument Response [au]')
        end  
        
        [path,~] = fileparts(Plaque_Cal_Filename);
        
        matfile_name = fullfile(path,['Hbb_Cal_Plaque_' datestr(cal.timestamp,'yyyymmdd_HHMMss') '.mat']);
        
        if isfile(matfile_name)
            promptMessage = sprintf('This calibration file already exists:\n\n%s\n\nDo you want to overwrite it?', matfile_name);
            titleBarCaption = 'Overwrite?';
            buttonText = questdlg(promptMessage, titleBarCaption, 'Yes', 'No', 'Yes');
            if strcmpi(buttonText, 'Yes')
                save(matfile_name,'cal');
                msgbox(['Calibration saved as ''' matfile_name ''''],'File Saved')
            end
        else
            save(matfile_name,'cal');
            msgbox(['Calibration saved as ''' matfile_name ''''],'File Saved')
        end
        saveas(gcf, fullfile(path, ['Cal_' datestr(cal.timestamp,'yyyymmdd_HHMMss') '_Fig4.jpg']))
        close gcf
        saveas(gcf, fullfile(path, ['Cal_' datestr(cal.timestamp,'yyyymmdd_HHMMss') '_Fig3.jpg']))
        close gcf
        saveas(gcf, fullfile(path, ['Cal_' datestr(cal.timestamp,'yyyymmdd_HHMMss') '_Fig2.jpg']))
        close gcf
        saveas(gcf, fullfile(path, ['Cal_' datestr(cal.timestamp,'yyyymmdd_HHMMss') '_Fig1.jpg']))
        close gcf
        
    end
end


function [pmtSigRef, pmtGamma] = FitPmtGainCurve(pmtGain, pmtSig, pmtGainRef)
% function [pmtSig0, gamma, resid, magn, rmse] = FitGainCurve(pmtGain, pmtSig, pmtGainRef)


%setting options for fmisearch
opts = optimset('fminsearch');      
opts = optimset(opts,'MaxIter',10000); 
opts = optimset(opts,'MaxFunEvals',1000);   % usually 100*number of params
opts = optimset(opts,'TolFun',1e-14);

% guess for paramters (PMT sig at reference gain, power law exponent)
x0 = [1 -10];

% minimization routine
x = fminsearch(@power_law, x0, opts, pmtGain, pmtSig, pmtGainRef);
% pmtSigHat = x(1).*(pmtGain./pmtGainRef).^x(2);

pmtSigRef = x(1);
pmtGamma = x(2);

return 


function err = power_law(x, pmtGain, pmtSigMeas, pmtGainRef)

    w = 1./pmtSigMeas.^2;
    pmtSigHat = x(1).*(pmtGain./pmtGainRef).^x(2);
    err = mean(w.*(pmtSigMeas-pmtSigHat).^2);

return

function y = nanmedian(x,dim)
% FORMAT: Y = NANMEDIAN(X,DIM)
% 
%    Median ignoring NaNs
%
%    This function enhances the functionality of NANMEDIAN as distributed
%    in the MATLAB Statistics Toolbox and is meant as a replacement (hence
%    the identical name).  
%
%    NANMEDIAN(X,DIM) calculates the mean along any dimension of the N-D
%    array X ignoring NaNs.  If DIM is omitted NANMEDIAN averages along the
%    first non-singleton dimension of X.
%
%    Similar replacements exist for NANMEAN, NANSTD, NANMIN, NANMAX, and
%    NANSUM which are all part of the NaN-suite.
%
%    See also MEDIAN

% -------------------------------------------------------------------------
%    author:      Jan Gläscher
%    affiliation: Neuroimage Nord, University of Hamburg, Germany
%    email:       glaescher@uke.uni-hamburg.de
%    
%    $Revision: 1.2 $ $Date: 2007/07/30 17:19:19 $

if isempty(x)
	y = [];
	return
end

if nargin < 2
	dim = min(find(size(x)~=1));
	if isempty(dim)
		dim = 1;
	end
end

siz  = size(x);
n    = size(x,dim);

% Permute and reshape so that DIM becomes the row dimension of a 2-D array
perm = [dim:max(length(size(x)),dim) 1:dim-1];
x = reshape(permute(x,perm),n,prod(siz)/n);


% force NaNs to bottom of each column
x = sort(x,1);

% identify and replace NaNs
nans = isnan(x);
x(isnan(x)) = 0;

% new dimension of x
[n m] = size(x);

% number of non-NaN element in each column
s = size(x,1) - sum(nans);
y = zeros(size(s));

% now calculate median for every element in y
% (does anybody know a more eefficient way than with a 'for'-loop?)
for i = 1:length(s)
	if rem(s(i),2) & s(i) > 0
		y(i) = x((s(i)+1)/2,i);
	elseif rem(s(i),2)==0 & s(i) > 0
		y(i) = (x(s(i)/2,i) + x((s(i)/2)+1,i))/2;
	end
end

% Protect against a column of NaNs
i = find(y==0);
y(i) = i + nan;

% permute and reshape back
siz(dim) = 1;
y = ipermute(reshape(y,siz(perm)),perm);

% $Id: nanmedian.m,v 1.2 2007/07/30 17:19:19 glaescher Exp glaescher $
