
function [Data,Results]=RunTimeCourse(GLMinput,Data,GLMname,DataB)

save_fig = 1; % if you want to save the figures
save_peak= 1; % if you want to save the time courses to later test for significance

addpath(fullfile(Data.Basepath,'SEA','timecourse_Scripts','fmt')) % where we have the gammapdf function (from fmt toolbox)

mfig_path=fullfile(Data.folders.timecoursesResults,'MatlabFigs');
text_path.p=fullfile(Data.folders.timecoursesResults,'DataFigures');

timefolder = fullfile(Data.folders.timecourses);   %folder where the extracted time courses are stored

subs = 1:length(Data.SubjectIDs); % participants for who to run the timecourse analys§is
text_path.ID = Data.SubjectIDs;

% list of different regions for analyses
regions={'32d_bilateral','accumbensBilateral','accumbensLeft','accumbensRight','dACCact','dACCLarge_bilateral','dACCNils','dACCpropsValActDesA2','dlPFC','lFP','lOFC','perigenual','preSMA','RCZa_bilateral','TPJRogier'};
regions={'32d_bilateral','accumbensBilateral','dACCLarge_bilateral','dlPFC','lFP','perigenual','preSMA','RCZa_bilateral','TPJRogier'};

regions={'cope3_cluster12_vmPFC_1','cope3_cluster12_vmPFC_6','cope3_cluster10_midCingulate_1','cope29_cluster7_dACC_1','cope29_cluster7_dACC_2','cope3_neg_cluster9_dACC_1','cope3_neg_cluster9_dACC_2',...
   'cope3_neg_cluster9_dACC_3','cope3_neg_cluster9_dACC_5','cope3_neg_cluster5_LFPl_1','cope8_cluster4_preSMA_1'...
   ,'cope8_cluster6_TPJ_3','cope3_neg_cluster11_IPS_5','cope29_cluster2_RpSTS_1','cope29_cluster7_dACC_1'}; %

regions={'cope3_cluster12_vmPFC_1','cope3_cluster12_vmPFC_6','cope3_cluster10_midCingulate_1','cope29_cluster7_dACC_1','cope29_cluster7_dACC_2','cope3_neg_cluster8_dACC_1','cope3_neg_cluster8_dACC_2',...
   'cope3_neg_cluster8_dACC_3','cope3_neg_cluster8_dACC_5','cope3_neg_cluster7_FPl_1','cope8_cluster6_preSMA_1'...
   ,'cope8_cluster6_TPJ_3','cope3_neg_cluster11_IPS_5','cope29_cluster2_RpSTS_1','cope29_cluster7_dACC_1','cope3_vmPFC','cope3_wholevmPFC'}; %


regions={'cope3_cluster12_vmPFC_1','cope3_cluster12_vmPFC_6','cope3_cluster10_midCingulate_1','cope29_cluster7_dACC_1','cope29_cluster7_dACC_2','cope3_neg_cluster8_dACC_1','cope3_neg_cluster8_dACC_2',...
   'cope3_neg_cluster8_dACC_3','cope3_neg_cluster8_dACC_5','cope3_neg_cluster7_FPl_1','cope8_cluster6_preSMA_1'...
   ,'cope29_cluster7_dACC_1','cope3_vmPFC','cope3_wholevmPFC'}; % ,'cope8_cluster6_TPJ_3','cope3_neg_cluster11_IPS_5','cope29_cluster2_RpSTS_1'


regions={'cope5_cluster1_preSMA','cope5_cluster2_rMFG','cope18_cluster1_M1','cope18_cluster2_TPJ','cope19_neg_cluster1_preSMA','cope19_neg_cluster2_rIPLE',...
    'cope19_neg_cluster3_precuneus','cope19_neg_cluster4_lIPLE','cope19_neg_cluster5_antinsula','cope21_cluster1_rIPLA','cope21_cluster2_RCZp',...
    'cope21_cluster3_lIPLB','cope21_cluster4_lOccipital','cope21_cluster5_rOccipital','cope21_cluster6_lOccipital','cope21_cluster7_rV5','cope22_cluster1_area25',...
    'cope22_cluster2_precuneus','cope30_cluster1_IPLE','cope30_neg_cluster1_11m','cope21_neg_cluster5_middle_temporal_gyrus','cope21_neg_dACC'};


regions={'cope5_cluster1_preSMA','cope5_cluster2_rMFG','cope18_cluster2_TPJ','cope19_neg_cluster1_preSMA','cope19_neg_cluster2_rIPLE',...
    'cope19_neg_cluster3_precuneus','cope19_neg_cluster4_lIPLE','cope19_neg_cluster5_antinsula','cope21_cluster1_rIPLA','cope21_cluster2_RCZp',...
    'cope21_cluster3_lIPLB','cope22_cluster1_area25','cope22_cluster2_precuneus','cope30_cluster1_IPLE','cope30_neg_cluster1_11m',... 
    'cope3_cluster12_vmPFC_1','cope3_cluster12_vmPFC_6','cope3_cluster10_midCingulate_1','cope29_cluster7_dACC_1','cope3_neg_cluster8_dACC_1','cope3_neg_cluster7_FPl_1','cope8_cluster6_preSMA_1'...
   ,'cope29_cluster7_dACC_1','cope3_vmPFC','cope3_wholevmPFC','cope21_neg_dACC'};



regions={'cope5_cluster1_preSMA','cope19_neg_cluster1_preSMA','cope29_cluster7_dACC_1','cope3_neg_cluster8_dACC_1','cope3_neg_cluster7_FPl_1','cope8_cluster6_preSMA_1'...
   ,'cope29_cluster7_dACC_1','cope21_neg_dACC'};

regions={'cope5_cluster1_preSMA','cope19_neg_cluster1_preSMA','cope29_cluster7_dACC_1','cope3_neg_cluster8_dACC_1','cope3_neg_cluster7_FPl_1','cope8_cluster6_preSMA_1'...
   ,'cope29_cluster7_dACC_1','cope21_neg_dACC','cope3_cluster12_vmPFC_1','cope3_cluster12_vmPFC_6','cope3_neg_cluster7_FPl_1','cope3_vmPFC','cope3_wholevmPFC'};


regions={'cope8_cluster6_preSMA_1','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1','cope3_vmPFC','cope3_wholevmPFC'};

regions={'cope5_cluster1_preSMA','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1','cope3_vmPFC','cope3_wholevmPFC'};

regions={'cope5_cluster1_preSMA','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1',...
  'cope8_cluster7_midParietal_1','cope8_cluster5_IPS_1','cope8_cluster6_TPJ_1','avgmask'};

%cope8_cluster5_IPS_1 Similar to dACC for value stuff
% cope8_cluster7_midParietal_1 strong signals for self but nothing for
% other choice (negative value signals)
% Nothing in TPJ

%'cope8_cluster7_midParietal_2','cope8_cluster7_midParietal_1','cope8_cluster6_TPJ_6','cope8_cluster6_TPJ_5','cope8_cluster6_TPJ_4',...
%    'cope8_cluster6_TPJ_3','cope8_cluster6_TPJ_2','cope8_cluster6_TPJ_1','cope8_cluster6_preSMA_2','cope8_cluster6_preSMA_1','cope8_cluster5_IPS_5',...
%    ]'cope8_cluster5_IPS_4','cope8_cluster5_IPS_3'.'cope8_cluster5_IPS_2','cope8_cluster5_IPS_1'



% XX=SUBJECTNR
% YY=COPENUMBER
% ZZ=CLUSTERNR
% A=SUBCLUSTERNR
% XX_copeYY_clusterZZ_AnatomicalLabel_A 
% XX_copeYY_neg_clusterZZ_AnatomicalLabel_A 

% % %
% % COPE3 positive, cluster 12, vmPFC
% % 4 40 -12 Good Center for vmPFC  c1
% % 2 18 -6 subcallosal area 25     c6
% COPE3 positive, cluster 10 midCingulate
% % -2 -8 42 fine midCingulate c1
% COPE29 positive, cluster 7, dACC
% % 4 30 40 rostral cingulate zonish upper bank of (para)cingulate gyrus c1
% % 8 24 50 preSMA or 8M c2
% COPE3 negative, cluster 9, dACC
% % -4 22 48 preSMA or 8M c1
% % 54 20 32 lateraldlPFC c2
% % -6 22 44 ACCposteriorrostalcingulatezone c3
% % -6 6 66 preSMA c5
% COPE3 negative, cluster 5, lFP
% % -40 60 6 goodlateralFP c1
% COPE8 positive, cluster 4, preSMA
% % -12 14 64 bit lateral c1


% cope3_neg_cluster8_PPC/IPS_2 [-30 -52 40]
% %cope8_cluster6_TPJ_3 [54 -56 16] for TPJp.
% %cope3_neg_cluster11_IPS_5 
% %cope29_cluster2_RpSTS_1 [58 -32 -4] pSTS.
% cope10_cluster10_IPS_1 [-32 -54 36] IPS; 
% cope12_cluster11_IPS_1  [44 -60 50] IPS; 
% %cope29_cluster7_dACC_1 [4 30 40] mid cingulate; 



