diff --git a/.gitattributes b/.gitattributes
new file mode 100644
index 0000000..dfe0770
--- /dev/null
+++ b/.gitattributes
@@ -0,0 +1,2 @@
+# Auto detect text files and perform LF normalization
+* text=auto
diff --git a/get_peaks.asv b/get_peaks.asv
new file mode 100644
index 0000000..c572281
--- /dev/null
+++ b/get_peaks.asv
@@ -0,0 +1,185 @@
+%after saving the data with new_data_set.m, we isolate the peaks for the
+%new analysis of the refined data
+
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+channelfilename = [newdatadir 'refined_dataset']; 
+data_file = load(channelfilename);
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ %% find peaks in order to analyze them on R and fit a LMER
+ 
+channum = 1: length(data_file.new_data);
+xabs = -199:1300;
+nyq = 500;
+mean_filtered_dSUA = struct();
+
+mean_blank_power = nan(length(channum),1);
+mean_fourhzpowerhighcontrast = nan(length(channum),1);
+
+all_locsdSUA_filtered = nan(1,length(channum));
+norm_data_peaks = struct();
+ 
+ for i = channum
+ data = squeeze(data_file.new_data(i).channel_data.hypo{1,2}.cont_su(401:1900,:,:));
+ 
+ filename = [data_file.new_data(i).channel_data.filename, f{2}];
+ filename = erase(filename, '.mat');
+
+ all_norm_pks = nan(4,length(data(1,:)));
+ 
+ 
+   blankcontrast = data_file.new_data(i).channel_data.contrast ==  0 & data_file.new_data(i).channel_data.fixedc ==  0;
+   highcontrast = data_file.new_data(i).channel_data.contrast >=  0.5 & data_file.new_data(i).channel_data.fixedc ==  0; 
+    
+    trialidx = 1:length(data_file.new_data(i).channel_data.sdftr_chan(1,:));
+    raw_bs = nan(length(xabs), length(trialidx));
+    filtered_dSUA = nan(length(xabs), length(trialidx));
+    all_norm_lpdSUA= nan(length(xabs),length(trialidx));
+
+    
+    powerstim = nan(length(trialidx),1025);
+    freqstim = nan(length(trialidx),1025);
+    fourhzpowerstim =nan(length(trialidx),1);
+    bsl = nan(1, length(trialidx));
+    mean_wnd1 = nan(1,length(trialidx));
+     for tridx = trialidx
+      
+    all_data = data_file.new_data(i).channel_data.sdftr_chan(401:1900,tridx);
+      
+    bsl(tridx) = mean(all_data(1:200));
+    raw_bs(:,tridx) = all_data(1:end)- bsl(tridx);
+   
+    lpc       = 4.5; %low pass cutoff
+    lWn       = lpc/nyq;
+    [bwb,bwa] = butter(4,lWn,'low');
+    lpdSUA      = filtfilt(bwb,bwa, raw_bs(:,tridx));
+  
+
+    filtered_dSUA(:,tridx) = lpdSUA;
+    all_norm_lpdSUA(:,tridx) = (lpdSUA - min(lpdSUA))/(max(lpdSUA)- min(lpdSUA));
+    mean_wnd1(tridx) = mean(lpdSUA(231:480));
+    
+    %%% power
+
+          
+    [powerstim(tridx,:), freqstim(tridx,:)] = calcFFT(all_data(200:1350));
+        
+    %find the index of the frequency vector closest to 4hz and point to the
+    %power value of this index for every trial, and store the value in
+    %fourhzpower
+    [val,index] = min(abs(4-freqstim(tridx,:)));
+    fourhzpowerstim(tridx,1) = powerstim(tridx,index);
+    
+     end
+   %power related variables
+   power0 = fourhzpowerstim(blankcontrast);
+   power5 = fourhzpowerstim(highcontrast);
+   mean_blank_power(i,1) = mean(power0);
+   mean_fourhzpowerhighcontrast(i,1) = mean(power5);
+   
+   %spiking activity related variables
+   mean_wnd1_5 =mean_wnd1(highcontrast);
+   filtered_dSUA_high = filtered_dSUA(:, highcontrast);
+   all_norm_lpdSUA_high = all_norm_lpdSUA(:, highcontrast);
+   %first peak location related variables 
+    sua_bsl =  mean(filtered_dSUA_high(1:200,:),1);
+    
+
+  %%%%%%%%%%% %reject trials below the 95%CI in the blank condition %%%%%%
+   for tr = 1:length(power5)
+       if power5(tr) > mean(power0)+1.96*std(power0) && mean_wnd1_5(tr) > mean(sua_bsl) + 1.96*std(sua_bsl)
+     
+           filtered_dSUA_high(:,tr) = filtered_dSUA_high(:,tr);
+           
+       else
+   
+           filtered_dSUA_high(:,tr) = nan(length(filtered_dSUA_high(:,tr)),1);
+       end
+   end
+   filtered_dSUA_high = filtered_dSUA_high(:,~all(isnan(filtered_dSUA_high))); % for nan - cols
+   
+   
+   all_locsdSUA_trials = nan(6,length(filtered_dSUA_high(1,:)));
+   if length(filtered_dSUA_high(1,:)) >=10
+   
+    %determine the first peak location for each trial of a given single
+    %unit
+    
+  
+    for trial = 1:length(filtered_dSUA_high(1,:))
+        
+         for ln = 30:550
+             if filtered_dSUA_high(200+ln,trial) < filtered_dSUA_high(200+ln+1,trial)
+              
+            
+             locsdSUA_trial = findpeaks(filtered_dSUA_high(200+ln:1350,trial));
+             %store 4 peaks locations 
+             all_locsdSUA_trials(1:length(locsdSUA_trial.loc),trial) = locsdSUA_trial.loc+200+ln;
+             break 
+             end 
+         end
+         
+        if length(locsdSUA_trial.loc) >= 4
+             %adjust location to the first data point of lpsu (+len),
+    lpsulocs = locsdSUA_trial.loc(1:4) + 200+ ln;
+    all_norm_pks(:,trial) = all_norm_lpdSUA_high(lpsulocs(1:4), trial);
+  
+         end 
+    end
+    elseif length(filtered_dSUA_high(1,:)) <=10  || length(locsdSUA_trial.loc) <= 4
+      all_locsdSUA_trials(:,:) = nan(length(all_locsdSUA_trials(:,1)),);
+      all_norm_pks(:,trial) = nan(length(all_norm_pks(:,1)),1);
+   end
+    namelist(1,1:length(sprintf('chan_%d',i))) = sprintf('chan_%d',i);
+ norm_data_peaks(i).namelist = all_norm_pks(:,~all(isnan(all_norm_pks)));
+channelfilename = [newdatadir filename];
+%save(strcat(channelfilename, '.mat'), 'all_pks');
+ end
+  allfilename = [newdatadir 'all_norm_data_peaks'];
+ save(strcat(allfilename, '.mat'), 'norm_data_peaks');
+ 
+    %{
+ %all_locs =nan(4,length(channel_data(1,1,:)));
+ for n =1:length(data(1,:))
+   lpc       = 4.5; %low pass cutoff
+   lWn       = lpc/nyq;
+   [bwb,bwa] = butter(4,lWn,'low');
+   lpsu      = filtfilt(bwb,bwa, data(:,n));
+   
+   norm_lpsu = (lpsu - min(lpsu))/(max(lpsu)- min(lpsu));
+   %all_lpsu(:,n) = lpsu;
+  % plot(lpsu)
+ if ~all(isnan(norm_lpsu)) 
+  for len = 30:550
+            if norm_lpsu(len) < norm_lpsu(len+1)
+   locs = findpeaks(norm_lpsu(len:1200));
+        break
+            end
+  end
+         
+         if length(locs.loc) >= 4
+             %adjust location to the first data point of lpsu (+len),
+    lpsulocs = locs.loc(1:4) + len;
+   all_norm_pks(:,n) = norm_lpsu(lpsulocs(1:4));
+  
+         end 
+ end
+% plot(locs.loc +len, lpsu(locs.loc +len))
+ end
+ 
+ namelist(1,1:length(sprintf('chan_%d',i))) = sprintf('chan_%d',i);
+ norm_data_peaks(i).namelist = all_norm_pks;
+channelfilename = [newdatadir filename];
+%save(strcat(channelfilename, '.mat'), 'all_pks');
+ end
+    %}
+ allfilename = [newdatadir 'all_norm_data_peaks'];
+ save(strcat(allfilename, '.mat'), 'norm_data_peaks');
\ No newline at end of file
diff --git a/get_peaks.m b/get_peaks.m
new file mode 100644
index 0000000..adb6d3a
--- /dev/null
+++ b/get_peaks.m
@@ -0,0 +1,186 @@
+%after saving the data with new_data_set.m, we isolate the peaks for the
+%new analysis of the refined data
+
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+channelfilename = [newdatadir 'refined_dataset']; 
+data_file = load(channelfilename);
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ %% find peaks in order to analyze them on R and fit a LMER
+ 
+channum = 1: length(data_file.new_data);
+xabs = -199:1300;
+nyq = 500;
+mean_filtered_dSUA = struct();
+
+mean_blank_power = nan(length(channum),1);
+mean_fourhzpowerhighcontrast = nan(length(channum),1);
+
+all_locsdSUA_filtered = nan(1,length(channum));
+norm_data_peaks = struct();
+ 
+ for i = channum
+ data = squeeze(data_file.new_data(i).channel_data.hypo{1,2}.cont_su(401:1900,:,:));
+ 
+ filename = [data_file.new_data(i).channel_data.filename, f{2}];
+ filename = erase(filename, '.mat');
+
+ all_norm_pks = nan(4,length(data(1,:)));
+ 
+ 
+   blankcontrast = data_file.new_data(i).channel_data.contrast ==  0 & data_file.new_data(i).channel_data.fixedc ==  0;
+   highcontrast = data_file.new_data(i).channel_data.contrast >=  0.5 & data_file.new_data(i).channel_data.fixedc ==  0; 
+    
+    trialidx = 1:length(data_file.new_data(i).channel_data.sdftr_chan(1,:));
+    raw_bs = nan(length(xabs), length(trialidx));
+    filtered_dSUA = nan(length(xabs), length(trialidx));
+    all_norm_lpdSUA= nan(length(xabs),length(trialidx));
+
+    
+    powerstim = nan(length(trialidx),1025);
+    freqstim = nan(length(trialidx),1025);
+    fourhzpowerstim =nan(length(trialidx),1);
+    bsl = nan(1, length(trialidx));
+    mean_wnd1 = nan(1,length(trialidx));
+     for tridx = trialidx
+      
+    all_data = data_file.new_data(i).channel_data.sdftr_chan(401:1900,tridx);
+      
+    bsl(tridx) = mean(all_data(1:200));
+    raw_bs(:,tridx) = all_data(1:end)- bsl(tridx);
+   
+    lpc       = 4.5; %low pass cutoff
+    lWn       = lpc/nyq;
+    [bwb,bwa] = butter(4,lWn,'low');
+    lpdSUA      = filtfilt(bwb,bwa, raw_bs(:,tridx));
+  
+
+    filtered_dSUA(:,tridx) = lpdSUA;
+    all_norm_lpdSUA(:,tridx) = (lpdSUA - min(lpdSUA))/(max(lpdSUA)- min(lpdSUA));
+    mean_wnd1(tridx) = mean(lpdSUA(231:480));
+    
+    %%% power
+
+          
+    [powerstim(tridx,:), freqstim(tridx,:)] = calcFFT(all_data(200:1350));
+        
+    %find the index of the frequency vector closest to 4hz and point to the
+    %power value of this index for every trial, and store the value in
+    %fourhzpower
+    [val,index] = min(abs(4-freqstim(tridx,:)));
+    fourhzpowerstim(tridx,1) = powerstim(tridx,index);
+    
+     end
+   %power related variables
+   power0 = fourhzpowerstim(blankcontrast);
+   power5 = fourhzpowerstim(highcontrast);
+   mean_blank_power(i,1) = mean(power0);
+   mean_fourhzpowerhighcontrast(i,1) = mean(power5);
+   
+   %spiking activity related variables
+   mean_wnd1_5 =mean_wnd1(highcontrast);
+   filtered_dSUA_high = filtered_dSUA(:, highcontrast);
+   all_norm_lpdSUA_high = all_norm_lpdSUA(:, highcontrast);
+   %first peak location related variables 
+    sua_bsl =  mean(filtered_dSUA_high(1:200,:),1);
+    
+
+  %%%%%%%%%%% %reject trials below the 95%CI in the blank condition %%%%%%
+   for tr = 1:length(power5)
+       if power5(tr) > mean(power0)+1.96*std(power0) && mean_wnd1_5(tr) > mean(sua_bsl) + 1.96*std(sua_bsl)
+     
+           filtered_dSUA_high(:,tr) = filtered_dSUA_high(:,tr);
+           
+       else
+   
+           filtered_dSUA_high(:,tr) = nan(length(filtered_dSUA_high(:,tr)),1);
+       end
+   end
+   filtered_dSUA_high = filtered_dSUA_high(:,~all(isnan(filtered_dSUA_high))); % for nan - cols
+   
+   
+   all_locsdSUA_trials = nan(6,length(filtered_dSUA_high(1,:)));
+   if length(filtered_dSUA_high(1,:)) >=10
+   
+    %determine the first peak location for each trial of a given single
+    %unit
+    
+  
+    for trial = 1:length(filtered_dSUA_high(1,:))
+        
+         for ln = 30:550
+             if filtered_dSUA_high(200+ln,trial) < filtered_dSUA_high(200+ln+1,trial)
+              
+            
+             locsdSUA_trial = findpeaks(filtered_dSUA_high(200+ln:1350,trial));
+             %store 4 peaks locations 
+             all_locsdSUA_trials(1:length(locsdSUA_trial.loc),trial) = locsdSUA_trial.loc+200+ln;
+             break 
+             end 
+         end
+         
+        if length(locsdSUA_trial.loc) >= 4
+             %adjust location to the first data point of lpsu (+len),
+    lpsulocs = locsdSUA_trial.loc(1:4) + 200+ ln;
+    all_norm_pks(:,trial) = all_norm_lpdSUA_high(lpsulocs(1:4), trial);
+  
+         end 
+    end
+    elseif length(filtered_dSUA_high(1,:)) <=10  || length(locsdSUA_trial.loc) <= 4
+      all_locsdSUA_trials(:,:) = nan(length(all_locsdSUA_trials(:,1)),length(all_locsdSUA_trials(1,:)));
+      all_norm_pks(:,:) = nan(length(all_norm_pks(:,1)),length(all_norm_pks(1,:)));
+   end
+ namelist(1,1:length(sprintf('chan_%d',i))) = sprintf('chan_%d',i);
+ norm_data_peaks(i).namelist = all_norm_pks(:,~all(isnan(all_norm_pks)));
+ all_norm_pks = all_norm_pks(:,~all(isnan(all_norm_pks)));
+channelfilename = [newdatadir 'su_peaks\' filename];
+save(strcat(channelfilename, '.mat'), 'all_norm_pks');
+ end
+  allfilename = [newdatadir 'all_norm_data_peaks'];
+ save(strcat(allfilename, '.mat'), 'norm_data_peaks');
+ 
+    %{
+ %all_locs =nan(4,length(channel_data(1,1,:)));
+ for n =1:length(data(1,:))
+   lpc       = 4.5; %low pass cutoff
+   lWn       = lpc/nyq;
+   [bwb,bwa] = butter(4,lWn,'low');
+   lpsu      = filtfilt(bwb,bwa, data(:,n));
+   
+   norm_lpsu = (lpsu - min(lpsu))/(max(lpsu)- min(lpsu));
+   %all_lpsu(:,n) = lpsu;
+  % plot(lpsu)
+ if ~all(isnan(norm_lpsu)) 
+  for len = 30:550
+            if norm_lpsu(len) < norm_lpsu(len+1)
+   locs = findpeaks(norm_lpsu(len:1200));
+        break
+            end
+  end
+         
+         if length(locs.loc) >= 4
+             %adjust location to the first data point of lpsu (+len),
+    lpsulocs = locs.loc(1:4) + len;
+   all_norm_pks(:,n) = norm_lpsu(lpsulocs(1:4));
+  
+         end 
+ end
+% plot(locs.loc +len, lpsu(locs.loc +len))
+ end
+ 
+ namelist(1,1:length(sprintf('chan_%d',i))) = sprintf('chan_%d',i);
+ norm_data_peaks(i).namelist = all_norm_pks;
+channelfilename = [newdatadir filename];
+%save(strcat(channelfilename, '.mat'), 'all_pks');
+ end
+    %}
+ allfilename = [newdatadir 'all_norm_data_peaks'];
+ save(strcat(allfilename, '.mat'), 'norm_data_peaks');
\ No newline at end of file
diff --git a/new_data_set.asv b/new_data_set.asv
new file mode 100644
index 0000000..0cc4878
--- /dev/null
+++ b/new_data_set.asv
@@ -0,0 +1,15 @@
+gooddatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\';
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+
+
+channelfilename = [gooddatadir 'good_single_units_data_4bmpmore']; 
+data_file = load(channelfilename);
+
+%Reject 160517 as no well triggered (1)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+new_data = data_file.good_data([2:45,47:54,56:end]);
+
+allfilename = [newdatadir 'refined_dataset'];
+save(strcat(allfilename, '.mat'), 'new_data', '-v7.3');
+ 
\ No newline at end of file
diff --git a/new_data_set.m b/new_data_set.m
new file mode 100644
index 0000000..196172d
--- /dev/null
+++ b/new_data_set.m
@@ -0,0 +1,15 @@
+gooddatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\';
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+
+
+oldchannelfilename = [gooddatadir 'good_single_units_data_4bmpmore']; 
+old_data_file = load(oldchannelfilename);
+
+%Reject 160517, K cell, as no well triggered (1)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+new_data = old_data_file.good_data([2:45,47:54,56:end]);
+
+allfilename = [newdatadir 'refined_dataset'];
+save(strcat(allfilename, '.mat'), 'new_data', '-v7.3');
+ 
\ No newline at end of file
diff --git a/single_units_analysis/align_clean_data.m b/single_units_analysis/align_clean_data.m
new file mode 100644
index 0000000..2821100
--- /dev/null
+++ b/single_units_analysis/align_clean_data.m
@@ -0,0 +1,207 @@
+%this script was developped in order to align the refined data(selected
+%channels after purification through get_clean_peaks_and_data.m) to the
+%first peak
+
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\su_peaks_02262020_2\all_units\';
+channelfilename = [newdatadir 'clean_SUA_sup_50']; 
+data_file = load(channelfilename);
+locsfilename = [newdatadir 'clean_SUA_locs'];
+all_locsdSUA = load(locsfilename);
+xabs = -199:1300;
+
+nyq = 500;
+channum = 1: length(data_file.clean_high_SUA);
+mean_filtered_dSUA = struct();
+suas_trials = struct();
+up_dist = nan(1, length(channum));
+all_locsdSUA_filtered = nan(1,length(channum));
+
+
+for i = channum  
+    if ~isempty(data_file.clean_high_SUA(i).namelist)
+    trialidx = 1:length(data_file.clean_high_SUA(i).namelist(1,:));
+   
+   filtered_dSUA = data_file.clean_high_SUA(i).namelist;
+  % filtered_dSUA = filtered_dSUA(:,~all(isnan(filtered_dSUA))); % for nan - cols
+   
+   
+          %store unaligned data
+    suas_trials(i).unaligned = filtered_dSUA;
+    %determine the first peak location for each trial of a given single
+    %unit
+    all_locsdSUA_trials = all_locsdSUA.peaks_locs(i).locs;
+    locs_peak1 = all_locsdSUA_trials(1, :);
+
+    up_dist_trials= length(xabs)- locs_peak1;
+    max_low_dist_unit = max(locs_peak1(1,:));
+    %create new matrix with the length(max(d)+max(xabs - d))
+    new_dist_unit = max_low_dist_unit + max(up_dist_trials); 
+    fp_locked_trials = nan(new_dist_unit,length(filtered_dSUA(1,:)));
+  
+     clear n
+     for n = 1:length(filtered_dSUA(1,:))
+     
+         lower_unit_bound =max_low_dist_unit-locs_peak1(n)+1;
+         upper_unit_bound =max_low_dist_unit-locs_peak1(n)+length(xabs);
+         fp_locked_trials(lower_unit_bound:upper_unit_bound,n) = filtered_dSUA(:,n);
+        
+             
+             
+     end
+     fp_locked_trials_out = fp_locked_trials(:,~all(isnan(fp_locked_trials))); % for nan - cols
+      %compute the mean single unit activity if more than 10 trials
+      mean_filtered_dSUA(i).mean_unit = mean(fp_locked_trials_out,2); % for nan - cols
+     %get the aligned data if it exists for the unit 
+     suas_trials(i).aligned = fp_locked_trials_out; 
+    end
+    
+   %{
+   figure();
+   x = 1:length(fp_locked_trials_out(:,1));
+   plot(x,fp_locked_trials)
+   hold on
+   plot(x, mean(fp_locked_trials_out,2),'LineWidth',1, 'Color', 'black')
+   
+   %mean_aligned = mean(fp_locked_trials_out,2);
+   %nanmean_aligned = nanmean(fp_locked_trials_out,2);
+   %}
+   %{
+   %%%%%%% align mean single units to first peak %%%%
+      %start finding the peaks at the first non NaN location 
+        for len1 =1:800
+            if ~isnan(mean_filtered_dSUA(i).mean_unit(len1))
+                break
+            end
+        end
+        for len2 = len1+30:1550
+            if ~all(isnan(mean_filtered_dSUA(i).mean_unit)) && mean_filtered_dSUA(i).mean_unit(len2) < mean_filtered_dSUA(i).mean_unit(len2+1)
+
+             locsdSUA_filtered = findpeaks(mean_filtered_dSUA(i).mean_unit(len2:end));
+             
+                if mean_filtered_dSUA(i).mean_unit(locsdSUA_filtered.loc(1)+len2) > 0.6*mean_filtered_dSUA(i).mean_unit(locsdSUA_filtered.loc(2)+len2)
+         %store first peak location 
+             all_locsdSUA_filtered(:,i) = locsdSUA_filtered.loc(1)+len2;
+                else
+             all_locsdSUA_filtered(:,i) = locsdSUA_filtered.loc(2)+len2;    
+            
+                end
+             break 
+            end 
+        end
+
+    %compute the distance between the first peak and the last datapoint and store  
+   %in a matrix
+   up_dist(:,i)= length(mean_filtered_dSUA(i).mean_unit)- all_locsdSUA_filtered(i);
+  %}
+end  
+
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+ pvalfilename = [pvaluesdir 'lmer_results_rawbs.csv'];
+ pvalues = dlmread(pvalfilename, ',', 1,1);
+ 
+%% Plots of peaks with pvalues 
+ 
+channum = 1: length(data_file.clean_high_SUA);
+
+for chan = channum
+    
+     mean_unaligned = mean(suas_trials(chan).unaligned,2);
+     unaligned = suas_trials(chan).unaligned;
+     mean_aligned = mean(suas_trials(chan).aligned,2);
+     aligned = suas_trials(chan).aligned;
+     if ~isempty(mean_unaligned) && ~isempty(mean_aligned)
+    h = figure;
+    xabs = -199:1300;
+    subplot(2, 1,1);
+    plot(xabs, unaligned)
+    hold on
+    plot(xabs, mean_unaligned,'LineWidth',1, 'Color', 'black')
+    hold on
+    plot([0 0], ylim,'k')
+    hold on
+    plot([1150 1150], ylim,'k')
+ 
+   
+    xalign = 1:length(suas_trials(chan).aligned(:,1));
+    subplot(2, 1,2);
+    plot(xalign, aligned)
+    hold on
+    plot(xalign, mean_aligned,'LineWidth',1, 'Color', 'black')
+  
+    ylh = ylabel({'\fontsize{9}Contacts','\fontsize{9}Spike Rate (spikes/s)'});
+   
+ 
+      set(gca, 'linewidth',2)
+      set(gca,'box','off')
+
+    sgtitle({sprintf('%s | %s', num2str(chan), char(layer(chan))), 'responses, p<0.05, associated to adaptation pvalues'}, 'Interpreter', 'none')
+    xlabel('Time from -50ms from stimulus onset (ms)')
+   set(gcf,'Units','inches') 
+   set(gcf,'position',[1 1 8.5 11])
+   % filename = strcat('C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\plots\',strcat( sprintf('final_bscorr_data_aligned_unaligned_P1040P2_1ms_stimonset_unit_%d', chan)));
+   % saveas(gcf, strcat(filename, '.png'));
+     end
+end
+
+%{
+
+%get the max distance between the first peak and the stimulus onset (or
+%first data point)
+ max_low_dist = max(all_locsdSUA_filtered(1,:));
+ %create new matrix with the length(max(d)+max(xabs - d))
+ new_dist = max_low_dist + max(up_dist);
+ clear fp_locked_data norm_fp_locked
+ fp_locked_data = nan(round(new_dist),length(channum));
+ norm_fp_locked = nan(round(new_dist),length(channum));
+ clear n
+ for n = 1:length(channum)
+    
+ lower_bound =round(max_low_dist)-round(all_locsdSUA_filtered(n))+1;
+ upper_bound =round(max_low_dist)-round(all_locsdSUA_filtered(n))+length(mean_filtered_dSUA(n).mean_unit);
+ % need to change this :
+       if ~all(isnan(mean_filtered_dSUA(n).mean_unit))
+    
+ fp_locked_data(lower_bound:upper_bound,n) = mean_filtered_dSUA(n).mean_unit;
+  %{
+ for len = 30:1500
+            if fp_locked_data(len,n) < fp_locked_data(len+1,n)
+locsdSUA_mean = findpeaks(fp_locked_data(len:end, n));
+        break
+            end 
+ end
+ %x1 = -locsdSUA_mean.loc(1)-len+1:length(fp_locked_data(:,n))-locsdSUA_mean.loc(1)-len;
+%}
+ 
+
+  norm_fp_locked(lower_bound:upper_bound,n) = (fp_locked_data(lower_bound:upper_bound,n)-min(fp_locked_data(lower_bound:upper_bound,n)))/(max(fp_locked_data(lower_bound:upper_bound,n))-min(fp_locked_data(lower_bound:upper_bound,n)));
+  
+        else
+       fp_locked_data(:,n) = nan(length(fp_locked_data(:,n)),1);
+       norm_fp_locked(:,n) = nan(length(fp_locked_data(:,n)),1);
+       
+       end
+        
+  
+ end
+
+ fp_locked_data_out = fp_locked_data(:,~all(isnan(fp_locked_data))); % for nan - cols
+ 
+
+figure();
+x = 1:length(fp_locked_data_out(:,1));
+ plot(x, fp_locked_data_out(:,:),'HandleVisibility','off')
+ hold on
+ plot(x, mean(fp_locked_data_out(:,:),2),'LineWidth',1, 'Color', 'black')
+ plot(xlim, [0 0],'k','HandleVisibility','off')
+
+  title('Trial-selected mean single units spiking activity contrast >.5' , 'Interpreter', 'none')
+  
+  %}
\ No newline at end of file
diff --git a/single_units_analysis/former_align_clean_data.m b/single_units_analysis/former_align_clean_data.m
new file mode 100644
index 0000000..ef32bfa
--- /dev/null
+++ b/single_units_analysis/former_align_clean_data.m
@@ -0,0 +1,185 @@
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\su_peaks_02252020_2\all_units\';
+channelfilename = [newdatadir 'clean_SUA_sup_50']; 
+data_file = load(channelfilename);
+
+xabs = -199:1300;
+
+nyq = 500;
+channum = 1: length(data_file.clean_high_SUA);
+mean_filtered_dSUA = struct();
+suas_trials = struct();
+up_dist = nan(1, length(channum));
+all_locsdSUA_filtered = nan(1,length(channum));
+
+
+for i = channum  
+    if ~isempty(data_file.clean_high_SUA(i).namelist)
+    trialidx = 1:length(data_file.clean_high_SUA(i).namelist(1,:));
+   
+   filtered_dSUA = data_file.clean_high_SUA(i).namelist;
+  % filtered_dSUA = filtered_dSUA(:,~all(isnan(filtered_dSUA))); % for nan - cols
+   
+   
+   bsl_peaks = nan(1, length(filtered_dSUA(1,:)));
+   for tr = trialidx
+      for loc = 1:200
+          if filtered_dSUA(loc,tr) < filtered_dSUA(loc+1,tr)
+      bsl_peak_locs = findpeaks(filtered_dSUA(loc:200,tr));
+      bsl_peaks(1,tr) = max(filtered_dSUA(bsl_peak_locs.loc+loc,tr));
+      break
+          end 
+      end
+   end
+   
+   
+   %if length(filtered_dSUA(1,:)) >=10
+       %store unaligned data
+    suas_trials(i).unaligned = filtered_dSUA;
+    %determine the first peak location for each trial of a given single
+    %unit
+    all_locsdSUA_trials = nan(1,length(filtered_dSUA(1,:)));
+    peak1 = nan(1, length(filtered_dSUA(1,:)));
+  
+    for trial = 1:length(filtered_dSUA(1,:))
+         for ln = 1:550
+             if ~all(isnan(filtered_dSUA(:,trial))) && filtered_dSUA(200+ln,trial) < filtered_dSUA(200+ln+1,trial)
+            
+                 locsdSUA_trial = findpeaks(filtered_dSUA(200+ln:1350,trial));
+             
+                if filtered_dSUA(locsdSUA_trial.loc(1)+200+ln,trial) > 0.4*filtered_dSUA(locsdSUA_trial.loc(2)+200+ln)
+         %store first peak location 
+             all_locsdSUA_trials(:,trial) = locsdSUA_trial.loc(1)+200+ln;
+                else
+             all_locsdSUA_trials(:,trial) = locsdSUA_trial.loc(2)+200+ln;    
+            
+                end
+                 %get distribution of peak 1 across trials
+             peak1(1,trial) = filtered_dSUA(all_locsdSUA_trials(1,trial),trial);
+             break 
+             end 
+         end
+         
+    end
+    up_dist_trials= length(xabs)- all_locsdSUA_trials;
+    max_low_dist_unit = max(all_locsdSUA_trials(1,:));
+    %create new matrix with the length(max(d)+max(xabs - d))
+    new_dist_unit = max_low_dist_unit + max(up_dist_trials); 
+    fp_locked_trials = nan(new_dist_unit,length(filtered_dSUA(1,:)));
+   
+    %find peak1 outliers
+    outlier = isoutlier(peak1);
+    %find bsl peaks outliers
+    out_bsl_peaks = isoutlier(bsl_peaks);
+     clear n
+     for n = 1:length(filtered_dSUA(1,:))
+        if ~isnan(all_locsdSUA_trials(n)) && outlier(n) == 0 %&& out_bsl_peaks(n) == 0
+         lower_unit_bound =max_low_dist_unit-all_locsdSUA_trials(n)+1;
+         upper_unit_bound =max_low_dist_unit-all_locsdSUA_trials(n)+length(xabs);
+     
+         fp_locked_trials(lower_unit_bound:upper_unit_bound,n) = filtered_dSUA(:,n);
+        end
+     end
+     fp_locked_trials_out = fp_locked_trials(:,~all(isnan(fp_locked_trials))); % for nan - cols
+      %compute the mean single unit activity if more than 10 trials
+      mean_filtered_dSUA(i).mean_unit = mean(fp_locked_trials_out,2); % for nan - cols
+     %get the aligned data if it exists for the unit 
+     suas_trials(i).aligned = fp_locked_trials_out; 
+  
+  % elseif length(filtered_dSUA(1,:)) <10
+  %    mean_filtered_dSUA(i).mean_unit = nan(length(fp_locked_trials(:,1)),1);  
+  %    suas_trials(i).unaligned = nan(1,1);
+  % end
+  end 
+   %{
+   figure();
+   x = 1:length(fp_locked_trials_out(:,1));
+   plot(x,fp_locked_trials)
+   hold on
+   plot(x, mean(fp_locked_trials_out,2),'LineWidth',1, 'Color', 'black')
+   
+   %mean_aligned = mean(fp_locked_trials_out,2);
+   %nanmean_aligned = nanmean(fp_locked_trials_out,2);
+   %}
+   %{
+   %%%%%%% align mean single units to first peak %%%%
+      %start finding the peaks at the first non NaN location 
+        for len1 =1:800
+            if ~isnan(mean_filtered_dSUA(i).mean_unit(len1))
+                break
+            end
+        end
+        for len2 = len1+30:1550
+            if ~all(isnan(mean_filtered_dSUA(i).mean_unit)) && mean_filtered_dSUA(i).mean_unit(len2) < mean_filtered_dSUA(i).mean_unit(len2+1)
+
+             locsdSUA_filtered = findpeaks(mean_filtered_dSUA(i).mean_unit(len2:end));
+             
+                if mean_filtered_dSUA(i).mean_unit(locsdSUA_filtered.loc(1)+len2) > 0.6*mean_filtered_dSUA(i).mean_unit(locsdSUA_filtered.loc(2)+len2)
+         %store first peak location 
+             all_locsdSUA_filtered(:,i) = locsdSUA_filtered.loc(1)+len2;
+                else
+             all_locsdSUA_filtered(:,i) = locsdSUA_filtered.loc(2)+len2;    
+            
+                end
+             break 
+            end 
+        end
+
+    %compute the distance between the first peak and the last datapoint and store  
+   %in a matrix
+   up_dist(:,i)= length(mean_filtered_dSUA(i).mean_unit)- all_locsdSUA_filtered(i);
+  %}
+end
+
+pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+ pvalfilename = [pvaluesdir 'lmer_results_rawbs.csv'];
+ pvalues = dlmread(pvalfilename, ',', 1,1);
+ 
+%% Plots of peaks with pvalues 
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = []; 
+channum = 1: length(data_file.clean_high_SUA);
+
+for chan = 1:15
+ 
+     mean_unaligned = mean(suas_trials(chan).unaligned,2);
+     unaligned = suas_trials(chan).unaligned;
+     mean_aligned = mean(suas_trials(chan).aligned,2);
+     aligned = suas_trials(chan).aligned;
+     if ~isempty(mean_unaligned) && ~isempty(mean_aligned)
+    h = figure;
+    xabs = -199:1300;
+    subplot(2, 1,1);
+    plot(xabs, unaligned)
+    hold on
+    plot(xabs, mean_unaligned,'LineWidth',1, 'Color', 'black')
+    hold on
+    plot([0 0], ylim,'k')
+    hold on
+    plot([1150 1150], ylim,'k')
+ 
+   
+    xalign = 1:length(suas_trials(chan).aligned(:,1));
+    subplot(2, 1,2);
+    plot(xalign, aligned)
+    hold on
+    plot(xalign, mean_aligned,'LineWidth',1, 'Color', 'black')
+  
+    ylh = ylabel({'\fontsize{9}Contacts','\fontsize{9}Spike Rate (spikes/s)'});
+   
+ 
+      set(gca, 'linewidth',2)
+      set(gca,'box','off')
+
+    sgtitle({sprintf('%s | %s', num2str(chan), char(layer(chan))), 'responses, p<0.05, associated to adaptation pvalues'}, 'Interpreter', 'none')
+    xlabel('Time from -50ms from stimulus onset (ms)')
+   set(gcf,'Units','inches') 
+   set(gcf,'position',[1 1 8.5 11])
+   % filename = strcat('C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\plots\',strcat( sprintf('bscorr_data_aligned_unaligned_P1040P2_1ms_stimonset_unit_%d_exclude4peaks', chan)));
+    %saveas(gcf, strcat(filename, '.png'));
+     end
+end
diff --git a/single_units_analysis/get_clean_norm_peaks_and_norm_data.m b/single_units_analysis/get_clean_norm_peaks_and_norm_data.m
new file mode 100644
index 0000000..97023cb
--- /dev/null
+++ b/single_units_analysis/get_clean_norm_peaks_and_norm_data.m
@@ -0,0 +1,147 @@
+%after saving the data with new_data_set.m, we isolate the peaks for the
+%new analysis of the refined data
+%we also save the data we want to plot (only the clean data)
+
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+channelfilename = [newdatadir 'refined_dataset']; 
+data_file = load(channelfilename);
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ %% find peaks in order to analyze them on R and fit a LMER
+ 
+channum = 1: length(data_file.new_data);
+xabs = -199:1300;
+nyq = 500;
+mean_filtered_dSUA = struct();
+
+data_peaks = struct();
+clean_high_SUA = struct();
+ 
+ for i = channum
+ 
+   filename = [data_file.new_data(i).channel_data.filename, f{2}];
+   filename = erase(filename, '.mat');
+
+   blankcontrast = data_file.new_data(i).channel_data.contrast ==  0 & data_file.new_data(i).channel_data.fixedc ==  0;
+   highcontrast = data_file.new_data(i).channel_data.contrast >=  0.5 & data_file.new_data(i).channel_data.fixedc ==  0; 
+    
+    trialidx = 1:length(data_file.new_data(i).channel_data.sdftr_chan(1,:));
+    raw_bs = nan(length(xabs), length(trialidx));
+    filtered_dSUA = nan(length(xabs), length(trialidx));
+    all_norm_lpdSUA= nan(length(xabs),length(trialidx));
+       
+    
+    powerstim = nan(length(trialidx),1025);
+    freqstim = nan(length(trialidx),1025);
+    fourhzpowerstim =nan(length(trialidx),1);
+    mean_wnd1 = nan(1,length(trialidx));
+    
+     all_norm_pks = nan(4,length(data_file.new_data(i).channel_data.sdftr_chan(1,highcontrast)));
+     
+     for tridx = trialidx
+      
+    all_data = data_file.new_data(i).channel_data.sdftr_chan(401:1900,tridx);
+ 
+    raw_bs(:,tridx) = all_data(1:end)- mean(all_data(1:200));
+   
+    lpc       = 4.5; %low pass cutoff
+    lWn       = lpc/nyq;
+    [bwb,bwa] = butter(4,lWn,'low');
+    lpdSUA      = filtfilt(bwb,bwa, raw_bs(:,tridx));
+  
+
+    filtered_dSUA(:,tridx) = lpdSUA;
+    all_norm_lpdSUA(:,tridx) = (lpdSUA - min(lpdSUA))/(max(lpdSUA)- min(lpdSUA));
+    mean_wnd1(tridx) = mean(lpdSUA(231:480));
+    
+    %%% power
+
+          
+    [powerstim(tridx,:), freqstim(tridx,:)] = calcFFT(all_data(200:1350));
+        
+    %find the index of the frequency vector closest to 4hz and point to the
+    %power value of this index for every trial, and store the value in
+    %fourhzpower
+    [val,index] = min(abs(4-freqstim(tridx,:)));
+    fourhzpowerstim(tridx,1) = powerstim(tridx,index);
+    
+     end
+     
+   %%%%%%%%%%% %reject trials below the 95%CI in the blank condition %%%%%%
+   %power related variables
+   power0 = fourhzpowerstim(blankcontrast);
+   power5 = fourhzpowerstim(highcontrast);
+  
+   %spiking activity related variables
+
+   filtered_dSUA_high = filtered_dSUA(:, highcontrast);
+   all_norm_lpdSUA_high = all_norm_lpdSUA(:, highcontrast);
+   %first peak location related variables
+   mean_wnd1_5 = mean_wnd1(highcontrast); 
+   sua_bsl =  mean(filtered_dSUA_high(1:200,:),1);
+    
+   for tr = 1:length(power5)
+       if power5(tr) > mean(power0)+1.96*std(power0) && mean_wnd1_5(tr) > mean(sua_bsl) + 1.96*std(sua_bsl)
+     
+           all_norm_lpdSUA_high(:,tr) = all_norm_lpdSUA_high(:,tr);
+           
+       else
+   
+           all_norm_lpdSUA_high(:,tr) = nan(length(all_norm_lpdSUA_high(:,tr)),1);
+       end
+   end
+   all_norm_lpdSUA_high = all_norm_lpdSUA_high(:,~all(isnan(all_norm_lpdSUA_high))); % for nan - cols
+   
+   
+   if length(all_norm_lpdSUA_high(1,:)) >=10
+   
+    %determine the first peak location for each trial of a given single
+    %unit
+    
+  
+    for trial = 1:length(all_norm_lpdSUA_high(1,:))
+        
+         for ln = 30:550
+             if all_norm_lpdSUA_high(200+ln,trial) < all_norm_lpdSUA_high(200+ln+1,trial)
+              
+            
+             locsdSUA_trial = findpeaks(all_norm_lpdSUA_high(200+ln:1350,trial));
+         
+             break 
+             end 
+         end
+         
+        if length(locsdSUA_trial.loc) >= 4
+             %adjust location to the first data point of lpsu (+ln),
+    lpsulocs = locsdSUA_trial.loc(1:4) + 200+ ln;
+    all_norm_pks(:,trial) = all_norm_lpdSUA_high(lpsulocs(1:4), trial);
+  
+        end 
+    end
+   namelist(1,1:length(sprintf('chan_%d',i))) = sprintf('chan_%d',i);
+   clean_high_SUA(i).namelist =  all_norm_lpdSUA_high;
+    
+    elseif length(all_norm_lpdSUA_high(1,:)) <10  || length(locsdSUA_trial.loc) <= 4
+       all_norm_pks(:,:) = nan(length(all_norm_pks(:,1)),length(all_norm_pks(1,:)));
+       clean_high_SUA(i).namelist =  nan(size(all_norm_lpdSUA_high));
+   end
+   
+ namelist(1,1:length(sprintf('chan_%d',i))) = sprintf('chan_%d',i);
+ data_peaks(i).namelist = all_norm_pks(:,~all(isnan(all_norm_pks)));
+ all_norm_pks = all_norm_pks(:,~all(isnan(all_norm_pks)));
+channelfilename = [newdatadir 'norm_su_peaks\' filename];
+%save(strcat(channelfilename, '.mat'), 'all_norm_pks');
+ end
+  allfilename = [newdatadir 'all_norm_data_peaks'];
+ %save(strcat(allfilename, '.mat'), 'data_peaks');
+ allfilename = [newdatadir 'clean_norm_SUA_sup_50'];
+ %save(strcat(allfilename, '.mat'), 'clean_high_SUA');
+ 
\ No newline at end of file
diff --git a/single_units_analysis/get_clean_peaks_and_data.m b/single_units_analysis/get_clean_peaks_and_data.m
new file mode 100644
index 0000000..067bd29
--- /dev/null
+++ b/single_units_analysis/get_clean_peaks_and_data.m
@@ -0,0 +1,236 @@
+%after saving the data with new_data_set.m, we isolate the peaks for the
+%new analysis of the refined data
+%we also save the data we want to plot (only the clean data)
+
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+channelfilename = [newdatadir 'refined_dataset']; 
+data_file = load(channelfilename);
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ %% find peaks in order to analyze them on R and fit a LMER
+ 
+channum = 1: length(data_file.new_data);
+xabs = -199:1300;
+nyq = 500;
+mean_filtered_dSUA = struct();
+
+clean_high_SUA = struct();
+data_peaks = struct();
+peaks_locs = struct();
+ 
+ for i = channum
+ 
+   filename = [data_file.new_data(i).channel_data.filename, f{2}];
+   filename = erase(filename, '.mat');
+ 
+   blankcontrast = data_file.new_data(i).channel_data.contrast ==  0 & data_file.new_data(i).channel_data.fixedc ==  0;
+   highcontrast = data_file.new_data(i).channel_data.contrast >=  0.5 & data_file.new_data(i).channel_data.fixedc ==  0; 
+    
+    trialidx = 1:length(data_file.new_data(i).channel_data.sdftr_chan(1,:));
+    raw_bs = nan(length(xabs), length(trialidx));
+    filtered_dSUA = nan(length(xabs), length(trialidx));
+    all_norm_lpdSUA= nan(length(xabs),length(trialidx));
+       
+    
+    powerstim = nan(length(trialidx),1025);
+    freqstim = nan(length(trialidx),1025);
+    fourhzpowerstim =nan(length(trialidx),1);
+    bsl = nan(1, length(trialidx));
+    mean_wnd1 = nan(1,length(trialidx));
+    
+    all_pks = nan(4,length(data_file.new_data(i).channel_data.sdftr_chan(1,highcontrast)));
+     
+     for tridx = trialidx
+      
+    all_data = data_file.new_data(i).channel_data.sdftr_chan(401:1900,tridx);
+    raw_bs(:,tridx) = all_data(1:end)- mean(all_data(1:200));
+   
+    lpc       = 4.5; %low pass cutoff
+    lWn       = lpc/nyq;
+    [bwb,bwa] = butter(4,lWn,'low');
+    lpdSUA      = filtfilt(bwb,bwa, raw_bs(:,tridx));
+  
+
+    filtered_dSUA(:,tridx) = lpdSUA;
+    %all_norm_lpdSUA(:,tridx) = (lpdSUA - min(lpdSUA))/(max(lpdSUA)- min(lpdSUA));
+    mean_wnd1(tridx) = mean(lpdSUA(201:480));
+    
+    %%% power
+
+          
+    [powerstim(tridx,:), freqstim(tridx,:)] = calcFFT(all_data(200:1350));
+        
+    %find the index of the frequency vector closest to 4hz and point to the
+    %power value of this index for every trial, and store the value in
+    %fourhzpower
+    [val,index] = min(abs(4-freqstim(tridx,:)));
+    fourhzpowerstim(tridx,1) = powerstim(tridx,index);
+    
+     end
+     
+  %%%%%%%%%%% %reject trials below the 95%CI in the blank condition %%%%%%
+   %power related variables
+   power0 = fourhzpowerstim(blankcontrast);
+   power5 = fourhzpowerstim(highcontrast);
+   
+   %spiking activity related variables
+   mean_wnd1_5 =mean_wnd1(highcontrast);
+   filtered_dSUA_high = filtered_dSUA(:, highcontrast);
+   %all_norm_lpdSUA_high = all_norm_lpdSUA(:, highcontrast);
+   %first peak location related variables 
+   sua_bsl =  mean(filtered_dSUA_high(1:200,:),1);
+  
+   for tr = 1:length(power5)
+      if mean_wnd1_5(tr) > mean(sua_bsl)+1.96*std(sua_bsl)  && power5(tr) > mean(power0)+1.96*std(power0) %
+     
+           filtered_dSUA_high(:,tr) = filtered_dSUA_high(:,tr);
+           
+       else
+   
+           filtered_dSUA_high(:,tr) = nan(length(filtered_dSUA_high(:,tr)),1);
+       end
+   end
+  
+      %determine the first peak location for each trial of a given single
+    %unit
+    all_locsdSUA_trials = nan(6,length(filtered_dSUA_high(1,:)));
+    %peak1 = nan(1, length(filtered_dSUA_high(1,:)));
+    for trial = 1:length(filtered_dSUA_high(1,:))
+        
+         for ln = 1:550
+             if filtered_dSUA_high(200+ln,trial) < filtered_dSUA_high(200+ln+1,trial) && ~all(isnan(filtered_dSUA_high(:,trial)))
+  
+             locsdSUA_trial = findpeaks(filtered_dSUA_high(200+ln:1350,trial));
+             %if peak1 is too small, peak2 becomes peak1
+                     if filtered_dSUA_high(locsdSUA_trial.loc(1)+200+ln,trial) >= 0.4*filtered_dSUA_high(locsdSUA_trial.loc(2)+200+ln)
+             %store first peak location 
+                     all_locsdSUA_trials(1:length(locsdSUA_trial.loc),trial) = locsdSUA_trial.loc(1:end)+200+ln;
+                     else %if filtered_dSUA_high(locsdSUA_trial.loc(1)+200+ln,trial) < 0.4*filtered_dSUA_high(locsdSUA_trial.loc(2)+200+ln) && filtered_dSUA_high(locsdSUA_trial.loc(2)+200+ln,trial) >= 0.4*filtered_dSUA_high(locsdSUA_trial.loc(3)+200+ln)
+                     all_locsdSUA_trials(1:length(locsdSUA_trial.loc(2:end)),trial) = locsdSUA_trial.loc(2:end)+200+ln; 
+                     
+                     end
+                     
+              break 
+             end 
+         end
+         
+        if nnz(~isnan(all_locsdSUA_trials(:,trial))) >= 4 && ~all(isnan(all_locsdSUA_trials(:,trial))) 
+             %adjust location to the first data point of lpsu (+ln),
+    %lpsulocs = ~isnan(all_locsdSUA_trials(:,trial));
+    all_pks(:,trial) = filtered_dSUA_high(all_locsdSUA_trials(1:4,trial), trial);
+    filtered_dSUA_high(:,trial) = filtered_dSUA_high(:,trial); 
+    all_locsdSUA_trials(:,trial) = all_locsdSUA_trials(:,trial);
+   % peak1(1,trial) = filtered_dSUA_high(all_locsdSUA_trials(1,trial), trial);
+        else
+    filtered_dSUA_high(:,trial) = nan(length(filtered_dSUA_high(:,trial)),1);
+    all_locsdSUA_trials(:,trial) = nan(size(all_locsdSUA_trials(:,trial)));
+        end 
+    end
+    %%% reject outlier peaks and the corresponding trials in
+    %%% filtered_dSUA_high
+   
+   
+   %reject if there is a peak 1 outlier, if the max peak value in the
+   %baseline is an outlier
+   
+    % First find peaks before stimulus onset 
+    
+    bsl_peaks = nan(1, length(filtered_dSUA_high(1,:)));
+    clear tr
+    for tr = 1:length(filtered_dSUA_high(1,:))
+        
+         for loc = 1:200
+            if filtered_dSUA_high(loc,tr) < filtered_dSUA_high(loc+1,tr) && ~all(isnan(filtered_dSUA_high(:,tr)))
+           bsl_peak_locs = findpeaks(filtered_dSUA_high(loc:200,tr));
+           bsl_peaks(1,tr) = max(filtered_dSUA_high(bsl_peak_locs.loc+loc,tr));
+          break
+            end 
+         end
+    end
+    
+    out_bsl_peaks = isoutlier(bsl_peaks);
+    
+    p1outliers = isoutlier(all_pks(1,:));
+    clear tr
+    for tr = 1:length(filtered_dSUA_high(1,:))
+        %exclude trials
+        if p1outliers(tr) == 0 && ~all(isnan(all_pks(:,tr))) && out_bsl_peaks(tr) ==0 
+            
+    filtered_dSUA_high(:,tr) = filtered_dSUA_high(:, tr);
+    all_pks(:, tr) = all_pks(:,tr);
+    all_locsdSUA_trials(:,tr) = all_locsdSUA_trials(:,tr);
+   % peak1(tr) = peak1(tr);
+        else 
+    filtered_dSUA_high(:,tr) = nan(length(filtered_dSUA_high(:,tr)),1);
+    all_pks(:,tr) = nan(length(all_pks(:,tr)),1);
+    all_locsdSUA_trials(:,tr) = nan(size(all_locsdSUA_trials(:,tr)));
+   % peak1(tr) = nan(1,1);
+        end
+    end
+    filtered_dSUA_high = filtered_dSUA_high(:,~all(isnan(filtered_dSUA_high))); % for nan - cols
+   all_locsdSUA_trials =  all_locsdSUA_trials(:,~all(isnan(all_locsdSUA_trials)));
+   all_pks = all_pks(:, ~all(isnan(all_pks)));
+   %peak1 = peak1(:,~isnan(peak1));
+   
+ %%% perform the peak1 outlier exclusion a second time
+ %{
+    p1outliers = isoutlier(all_pks(1,:));
+    clear tr
+    for tr = 1:length(filtered_dSUA_high(1,:))
+        %exclude trials
+        if  ~all(isnan(all_pks(:,tr))) && p1outliers(tr) == 0 %&& out_bsl_peaks(tr) ==0 
+            
+    filtered_dSUA_high(:,tr) = filtered_dSUA_high(:, tr);
+    all_pks(:, tr) = all_pks(:,tr);
+    all_locsdSUA_trials(:,tr) = all_locsdSUA_trials(:,tr);
+        else 
+    filtered_dSUA_high(:,tr) = nan(length(filtered_dSUA_high(:,tr)),1);
+    all_pks(:,tr) = nan(length(all_pks(:,tr)),1);
+    all_locsdSUA_trials(:,tr) = nan(size(all_locsdSUA_trials(:,tr)));
+        end
+    end
+ 
+    
+    filtered_dSUA_high = filtered_dSUA_high(:,~all(isnan(filtered_dSUA_high))); % for nan - cols
+   all_locsdSUA_trials =  all_locsdSUA_trials(:,~all(isnan(all_locsdSUA_trials)));
+   all_pks = all_pks(:, ~all(isnan(all_pks)));
+ %}
+   %%%%
+  
+   if length(filtered_dSUA_high(1,:)) >=10
+   clean_high_SUA(i).namelist =  filtered_dSUA_high;
+   peaks_locs(i).locs = all_locsdSUA_trials; 
+   elseif length(filtered_dSUA_high(1,:)) <10  
+    all_pks(:,:) = nan(length(all_pks(:,1)),length(all_pks(1,:)));
+    clean_high_SUA(i).namelist =  [];
+    peaks_locs(i).locs = []; 
+   end
+   
+ data_peaks(i).namelist = all_pks(:,~all(isnan(all_pks)));
+ all_pks = all_pks(:,~all(isnan(all_pks)));
+channelfilename = [newdatadir 'su_peaks_02262020_2\' filename];
+save(strcat(channelfilename, '.mat'), 'all_pks');
+ end
+ allfilename = [newdatadir 'su_peaks_02262020_2\all_units\all_data_peaks'];
+ save(strcat(allfilename, '.mat'), 'data_peaks');
+ allfilename = [newdatadir 'su_peaks_02262020_2\all_units\clean_SUA_sup_50'];
+ save(strcat(allfilename, '.mat'), 'clean_high_SUA');
+ allfilename = [newdatadir 'su_peaks_02262020_2\all_units\clean_SUA_locs'];
+ save(strcat(allfilename, '.mat'), 'peaks_locs');
+ 
+   cnt =0;
+   for i =1:length(data_peaks)
+       if ~isempty(data_peaks(i).namelist)
+           cnt = cnt+1;
+       end
+        
+   end
+   
\ No newline at end of file
diff --git a/single_units_analysis/lmer_results_anal.asv b/single_units_analysis/lmer_results_anal.asv
new file mode 100644
index 0000000..92d28de
--- /dev/null
+++ b/single_units_analysis/lmer_results_anal.asv
@@ -0,0 +1,182 @@
+%this script was written after get_clean_peaks_and_data.m in order to
+%analyze the lmer results and plot the clean data.
+
+%loading the clean data
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\su_peaks_02262020_2\all_units\';
+channelfilename = [newdatadir 'clean_SUA_sup_50']; 
+data_file = load(channelfilename);
+
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ 
+ pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+ pvalfilename = [pvaluesdir 'lmer_results_02252020_2.csv'];
+ pvalues = dlmread(pvalfilename, ',', 1,1);
+ 
+%% Rough plots of peaks with pvalues (not very representative, as mean unit activity)
+  channum = 1: length(data_file.clean_high_SUA);
+
+ %for n = 1:3 
+for chan = 1:12:length(channum)
+h = figure;
+xabs = -199:1300;
+idx = [1 3 5 7 9 11 2 4 6 8 10 12];
+nyq = 500;
+all_mean_data = nan(length(xabs), length(1:12));
+clear i ;
+ for i = 1:12
+ 
+   mean_data = mean(data_file.clean_high_SUA(chan+i-1).namelist(1:1500,:),2);
+   
+  
+   lpc       = 4.5; %low pass cutoff
+   lWn       = lpc/nyq;
+   [bwb,bwa] = butter(4,lWn,'low');
+   
+  if ~all(isnan(mean_data))
+   lpsu      = filtfilt(bwb,bwa, mean_data);
+   
+    sp = subplot(length(1:6), 2, idx(i));
+    plot(xabs, lpsu)
+    hold on
+    plot([0 0], ylim,'k')
+    hold on
+    plot([1150 1150], ylim,'k')
+
+    if i == length(6)/2
+        ylh = ylabel({'\fontsize{9}Contacts','\fontsize{9}Spike Rate (spikes/s)'});
+    end
+   if i < 6 || (i >= 7)&&(i < 12)
+        set(sp, 'XTick', [])
+   end
+      ylabelh = text(max(xabs), mean(lpsu,1), strcat(num2str(chan+i-1),' | ', layer(chan+i-1)),'HorizontalAlignment','left','FontName', 'Arial','FontSize', 10);
+   for npeak = 1:4
+         for len = 231:480 %from 200 + 30 (lgn response onset) 
+            if lpsu(len) < lpsu(len+1)
+   locs = findpeaks(lpsu(len:1450));
+        break
+            end
+         end
+         
+         if length(locs.loc) >= 4
+             %adjust location to the first data point of lpsu (+len), then adjust
+             %to xabs (-200)
+   xlocation = locs.loc(npeak)+len-200;
+         end 
+            
+   text(xlocation, mean(lpsu,1), strcat(num2str(sprintf('%.2f', pvalues(chan+i-1,npeak)))),'HorizontalAlignment','center','FontName', 'Arial','FontSize', 7);
+   end 
+
+ end
+ end
+      set(gca, 'linewidth',2)
+      set(gca,'box','off')
+
+    sgtitle({f{2}, 'all good responses, p<0.05, associated to adaptation pvalues'}, 'Interpreter', 'none')
+    xlabel('Time from -50ms from stimulus onset (ms)')
+   set(gcf,'Units','inches') 
+   set(gcf,'position',[1 1 8.5 11])
+    %filename = strcat('C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\plots\',strcat(f{2}, sprintf('x_%d_better_raw_data_peakspvalues_2dec', chan+i-1)));
+    %saveas(gcf, strcat(filename, '.png'));
+
+end
+
+%%
+
+ %% compute proportion of significant adaptation per peak and proportion of neurons adapting for a certain amount of 
+%peak from peak 2 to 4
+
+channeldir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\su_peaks_02262020_2\all_units\';
+pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+pvalfilename = [pvaluesdir 'lmer_results_02262020_2.csv'];
+pvalues = dlmread(pvalfilename, ',', 1,1);
+
+peakvals = load([channeldir 'all_data_peaks']);
+
+ 
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+
+ layer_idx = find(strcmp(layer, 'K'));
+
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+
+ all_locs = nan(4,length(layer_idx));
+ all_pks = nan(4,length(layer_idx));
+ all_mean_data = nan(4, length(layer_idx));
+
+  
+ cntpk2 = 0;
+ cntpk3 = 0;
+ cntpk4 = 0;
+ cntpk2pk3 = 0;
+ cntpk2pk3pk4 = 0;
+ cntpk3pk4 =0;
+ cntincpk4 =0;
+ 
+ cntnspk3 =0;
+ cntnspk4=0;
+ 
+  for nunit = 1:length(layer_idx)
+ mean_data = nanmean(peakvals.data_peaks(layer_idx(nunit)).namelist,2);
+   
+   all_mean_data(:,nunit) = mean_data;
+
+     if all_mean_data(2,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),2) < .05
+         cntpk2 = cntpk2 +1;
+     end
+     if all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05
+         cntpk3 = cntpk3 +1;
+     end
+     if all_mean_data(4,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05
+         cntpk4 = cntpk4 +1;
+     end
+     if all_mean_data(2,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),2) < .05 && ...
+             all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05
+         cntpk2pk3 = cntpk2pk3 +1;
+     end
+     if all_mean_data(2,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),2) < .05 && ...
+             all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05 ...
+             && all_mean_data(4,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05
+         cntpk2pk3pk4 = cntpk2pk3pk4 +1;
+     end
+       if all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05 ...
+             && all_mean_data(4,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05 
+         cntpk3pk4 = cntpk3pk4 +1;
+       end
+     if all_mean_data(4,nunit) > all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05
+         cntincpk4 = cntincpk4 +1;
+     end
+     if pvalues(layer_idx(nunit),3) > .05
+         cntnspk3 = cntnspk3 +1;
+     end
+    if pvalues(layer_idx(nunit),4) > .05
+         cntnspk4 = cntnspk4 +1;
+     end
+  end
+ all_mean_data = all_mean_data(:, ~all(isnan(all_mean_data)));
+ 
+ percentpk2 = cntpk2*100/length(all_mean_data(1,:));
+ percentpk3 = cntpk3*100/length(all_mean_data(1,:)); 
+ percentpk4 = cntpk4*100/length(all_mean_data(1,:));
+ percentpk2pk3 = cntpk2pk3*100/length(all_mean_data(1,:));
+ percentpk2pk3pk4 = cntpk2pk3pk4*100/length(all_mean_data(1,:));
+ percentpk3pk4 = cntpk3pk4*100/length(all_mean_data(1,:));
+ percentincpk4 = cntincpk4*100/length(all_mean_data(1,:));
+ percentnspk3 = cntnspk3*100/length(all_mean_data(1,:));
+ percentnspk4 = cntnspk4*100/length(all_mean_data(1,:));
+ 
+ ncells =length(all_mean_data(1,:));
\ No newline at end of file
diff --git a/single_units_analysis/lmer_results_anal.m b/single_units_analysis/lmer_results_anal.m
new file mode 100644
index 0000000..8d095ce
--- /dev/null
+++ b/single_units_analysis/lmer_results_anal.m
@@ -0,0 +1,182 @@
+%this script was written after get_clean_peaks_and_data.m in order to
+%analyze the lmer results and plot the clean data.
+
+%loading the clean data
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\su_peaks_02262020_2\all_units\';
+channelfilename = [newdatadir 'clean_SUA_sup_50']; 
+data_file = load(channelfilename);
+
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ 
+ pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+ pvalfilename = [pvaluesdir 'lmer_results_02252020_2.csv'];
+ pvalues = dlmread(pvalfilename, ',', 1,1);
+ 
+%% Rough plots of peaks with pvalues (not very representative, as mean unit activity)
+  channum = 1: length(data_file.clean_high_SUA);
+
+ %for n = 1:3 
+for chan = 1:12:length(channum)
+h = figure;
+xabs = -199:1300;
+idx = [1 3 5 7 9 11 2 4 6 8 10 12];
+nyq = 500;
+all_mean_data = nan(length(xabs), length(1:12));
+clear i ;
+ for i = 1:12
+ 
+   mean_data = mean(data_file.clean_high_SUA(chan+i-1).namelist(1:1500,:),2);
+   
+  
+   lpc       = 4.5; %low pass cutoff
+   lWn       = lpc/nyq;
+   [bwb,bwa] = butter(4,lWn,'low');
+   
+  if ~all(isnan(mean_data))
+   lpsu      = filtfilt(bwb,bwa, mean_data);
+   
+    sp = subplot(length(1:6), 2, idx(i));
+    plot(xabs, lpsu)
+    hold on
+    plot([0 0], ylim,'k')
+    hold on
+    plot([1150 1150], ylim,'k')
+
+    if i == length(6)/2
+        ylh = ylabel({'\fontsize{9}Contacts','\fontsize{9}Spike Rate (spikes/s)'});
+    end
+   if i < 6 || (i >= 7)&&(i < 12)
+        set(sp, 'XTick', [])
+   end
+      ylabelh = text(max(xabs), mean(lpsu,1), strcat(num2str(chan+i-1),' | ', layer(chan+i-1)),'HorizontalAlignment','left','FontName', 'Arial','FontSize', 10);
+   for npeak = 1:4
+         for len = 231:480 %from 200 + 30 (lgn response onset) 
+            if lpsu(len) < lpsu(len+1)
+   locs = findpeaks(lpsu(len:1450));
+        break
+            end
+         end
+         
+         if length(locs.loc) >= 4
+             %adjust location to the first data point of lpsu (+len), then adjust
+             %to xabs (-200)
+   xlocation = locs.loc(npeak)+len-200;
+         end 
+            
+   text(xlocation, mean(lpsu,1), strcat(num2str(sprintf('%.2f', pvalues(chan+i-1,npeak)))),'HorizontalAlignment','center','FontName', 'Arial','FontSize', 7);
+   end 
+
+ end
+ end
+      set(gca, 'linewidth',2)
+      set(gca,'box','off')
+
+    sgtitle({f{2}, 'all good responses, p<0.05, associated to adaptation pvalues'}, 'Interpreter', 'none')
+    xlabel('Time from -50ms from stimulus onset (ms)')
+   set(gcf,'Units','inches') 
+   set(gcf,'position',[1 1 8.5 11])
+    %filename = strcat('C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\plots\',strcat(f{2}, sprintf('x_%d_better_raw_data_peakspvalues_2dec', chan+i-1)));
+    %saveas(gcf, strcat(filename, '.png'));
+
+end
+
+%%
+
+ %% compute proportion of significant adaptation per peak and proportion of neurons adapting for a certain amount of 
+%peak from peak 2 to 4
+
+channeldir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\su_peaks_02262020_2\all_units\';
+pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+pvalfilename = [pvaluesdir 'lmer_results_02262020_2.csv'];
+pvalues = dlmread(pvalfilename, ',', 1,1);
+
+peakvals = load([channeldir 'all_data_peaks']);
+
+ 
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+
+ layer_idx = find(strcmp(layer, 'M'));
+
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+
+ all_locs = nan(4,length(layer_idx));
+ all_pks = nan(4,length(layer_idx));
+ all_mean_data = nan(4, length(layer_idx));
+
+  
+ cntpk2 = 0;
+ cntpk3 = 0;
+ cntpk4 = 0;
+ cntpk2pk3 = 0;
+ cntpk2pk3pk4 = 0;
+ cntpk3pk4 =0;
+ cntincpk4 =0;
+ 
+ cntnspk3 =0;
+ cntnspk4=0;
+ 
+  for nunit = 1:length(layer_idx)
+ mean_data = nanmean(peakvals.data_peaks(layer_idx(nunit)).namelist,2);
+   
+   all_mean_data(:,nunit) = mean_data;
+
+     if all_mean_data(2,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),2) < .05
+         cntpk2 = cntpk2 +1;
+     end
+     if all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05
+         cntpk3 = cntpk3 +1;
+     end
+     if all_mean_data(4,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05
+         cntpk4 = cntpk4 +1;
+     end
+     if all_mean_data(2,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),2) < .05 && ...
+             all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05
+         cntpk2pk3 = cntpk2pk3 +1;
+     end
+     if all_mean_data(2,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),2) < .05 && ...
+             all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05 ...
+             && all_mean_data(4,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05
+         cntpk2pk3pk4 = cntpk2pk3pk4 +1;
+     end
+       if all_mean_data(3,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),3) < .05 ...
+             && all_mean_data(4,nunit) < all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05 
+         cntpk3pk4 = cntpk3pk4 +1;
+       end
+     if all_mean_data(4,nunit) > all_mean_data(1,nunit) && pvalues(layer_idx(nunit),4) < .05
+         cntincpk4 = cntincpk4 +1;
+     end
+     if pvalues(layer_idx(nunit),3) > .05
+         cntnspk3 = cntnspk3 +1;
+     end
+    if pvalues(layer_idx(nunit),4) > .05
+         cntnspk4 = cntnspk4 +1;
+     end
+  end
+ all_mean_data = all_mean_data(:, ~all(isnan(all_mean_data)));
+ 
+ percentpk2 = cntpk2*100/length(all_mean_data(1,:));
+ percentpk3 = cntpk3*100/length(all_mean_data(1,:)); 
+ percentpk4 = cntpk4*100/length(all_mean_data(1,:));
+ percentpk2pk3 = cntpk2pk3*100/length(all_mean_data(1,:));
+ percentpk2pk3pk4 = cntpk2pk3pk4*100/length(all_mean_data(1,:));
+ percentpk3pk4 = cntpk3pk4*100/length(all_mean_data(1,:));
+ percentincpk4 = cntincpk4*100/length(all_mean_data(1,:));
+ percentnspk3 = cntnspk3*100/length(all_mean_data(1,:));
+ percentnspk4 = cntnspk4*100/length(all_mean_data(1,:));
+ 
+ ncells =length(all_mean_data(1,:));
\ No newline at end of file
diff --git a/single_units_analysis/new_data_set.m b/single_units_analysis/new_data_set.m
new file mode 100644
index 0000000..196172d
--- /dev/null
+++ b/single_units_analysis/new_data_set.m
@@ -0,0 +1,15 @@
+gooddatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\';
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+
+
+oldchannelfilename = [gooddatadir 'good_single_units_data_4bmpmore']; 
+old_data_file = load(oldchannelfilename);
+
+%Reject 160517, K cell, as no well triggered (1)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+new_data = old_data_file.good_data([2:45,47:54,56:end]);
+
+allfilename = [newdatadir 'refined_dataset'];
+save(strcat(allfilename, '.mat'), 'new_data', '-v7.3');
+ 
\ No newline at end of file
diff --git a/single_units_analysis/plot_units_peak_aligned.m b/single_units_analysis/plot_units_peak_aligned.m
new file mode 100644
index 0000000..04a5ac7
--- /dev/null
+++ b/single_units_analysis/plot_units_peak_aligned.m
@@ -0,0 +1,90 @@
+%this script is intended to plot mean single units activity with both
+%unaligned data and aligned data (to the first peak)
+
+%loading the clean data
+newdatadir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\';
+channelfilename = [newdatadir 'clean_norm_SUA_sup_50']; 
+data_file = load(channelfilename);
+
+
+%exclude 160517, (first unit, left empty, it is a K neuron)
+%Reject 180806 p1 uclust17, M cell, as doesn't seem well triggered (46)
+%Reject 181207 (B) uclust22, M cell, as doesn't seem well triggered (55)
+ layer = {'K','M','P','K','K','K','M','P','P','','M','M','','','M','','','P','','M','','M','M','','P','M','','P', ...
+'P','','','K','P','M','M','M','P','','P','K','P','P','','P','P','M','','P','M','P','M','P','','P','M','M','P','','M','M','P','M', ...
+'','','M','M','M','P','M','M','M','M','P','P'};
+layer([1,46,55]) = [];
+ f = {'DE0_NDE50','DE50_NDE0','DE50_NDE50'};
+ 
+ 
+ pvaluesdir = 'C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\lmer_results_peaks\';
+ pvalfilename = [pvaluesdir 'lmer_results_norm.csv'];
+ pvalues = dlmread(pvalfilename, ',', 1,1);
+ 
+%% Rough plots of peaks with pvalues (not very representative, as mean unit activity)
+  channum = 1: length(data_file.clean_high_SUA);
+
+ %for n = 1:3 
+for chan = 1:12:length(channum)
+h = figure;
+xabs = -199:1300;
+idx = [1 3 5 7 9 11 2 4 6 8 10 12];
+nyq = 500;
+all_mean_data = nan(length(xabs), length(1:12));
+clear i ;
+ for i = 1:12
+ 
+   mean_data = mean(data_file.clean_high_SUA(chan+i-1).namelist(1:1500,:),2);
+   
+  
+   lpc       = 4.5; %low pass cutoff
+   lWn       = lpc/nyq;
+   [bwb,bwa] = butter(4,lWn,'low');
+   
+  if ~all(isnan(mean_data))
+   lpsu      = filtfilt(bwb,bwa, mean_data);
+   
+    sp = subplot(length(1:6), 2, idx(i));
+    plot(xabs, lpsu)
+    hold on
+    plot([0 0], ylim,'k')
+    hold on
+    plot([1150 1150], ylim,'k')
+
+    if i == length(6)/2
+        ylh = ylabel({'\fontsize{9}Contacts','\fontsize{9}Spike Rate (spikes/s)'});
+    end
+   if i < 6 || (i >= 7)&&(i < 12)
+        set(sp, 'XTick', [])
+   end
+      ylabelh = text(max(xabs), mean(lpsu,1), strcat(num2str(chan+i-1),' | ', layer(chan+i-1)),'HorizontalAlignment','left','FontName', 'Arial','FontSize', 10);
+   for npeak = 1:4
+         for len = 231:480 %from 200 + 30 (lgn response onset) 
+            if lpsu(len) < lpsu(len+1)
+   locs = findpeaks(lpsu(len:1450));
+        break
+            end
+         end
+         
+         if length(locs.loc) >= 4
+             %adjust location to the first data point of lpsu (+len), then adjust
+             %to xabs (-200)
+   xlocation = locs.loc(npeak)+len-200;
+         end 
+            
+   text(xlocation, mean(lpsu,1), strcat(num2str(sprintf('%.2f', pvalues(chan+i-1,npeak)))),'HorizontalAlignment','center','FontName', 'Arial','FontSize', 7);
+   end 
+
+ end
+ end
+      set(gca, 'linewidth',2)
+      set(gca,'box','off')
+
+    sgtitle({f{2}, 'all good responses, p<0.05, associated to adaptation pvalues'}, 'Interpreter', 'none')
+    xlabel('Time from -50ms from stimulus onset (ms)')
+   set(gcf,'Units','inches') 
+   set(gcf,'position',[1 1 8.5 11])
+    %filename = strcat('C:\Users\maier\Documents\LGN_data\single_units\inverted_power_channels\good_single_units_data_4bumps_more\new_peak_alignment_anal\plots\',strcat(f{2}, sprintf('x_%d_better_raw_data_peakspvalues_2dec', chan+i-1)));
+    %saveas(gcf, strcat(filename, '.png'));
+
+end
\ No newline at end of file