fn_photo = ['\Users\Wassum Lab\Documents\Fiber Photometry\Raw Data\AL Pav Ph 22 SSS deval test 2 suc deval NPM2021-01-20T13_15_41.csv']; %path to file containing photometry data
fn_behavior = ['\Users\Wassum Lab\Documents\Fiber Photometry\Raw Data\AL Pav Ph 22 SSS deval test 2 suc deval MPC2021-01-20T13_15_41.csv']; %path to file containing behavioral data

%load imaging & behavioral data
photoData_raw = load(fn_photo); %load GCaMP data file
behaviorData_raw = load(fn_behavior); %load the Med TTLs

%convert timestamps to seconds 
photoData_raw(:,1) = photoData_raw(:,1)./1000;
behaviorData_raw(:,1) = behaviorData_raw(:,1)./1000;

%plot raw gcamp
plot(photoData_raw(:,1), photoData_raw(:,2))

%synchronize timestamps
imageStart = photoData_raw(1,1);
photoData_raw(:,1) = photoData_raw(:,1)-imageStart;
behaviorData_raw(:,1) = behaviorData_raw(:,1)-imageStart;
behaviorData_raw(find(behaviorData_raw(:,1)>photoData_raw(end,1)),:)=[];%remove behavior data time points that are past photometry time points
plot(photoData_raw(:,1), photoData_raw(:,2))

subplot(3,1,1)
plot(behaviorData_raw(:,1),behaviorData_raw(:,3));
xlim([0 photoData_raw(end,1)])
subplot(3,1,2)
plot(behaviorData_raw(:,1),behaviorData_raw(:,4));
xlim([0 photoData_raw(end,1)])
subplot(3,1,3)
plot(photoData_raw(:,1),photoData_raw(:,2))
xlim([0 photoData_raw(end,1)])

%run this if animal became detached and signal drops
thresh=2600;
time = [80 96];

if size(time,2)==2
    weird_frames=find(photoData_raw(:,1)>time(1) & photoData_raw(:,1)<time(2)); %set the limits to take out weird frames and set time parameters you could like to remove the weird frames from
    time_points_delete=photoData_raw(weird_frames,1); %determine the timepoints you are removing 
    delete_frames=find(photoData_raw(:,1) <= time_points_delete(end) & photoData_raw(:,1) >= time_points_delete(1)); %choose which frames to delete
    delete_frames_medPC=(find(behaviorData_raw(:,1) >= time_points_delete(1)): (find(behaviorData_raw(:,1) >= time_points_delete(end),1)-1))'; %finds same time period to delete from the medPC file
    gcamp_time_subtract= photoData_raw(delete_frames(end)+1,1)-photoData_raw(delete_frames(1),1);
    medPC_time_subtract= behaviorData_raw(delete_frames_medPC(end)+1,1)-behaviorData_raw(delete_frames_medPC(1),1);
    photoData_raw(delete_frames,:)=[];
    behaviorData_raw(delete_frames_medPC,:)=[];
    photoData_raw(delete_frames(1):end,1)=photoData_raw(delete_frames(1):end,1)-gcamp_time_subtract;
    behaviorData_raw(delete_frames_medPC(1):end,1)=behaviorData_raw(delete_frames_medPC(1):end,1)-medPC_time_subtract;
else 
    weird_frames=find(photoData_raw(:,2)<thresh & photoData_raw(:,1)>time(1)); %set the limits to take out weird frames and set time parameters you could like to remove the weird frames from
    time_points_delete=photoData_raw(weird_frames,1); %determine the timepoints you are removing 
    delete_frames=find(photoData_raw(:,1)>= time_points_delete(1)); %choose which frames to delete
    photoData_raw(delete_frames:end,:)=[];
    behaviorData_raw(find(behaviorData_raw(:,1)>photoData_raw(end,1)),:)=[];%remove behavior data time points that are past photometry time points
end

plot(photoData_raw(:,1), photoData_raw(:,2))

%check for negative timestamps & add a buffer if necessary
if behaviorData_raw(1,1)<0;
    pad = abs(behaviorData_raw(1,1))+0.5;
    photoData_raw(:,1) = photoData_raw(:,1) + pad;
    behaviorData_raw(:,1) = behaviorData_raw(:,1) + pad;
end


%look at de-interleaved data and if it looks off, remove any detected dropped frames
evenData = photoData_raw(2:2:end,:);
oddData = photoData_raw(1:2:end,:);

plot(evenData(:,1),evenData(:,2))
hold on
plot(oddData(:,1),oddData(:,2))