%cope 29: [38 50 18] 9/46V; [-36 30 28] 9/46V; 




%% Load and upsample the timecourses
out=UpsampleTimecourse(timefolder,GLMinput,regions,Data);
SelfDecTC              = out.SelfDec;
OtherDecTC                = out.OtherDec;
SamplingTC              = out.Sampling;
OutcomeTC                = out.Outcome;

settings(1).pseudosamples  = out.SelfDec.pseudosamples; %how data was upsampled, needed for plotting
settings(1).window         = out.SelfDec.window;
settings(1).TR             = out.SelfDec.TR;
settings(1).upsample       = out.SelfDec.upsample;

settings(2).pseudosamples  = out.OtherDec.pseudosamples; %how data was upsampled, needed for plotting
settings(2).window         = out.OtherDec.window;
settings(2).TR             = out.OtherDec.TR;
settings(2).upsample       = out.OtherDec.upsample;
settings(2).TCShift        = out.OtherDec.TCShift;

settings(3).pseudosamples  = out.Sampling.pseudosamples; %how data was upsampled, needed for plotting
settings(3).window         = out.Sampling.window;
settings(3).TR             = out.Sampling.TR;
settings(3).upsample       = out.Sampling.upsample;

settings(4).pseudosamples  = out.Outcome.pseudosamples; %how data was upsampled, needed for plotting
settings(4).window         = out.Outcome.window;
settings(4).TR             = out.Outcome.TR;
settings(4).upsample       = out.Outcome.upsample;


      
%% Run the regressions for the timecourse
    % For decision phase (time-locked to choice possible)
    
    tlabs=[];
    for l=(length(GLMinput.SelfDec.Labels)+1):-1:2;tlabs{l}=GLMinput.SelfDec.Labels{l-1};end
    tlabs{1}='SDCon';
    GLMinput.SelfDec.Labels=tlabs;
    outT=TimeCourseRunRegression(GLMinput.SelfDec,SelfDecTC,Data);
    Results(1).design(1).time  = 'SelfDecision';
    Results(1).design(1).inp   = GLMinput.SelfDec.DesignMatrix;
    Results(1).design(1).out   = outT;
    
    Results(1).design(1).labels= outT.labels;
    Results(1).design(1).designName=GLMname;
    
        tlabs=[];
    for l=(length(GLMinput.OtherDec.Labels)+1):-1:2;tlabs{l}=GLMinput.OtherDec.Labels{l-1};end
    tlabs{1}='ODCon';
    GLMinput.OtherDec.Labels=tlabs;
    
    outT=TimeCourseRunRegression(GLMinput.OtherDec,OtherDecTC,Data);
    Results(2).design(1).time  = 'OtherDecision';
    Results(2).design(1).inp   = GLMinput.OtherDec.DesignMatrix;
    Results(2).design(1).out   = outT;
     
    Results(2).design(1).labels= outT.labels;
     Results(2).design(1).designName=GLMname;
    
     
    tlabs=[];
    for l=(length(GLMinput.Sampling.Labels)+1):-1:2;tlabs{l}=GLMinput.Sampling.Labels{l-1};end
    tlabs{1}='SamplingCon';
    GLMinput.Sampling.Labels=tlabs;
    outT=TimeCourseRunRegression(GLMinput.Sampling,SamplingTC,Data);
    Results(3).design(1).time  = 'Sampling';
    Results(3).design(1).inp   = GLMinput.Sampling.DesignMatrix;
    Results(3).design(1).out   = outT;     
    Results(3).design(1).labels= outT.labels;    
    Results(3).design(1).designName=GLMname;
    
    tlabs=[];
    for l=(length(GLMinput.Outcome.Labels)+1):-1:2;tlabs{l}=GLMinput.Outcome.Labels{l-1};end
    tlabs{1}='OutcomeCon';
    GLMinput.Outcome.Labels=tlabs;
    outT=TimeCourseRunRegression(GLMinput.Outcome,OutcomeTC,Data);
    Results(4).design(1).time  = 'Outcome';
    Results(4).design(1).inp   = GLMinput.Outcome.DesignMatrix;
    Results(4).design(1).out   = outT;
    Results(4).design(1).labels=outT.labels;
    Results(4).design(1).designName=GLMname;
   
