function [Data,Results]=RunTimeCourse(GLMinput,Data,GLMname,DataB)
%% This the function to run timecourses with

% written originally by Jacqueline Scholl, adapted by Nils Kolling
%%
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, 'git_repository', 'Matlab_Scripts')) % 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
% These timecourses are labelled according to the cope used to identify
% them (Sampling_ChosenEvid analysis), then the region label, except for
% avgmask which is the average global timecourse needed for PPI analyses
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'};

% The naming convention for the timecourse texts is :
% XX=SUBJECTNR
% YY=COPENUMBER
% ZZ=CLUSTERNR
% A=SUBCLUSTERNR
% XX_copeYY_clusterZZ_AnatomicalLabel_A
% XX_copeYY_neg_clusterZZ_AnatomicalLabel_A


%% Load and upsample the timecourses %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%(Also time-locking to the respective event of interest, except for Other-decision where we cut out the first 3 seconds, because of previous event and HRF lag means that data is not relevant anyway)
% in this script timewindows are specified manually (SD=12, OD=6 but 3 moved, so from 3 to 9, Sampling=12, Outcome=8)
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 Self-decision phase
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;


% For Other-decision phase
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;

% For Sampling phase
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;



%% Plot the results of the timecourse analyses for the paper: only selected areas and regressors %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(GLMname,'Sampling_ChosenEvid')
    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
    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
    
    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
    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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    ParameterToTest=[2 3];
    RegionsToTest=[1 2 3 4 5 6 7];
    
    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
                dataforAverage=[]; datape=[];

    for R=1:length(RegionsToTest)
        
        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));
            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(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));
            plot(repmat(CPIndex(R)+1,size(adjustedConvolvedSignal,1),1),adjustedConvolvedSignal(:,R,P),'*k','Markersize',2)

            dataforAverage(:,CPIndex(R)+1,R)=adjustedConvolvedSignal(:,R,P);
            datape(CPIndex(R)+1,R)=PaperHRFTestsFig_ps(R,CPIndex(R)+1);
            
            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;
        end
    end
    
    ParameterToTest=[2 3];
    RegionsToTest=[1 2 3 4 5 6 7];
    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));
    
    
    for R=1:length(RegionsToTest)
        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));
            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(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));
            plot(repmat(CPIndex(R)+1,size(adjustedConvolvedSignal,1),1),adjustedConvolvedSignal(:,R,P),'*k','Markersize',2)
            dataforAverage(:,CPIndex(R)+1,R)=adjustedConvolvedSignal(:,R,P);
            datape(CPIndex(R)+1,R)=PaperHRFTestsFig_ps(R,CPIndex(R)+1);
            
            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;
        end
    end
    
    % Make Violin Plots
    
    for R=1:length(RegionsToTest)
        
        %dataforAverage(:,CPIndex(R)+1,R)=adjustedConvolvedSignal(:,R,P);
        %datape(CPIndex(R)+1,R)=PaperHRFTestsFig_ps(R,CPIndex(R)+1);
        
        
        figure('name',['Region:' Results(2).design(1).out.regions{RegionsToTest(R)} ' HRF Violins'],'color',[1 1 1]);
        subplot(1,1,1,'XTick',1:length(PaperHRFTestsLabels),'XTickLabel',PaperHRFTestsLabels,'XTicklabelrotation',45,'fontsize',20);hold on;

        %violinplot(dataforAverage(:,:,R));
        hold on;
        for pe=1:length(datape(:,R))
            significant=nan ;
            supersignificant=nan ;
            superdupersignificant=nan ;
            
            significant(datape(pe,R)<0.05)=true;
            supersignificant(datape(pe,R)<0.01)=true;
            superdupersignificant(datape(pe,R)<0.001)=true;
            
            plot(pe+0.1.*significant,[mean(dataforAverage(:,pe,R))+(std(dataforAverage(:,pe,R))./length(dataforAverage(:,pe,R)).^.5)./2].*significant,'*k','markersize',15)
            plot(pe+0.2.*supersignificant,[mean(dataforAverage(:,pe,R))+(std(dataforAverage(:,pe,R))./length(dataforAverage(:,pe,R)).^.5)./2].*supersignificant,'*k','markersize',15)
            plot(pe+0.3.*superdupersignificant,[mean(dataforAverage(:,pe,R))+(std(dataforAverage(:,pe,R))./length(dataforAverage(:,pe,R)).^.5)./2].*superdupersignificant,'*k','markersize',15)
        end
        
        plot([0 size(dataforAverage,2)+1],[0 0],'k','Linewidth',2)

        
      saveas(gcf, fullfile(mfig_path,ttestregions{R},['HRFConvolv_Violin.jpg'])); %save the figure
      savefig(gcf, fullfile(mfig_path,ttestregions{R},['HRFConvolv_Violin'])); %save the figure
        
    end
    
    
    %% Save lots of timecourse figures and underlying data
    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
    
end


