%% Index the cells in Suite2P-analzying data that had calcium transients
% takes the first column from the matrix 'iscell' 
%   (the 0s for 'not a cell' and 1s for 'this is a cell'; we have 17 '1s')
iscell_firstcolumn = iscell(:,1);
% returns the indices of the 'first column of iscell' that point to nonzero elements
%   so now you know which row/cells #s had '1s' indicating 'this is a cell'
iscell_id = find(iscell_firstcolumn);
% in 'F' (the raw fluorescence data for cells) take the data for the
%   17 cells that had '1s'
Fcell = F(iscell_id,:);
% in 'Fneu' (the raw fluorescence data for neuropil) take the data for the
%   17 cells that had '1s'
Fneu2 = Fneu(iscell_id,:);


%% Neuropil correction
% constant used for neuropil correction
k = 0.7;

% neuropil subtraction equation (correct for neuropil contamination)
%   subtract a neuropil signal from the somatic signal; 
%   define a neuropil region around the ROI (soma) and then subtract the
%   spatially averaged (local) neuropil signal from the somatic signal;
%   however, need to decide how much weight to give to the subtracted
%   neuropil signal (therefore, we use the constant k)
Fcorr = Fcell - k*Fneu2;

% takes the 1st of the 17 cells
% can simply change this number to get a different cell
nr =1;
%I think that this is just to make sure these figures have a distinct #
figure(1000+nr)
 % plots the raw fluo data for the 1st cell with signal, time vs F value
plot(Fcell(nr,:))
% retains the current plot and certain axes properties so that subsequent graphing commands add to the existing graph
hold on
% plots the raw fluo data for the neuropil of the 1st cell with signal, time vs F value 
plot(Fneu2(nr,:))
% plots the signal, now corrected for neuropil contamination
plot(Fcorr(nr,:))
hold off


%% Add variables to the workspace
% counts how many cells we have with '1s' (17) and add it to the workspace
nCells = size(Fcorr,1);
% we create a variable saying how many trials we ran
nTrials = 4;
% we create a variable saying how many drifting gratings we used
nDir = 8;

% we create a variable saying how much time preceded the stimulus (2s)
preStim = 2; 
% we create a variable saying how long the stimulus lasted (3s)
stim = 3;
% we create a variable saying how much time followed/succeeded the stimulus (2s)
poststim = 1.5;
% we create a variable defining how much time we will use for our baseline
% (1s, averaged, before the stimulus starts)
baseline = 1;

% define a frame rate (in Hz)
frame_rate = 30.94;
% define a frame time (in sec)
frame_time = 1/frame_rate; %0.03s %30ms

% these are rounded, because can't use frame # if #s aren't rounded...
% could also put it in time...
% 'time preceding the stimulus (2s)' x 'frame rate' = 
%   how many frames before the stimulus starts/stimulus onset
preStim_fr  = round( preStim*frame_rate );
% 'how long the stimulus lasted (3s)' x 'frame rate' = 
%   how many frames during stimulus presentation
stim_fr     = round( stim*frame_rate );
% 'time following/succeeding the stimulus (2s)' x 'frame rate' = 
%   how many frames after the stimulus ends & before the next stimulus onset 
poststim_fr = round( poststim*frame_rate );
% 'baseline (1s, averaged, before the stimulus)' x 'frame rate' = 
%   how many frames during the baseline
baseline_fr = round( baseline*frame_rate );
% 'how many frames during the baseline' + 
% 'how many frames during stimulus presentation'
% = 'how many frames during baselines AND stimulus presentation (1 trial)'
trial_fr = baseline_fr + stim_fr;


%% Create a matrix with the 4 dimensions of my experiment
% nan = not a number
% make a matrix with 4 dimensions:
%   'frames in 1 trial', '# of drifting gratings used', 
%   'how many trials we ran', 'how many cells we have with '1s''
Ftrial = nan(trial_fr, nDir, nTrials, nCells);


%% Isolate traces corresponding to particular grating stimuli for each cell across trials
% run the loop for all cells, 1-17
for c = 1 : nCells
    % run the loop for trials, 1-4
    for r = 1 : nTrials
    
        %run the loop for all drifting gratings used, 1-8
        for d = 1 : nDir
            
            %prestim +1 gives you the first stim onset in frames (63)
            %prestim + stim gives you the first stim offset in frames
            %loop through all 8 gratings for the first trial
            %then 2nd trial, and go through the next 8 gratings
            %do for all 4 trials
            
           ix_onset  = preStim_fr + (r-1)*(nDir*(stim_fr+poststim_fr)) + (d-1)*(stim_fr+poststim_fr) + 1;
           ix_offset = preStim_fr + (r-1)*(nDir*(stim_fr+poststim_fr)) + (d-1)*(stim_fr+poststim_fr) + stim_fr;
                      
           %Fcorr has the traces for all 17 cells
           %for cell 1, take trace from onset to offset
            Fstim_tmp     = Fcorr(c,ix_onset:ix_offset);
           %for cell 1, take trace from start of baseline 
           %until end of baseline
            Fbaseline_tmp = Fcorr(c,ix_onset-baseline_fr:ix_onset-1);
            %concatenate baseline and stim time into a vector
            % (for 1 cell, grating, trial...)
            Ftrial_tmp = [Fbaseline_tmp, Fstim_tmp];
            
            %calculate the deltaF/F
            Fbaseline_mean = mean(Fbaseline_tmp);
            dFoF_tmp = (Fstim_tmp - Fbaseline_mean)/(Fbaseline_mean);
            
            %All of the data is in here
            Ftrial(:,d,r,c) = Ftrial_tmp;
            %All of the dFoF values
            dFoFtrial(:,d,r,c) = dFoF_tmp;   
            %dFoF_avg = mean(dFoF_tmp)
                        
        end
    end
    
end


%% For a single cell make 8 plots, one for each grating, comparing time with F
%for 1 cell, 8 plots for each grating; change the #, get a different cell
nr = 1;

figure(2000+nr); clf

%loop through the 1-8 gratings Ftrial
for d = 1 : nDir
    %make 8 separate plots with the 4 trials in each
    %all 4 belong to 1 cell
    subplot(1,nDir,d);
    hold on
        for r = 1 : nTrials
        plot( Ftrial(:,d,r,nr) )
    end
end


%% For a single cell make 8 plots, one for each grating, comparing time with dF/F
%for 1 cell, 8 plots for each grating; change the #, get a different cell
nr = 1;

figure(3000+nr); clf

%loop through the 1-8 gratings dFoF
for d = 1 : nDir
    %make 8 separate plots with the 4 trials in each
    %all 4 belong to 1 cell
    subplot(1,nDir,d);
    hold on
        for r = 1 : nTrials
        plot( dFoFtrial(:,d,r,c) )
        % plot (dFoF_avg, 'LineWidth', 100)
        % how to add the mean to each plot...
    end
end


%% Plot mean dF/F for different stimulus gratings