%% Plot the results of time course analyses - all areas together
% Each figure: one area, subplots: different regressors of the GLM
%plotlist.ploteverything = true; % if this is true, plots all designs and all regressors for each region and does not check the other choices in plotlist
clear plotl
plotl.designList=1;
plotl.plotWhat = 'SortByRegressor';
PlotTC(Results(1),plotl,settings(1)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting

%% Sorth by region
%plotl.plotWhat = 'SortByRegion';
%PlotTC(Results(1),plotl,settings(1));

plotl.plotWhat = 'SortByRegressor';
PlotTC(Results(2),plotl,settings(2)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
%% Sorth by region
%plotl.plotWhat = 'SortByRegion';
%PlotTC(Results(2),plotl,settings(2));

plotl.plotWhat = 'SortByRegressor';
PlotTC(Results(3),plotl,settings(3)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
% %% Sorth by region
% %plotl.plotWhat = 'SortByRegion';
% %PlotTC(Results(3),plotl,settings(3));
% 
plotl.plotWhat = 'SortByRegressor';
PlotTC(Results(4),plotl,settings(4)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
%% Sorth by region
%plotl.plotWhat = 'SortByRegion';
%PlotTC(Results(4),plotl,settings(4));


%keyboard;
%close all
%% Plot the results of the timecourse analyses for the paper: only selected areas and regressors

clear plotl
plotl.plotWhat='TimeCoursePaperFigures';
%
%regions={'cope5_cluster1_preSMA','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1',...
%  'cope8_cluster7_midParietal_1','cope8_cluster5_IPS_1','cope8_cluster6_TPJ_1'};

 plotl.DesRegionsRegsFigure =[1 1 2 1 2 1 1;   % columns: design, regions, regressor, figure ID, number of rows, number of columns, which subplot
     1 1 3 1 2 1 1;
     1 2 2 1 2 1 2;    
     1 2 3 1 2 1 2;
     1 3 2 2 2 1 1;    
     1 3 3 2 2 1 1;
     1 4 2 2 2 1 2;    
     1 4 3 2 2 1 2;
     1 5 2 3 2 1 1;    
     1 5 3 3 2 1 1;
     1 6 2 3 2 1 2;    
     1 6 3 3 2 1 2;];

 
plotl.Colors=[51 102 0;
    204 0 0;
    51 102 0;
    204 0 0;
    51 102 0;
    204 0 0;
     51 102 0;
    204 0 0; 
    51 102 0;
    204 0 0; 
    51 102 0;
    204 0 0];
% 
plotl.TimeLimit =settings(1).window*ones(size(plotl.DesRegionsRegsFigure,1),3);
plotl.ylims=[-0.16 0.06;
    -0.16 0.06;
    -0.06 0.15;
    -0.06 0.15;
    -0.21 0.06;
    -0.21 0.06;
     -0.21 0.06;
    -0.21 0.06;
    -.12 .12;
    -.12 .12;
    -.12 .12;
    -.12 .12];
% plot the figures
[figs,plotData,timesSec]=PlotTC(Results(1),plotl,settings(1));
% save the figures and text files
% SaveTxtFiles(plotData{1},timesSec,text_path,'TC_dACC_prospV') % store the text file with the individual subject data

saveas(figs(1), fullfile(mfig_path,['TC_preSMA_vmPFC_SD.jpg'])); %save the figure
 saveas(figs(2), fullfile(mfig_path,['TC_FPl_dACC_SD.jpg'])); %save the figure
 saveas(figs(3), fullfile(mfig_path,['TC_midParietal_IPS_SD.jpg'])); %save the figure
savefig(figs(1), fullfile(mfig_path,['TC_preSMA_vmPFC_SD'])); %save the figure
 savefig(figs(2), fullfile(mfig_path,['TC_FPl_dACC_SD'])); %save the figure
 savefig(figs(3), fullfile(mfig_path,['TC_midParietal_IPS_SD'])); %save the figure

% keyboard
 
 clear plotl
 plotl.plotWhat='TimeCoursePaperFigures';
 
   plotl.DesRegionsRegsFigure =[1 1 2 1 2 1 1;   % columns: design, regions, regressor, figure ID, number of rows, number of columns
     1 1 3 1 2 1 1;
     1 2 2 1 2 1 2;    
     1 2 3 1 2 1 2;
     1 3 2 2 2 1 1;    
     1 3 3 2 2 1 1;
     1 4 2 2 2 1 2;    
     1 4 3 2 2 1 2;
     1 5 2 3 2 1 1;    
     1 5 3 3 2 1 1;
     1 6 2 3 2 1 2;    
     1 6 3 3 2 1 2;];
% 
plotl.Colors=[51 102 0;
    204 0 0;
    51 102 0;
    204 0 0;
    51 102 0;
    204 0 0;
     51 102 0;
    204 0 0; 
    51 102 0;
    204 0 0; 
    51 102 0;
    204 0 0];
% 
plotl.TimeLimit =settings(1).window*ones(size(plotl.DesRegionsRegsFigure,1),3);
plotl.ylims=[-0.1 0.05;
    -0.1 0.05;
    0 0.1;
    0 0.1;
    -0.1 0.06;
    -0.1 0.06;
     -0.1 0.06;
    -0.1 0.06;
    -.1 .02;
    -.1 .02;
    -.1 .02;
    -.1 .02];
% plot the figures
[figs,plotData,timesSec]=PlotTC(Results(2),plotl,settings(2));
% save the figures and text files
% SaveTxtFiles(plotData{2},timesSec,text_path,'TC_dlPFC_prospV')
% SaveTxtFiles(plotData{3},timesSec,text_path,'TC_dACC_myoVal')
% SaveTxtFiles(plotData{4},timesSec,text_path,'TC_dlPFC_myoVal')
% SaveTxtFiles(plotData{5},timesSec,text_path,'TC_dACC_difficulty')
% SaveTxtFiles(plotData{6},timesSec,text_path,'TC_RCZ_prospV')
% SaveTxtFiles(plotData{7},timesSec,text_path,'TC_RCZ_prospVWithFutureChangeCovar')

regions={'cope5_cluster1_preSMA','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1',...
  'cope8_cluster7_midParietal_1','cope8_cluster5_IPS_1','cope8_cluster6_TPJ_1','avgmask'};

 saveas(figs(1), fullfile(mfig_path,['TC_preSMA_vmPFC_OD.jpg'])); %save the figure
 saveas(figs(2), fullfile(mfig_path,['TC_FPl_dACC_OD.jpg'])); %save the figure
 saveas(figs(3), fullfile(mfig_path,['TC_midParietal_IPS_OD.jpg'])); %save the figure
savefig(figs(1), fullfile(mfig_path,['TC_preSMA_vmPFC_OD'])); %save the figure
 savefig(figs(2), fullfile(mfig_path,['TC_FPl_dACC_OD'])); %save the figure
 savefig(figs(3), fullfile(mfig_path,['TC_midParietal_IPS_OD'])); %save the figure

 
clear plotl
 plotl.plotWhat='TimeCoursePaperFigures';
 
   plotl.DesRegionsRegsFigure =[1 1 2 1 2 3;   % columns: design, regions, regressor, figure ID, number of rows, number of columns
    1 4 2 1 2 3;
    1 5 2 1 2 3;
    1 6 2 1 2 3;
    1 7 2 1 2 3;];
% 
plotl.Colors=[0 204 0;
    0 204 0;
    0 204 0;
    0 204 0;
    0 204 0];
% 
plotl.TimeLimit =settings(1).window*ones(size(plotl.DesRegionsRegsFigure,1),3);
plotl.ylims=[-.04 .04;
    -.04 .04;
    -.04 .04;
    -.04 .04;
    -.04 .04];
% plot the figures
[figs,plotData,timesSec]=PlotTC(Results(3),plotl,settings(3));
 
saveas(figs(1), fullfile(mfig_path,['TC_SamplingOcc_SharedDiff.jpg'])); %save the figure
savefig(figs(1), fullfile(mfig_path,['TC_SamplingOcc_SharedDiff'])); %save the figure

 
 
% 
% %% Testing timecourse significance using leave one subject out peak selection
% %keyboard;
% %close all
% ides=5;

 
ParameterToTest=[2 3];
RegionsToTest=[1 2 3 4 5 6 7];
BDesignNr=18; behavregressors=[2 1];

%keyboard
peakDetectLimit = [4 settings(1).window];
% % HRF settings - as in other paper, do not modify (unless there is a good reason)
sigma1=3;%old here:3;
my1=7;   % old here 6
alpha1=my1^2/sigma1^2;
beta1=my1/sigma1^2;
hrf_delay=0;
% Code the HRF
t=0:(settings(1).TR /settings(1).upsample):(settings(1).window-hrf_delay);
t=t(1:end-1);
hrf = gammapdf(alpha1,beta1,t);
simplehrf_peak_point = find(hrf==max(hrf));
 % delayhrf=[zeros(1,ceil(hrf_delay/(settings(1).TR  /settings(1).upsample)))  hrf];
% %delayhrf=[zeros(1,pseudosamples-length(hrf))  hrf];
% hrf_peak_time = find(delayhrf==max(delayhrf));
CPIndex=[]; clear PaperHRFTestsFigData;
for R=1:length(RegionsToTest)
    PaperHRFTestsFig(R)=figure('name',['Region:' Results(2).design(1).out.regions{RegionsToTest(R)} ' HRF t-tests'],'color',[1 1 1]);hold on;subplot(1,1,1)
    PaperHRFTestsFigData{R}{1}=[];
    PaperHRFTestsLabels=[];
    CPIndex(R)=0;
   PaperHRFTestsFig_ps=[];
   PaperHRFTestsFig_stats{R}{1}={};
   PaperHRFTestsFigDataTC{R}{1}=[];
end

 for R=1:length(RegionsToTest)
    Ttestfigs(R)=figure('name',['Region:' Results(1).design(1).out.regions{RegionsToTest(R)} ' peak and HRF t-tests'],'color',[1 1 1]);hold on;
    %CorrBfigs(R)=figure('name',['Region:' Results(1).design(1).out.regions{RegionsToTest(R)} ' peak and HRF behavioural corrs']);hold on;

    for P=1:length(ParameterToTest)
        TCData=squeeze(Results(1).design(1).out.betas(:,RegionsToTest(R),:,ParameterToTest(P)));
        ttestLabels{P}=Results(1).design(1).out.labels{ParameterToTest(P)};
        ttestregions{R}=Results(1).design(1).out.regions{RegionsToTest(R)};
        % First Find the Peak
        dataT.v=TCData;
        [ParameterPeaks(:,P,R),ParameterPeakpoint(:,P,R)] = ComputePeakValue(dataT,peakDetectLimit,settings(1));
        title([ ttestregions{R} '  ' ttestLabels{P} ])
        PeakAdjust=floor(settings(1).pseudosamples/settings(1).window*peakDetectLimit(1))-1+ParameterPeakpoint(:,P,R)-simplehrf_peak_point;
        for adj=1:length(PeakAdjust)
            if PeakAdjust(adj)>0
        adjustedhrfs(adj,:)=[zeros(1,PeakAdjust(adj))  hrf(1:end-PeakAdjust(adj))];
            elseif PeakAdjust(adj)<0
                    adjustedhrfs(adj,:)=[hrf(1-PeakAdjust(adj):end) zeros(1,-PeakAdjust(adj))];    
            else
              adjustedhrfs(adj,:)=  hrf;
            end
        end
        adjustedConvolvedSignal(:,R,P)= sum(TCData.*adjustedhrfs,2);
        [H,CadjPValue,CI,STATS]=ttest(adjustedConvolvedSignal(:,R,P));
        ConvadjPValue(R,P)=CadjPValue;%PeakStats(P,R).STATS=STATS;
        nonParaadjConvPValue(R,P)= signrank(adjustedConvolvedSignal(:,R,P));
        [h,ParameterPeakPValues(R,P),ci,stats]=ttest(ParameterPeaks(:,P,R));
        figure(Ttestfigs(R))      
        subplot(2,2,P);hold on;title([ ttestregions{R} '  ' ttestLabels{P} ])
        bar(1,mean(adjustedConvolvedSignal(:,R,P)));errorbar(1,mean(adjustedConvolvedSignal(:,R,P)),std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5))
        bar(2,mean(ParameterPeaks(:,P,R)));errorbar(2,mean(ParameterPeaks(:,P,R)),std(ParameterPeaks(:,P,R))./(length(ParameterPeaks(:,P,R)).^0.5))
        
        figure(PaperHRFTestsFig(R));hold on
        bar(CPIndex(R)+1,mean(adjustedConvolvedSignal(:,R,P)));
        errorbar(CPIndex(R)+1,mean(adjustedConvolvedSignal(:,R,P)),...
            std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5),'.k','Linewidth',4)
        [h,PaperHRFTestsFig_ps(R,CPIndex(R)+1),ci,PaperHRFTestsFig_stats{R}{CPIndex(R)+1}]=ttest(adjustedConvolvedSignal(:,R,P));
        PaperHRFTestsLabels{CPIndex(R)+1}=ttestLabels{P};
        PaperHRFTestsFigData{R}{CPIndex(R)+1}=adjustedConvolvedSignal(:,R,P);
        PaperHRFTestsFigDataTC{R}{CPIndex(R)+1}=TCData;
        significant=nan ;
        supersignificant=nan ;
        superdupersignificant=nan ;
        
        significant(PaperHRFTestsFig_ps(R,CPIndex(R)+1)<0.05)=true;
        supersignificant(PaperHRFTestsFig_ps(R,CPIndex(R)+1)<0.01)=true;
        superdupersignificant(PaperHRFTestsFig_ps(R,CPIndex(R)+1)<0.001)=true;
        
        
        plot(CPIndex(R)+1.1.*significant,[mean(adjustedConvolvedSignal(:,R,P))+(std(adjustedConvolvedSignal(:,R,P))./length(adjustedConvolvedSignal(:,R,P)).^.5)./2].*significant,'*k','markersize',15)
        plot(CPIndex(R)+1.2.*supersignificant,[mean(adjustedConvolvedSignal(:,R,P))+(std(adjustedConvolvedSignal(:,R,P))./length(adjustedConvolvedSignal(:,R,P)).^.5)./2].*supersignificant,'*k','markersize',15)
        plot(CPIndex(R)+1.3.*superdupersignificant,[mean(adjustedConvolvedSignal(:,R,P))+(std(adjustedConvolvedSignal(:,R,P))./length(adjustedConvolvedSignal(:,R,P)).^.5)./2].*superdupersignificant,'*k','markersize',15)

        HRF_Figstats.RegionName{R}=Results(1).design(1).out.regions{RegionsToTest(R)};
        HRF_Figstats.Region{R}.ParameterLabels{CPIndex(R)+1}=ttestLabels{P};
        HRF_Figstats.Region{R}.pvalue(CPIndex(R)+1)=PaperHRFTestsFig_ps(R,CPIndex(R)+1);         
        HRF_Figstats.Region{R}.tvalue(CPIndex(R)+1)=PaperHRFTestsFig_stats{R}{CPIndex(R)+1}.tstat;        
        HRF_Figstats.Region{R}.df(CPIndex(R)+1)=PaperHRFTestsFig_stats{R}{CPIndex(R)+1}.df;

        CPIndex(R)=CPIndex(R)+1;

        %between-subject correlation with behaviour
%         figure(CorrBfigs(R));
%         subplot(2,2,P);hold on;
%         plot(adjustedConvolvedSignal(:,R,P),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:),'.b','Markersize',5);
%         
%         [temp1,temp2]=corrcoef(adjustedConvolvedSignal(:,R,P),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:));
%         rvalcon=temp1(end,1);
%         pvalcon=temp2(end,1);
% 
%         title([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rcon=' num2str(rvalcon) ',p=' num2str(pvalcon) ]);
%                 
%         [temp1,temp2]=corrcoef(ParameterPeaks(:,P,R),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:));
%         rvalpeak=temp1(end,1);
%         pvalpeak=temp2(end,1);
%         
%         subplot(2,2,length(ParameterToTest)+P);hold on;
%         plot(ParameterPeaks(:,P,R),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:),'.r','Markersize',5);        
%         title([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rpeak=' num2str(rvalpeak) ',p=' num2str(pvalpeak)  ]);
%        
%      disp([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rcon=' num2str(rvalcon) ',p=' num2str(pvalcon) ]);        
%      disp([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rpeak=' num2str(rvalpeak) ',p=' num2str(pvalpeak)  ]);
       
    end
    %textsave=squeeze(adjustedConvolvedSignal(:,R,:));
    %dlmwrite(fullfile(mfig_path,['supplFig9_convolvedHRF_Values_' ttestregions{R} '.txt']),textsave)   
 end

 %keyboard;
 
 
 ParameterToTest=[2 3];
RegionsToTest=[1 2 3 4 5 6 7];
BDesignNr=18; behavregressors=[4 3];

%keyboard
peakDetectLimit = [1 settings(2).window];
% % HRF settings - as in other paper, do not modify (unless there is a good reason)
sigma1=3;%old here:3;
my1=7;   % old here 6
alpha1=my1^2/sigma1^2;
beta1=my1/sigma1^2;
hrf_delay=0;
% Code the HRF
t=0:(settings(2).TR /settings(2).upsample):(settings(2).window-hrf_delay);
t=t(1:end-1);
hrf = gammapdf(alpha1,beta1,t);
simplehrf_peak_point = find(hrf==max(hrf));
 % delayhrf=[zeros(1,ceil(hrf_delay/(settings(2).TR  /settings(2).upsample)))  hrf];
% %delayhrf=[zeros(1,pseudosamples-length(hrf))  hrf];
% hrf_peak_time = find(delayhrf==max(delayhrf));
 for R=1:length(RegionsToTest)
    Ttestfigs(R)=figure('name',['Region:' Results(2).design(1).out.regions{RegionsToTest(R)} ' peak and HRF t-tests'],'color',[1 1 1]);hold on;
    %CorrBfigs(R)=figure('name',['Region:' Results(2).design(1).out.regions{RegionsToTest(R)} ' peak and HRF behavioural corrs']);hold on;

    for P=1:length(ParameterToTest)
        TCData=squeeze(Results(2).design(1).out.betas(:,RegionsToTest(R),:,ParameterToTest(P)));
        ttestLabels{P}=Results(2).design(1).out.labels{ParameterToTest(P)};
        ttestregions{R}=Results(2).design(1).out.regions{RegionsToTest(R)};
        % First Find the Peak
        dataT.v=TCData;
        [ParameterPeaks(:,P,R),ParameterPeakpoint(:,P,R)] = ComputePeakValue(dataT,peakDetectLimit,settings(2));
        title([ ttestregions{R} '  ' ttestLabels{P} ])
        PeakAdjust=floor(settings(2).pseudosamples/settings(2).window*peakDetectLimit(1))-1+ParameterPeakpoint(:,P,R)-simplehrf_peak_point;
        adjustedhrfs=[];
        for adj=1:length(PeakAdjust)
            if PeakAdjust(adj)>0
        adjustedhrfs(adj,:)=[zeros(1,PeakAdjust(adj))  hrf(1:end-PeakAdjust(adj))];
            elseif PeakAdjust(adj)<0
                    adjustedhrfs(adj,:)=[hrf(1-PeakAdjust(adj):end) zeros(1,-PeakAdjust(adj))];    
            else
              adjustedhrfs(adj,:)=  hrf;
            end
        end
        adjustedConvolvedSignal(:,R,P)= sum(TCData.*adjustedhrfs,2);
        [H,CadjPValue,CI,STATS]=ttest(adjustedConvolvedSignal(:,R,P));
        ConvadjPValue(R,P)=CadjPValue;%PeakStats(P,R).STATS=STATS;
        nonParaadjConvPValue(R,P)= signrank(adjustedConvolvedSignal(:,R,P));
        [h,ParameterPeakPValues(R,P),ci,stats]=ttest(ParameterPeaks(:,P,R));
        figure(Ttestfigs(R))      
        subplot(2,2,P);hold on;title([ ttestregions{R} '  ' ttestLabels{P} ])
        bar(1,mean(adjustedConvolvedSignal(:,R,P)));errorbar(1,mean(adjustedConvolvedSignal(:,R,P)),std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5))
        bar(2,mean(ParameterPeaks(:,P,R)));errorbar(2,mean(ParameterPeaks(:,P,R)),std(ParameterPeaks(:,P,R))./(length(ParameterPeaks(:,P,R)).^0.5))
        
        figure(PaperHRFTestsFig(R));hold on
        bar(CPIndex(R)+1,mean(adjustedConvolvedSignal(:,R,P)));
        errorbar(CPIndex(R)+1,mean(adjustedConvolvedSignal(:,R,P)),...
            std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5),'.k','Linewidth',4)
        [h,PaperHRFTestsFig_ps(R,CPIndex(R)+1),ci,PaperHRFTestsFig_stats{R}{CPIndex(R)+1}]=ttest(adjustedConvolvedSignal(:,R,P));
        PaperHRFTestsLabels{CPIndex(R)+1}=ttestLabels{P};
        PaperHRFTestsFigData{R}{CPIndex(R)+1}=adjustedConvolvedSignal(:,R,P);
        PaperHRFTestsFigDataTC{R}{CPIndex(R)+1}=TCData;
        significant=nan ;
        supersignificant=nan ;
        superdupersignificant=nan ;
        
        significant(PaperHRFTestsFig_ps(R,CPIndex(R)+1)<0.05)=true;
        supersignificant(PaperHRFTestsFig_ps(R,CPIndex(R)+1)<0.01)=true;
        superdupersignificant(PaperHRFTestsFig_ps(R,CPIndex(R)+1)<0.001)=true;
        
        
        plot(CPIndex(R)+1.1.*significant,[mean(adjustedConvolvedSignal(:,R,P))+(std(adjustedConvolvedSignal(:,R,P))./length(adjustedConvolvedSignal(:,R,P)).^.5)./2].*significant,'*k','markersize',15)
        plot(CPIndex(R)+1.2.*supersignificant,[mean(adjustedConvolvedSignal(:,R,P))+(std(adjustedConvolvedSignal(:,R,P))./length(adjustedConvolvedSignal(:,R,P)).^.5)./2].*supersignificant,'*k','markersize',15)
        plot(CPIndex(R)+1.3.*superdupersignificant,[mean(adjustedConvolvedSignal(:,R,P))+(std(adjustedConvolvedSignal(:,R,P))./length(adjustedConvolvedSignal(:,R,P)).^.5)./2].*superdupersignificant,'*k','markersize',15)

        HRF_Figstats.RegionName{R}=Results(1).design(1).out.regions{RegionsToTest(R)};
        HRF_Figstats.Region{R}.ParameterLabels{CPIndex(R)+1}=ttestLabels{P};
        HRF_Figstats.Region{R}.pvalue(CPIndex(R)+1)=PaperHRFTestsFig_ps(R,CPIndex(R)+1);         
        HRF_Figstats.Region{R}.tvalue(CPIndex(R)+1)=PaperHRFTestsFig_stats{R}{CPIndex(R)+1}.tstat;        
        HRF_Figstats.Region{R}.df(CPIndex(R)+1)=PaperHRFTestsFig_stats{R}{CPIndex(R)+1}.df;
        
        CPIndex(R)=CPIndex(R)+1;
        
        
          %between-subject correlation with behaviour
%         figure(CorrBfigs(R));
%         subplot(2,2,P);hold on;
%         plot(adjustedConvolvedSignal(:,R,P),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:),'.b','Markersize',5);
%         
%         [temp1,temp2]=corrcoef(adjustedConvolvedSignal(:,R,P),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:));
%         rvalcon=temp1(end,1);
%         pvalcon=temp2(end,1);
% 
%         title([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rcon=' num2str(rvalcon) ',p=' num2str(pvalcon) ]);
%                 
%         [temp1,temp2]=corrcoef(ParameterPeaks(:,P,R),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:));
%         rvalpeak=temp1(end,1);
%         pvalpeak=temp2(end,1);
%         
%         subplot(2,2,length(ParameterToTest)+P);hold on;
%         plot(ParameterPeaks(:,P,R),DataB.GLMs(BDesignNr).betaweights(1+behavregressors(P),:),'.r','Markersize',5);        
%         title([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rpeak=' num2str(rvalpeak) ',p=' num2str(pvalpeak)  ]);
%        
%      disp([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rcon=' num2str(rvalcon) ',p=' num2str(pvalcon) ]);        
%      disp([ ttestregions{R} ' ' ttestLabels{P}  ' ' DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressors(P)) 'rpeak=' num2str(rvalpeak) ',p=' num2str(pvalpeak)  ]);
    end
    %textsave=squeeze(adjustedConvolvedSignal(:,R,:));
    %dlmwrite(fullfile(mfig_path,['supplFig9_convolvedHRF_Values_' ttestregions{R} '.txt']),textsave)   
 end
 
 for R=1:length(RegionsToTest)
     
     figure(PaperHRFTestsFig(R))
     subplot(1,1,1,'XTick',1:length(PaperHRFTestsLabels),'XTickLabel',PaperHRFTestsLabels,'XTicklabelrotation',45,'fontsize',20)
    if ~exist(fullfile(mfig_path,ttestregions{R}))
        mkdir(fullfile(mfig_path,ttestregions{R}))
    end
    saveas(PaperHRFTestsFig(R), fullfile(mfig_path,ttestregions{R},['HRFConvolvBarPlot.jpg'])); %save the figure
    savefig(PaperHRFTestsFig(R), fullfile(mfig_path,ttestregions{R},['HRFConvolvBarPlot'])); %save the figure
    for P=1:length(PaperHRFTestsLabels)
    textsave=PaperHRFTestsFigDataTC{R}{P};
    dlmwrite(fullfile(mfig_path,ttestregions{R},[PaperHRFTestsLabels{P} '_TC.txt']),textsave,'delimiter','\t')
    textsave=PaperHRFTestsFigData{R}{P};
    dlmwrite(fullfile(mfig_path,ttestregions{R},[PaperHRFTestsLabels{P} '_BarData.txt']),textsave,'delimiter','\t')
    end
 end
%keyboard

 
 
 %% Correlation between behaviour and timecourse (or peak/hrf signal) [can also implement checking for other b regressor to neural relationship]
 %DataB
 %keyboard
 
 BDesignNr=18;
  RegionLabelsToTest={'cope5_cluster1_preSMA','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1',...
  'cope8_cluster7_midParietal_1','cope8_cluster5_IPS_1'};
 for r=1:length(RegionLabelsToTest)
     RegionsToTest(R)=find(strcmp(Results(1).design(1).out.regions,RegionLabelsToTest(r)));    
 end
 %keyboard
 disp(DataB.GLMs(BDesignNr).DesignName)

 
 for R=1:length(RegionLabelsToTest)
     
               figure('name',['region:' RegionLabelsToTest{R}],'color',[1 1 1])
     
     behavregressor=1;
          cLabel=DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressor);
     behavbetas=DataB.GLMs(BDesignNr).betaweights(1+behavregressor,:); % to skip constant
     neuralbetas= squeeze(Results(1).design(1).out.betas(:,RegionsToTest(R),:,3));
     for t=1:size(neuralbetas,2)
         [temp,temp2]=corrcoef(behavbetas',neuralbetas(:,t));
         corrTimecourse{1}(t)=temp(end,1);
         corrPs{1}(t)=temp2(end,1);
     end
     
     subplot(2,2,1);plot(corrTimecourse{1});title(['SD shared Evidence ' cLabel]);hold on;
     sig=nan(length(corrTimecourse{1}),1);
     sig(corrPs{1}<0.05)=zeros;     
     plot(sig,'k','Linewidth',5)
     
     behavregressor=2;     
          cLabel=DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressor);
     behavbetas=DataB.GLMs(BDesignNr).betaweights(1+behavregressor,:); % to skip constant
     neuralbetas= squeeze(Results(1).design(1).out.betas(:,RegionsToTest(R),:,2));
     for t=1:size(neuralbetas,2)
         [temp,temp2]=corrcoef(behavbetas',neuralbetas(:,t));
         corrTimecourse{2}(t)=temp(end,1);
         corrPs{2}(t)=temp2(end,1);
     end
subplot(2,2,2);plot(corrTimecourse{2});title(['SD occl Evidence ' cLabel]);hold on;
     sig=nan(length(corrTimecourse{2}),1);
     sig(corrPs{2}<0.05)=zeros;     
     plot(sig,'k','Linewidth',5)
     
     behavregressor=3;
          cLabel=DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressor);
     behavbetas=DataB.GLMs(BDesignNr).betaweights(1+behavregressor,:); % to skip constant
     neuralbetas= squeeze(Results(2).design(1).out.betas(:,RegionsToTest(R),:,3));
     for t=1:size(neuralbetas,2)
         [temp,temp2]=corrcoef(behavbetas',neuralbetas(:,t));
         corrTimecourse{3}(t)=temp(end,1);
         corrPs{3}(t)=temp2(end,1);
     end
     
     subplot(2,2,3);plot(corrTimecourse{3});title(['OD shared Evidence ' cLabel]);hold on;
     sig=nan(length(corrTimecourse{3}),1);
     sig(corrPs{3}<0.05)=zeros;     
     plot(sig,'k','Linewidth',5)
     
     
     behavregressor=4;
     cLabel=DataB.GLMs(BDesignNr).DesignMatrixLabels(behavregressor);
     behavbetas=DataB.GLMs(BDesignNr).betaweights(1+behavregressor,:); % to skip constant
     neuralbetas= squeeze(Results(2).design(1).out.betas(:,RegionsToTest(R),:,2));
     for t=1:size(neuralbetas,2)
         [temp,temp2]=corrcoef(behavbetas',neuralbetas(:,t));
         corrTimecourse{4}(t)=temp(end,1);
         corrPs{4}(t)=temp2(end,1);
     end
     
     subplot(2,2,4);plot(corrTimecourse{4});title(['OD occluded Evidence' cLabel]);hold on;
     sig=nan(length(corrTimecourse{4}),1);
     sig(corrPs{4}<0.05)=zeros;
     plot(sig,'k','Linewidth',5)
     
     
     
 end
 
 % Send all timecourses, data (txt) and stats
 
 
 
% ConvadjPValue
% nonParaadjConvPValue
% %figure;
% %plot



%keyboard;



%% PPI settings 

%keyboard
%regions={'cope5_cluster1_preSMA','cope3_cluster12_vmPFC_1','cope3_neg_cluster7_FPl_1','cope3_neg_cluster8_dACC_1',...
%  'cope8_cluster7_midParietal_1','cope8_cluster5_IPS_1','cope8_cluster6_TPJ_1'};

         PPIcombos = {'cope3_neg_cluster7_FPl_1','cope8_cluster6_TPJ_1','SD_ChOccEvid','avgmask',1};
 

ppiSettings(1).pseudosamples  = out.SelfDec.pseudosamples; %how data was upsampled, needed for plotting
ppiSettings(1).window         = out.SelfDec.window;
ppiSettings(1).TR             = out.SelfDec.TR;
ppiSettings(1).upsample       = out.SelfDec.upsample;

         for ippi = 1:size(PPIcombos,1)
             % Select seed region, target region, and psychological regressor
             ppiSettings(1).seedRegion   = PPIcombos{ippi,1};
             ppiSettings(1).targetRegion = PPIcombos{ippi,2};
             ppiSettings(1).psychReg     = PPIcombos{ippi,3};
             ppiSettings(1).globalTC     = PPIcombos{ippi,4};
             ppiSettings(1).split        =false;
             % For decision phase (time-locked to choice possible)
             PPIGLMinput.SelfDec=GLMinput.SelfDec;
             PPIGLMinput.SelfDec.Labels =[];
             for l=1:length(GLMinput.SelfDec.Labels)-1
                 PPIGLMinput.SelfDec.Labels{l}=GLMinput.SelfDec.Labels{l+1};
             end
             outT=PPIRunRegression(PPIGLMinput.SelfDec,SelfDecTC,ppiSettings(1));
             PPIResults(1).time = 'SelfDecision';
             PPIResults(1).design(ippi).inp = PPIGLMinput.SelfDec.DesignMatrix;
             PPIResults(1).design(ippi).out = outT;
         end
         
         %% Plot the results of the PPI - plot all time courses
         % if I want to, can specify colors
         plotl=[];
         PlotPPI(PPIResults(1),plotl,ppiSettings(1)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting

         
         PPIcombos = {'cope3_neg_cluster7_FPl_1','cope8_cluster6_TPJ_1','OD_ChOccEvid','avgmask',1;...
             ... 'cope8_cluster6_TPJ_1','cope3_neg_cluster7_FPl_1','OD_ChOccEvid','avgmask',1;
             };%'cope3_neg_cluster7_FPl_1','cope5_cluster1_preSMA','OD_ChOccEvid','avgmask',1;...
             %'cope5_cluster1_preSMA','cope3_neg_cluster7_FPl_1','OD_ChOccEvid','avgmask',1

 
             
             ppiSettings(2).pseudosamples  = out.OtherDec.pseudosamples; %how data was upsampled, needed for plotting
             ppiSettings(2).window         = out.OtherDec.window;
             ppiSettings(2).TR             = out.OtherDec.TR;
             ppiSettings(2).upsample       = out.OtherDec.upsample;

         for ippi = 1:size(PPIcombos,1)
             % Select seed region, target region, and psychological regressor
             ppiSettings(2).seedRegion   = PPIcombos{ippi,1};
             ppiSettings(2).targetRegion = PPIcombos{ippi,2};
             ppiSettings(2).psychReg     = PPIcombos{ippi,3};
             ppiSettings(2).globalTC     = PPIcombos{ippi,4};
             ppiSettings(2).split        =false;
             PPIGLMinput.OtherDec=GLMinput.OtherDec;
             PPIGLMinput.OtherDec.Labels =[];
             for l=1:length(GLMinput.OtherDec.Labels)-1
                 PPIGLMinput.OtherDec.Labels{l}=GLMinput.OtherDec.Labels{l+1};
             end
             % For decision phase (time-locked to choice possible)
             outT=PPIRunRegression(PPIGLMinput.OtherDec,OtherDecTC,ppiSettings(2));
             PPIResults(2).time = 'OtherDecision';
             PPIResults(2).design(ippi).inp = PPIGLMinput.OtherDec.DesignMatrix;
             PPIResults(2).design(ippi).out = outT;
         end
         %% Plot the results of the PPI - plot all time courses
         % if I want to, can specify colors
         plotl=[];
         PlotPPI(PPIResults(2),plotl,ppiSettings(2)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting

        
      %   keyboard
         
          PPIcombos = {'cope3_neg_cluster7_FPl_1','cope8_cluster6_TPJ_1','occludedvssharedMainDiff','avgmask',1}; %'cope8_cluster6_TPJ_1','cope3_neg_cluster7_FPl_1','occludedvssharedMainDiff','avgmask',1;...

 
 
% PPIcombos = {'dACCact','dlPFC','DecProspV','avgmask',1;...
%              'dACCact','dlPFC','DecV1sum', 'avgmask',1;
%              'dlPFC','dACCact','DecProspV','avgmask',1;...
%              'dlPFC','dACCact','DecV1sum', 'avgmask',1;...
%              'perigAnatomicalBinar','accumbensBilateral','DecCosts','avgmask',1;...
%              'accumbensBilateral','perigAnatomicalBinar','DecCosts','avgmask',1};
      

ppiSettings(3).pseudosamples  = out.Sampling.pseudosamples; %how data was upsampled, needed for plotting
ppiSettings(3).window         = out.Sampling.window;
ppiSettings(3).TR             = out.Sampling.TR;
ppiSettings(3).upsample       = out.Sampling.upsample;

         for ippi = 1:size(PPIcombos,1)
             % Select seed region, target region, and psychological regressor
             ppiSettings(3).seedRegion   = PPIcombos{ippi,1};
             ppiSettings(3).targetRegion = PPIcombos{ippi,2};
             ppiSettings(3).psychReg     = PPIcombos{ippi,3};
             ppiSettings(3).globalTC     = PPIcombos{ippi,4};
             ppiSettings(3).split        =false;
              PPIGLMinput.Sampling=GLMinput.Sampling;
             PPIGLMinput.Sampling.Labels =[];
             for l=1:length(GLMinput.Sampling.Labels)-1
                 PPIGLMinput.Sampling.Labels{l}=GLMinput.Sampling.Labels{l+1};
             end
             % For decision phase (time-locked to choice possible)
             outT=PPIRunRegression(PPIGLMinput.Sampling,SamplingTC,ppiSettings(3));
             PPIResults(3).time = 'Sampling';
             PPIResults(3).design(ippi).inp = PPIGLMinput.Sampling.DesignMatrix;
             PPIResults(3).design(ippi).out = outT;
         end
         
         %% Plot the results of the PPI - plot all time courses
         % if I want to, can specify colors
         plotl=[];
         PlotPPI(PPIResults(3),plotl,ppiSettings(3)) % inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
   
         Paper_Peak_PPIT_testFig=figure('name',['PPI_Peak_ttests'],'color',[1 1 1]);hold on;
         Ccount=0;
         clear PaperPPI_Peak_TestsFigData
         clear PaperPPI_Peak_TestsLabels
         PaperPPI_Peak_TestsFigData{1}=[];
         PaperPPI_Peak_TestsLabels{1}=[];
         PaperPPI_Peak_TestsFig_ps=[];
         
         PaperPPI_Peak_TestsFig_stats={};
         PaperPPI_Peak_TestsFigDataTC={};
         
         
  ParameterToTest=[6];
RegionsToTest=[1];
%keyboard
peakDetectLimit = [1 6];
% % HRF settings - as in other paper, do not modify (unless there is a good reason)
sigma1=3;%old here:3;
my1=7;   % old here 6
alpha1=my1^2/sigma1^2;
beta1=my1/sigma1^2;
hrf_delay=0;
% Code the HRF
t=0:(settings(2).TR /ppiSettings(1).upsample):(ppiSettings(1).window-hrf_delay);
t=t(1:end-1);
hrf = gammapdf(alpha1,beta1,t);
simplehrf_peak_point = find(hrf==max(hrf));
 % delayhrf=[zeros(1,ceil(hrf_delay/(settings(2).TR  /settings(2).upsample)))  hrf];
% %delayhrf=[zeros(1,pseudosamples-length(hrf))  hrf];
% hrf_peak_time = find(delayhrf==max(delayhrf));
 for R=1:length(RegionsToTest)
    PPITtestfigs(R)=figure('name',['Region:' PPIResults(1).design(1).out.regions{RegionsToTest(R)} ' peak and HRF t-tests'],'color',[1 1 1]);hold on;
    

    for P=1:length(ParameterToTest)
        TCData=squeeze(PPIResults(1).design(1).out.betas2(:,RegionsToTest(R),:,ParameterToTest(P)));
        ttestLabels{P}=PPIResults(1).design(1).out.labels2{ParameterToTest(P)};
        ttestregions{R}=PPIResults(1).design(1).out.regions{RegionsToTest(R)};
        % First Find the Peak
        dataT.v=TCData;
        [ParameterPeaks(:,P,R),ParameterPeakpoint(:,P,R)] = ComputePeakValue(dataT,peakDetectLimit,ppiSettings(1));
        title([ ttestregions{R} '  ' ttestLabels{P} ])
        PeakAdjust=floor(ppiSettings(1).pseudosamples/ppiSettings(1).window*peakDetectLimit(1))-1+ParameterPeakpoint(:,P,R)-simplehrf_peak_point;
        adjustedhrfs=[];
        for adj=1:length(PeakAdjust)
            if PeakAdjust(adj)>0
        adjustedhrfs(adj,:)=[zeros(1,PeakAdjust(adj))  hrf(1:end-PeakAdjust(adj))];
            elseif PeakAdjust(adj)<0
                    adjustedhrfs(adj,:)=[hrf(1-PeakAdjust(adj):end) zeros(1,-PeakAdjust(adj))];    
            else
              adjustedhrfs(adj,:)=  hrf;
            end
        end
        adjustedConvolvedSignal(:,R,P)= sum(TCData.*adjustedhrfs,2);
        [H,CadjPValue,CI,STATS]=ttest(adjustedConvolvedSignal(:,R,P));
        ConvadjPValue(R,P)=CadjPValue;%PeakStats(P,R).STATS=STATS;
        nonParaadjConvPValue(R,P)= signrank(adjustedConvolvedSignal(:,R,P));
        [h,ParameterPeakPValues(R,P),ci,stats]=ttest(ParameterPeaks(:,P,R));
        figure(PPITtestfigs(R))      
        subplot(2,2,P);hold on;title([ ttestregions{R} '  ' ttestLabels{P} ])
        bar(1,mean(adjustedConvolvedSignal(:,R,P)));errorbar(1,mean(adjustedConvolvedSignal(:,R,P)),std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5))
        bar(2,mean(ParameterPeaks(:,P,R)));errorbar(2,mean(ParameterPeaks(:,P,R)),std(ParameterPeaks(:,P,R))./(length(ParameterPeaks(:,P,R)).^0.5))
        
        
        
        figure(Paper_Peak_PPIT_testFig);hold on
        bar(Ccount+1,mean(ParameterPeaks(:,R,P)));
        errorbar(Ccount+1,mean(ParameterPeaks(:,R,P)),...
            std(ParameterPeaks(:,R,P))./(length(ParameterPeaks(:,R,P)).^0.5),'.k','Linewidth',4)
        [h,PaperPPI_Peak_TestsFig_ps(Ccount+1),ci,PaperPPI_Peak_TestsFig_stats{Ccount+1}]=ttest(ParameterPeaks(:,R,P));
        PaperPPI_Peak_TestsLabels{Ccount+1}=ttestLabels{P};
        PaperPPI_Peak_TestsFigData{Ccount+1}=ParameterPeaks(:,R,P);
        PaperPPI_Peak_TestsFigDataTC{Ccount+1}=TCData;
        significant=nan ;
        supersignificant=nan ;
        superdupersignificant=nan ;
        
        significant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.05)=true;
        supersignificant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.01)=true;
        superdupersignificant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.001)=true;
                
        plot(Ccount+1.1.*significant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*significant,'*k','markersize',15)
        plot(Ccount+1.2.*supersignificant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*supersignificant,'*k','markersize',15)
        plot(Ccount+1.3.*superdupersignificant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*superdupersignificant,'*k','markersize',15)
         PPI_Peak_Figstats.seedRegions{Ccount+1}=ppiSettings(3).seedRegion  ;
         PPI_Peak_Figstats.targetRegins{Ccount+1}=ppiSettings(3).targetRegion ; 
        
        PPI_Peak_Figstats.Region{R}.ParameterLabels{Ccount+1}=ttestLabels{P};
        PPI_Peak_Figstats.Region{R}.pvalue(Ccount+1)=PaperPPI_Peak_TestsFig_ps(Ccount+1);         
        PPI_Peak_Figstats.Region{R}.tvalue(Ccount+1)=PaperPPI_Peak_TestsFig_stats{Ccount+1}.tstat;        
        PPI_Peak_Figstats.Region{R}.df(Ccount+1)=PaperPPI_Peak_TestsFig_stats{Ccount+1}.df;
        Ccount=Ccount+1;
        
    end
    %textsave=squeeze(adjustedConvolvedSignal(:,R,:));
    %dlmwrite(fullfile(mfig_path,['supplFig9_convolvedHRF_Values_' ttestregions{R} '.txt']),textsave)   
 end

  ParameterToTest=[6];
RegionsToTest=[1];
%keyboard
peakDetectLimit = [1 6];
% % HRF settings - as in other paper, do not modify (unless there is a good reason)
sigma1=3;%old here:3;
my1=7;   % old here 6
alpha1=my1^2/sigma1^2;
beta1=my1/sigma1^2;
hrf_delay=0;
% Code the HRF
t=0:(settings(2).TR /ppiSettings(2).upsample):(ppiSettings(2).window-hrf_delay);
t=t(1:end-1);
hrf = gammapdf(alpha1,beta1,t);
simplehrf_peak_point = find(hrf==max(hrf));
 % delayhrf=[zeros(1,ceil(hrf_delay/(settings(2).TR  /settings(2).upsample)))  hrf];
% %delayhrf=[zeros(1,pseudosamples-length(hrf))  hrf];
% hrf_peak_time = find(delayhrf==max(delayhrf));
 for R=1:length(RegionsToTest)
    Ttestfigs(R)=figure('name',['Region:' PPIResults(2).design(1).out.regions{RegionsToTest(R)} ' peak and HRF t-tests'],'color',[1 1 1]);hold on;
    for P=1:length(ParameterToTest)
        TCData=squeeze(PPIResults(2).design(1).out.betas2(:,RegionsToTest(R),:,ParameterToTest(P)));
        ttestLabels{P}=PPIResults(2).design(1).out.labels2{ParameterToTest(P)};
        ttestregions{R}=PPIResults(2).design(1).out.regions{RegionsToTest(R)};
        % First Find the Peak
        dataT.v=TCData;
        [ParameterPeaks(:,P,R),ParameterPeakpoint(:,P,R)] = ComputePeakValue(dataT,peakDetectLimit,ppiSettings(2));
        title([ ttestregions{R} '  ' ttestLabels{P} ])
        PeakAdjust=floor(ppiSettings(2).pseudosamples/ppiSettings(2).window*peakDetectLimit(1))-1+ParameterPeakpoint(:,P,R)-simplehrf_peak_point;
        adjustedhrfs=[];
        for adj=1:length(PeakAdjust)
            if PeakAdjust(adj)>0
        adjustedhrfs(adj,:)=[zeros(1,PeakAdjust(adj))  hrf(1:end-PeakAdjust(adj))];
            elseif PeakAdjust(adj)<0
                    adjustedhrfs(adj,:)=[hrf(1-PeakAdjust(adj):end) zeros(1,-PeakAdjust(adj))];    
            else
              adjustedhrfs(adj,:)=  hrf;
            end
        end
        adjustedConvolvedSignal(:,R,P)= sum(TCData.*adjustedhrfs,2);
        [H,CadjPValue,CI,STATS]=ttest(adjustedConvolvedSignal(:,R,P));
        ConvadjPValue(R,P)=CadjPValue;%PeakStats(P,R).STATS=STATS;
        nonParaadjConvPValue(R,P)= signrank(adjustedConvolvedSignal(:,R,P));
        [h,ParameterPeakPValues(R,P),ci,stats]=ttest(ParameterPeaks(:,P,R));
        figure(Ttestfigs(R))      
        subplot(2,2,P);hold on;title([ ttestregions{R} '  ' ttestLabels{P} ])
        bar(1,mean(adjustedConvolvedSignal(:,R,P)));errorbar(1,mean(adjustedConvolvedSignal(:,R,P)),std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5))
        bar(2,mean(ParameterPeaks(:,P,R)));errorbar(2,mean(ParameterPeaks(:,P,R)),std(ParameterPeaks(:,P,R))./(length(ParameterPeaks(:,P,R)).^0.5))
        
        figure(Paper_Peak_PPIT_testFig);hold on
        bar(Ccount+1,mean(ParameterPeaks(:,R,P)));
        errorbar(Ccount+1,mean(ParameterPeaks(:,R,P)),...
            std(ParameterPeaks(:,R,P))./(length(ParameterPeaks(:,R,P)).^0.5),'.k','Linewidth',4)
        [h,PaperPPI_Peak_TestsFig_ps(Ccount+1),ci,PaperPPI_Peak_TestsFig_stats{Ccount+1}]=ttest(ParameterPeaks(:,R,P));
        PaperPPI_Peak_TestsLabels{Ccount+1}=ttestLabels{P};
        PaperPPI_Peak_TestsFigData{Ccount+1}=ParameterPeaks(:,R,P);
        PaperPPI_Peak_TestsFigDataTC{Ccount+1}=TCData;
        significant=nan ;
        supersignificant=nan ;
        superdupersignificant=nan ;
        
        significant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.05)=true;
        supersignificant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.01)=true;
        superdupersignificant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.001)=true;
        
        plot(Ccount+1.1.*significant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*significant,'*k','markersize',15)
        plot(Ccount+1.2.*supersignificant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*supersignificant,'*k','markersize',15)
        plot(Ccount+1.3.*superdupersignificant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*superdupersignificant,'*k','markersize',15)
         PPI_Peak_Figstats.seedRegions{Ccount+1}=ppiSettings(3).seedRegion  ;
         PPI_Peak_Figstats.targetRegins{Ccount+1}=ppiSettings(3).targetRegion ; 
        
        PPI_Peak_Figstats.Region{R}.ParameterLabels{Ccount+1}=ttestLabels{P};
        PPI_Peak_Figstats.Region{R}.pvalue(Ccount+1)=PaperPPI_Peak_TestsFig_ps(Ccount+1);         
        PPI_Peak_Figstats.Region{R}.tvalue(Ccount+1)=PaperPPI_Peak_TestsFig_stats{Ccount+1}.tstat;        
        PPI_Peak_Figstats.Region{R}.df(Ccount+1)=PaperPPI_Peak_TestsFig_stats{Ccount+1}.df;
        Ccount=Ccount+1;
    
    end
    %textsave=squeeze(adjustedConvolvedSignal(:,R,:));
    %dlmwrite(fullfile(mfig_path,['supplFig9_convolvedHRF_Values_' ttestregions{R} '.txt']),textsave)   
 end        
         
         
ParameterToTest=[9];
RegionsToTest=[1];
%keyboard
peakDetectLimit = [1 6];
% % HRF settings - as in other paper, do not modify (unless there is a good reason)
sigma1=3;%old here:3;
my1=7;   % old here 6
alpha1=my1^2/sigma1^2;
beta1=my1/sigma1^2;
hrf_delay=0;
% Code the HRF
t=0:(settings(2).TR /ppiSettings(3).upsample):(ppiSettings(3).window-hrf_delay);
t=t(1:end-1);
hrf = gammapdf(alpha1,beta1,t);
simplehrf_peak_point = find(hrf==max(hrf));
 % delayhrf=[zeros(1,ceil(hrf_delay/(settings(2).TR  /settings(2).upsample)))  hrf];
% %delayhrf=[zeros(1,pseudosamples-length(hrf))  hrf];
% hrf_peak_time = find(delayhrf==max(delayhrf));
 for R=1:length(RegionsToTest)
    Ttestfigs(R)=figure('name',['Region:' PPIResults(3).design(1).out.regions{RegionsToTest(R)} ' peak and HRF t-tests'],'color',[1 1 1]);hold on;
    for P=1:length(ParameterToTest)
        TCData=squeeze(PPIResults(3).design(1).out.betas2(:,RegionsToTest(R),:,ParameterToTest(P)));
        ttestLabels{P}=PPIResults(3).design(1).out.labels2{ParameterToTest(P)};
        ttestregions{R}=PPIResults(3).design(1).out.regions{RegionsToTest(R)};
        % First Find the Peak
        dataT.v=TCData;
        [ParameterPeaks(:,P,R),ParameterPeakpoint(:,P,R)] = ComputePeakValue(dataT,peakDetectLimit,ppiSettings(3));
        title([ ttestregions{R} '  ' ttestLabels{P} ])
        PeakAdjust=floor(ppiSettings(3).pseudosamples/ppiSettings(3).window*peakDetectLimit(1))-1+ParameterPeakpoint(:,P,R)-simplehrf_peak_point;
        adjustedhrfs=[];
        for adj=1:length(PeakAdjust)
            if PeakAdjust(adj)>0
        adjustedhrfs(adj,:)=[zeros(1,PeakAdjust(adj))  hrf(1:end-PeakAdjust(adj))];
            elseif PeakAdjust(adj)<0
                    adjustedhrfs(adj,:)=[hrf(1-PeakAdjust(adj):end) zeros(1,-PeakAdjust(adj))];    
            else
              adjustedhrfs(adj,:)=  hrf;
            end
        end
        adjustedConvolvedSignal(:,R,P)= sum(TCData.*adjustedhrfs,2);
        [H,CadjPValue,CI,STATS]=ttest(adjustedConvolvedSignal(:,R,P));
        ConvadjPValue(R,P)=CadjPValue;%PeakStats(P,R).STATS=STATS;
        nonParaadjConvPValue(R,P)= signrank(adjustedConvolvedSignal(:,R,P));
        [h,ParameterPeakPValues(R,P),ci,stats]=ttest(ParameterPeaks(:,P,R));
        figure(Ttestfigs(R))      
        subplot(2,2,P);hold on;title([ ttestregions{R} '  ' ttestLabels{P} ])
        bar(1,mean(adjustedConvolvedSignal(:,R,P)));errorbar(1,mean(adjustedConvolvedSignal(:,R,P)),std(adjustedConvolvedSignal(:,R,P))./(length(adjustedConvolvedSignal(:,R,P)).^0.5))
        bar(2,mean(ParameterPeaks(:,P,R)));errorbar(2,mean(ParameterPeaks(:,P,R)),std(ParameterPeaks(:,P,R))./(length(ParameterPeaks(:,P,R)).^0.5))
        
        figure(Paper_Peak_PPIT_testFig);hold on
        bar(Ccount+1,mean(ParameterPeaks(:,R,P)));
        errorbar(Ccount+1,mean(ParameterPeaks(:,R,P)),...
            std(ParameterPeaks(:,R,P))./(length(ParameterPeaks(:,R,P)).^0.5),'.k','Linewidth',4)
        [h,PaperPPI_Peak_TestsFig_ps(Ccount+1),ci,PaperPPI_Peak_TestsFig_stats{Ccount+1}]=ttest(ParameterPeaks(:,R,P));
        PaperPPI_Peak_TestsLabels{Ccount+1}=ttestLabels{P};
        PaperPPI_Peak_TestsFigData{Ccount+1}=ParameterPeaks(:,R,P);
        PaperPPI_Peak_TestsFigDataTC{Ccount+1}=TCData;
        significant=nan ;
        supersignificant=nan ;
        superdupersignificant=nan ;
        
        significant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.05)=true;
        supersignificant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.01)=true;
        superdupersignificant(PaperPPI_Peak_TestsFig_ps(Ccount+1)<0.001)=true;
                
        plot(Ccount+1.1.*significant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*significant,'*k','markersize',15)
        plot(Ccount+1.2.*supersignificant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*supersignificant,'*k','markersize',15)
        plot(Ccount+1.3.*superdupersignificant,[mean(ParameterPeaks(:,R,P))+(std(ParameterPeaks(:,R,P))./length(ParameterPeaks(:,R,P)).^.5)./2].*superdupersignificant,'*k','markersize',15)
        
        
        PPI_Peak_Figstats.seedRegions{Ccount+1}=ppiSettings(3).seedRegion  ;
         PPI_Peak_Figstats.targetRegins{Ccount+1}=ppiSettings(3).targetRegion ; 
        
        PPI_Peak_Figstats.Region{R}.ParameterLabels{Ccount+1}=ttestLabels{P};
        PPI_Peak_Figstats.Region{R}.pvalue(Ccount+1)=PaperPPI_Peak_TestsFig_ps(Ccount+1);         
        PPI_Peak_Figstats.Region{R}.tvalue(Ccount+1)=PaperPPI_Peak_TestsFig_stats{Ccount+1}.tstat;        
        PPI_Peak_Figstats.Region{R}.df(Ccount+1)=PaperPPI_Peak_TestsFig_stats{Ccount+1}.df;
        
        Ccount=Ccount+1;
    end
    %textsave=squeeze(adjustedConvolvedSignal(:,R,:));
    %dlmwrite(fullfile(mfig_path,['supplFig9_convolvedHRF_Values_' ttestregions{R} '.txt']),textsave)   
 end
 
 save(fullfile(mfig_path,'HRF_FigStats.mat'),'HRF_Figstats')
  save(fullfile(mfig_path,'PPI_Peak_FigStats.mat'),'PPI_Peak_Figstats')

figure(Paper_Peak_PPIT_testFig);
subplot(1,1,1,'XTick',1:length(PaperPPI_Peak_TestsLabels),'XTickLabel',PaperPPI_Peak_TestsLabels,'XTicklabelrotation',45,'fontsize',20)
     
    saveas(Paper_Peak_PPIT_testFig, fullfile(mfig_path,['Peak_PPI_BarPlot.jpg'])); %save the figure
    savefig(Paper_Peak_PPIT_testFig, fullfile(mfig_path,['Peak_PPI_BarPlot'])); %save the figure
    for P=1:length(PaperPPI_Peak_TestsLabels)
    textsave=PaperPPI_Peak_TestsFigDataTC{P};
    dlmwrite(fullfile(mfig_path,[PaperPPI_Peak_TestsLabels{P} '_TC.txt']),textsave,'delimiter','\t')
    textsave=PaperPPI_Peak_TestsFigData{P};
    dlmwrite(fullfile(mfig_path,[PaperPPI_Peak_TestsLabels{P} '_BarData.txt']),textsave,'delimiter','\t')
    end
 %% For PPI's statistics and plots check whether Fpole for TPJ always or change direction or sampling. Maybe we only show Peaks or HRF convolved results



%keyboard;
