%Finding RPE neurons
clear all;
load ('R_2R.mat');

Suc=R_2R.SucN; %neurons with greater firing for sucrose

%find cue excited neurons
Cue=strcmp('Cue', R_2R.Erefnames);
RD1=strcmp('RD1', R_2R.Erefnames);
RD2=strcmp('RD2', R_2R.Erefnames);
CueEx=R_2R.Ev(Cue).RespDir==1;

%neurons that are cue excited and sucrose preferring
Sel=CueEx&Suc;

sucrose=[1  0.6  0.1];
maltodextrin=[.9  0.3  .9];
water=[0.00  0.75  0.75];
total=[0.3  0.1  0.8];
exc=[0 113/255 188/255];
inh=[240/255 0 50/255];
cued=[0.5  0.1  0.8];
uncued=[0.3  0.7  0.1];

%divide neurons up by region
NAneurons=strcmp(R_2R.Ninfo(:,3),'NA');
VPneurons=strcmp(R_2R.Ninfo(:,3),'VP');

%% cue response
Xaxis=[-0.5 1];
inttime=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
paramtime=R_2R.Param.Tm(inttime);

%NAc
subplot(2,4,1);
hold on;

%plot sucrose preferring and cue excited neurons to cue
psth1=nanmean(R_2R.Ev(Cue).PSTHz(Sel&NAneurons,inttime),1); 
sem1=nanste(R_2R.Ev(Cue).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

plot(paramtime,psth1,'Color',total,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],total,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-1 10],':','color','k','linewidth',0.75);
axis([Xaxis(1) Xaxis(2) -0.6 6.5]);
ylabel('Mean firing (z-score)');

%VP
subplot(2,4,5);
hold on;

%plot sucrose preferring and cue excited neurons to cue
psth1=nanmean(R_2R.Ev(Cue).PSTHz(Sel&VPneurons,inttime),1); 
sem1=nanste(R_2R.Ev(Cue).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

plot(paramtime,psth1,'Color',total,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],total,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-1 10],':','color','k','linewidth',0.75);
axis([Xaxis(1) Xaxis(2) -0.6 8]);
ylabel('Mean firing (z-score)');
xlabel('Seconds from cue');

%% reward response
Xaxis=[-1 3.5];
inttime=find(R_2R.Param.Tm>=Xaxis(1) & R_2R.Param.Tm<=Xaxis(2));
paramtime=R_2R.Param.Tm(inttime);

%NAc
subplot(2,4,2:4);
hold on;

%plot sucrose preferring and cue excited neurons to sucrose
psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel&NAneurons,inttime),1); 
sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

%plot sucrose preferring and cue extied neurons to malt
psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel&NAneurons,inttime),1); 
sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel&NAneurons,inttime),1); %calculate standard error of the mean
up2=psth2+sem2;
down2=psth2-sem2;

plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-1 10],':','color','k','linewidth',0.75);
axis([Xaxis(1) Xaxis(2) -0.8 6.5]);
legend('Suc trials','Mal trials');
set(gca,'ytick',[]);
title(['NAc (n=' num2str(sum(Sel&NAneurons)) ', ' num2str(round(100*sum(Sel&NAneurons)/sum(NAneurons))) '%)']);

%VP
subplot(2,4,6:8);
hold on;

%plot sucrose preferring and cue excited neurons to sucrose
psth1=nanmean(R_2R.Ev(RD1).PSTHz(Sel&VPneurons,inttime),1); 
sem1=nanste(R_2R.Ev(RD1).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
up1=psth1+sem1;
down1=psth1-sem1;

%plot sucrose preferring and cue extied neurons to malt
psth2=nanmean(R_2R.Ev(RD2).PSTHz(Sel&VPneurons,inttime),1); 
sem2=nanste(R_2R.Ev(RD2).PSTHz(Sel&VPneurons,inttime),1); %calculate standard error of the mean
up2=psth2+sem2;
down2=psth2-sem2;

plot(paramtime,psth1,'Color',sucrose,'linewidth',1);
plot(paramtime,psth2,'Color',maltodextrin,'linewidth',1);
patch([paramtime,paramtime(end:-1:1)],[up1,down1(end:-1:1)],sucrose,'EdgeColor','none');alpha(0.5);
patch([paramtime,paramtime(end:-1:1)],[up2,down2(end:-1:1)],maltodextrin,'EdgeColor','none');alpha(0.5);
plot([-2 5],[0 0],':','color','k','linewidth',0.75);
plot([0 0],[-1 10],':','color','k','linewidth',0.75);
axis([Xaxis(1) Xaxis(2) -0.8 8]);
legend('Suc trials','Mal trials');
set(gca,'ytick',[]);
xlabel('Seconds from RD');
title(['VP (n=' num2str(sum(Sel&VPneurons)) ', ' num2str(round(100*sum(Sel&VPneurons)/sum(VPneurons))) '%)']);