%% PPI analyses %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(GLMname,'Sampling_ChosenEvid')
    
 
    
   
    
    
    %% PPI in the Sampling phase %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    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;...
    ppiSettings(1).pseudosamples  = out.Sampling.pseudosamples; %how data was upsampled, needed for plotting
    ppiSettings(1).window         = out.Sampling.window;
    ppiSettings(1).TR             = out.Sampling.TR;
    ppiSettings(1).upsample       = out.Sampling.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;
        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(1));
        PPIResults(1).time = 'Sampling';
        PPIResults(1).design(ippi).inp = PPIGLMinput.Sampling.DesignMatrix;
        PPIResults(1).design(ippi).out = outT;
    end
    
     %% PPI in the other decision phase %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    PPIcombos = {'cope3_neg_cluster7_FPl_1','cope8_cluster6_TPJ_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
        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
    
    %% now run the peak extraction and plotting for the PPI results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Paper_Peak_PPIT_testFig=figure('name',['PPI_Peak_ttests'],'color',[1 1 1]);hold on;
    Ccount=0;dataforAverage=[];datap=[];
    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=[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(1).upsample):(ppiSettings(1).window-hrf_delay);
    t=t(1:end-1);
    hrf = gammapdf(alpha1,beta1,t);
    simplehrf_peak_point = find(hrf==max(hrf));
    for R=1:length(RegionsToTest)
        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));
            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(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)
            plot(repmat(Ccount+1,size(ParameterPeaks,1),1),ParameterPeaks(:,R,P),'*k','Markersize',2)

            dataforAverage(:,Ccount+1)=ParameterPeaks(:,R,P);
            
            [h,PaperPPI_Peak_TestsFig_ps(Ccount+1),ci,PaperPPI_Peak_TestsFig_stats{Ccount+1}]=ttest(ParameterPeaks(:,R,P));
            datap(Ccount+1)=PaperPPI_Peak_TestsFig_ps(Ccount+1);

            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(2).seedRegion  ;
            PPI_Peak_Figstats.targetRegins{Ccount+1}=ppiSettings(2).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
    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));
    for R=1:length(RegionsToTest)
        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));
            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(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)
            plot(repmat(Ccount+1,size(ParameterPeaks,1),1),ParameterPeaks(:,R,P),'*k','Markersize',2)
            dataforAverage(:,Ccount+1)=ParameterPeaks(:,R,P);
            [h,PaperPPI_Peak_TestsFig_ps(Ccount+1),ci,PaperPPI_Peak_TestsFig_stats{Ccount+1}]=ttest(ParameterPeaks(:,R,P));
                        datap(Ccount+1)=PaperPPI_Peak_TestsFig_ps(Ccount+1);
            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(2).seedRegion  ;
            PPI_Peak_Figstats.targetRegins{Ccount+1}=ppiSettings(2).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
    end
    figure;
        subplot(1,1,1,'XTick',1:length(PaperPPI_Peak_TestsLabels),'XTickLabel',PaperPPI_Peak_TestsLabels,'XTicklabelrotation',45,'fontsize',20);hold on;
%    violinplot(dataforAverage);
    hold on;
    for pe=1:length(datap)
    significant=nan ;
            supersignificant=nan ;
            superdupersignificant=nan ;
            
            significant(datap(pe)<0.05)=true;
            supersignificant(datap(pe)<0.01)=true;
            superdupersignificant(datap(pe)<0.001)=true;
            
            plot(pe+0.1.*significant,[mean(dataforAverage(:,pe))+(std(dataforAverage(:,pe))./length(dataforAverage(:,pe)).^.5)./2].*significant,'*k','markersize',15)
            plot(pe+0.2.*supersignificant,[mean(dataforAverage(:,pe))+(std(dataforAverage(:,pe))./length(dataforAverage(:,pe)).^.5)./2].*supersignificant,'*k','markersize',15)
            plot(pe+0.3.*superdupersignificant,[mean(dataforAverage(:,pe))+(std(dataforAverage(:,pe))./length(dataforAverage(:,pe)).^.5)./2].*superdupersignificant,'*k','markersize',15)
    end
            plot([0 size(dataforAverage,2)+1],[0 0],'k','Linewidth',2)
 saveas(gcf, fullfile(mfig_path,['Peak_PPI_Violin.jpg'])); %save the figure
    savefig(gcf, fullfile(mfig_path,['Peak_PPI_Violin'])); %save the figure
    
    
    %% now save figure and data underlying it
    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
    
end

%% The end of the MS script. below some more potentially useful information

%% Additional info: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strcmp(GLMname,'ChoiceDivergenceAnalysis')
    %% Plot the results of time course analyses - all areas together and all regressors
    % 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));
    
    
    %% Possible additional analysis of OutcomePhase:
    % We are not using outcome timecourse analysis, so could remove this
    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;
    
    clear plotl
    plotl.designList=1;
    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
end

if 0 % change to 1 if you want to see all PPI results. might not work with all glms because PPI are regressor specific
    %% Plot all results of PPI
    plotl=[];% inputs: Data that has the regression weights, regions to plot, regressors to plot, designs to plot and settings for plotting
    PlotPPI(PPIResults(1),plotl,ppiSettings(1)) % Other-dec
    PlotPPI(PPIResults(2),plotl,ppiSettings(2)) % Sampling period
end