IFI = diff(photoData_raw(:,1)); %finds inter-frame interval
plot(IFI)
[~,droppedFrames] = findpeaks(IFI,'MinPeakProminence',0.05); %find any extra-long gaps between successive frames, change 0.05 threshold as needed

if ~isempty(droppedFrames)
    photoData_raw(droppedFrames,:) = [];
else
    droppedFrames = 0;
end

%assign data as gcamp or isos - Haven't figured how to separate
%these...have to decide if even frames or odd frames are the gcamp
evenData = photoData_raw(2:2:end,:);
oddData = photoData_raw(1:2:end,:);

plot(evenData(:,1),evenData(:,2))
hold on
plot(oddData(:,1),oddData(:,2))

gcampData=oddData; %change based on whether gcamp is the odd or even frames
isosData=evenData;

if length(gcampData(:,1))>length(isosData(:,1));
    gcampData(end,:) = [];
elseif length(isosData(:,1))>length(gcampData(:,1));
    isosData(end,:) = [];
end

%package uncorrected data for reference
gcamp_raw = tsd(gcampData(:,1),gcampData(:,2)); %this converts the gcamp data from a matrix to a tsd format where gcamp_raw.range is the timestamp and gcamp_raw.data is the gcamp values
isos_raw = tsd(isosData(:,1),isosData(:,2));


%% calculate df/f (470-fitted 415)/fitted 415

temp_fit=fit(isosData(:,2),gcampData(:,2),fittype('poly1'));%,'Robust','on');  %fit isos to gcamp
isos_fit=tsd(isosData(:,1), temp_fit(isosData(:,2))); %scale isos to gcamp based on fit 
plot(gcampData(:,1),gcampData(:,2),isos_fit.range,isos_fit.data) %plot the gcamp and fitted isos trace on same graph
df_F=(gcampData(:,2)-isos_fit.data)./isos_fit.data; %subtract isos from gcamp and normalize to isos
plot (gcampData(:,1),df_F(:,1)) %plot subtracted/normalized gcamp data
gcamp_minus_isos_df_F=tsd(gcampData(:,1), df_F); %package df_f as tsd


z_df_F=(gcamp_minus_isos_df_F.data-mean(gcamp_minus_isos_df_F.data))/std(gcamp_minus_isos_df_F.data); %z-score df/F across entire session
z_df_F =tsd(gcamp_minus_isos_df_F.range, z_df_F); %package z_df_F as tsd

%% find the behavioral timestamp corresponding to the TTL onsets

%Usually the first TTL on the MedPC interface for box 2 in photometry room
%has a smaller range of units than the other TTLs. This scales the values
%of TTL1 so the same threshold can be applied to all TTLs
if max(behaviorData_raw(:,2))<1
    behaviorData_raw(:,2)=behaviorData_raw(:,2)*1000;
end

%This is going to loop through all of the columns in corresponding to
%unique TTLs (ie event types) in the behaviorData variable and will find
%the onset of each trial for that event type by finding all timepoints in
%which the values in the TTL column exceed a set threshold
nVars = size(behaviorData_raw,2);
medVars = cell(nVars-1,1);
for iV = 2:nVars;
    clear signal tRise tFall
    signal = tsd(behaviorData_raw(:,1), behaviorData_raw(:,iV));
    [tRise, tFall]  = ttlEdges(signal);
    
    medVars{iV-1} = ts(tRise.data);
end

%% Align raw gcamp and fitted isos to medPC TTLs
%loop through medPC TTLs and extract gcamp data that centers around the
%onset of that TTL (ie TTL onset = time 0) within a given window
%this section will give you the aligned raw gcamp data and the aligned
%fitted isos data

%titles={'CS1','CS2','US1','US2','US1col','US2col','pN','discard1','Np','Tp','discard2','CSavg','USavg','UScolavg'};
titles={'start','CS1','CS2','US1','US2','US1col','US2col','pN','pT','Np','Tp'};%,'CSavg','USavg','UScolavg'};

input1=gcamp_raw;
input2= isos_fit;

[aligned_gcamp_raw,aligned_isos_fit]=align_gcampRaw_isosFit_medPC(titles,input1,input2,medVars,0,'CS',[-5 125],[-5 30]); %Inside the parentheses input titles,raw_gcamp input variable, fitted isos input variable, medVars, number of vars to average across, string for variables with specific time window, time window for specific variables, defaulttime window to cut around

%% Align gcamp to medPC TTLs - gcamp minus isos df_F input, 1sec baseline Z output
%loop through medPC TTLs and extract gcamp data that centers around the
%onset of that TTL (ie TTL onset = time 0) within a given window 
%this section will give you aligned df_F (gcamp vs isos) and aligned df_F 
%that has been z-scored relative to baseline window


input3=gcamp_minus_isos_df_F; %change this to whichever version of the processed gcamp data you want to use for behavioral alignment

[aligned_gcamp_minus_isos,aligned_gcamp_minus_isos_zBL]=align_df_F_zBL_medPC(titles,input3, medVars,0,'CS',[-5 125],[-5 30],[-4 -5]); %Inside the parentheses input titles,df_F_input, medVars,
%number of vars to average across, string for variables with specific time window, time window for specific variables, defaulttime window to cut around, time window
%of your baseline IN THAT ORDER

%% Align whole session z-scored df_F (gcamp vs isos) to medPC TTLs
%loop through medPC TTLs and extract gcamp data that centers around the
%onset of that TTL (ie TTL onset = time 0) within a given window
%this section will give you aligned df_F z-scored relative to the whole session avg and SD

input4=z_df_F; %change this to whichever version of the processed gcamp data you want to use for behavioral alignment

[aligned_whole_session_z_df_F]=align_z_df_F_medPC(titles,input4, medVars,0,'CS',[-5 125],[-5 30]); %Inside the parentheses input titles,df_F_input, medVars,
%number of vars to average across, string for variables with specific time window, time window for specific variables, defaulttime window to cut around

%% export aligned data to excel

for i=1:size(aligned_gcamp_raw,2);
    temp_table=array2table(aligned_gcamp_raw{2,i});
    sheetname=strcat(char(aligned_gcamp_raw{1,i}),'_gcampRaw');
    writetable(temp_table,'AS03_25L_Pav1pt2_gcampRaw_isosFit2gcamp.xlsx','Sheet',sheetname);
    temp_table=array2table(aligned_isos_fit{2,i});
    sheetname=strcat(char(aligned_isos_fit{1,i}),'_isosFit');
    writetable(temp_table,'AS03_25L_Pav1pt2_gcampRaw_isosFit2gcamp.xlsx','Sheet',sheetname);
end


for i=1:size(aligned_gcamp_minus_isos,2);
    temp_table=array2table(aligned_gcamp_minus_isos{1,i});
    sheetname=char(titles{i});
    writetable(temp_table,'AS03_25L_Pav1pt2_gcamp_minus_isos.xlsx','Sheet',sheetname);
end

% for i=1:size(aligned_gcamp_minus_isos_zBL,2);
%     temp_table=array2table(aligned_gcamp_minus_isos_zBL{2,i});
%     sheetname=char(aligned_gcamp_minus_isos_zBL{1,i});
%     writetable(temp_table,'AS03_10L_Pav2_gcamp_minus_isos_zBL.xlsx','Sheet',sheetname);
% end

for i=1:size(aligned_whole_session_z_df_F,2);
    temp_table=array2table(aligned_whole_session_z_df_F{2,i});
    sheetname=char(aligned_whole_session_z_df_F{1,i});
    writetable(temp_table,'AS03_25L_Pav1pt2_whole_session_z_df_F.xlsx','Sheet',sheetname);
end


subplot(12,1,1)
plot(behaviorData_raw(:,1),behaviorData_raw(:,2))
subplot(12,1,2)
plot(behaviorData_raw(:,1),behaviorData_raw(:,3))
subplot(12,1,3)
plot(behaviorData_raw(:,1),behaviorData_raw(:,4))
subplot(12,1,4)
plot(behaviorData_raw(:,1),behaviorData_raw(:,5))
subplot(12,1,5)
plot(behaviorData_raw(:,1),behaviorData_raw(:,6))
subplot(12,1,6)
plot(behaviorData_raw(:,1),behaviorData_raw(:,7))
subplot(12,1,7)
plot(behaviorData_raw(:,1),behaviorData_raw(:,8))
subplot(12,1,8)
plot(behaviorData_raw(:,1),behaviorData_raw(:,9))
subplot(12,1,9)
plot(behaviorData_raw(:,1),behaviorData_raw(:,10))
subplot(12,1,10)
plot(behaviorData_raw(:,1),behaviorData_raw(:,11))
subplot(12,1,11)
plot(behaviorData_raw(:,1),behaviorData_raw(:,12))
subplot(12,1,12)
plot(photoData_raw(:,1),photoData_raw(:,2))