function tutorial_visual_group(ProtocolName, reports_dir)
% TUTORIAL_VISUAL_GROUP: Runs the Brainstorm/SPM group analysis pipeline (group analysis, BIDS VERSION)
%
% ONLINE TUTORIALS: 
%    - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualSingle
%    - https://neuroimage.usc.edu/brainstorm/Tutorials/VisualGroup
%
% INPUTS:
%    - ProtocolName : Name of the protocol in which the recordings for the 19 subjects have been imported (TutorialVisual or TutorialGroup)
%                     This folder must include: the run-level averages (recordings and time-frequency) and the source kernels for each run
%    - reports_dir  : If defined, exports all the reports as HTML to this folder

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
% 
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Author: Francois Tadel, Elizabeth Bock, 2016-2018

% ===== CHECK PROTOCOL =====
% Start Brainstorm without the GUI
if ~brainstorm('status')
    brainstorm nogui
end
% Output folder for reports
if (nargin < 2) || isempty(reports_dir) || ~isdir(reports_dir)
    reports_dir = [];
end
% You have to specify the folder in which the tutorial dataset is unzipped
if (nargin < 1) || isempty(ProtocolName)
    ProtocolName = 'TutorialGroup';
end
% Select current protocol
iProtocol = bst_get('Protocol', ProtocolName);
if isempty(iProtocol)
    error(['Unknown protocol: ' ProtocolName]);
end
if (iProtocol ~= bst_get('iProtocol'))
    gui_brainstorm('SetCurrentProtocol', iProtocol);
end
% Process: Select results files in: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch
sFiles = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname', 'sub-01', ...
    'condition',   'sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch');
if isempty(sFiles)
    error(['No source files available in folder: sub-01/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg_notch.' 10 ...
           'You should run tutorial_visual_single first, or download the ' 10 ...
           'protocol TutorialGroup.zip from the Brainstorm website.']);
end
% Process: Select file comments with tag: Unconstr
sFiles = bst_process('CallProcess', 'process_select_tag', sFiles, [], 'tag', 'Unconstr');
if ~isempty(sFiles)
    isUnconstrained = 1;
else
    isUnconstrained = 0;
end
% Start a new report
bst_report('Start');


% ===== DELETE ALL THE EXISTING FILES =====
GroupSubjectName = bst_get('NormalizedSubjectName');
DirIntra = bst_get('DirAnalysisIntra');
% Delete group analysis
[sSubject, iSubject] = bst_get('Subject', GroupSubjectName);
if ~isempty(sSubject)
    db_delete_subjects(iSubject);
end
% Get all the files in the intra subjects
sAvgSrc = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% If there are existing files: delete the contents of the folders
if ~isempty(sAvgSrc)
    % Get the base folder of the protocol
    ProtocolInfo = bst_get('ProtocolInfo', iProtocol);
    % Get all the folders
    AllIntraDir = unique(cellfun(@bst_fileparts, {sAvgSrc.FileName}, 'UniformOutput', 0));
    % Loop over the intra folders
    for iSubj = 1:length(AllIntraDir)
        % Get all the files
        subjIntraDir = dir(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, '*.mat'));
        subjIntraFiles = setdiff({subjIntraDir.name}, 'brainstormstudy.mat');
        % Delete all the files
        for iFile = 1:length(subjIntraFiles)
            delete(fullfile(ProtocolInfo.STUDIES, AllIntraDir{iSubj}, subjIntraFiles{iFile}));
        end
    end
    % Reload protocol
    db_reload_database(iProtocol);
end


%% ===== SUBJECT AVERAGES: FAMOUS, UNFAMILIAR, SCRAMBLED ==============================
%  ====================================================================================

% ===== SUBJECT AVERAGE: MEG+EEG =====
% Process: Select data files in: */*/Avg
sAvgRun = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname', 'All', ...
    'tag',         'Avg');
% Process: Weighted Average: By trial group (subject average)
sAvgSubj = bst_process('CallProcess', 'process_average', sAvgRun, [], ...
    'avgtype',    6, ...  % By trial group (subject average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   1, ...
    'keepevents', 1);
% Process: Average: By trial group (grand average)
sAvgGroup = bst_process('CallProcess', 'process_average', sAvgSubj, [], ...
    'avgtype',    7, ...  % By trial group (grand average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== SUBJECT AVERAGES: SOURCES (EEG) =====
% Process: Select source files in: */*/EEG
sAvgRunSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'EEG');
% Process: Weighted Average: By trial group (subject average) - EEG
sAvgSubjSrcEeg = bst_process('CallProcess', 'process_average', sAvgRunSrcEeg, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 0);
% Process: Add tag: EEG
sAvgSubjSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcEeg, [], ...
    'tag',    'EEG', ...
    'output', 1);  % Add to comment

% ===== SUBJECT AVERAGES: SOURCES (MEG) =====
% Process: Select source files in: */*/MEG
sAvgRunSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], 'tag', 'MEG');
% Process: Weighted Average: By trial group (subject average) - MEG
sAvgSubjSrcMeg = bst_process('CallProcess', 'process_average', sAvgRunSrcMeg, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: MEG
sAvgSubjSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjSrcMeg, [], ...
    'tag',    'MEG', ...
    'output', 1);  % Add to comment

% ===== SUBJECT AVERAGES: TIMEFREQ =====
% Process: Select time-frequency files in: */*
sAvgRunTf = bst_process('CallProcess', 'process_select_files_timefreq', [], []);
% Process: Weighted Average: By trial group (subject average)
sAvgSubjTf = bst_process('CallProcess', 'process_average', sAvgRunTf, [], ...
    'avgtype',         6, ...  % By trial group (subject average)
    'avg_func',        1, ...  % Arithmetic average:  mean(x)
    'weighted',        1, ...
    'scalenormalized', 1);



%% ===== SUBJECT AVERAGES: FACES ======================================================
%  ====================================================================================

% ===== FACES AVERAGE: MEG+EEG ======
% Process: Select data files in: */*
sAllData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous (6 files)"
sAvgSubjFamous = bst_process('CallProcess', 'process_select_tag', sAllData, [], 'tag', 'WAvg: Avg: Famous (6 files)');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar (6 files)"
sAvgSubjUnfamiliar = bst_process('CallProcess', 'process_select_tag', sAllData, [],  'tag', 'WAvg: Avg: Unfamiliar (6 files)');
% Process: Weighted Average A&B
sAvgSubjFaces = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamous, sAvgSubjUnfamiliar, ...
    'weighted', 1);
% Process: Add tag: Faces
sAvgSubjFaces = bst_process('CallProcess', 'process_add_tag', sAvgSubjFaces, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFaces = bst_process('CallProcess', 'process_set_comment', sAvgSubjFaces, [], ...
    'tag',     'WAvg: Avg: Faces', ...
    'isindex', 0);
% Process: Average: By trial group (grand average)
sAvgGroupFaces = bst_process('CallProcess', 'process_average', sAvgSubjFaces, [], ...
    'avgtype',    7, ...  % By trial group (grand average)
    'avg_func',   1, ...  % Arithmetic average:  mean(x)
    'weighted',   0, ...
    'keepevents', 0);

% ===== FACES AVERAGE: SOURCES (EEG) ======
% Process: Select source files in: */*/EEG
sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'EEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ...
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: Faces
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcEeg, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcEeg, [], ...
    'tag',     'WAvg: Avg: Faces | EEG', ...
    'isindex', 0);

% ===== FACES AVERAGE: SOURCES (MEG) ======
% Process: Select source files in: */*
sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'MEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: Famous
sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Famous');
% Process: Select file comments with tag: Unfamiliar
sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ...
    'weighted',        1, ...
    'scalenormalized', 1);
% Process: Add tag: Faces
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesSrcMeg, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesSrcMeg, [], ...
    'tag',     'WAvg: Avg: Faces | MEG', ...
    'isindex', 0);

% ===== FACES AVERAGE: TIMEFREQ ======
% Process: Select time-frequency files in: */*
sAllTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarTf = bst_process('CallProcess', 'process_select_tag', sAllTf, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Weighted Average A&B
sAvgSubjFacesTf = bst_process('CallProcess', 'process_average_ab', sAvgSubjFamousTf, sAvgSubjUnfamiliarTf, ...
    'weighted', 1);
% Process: Add tag: Faces
sAvgSubjFacesTf = bst_process('CallProcess', 'process_add_tag', sAvgSubjFacesTf, [], ...
    'tag',    'Faces', ...
    'output', 2);  % Add to file name
% Process: Set comment
sAvgSubjFacesTf = bst_process('CallProcess', 'process_set_comment', sAvgSubjFacesTf, [], ...
    'tag',     'WAvg: Avg: Faces', ...
    'isindex', 0);



%% ===== SUBJECT AVERAGES: WITHIN-SUBJECT DIFFERENCES =================================
%  ====================================================================================

% ===== FACES - SCRAMBLED: SOURCES EEG =====
% Process: Select source files in: */*/EEG
sSrcEeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'EEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Faces"
sAvgSubjFacesSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Faces');
% Process: Select file comments with tag: "WAvg: Avg: Scrambled"
sAvgSubjScrambledSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Scrambled');
% Process: Difference: A-B
sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcEeg, sAvgSubjScrambledSrcEeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFacesSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcEeg, [], ...
    'tag',     'Faces - Scrambled | EEG', ...
    'isindex', 0);

% ===== FACES - SCRAMBLED: SOURCES MEG =====
% Process: Select source files in: */*/MEG
sSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'tag',           'MEG', ...
    'includeintra',  1);
% Process: Select file comments with tag: "WAvg: Avg: Faces"
sAvgSubjFacesSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Faces');
% Process: Select file comments with tag: "WAvg: Avg: Scrambled"
sAvgSubjScrambledSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Scrambled');
% Process: Difference: A-B
sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFacesSrcMeg, sAvgSubjScrambledSrcMeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFacesSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFacesSubjSrcMeg, [], ...
    'tag',     'Faces - Scrambled | MEG', ...
    'isindex', 0);

% ===== FAMOUS - UNFAMILIAR: SOURCE EEG =====
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcEeg = bst_process('CallProcess', 'process_select_tag', sSrcEeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Difference: A-B
sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcEeg, sAvgSubjUnfamiliarSrcEeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFamousSubjSrcEeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcEeg, [], ...
    'tag',     'Famous - Unfamiliar | EEG', ...
    'isindex', 0);

% ===== FAMOUS - UNFAMILIAR: SOURCE MEG =====
% Process: Select file comments with tag: "WAvg: Avg: Famous"
sAvgSubjFamousSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Famous');
% Process: Select file comments with tag: "WAvg: Avg: Unfamiliar"
sAvgSubjUnfamiliarSrcMeg = bst_process('CallProcess', 'process_select_tag', sSrcMeg, [], 'tag', 'WAvg: Avg: Unfamiliar');
% Process: Difference: A-B
sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_diff_ab', sAvgSubjFamousSrcMeg, sAvgSubjUnfamiliarSrcMeg, ...
    'source_abs', 0);
% Process: Set comment
sDiffFamousSubjSrcMeg = bst_process('CallProcess', 'process_set_comment', sDiffFamousSubjSrcMeg, [], ...
    'tag',     'Famous - Unfamiliar | MEG', ...
    'isindex', 0);



%% ===== SUBJECT AVERAGES: FILTER AND NORMALIZE =======================================
%  ====================================================================================

% ===== FILTER: MEG+EEG =====
% Process: Select data files in: */Intra
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Low-pass:32Hz
sAvgData = bst_process('CallProcess', 'process_bandpass', sAvgData, [], ...
    'sensortypes', 'MEG, EEG', ...
    'highpass',    0, ...
    'lowpass',     32, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Extract time: [-202ms,900ms]
sAvgData = bst_process('CallProcess', 'process_extract_time', sAvgData, [], ...
    'timewindow', [-0.2018, 0.9], ...
    'overwrite',  1);

% ===== FILTER AND NORMALIZE: SOURCES =====
% Process: Select source files in: */Intra
sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Low-pass:32Hz
sAvgSrc = bst_process('CallProcess', 'process_bandpass', sAvgSrc, [], ...
    'highpass',    0, ...
    'lowpass',     32, ...
    'attenuation', 'strict', ...  % 60dB
    'mirror',      0, ...
    'overwrite',   1);
% Process: Extract time: [-202ms,900ms]
sAvgSrc = bst_process('CallProcess', 'process_extract_time', sAvgSrc, [], ...
    'timewindow', [-0.2018, 0.9], ...
    'overwrite',  1);
% Process: Z-score transformation: [-202ms,-5ms]
sAvgSrc = bst_process('CallProcess', 'process_baseline_norm', sAvgSrc, [], ...
    'baseline',   [-0.2018, -0.005], ...
    'source_abs', 0, ...
    'method',     'zscore', ...  % Z-score transformation:    x_std = (x - &mu;) / &sigma;
    'overwrite',  1);

% ===== NORMALIZE: TIMEFREQ =====
% Process: Select time-frequency files in: */Intra
sAvgTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Event-related perturbation (ERS/ERD): [-200ms,-4ms]
sAvgTf = bst_process('CallProcess', 'process_baseline_norm', sAvgTf, [], ...
    'baseline',  [-0.2, -0.004], ...
    'method',    'ersd', ...  % Event-related perturbation (ERS/ERD):    x_std = (x - &mu;) / &mu; * 100
    'overwrite', 1);

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_1prepare.html']));
end



%% ===== SUBJECT AVERAGES: SCREEN CAPTURES ============================================
%  ====================================================================================
% Start a new report
bst_report('Start');
% Set colormap: local color scale
bst_colormaps('SetMaxMode', 'stat2', 'local');

% ===== FACES: RECORDINGS ======
% Process: Select data files in: */*/WAvg: Avg: Faces |
sAvgFaces = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'tag',           'WAvg: Avg: Faces |', ...
    'includeintra',  1);
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ...
    'target',         1, ...  % Sensor/MRI registration
    'modality',       1);     % MEG
% Process: Snapshot: Recordings time series
bst_process('CallProcess', 'process_snapshot', sAvgFaces, [], ...
    'target',         5, ...  % Recordings time series
    'modality',       4, ...  % EEG
    'time',           0.1055);

% ===== FACES: SOURCES MEG ======
% Process: Select source files in: */*/WAvg: Avg: Faces | MEG
sAvgFacesSrcMeg = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Faces | MEG', ...
    'includeintra',  1);
% Process: Snapshot: Sources (one time)
bst_process('CallProcess', 'process_snapshot', sAvgFacesSrcMeg, [], ...
    'target',         8, ...  % Sources (one time)
    'orient',         6, ...  % back
    'time',           0.1055, ...
    'threshold',      10);

% ===== FACES: TIME-FREQ EEG ======
% Process: Select time-frequency files in: */*/WAvg: Avg: Faces
sAvgFacesTf = bst_process('CallProcess', 'process_select_files_timefreq', [], [], ...
    'subjectname',   'All', ...
    'tag',           'WAvg: Avg: Faces', ...
    'includeintra',  1);
% Process: Snapshot: Time-frequency maps
bst_process('CallProcess', 'process_snapshot', sAvgFacesTf, [], ...
    'target',         14, ...  % Time-frequency maps
    'time',           0.1055, ...
    'rowname',        'EEG070');

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_2snapshots.html']));
end



%% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG =====================================
%  ====================================================================================
% Start a new report
bst_report('Start');
% Set colormap: global color scale
bst_colormaps('SetMaxMode', 'meg', 'global');
bst_colormaps('SetMaxMode', 'eeg', 'global');
% Set display properties for statistics: p<0.05, FDR-corrected
StatThreshOptions = bst_get('StatThreshOptions');
StatThreshOptions.pThreshold = 0.05;
StatThreshOptions.Correction = 'fdr';
StatThreshOptions.Control    = [1 2 3];
bst_set('StatThreshOptions', StatThreshOptions);

% ===== GRAND AVERAGES =====
% Process: Select data files in: Group_analysis/Intra
sAvgData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   GroupSubjectName, ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Take screen captures
DataScreenCapture(sAvgData(1), 1, 1, '');
DataScreenCapture(sAvgData(2), 1, 1, '');
DataScreenCapture(sAvgData(3), 1, 1, '');
DataScreenCapture(sAvgData(4), 1, 1, '');
% Define function for screen captures
function DataScreenCapture(sFiles, isMEG, isEEG, Comment)
    if isMEG
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         5, ...  % Recordings time series
            'modality',       3, ...  % MEG (Magnetometers)
            'time',           0.1055, ...
            'Comment',        ['MEG MAG: ' Comment]);
        % Process: Snapshot: Recordings topography (contact sheet)
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         7, ...  % Recordings topography (contact sheet)
            'modality',       3, ...  % MEG (Magnetometers)
            'contact_time',   [0.050, 0.300], ...
            'contact_nimage', 6, ...
            'Comment',        ['MEG MAG: ' Comment]);
    end
    if isEEG
        % Process: Snapshot: Recordings time series
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         5, ...  % Recordings time series
            'modality',       4, ...  % EEG
            'time',           0.1055, ...
            'Comment',        ['EEG: ' Comment]);
        % Process: Snapshot: Recordings topography (contact sheet)
        bst_process('CallProcess', 'process_snapshot', sFiles, [], ...
            'target',         7, ...  % Recordings topography (contact sheet)
            'modality',       4, ...  % EEG
            'contact_time',   [0.050, 0.300], ...
            'contact_nimage', 6, ...
            'Comment',        ['EEG: ' Comment]);
    end
end

% ===== FACES-SCRAMBLED: MEAN DIFF =====
% Process: Select data files in: */DirIntra/"Avg: Faces | "
sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Faces | ', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)"
sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Scrambled (6 files)', ...
    'includeintra',  1);
% Process: Difference of means [mean]
sDiffFacesData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Set comment
sDiffFacesData = bst_process('CallProcess', 'process_set_comment', sDiffFacesData, [], ...
    'tag',     'Faces - Scrambled', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sDiffFacesData, 1, 1, sDiffFacesData.Comment);

% ===== FACES-SCRAMBLED: PARAMETRIC T-TEST =====
% Process: Select data files in: */DirIntra/"Avg: Faces | "
sSubjAvgFacesData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Faces | ', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Scrambled (6 files)"
sSubjAvgScrambledData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Scrambled (6 files)', ...
    'includeintra',  1);
% Process: t-test paired [ALL]          H0:(A=B), H1:(A<>B)
sStatParamFacesData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'timewindow',    [], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
    'tail',          'two');  % Two-tailed
% Process: Set comment
sStatParamFacesData = bst_process('CallProcess', 'process_set_comment', sStatParamFacesData, [], ...
    'tag',     'Faces - Scrambled: Parametric t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatParamFacesData, 1, 1, sStatParamFacesData.Comment);

% ===== FACES-SCRAMBLED: NON-PARAMETRIC T-TEST =====
% Process: Perm t-test paired [All]          H0:(A=B), H1:(A<>B)
sStatPermFacesData = bst_process('CallProcess', 'process_test_permutation2p', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'timewindow',     [], ...
    'sensortypes',    '', ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgrow',         0, ...
    'iszerobad',      0, ...
    'Comment',        '', ...
    'test_type',      'ttest_paired', ...  % Paired Student's t-test T = mean(A-B) / std(A-B) * sqrt(n)
    'randomizations', 1000, ...
    'tail',           'two');  % Two-tailed
% Process: Set comment
sStatPermFacesData = bst_process('CallProcess', 'process_set_comment', sStatPermFacesData, [], ...
    'tag',     'Faces - Scrambled: Permutation t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatPermFacesData, 1, 1, sStatPermFacesData.Comment);

% ===== FACES-SCRAMBLED: CLUSTER T-TEST =====
% Process: FT t-test paired cluster [all EEG]          H0:(A=B), H1:(A<>B)
sStatClustFacesData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFacesData, sSubjAvgScrambledData, ...
    'sensortypes',    'EEG', ...
    'timewindow',     [], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  2, ...  % Paired t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...  % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% If FieldTrip is available and results were returned
if ~isempty(sStatClustFacesData)
    % Process: Set comment
    sStatClustFacesData = bst_process('CallProcess', 'process_set_comment', sStatClustFacesData, [], ...
        'tag',     'Faces - Scrambled: Cluster t-test EEG', ...
        'isindex', 0);
    % Take screen captures
    DataScreenCapture(sStatClustFacesData, 0, 1, sStatClustFacesData.Comment);
end


%% ===== GROUP ANALYSIS: FACES-SCRAMBLED: MEG/EEG =====================================
%  ====================================================================================
% ===== FAMOUS-UNFAMILIAR: MEAN DIFF =====
% Process: Select data files in: */DirIntra/"Avg: Famous | "
sSubjAvgFamousData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'Avg: Famous (6 files)', ...
    'includeintra',  1);
% Process: Select data files in: */DirIntra/"WAvg: Avg: Unfamiliar (6 files)"
sSubjAvgUnfamiliarData = bst_process('CallProcess', 'process_select_files_data', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'tag',           'WAvg: Avg: Unfamiliar (6 files)', ...
    'includeintra',  1);
% Process: Difference of means [mean]
sDiffFamousData = bst_process('CallProcess', 'process_diff_mean', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'avg_func', 1, ...  % Arithmetic average mean(A) - mean(B)
    'weighted', 0);
% Process: Set comment
sDiffFamousData = bst_process('CallProcess', 'process_set_comment', sDiffFamousData, [], ...
    'tag',     'Famous - Unfamiliar', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sDiffFamousData, 1, 1, sDiffFamousData.Comment);

% ===== FAMOUS-UNFAMILIAR: PARAMETRIC T-TEST =====
% Process: t-test paired [All]          H0:(A=B), H1:(A<>B)
sStatParamFamousData = bst_process('CallProcess', 'process_test_parametric2p', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'timewindow',    [], ...
    'sensortypes',   '', ...
    'isabs',         0, ...
    'avgtime',       0, ...
    'avgrow',        0, ...
    'Comment',       '', ...
    'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
    'tail',          'two');  % Two-tailed
% Process: Set comment
sStatParamFamousData = bst_process('CallProcess', 'process_set_comment', sStatParamFamousData, [], ...
    'tag',     'Famous - Unfamiliar: Parametric t-test', ...
    'isindex', 0);
% Take screen captures
DataScreenCapture(sStatParamFamousData, 1, 1, sStatParamFamousData.Comment);

% ===== FAMOUS-UNFAMILIAR: CLUSTER T-TEST =====
% Process: FT t-test paired cluster [all EEG]          H0:(A=B), H1:(A<>B)
sStatClustFamousData = bst_process('CallProcess', 'process_ft_timelockstatistics', sSubjAvgFamousData, sSubjAvgUnfamiliarData, ...
    'sensortypes',    'EEG', ...
    'timewindow',     [], ...
    'isabs',          0, ...
    'avgtime',        0, ...
    'avgchan',        0, ...
    'randomizations', 1000, ...
    'statistictype',  2, ...  % Paired t-test
    'tail',           'two', ...  % Two-tailed
    'correctiontype', 2, ...  % cluster
    'minnbchan',      0, ...
    'clusteralpha',   0.05);
% If FieldTrip is available and results were returned
if ~isempty(sStatClustFacesData)
    % Process: Set comment
    sStatClustFamousData = bst_process('CallProcess', 'process_set_comment', sStatClustFamousData, [], ...
        'tag',     'Famous - Unfamiliar: Cluster t-test EEG', ...
        'isindex', 0);
    % Take screen captures
    DataScreenCapture(sStatClustFamousData, 0, 1, sStatClustFamousData.Comment);
end

% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_3meeg.html']));
end



%% ===== FACES-SCRAMBLED: SOURCES (MEG) ===============================================
%  ====================================================================================
% Start a new report
bst_report('Start');

% ===== PROJECT SOURCES =====
% Process: Select source files in: */Intra
sAvgSrc = bst_process('CallProcess', 'process_select_files_results', [], [], ...
    'subjectname',   'All', ...
    'condition',     DirIntra, ...
    'includeintra',  1);
% Process: Absolute values / Norm for unconstrained
if isUnconstrained
    % Process: Unconstrained to flat map
    sAvgSrc = bst_process('CallProcess', 'process_source_flat', sAvgSrc, [], ...
        'method', 1);  % Norm: sqrt(x^2+y^2+z^2)
else
    sAvgSrc = bst_process('CallProcess', 'process_absolute', sAvgSrc, [], ...
        'overwrite', 1);
end
% Process: Project on default anatomy
sProjSrc = bst_process('CallProcess', 'process_project_sources', sAvgSrc, [], ...
    'headmodeltype', 'surface');  % Cortex surface
% Process: Spatial smoothing (3.00,abs)
sProjSrc = bst_process('CallProcess', 'process_ssmooth_surfstat', sProjSrc, [], ...
    'fwhm',       3, ...
    'overwrite',  1, ...
    'source_abs', 1);

% ===== MOVE FILES =====
% File tags / destination folder
ListTags = {'Avg: Faces | EEG',                 'Faces_EEG'; ...
            'Avg: Faces | MEG',                 'Faces_MEG'; ...
            'Avg: Famous (6 files) | EEG',      'Famous_EEG'; ...
            'Avg: Famous (6 files) | MEG',      'Famous_MEG'; ...
            'Avg: Unfamiliar (6 files) | EEG',  'Unfamiliar_EEG'; ...
            'Avg: Unfamiliar (6 files) | MEG',  'Unfamiliar_MEG'; ...
            'Avg: Scrambled (6 files) | EEG',   'Scrambled_EEG'; ...
            'Avg: Scrambled (6 files) | MEG',   'Scrambled_MEG'; ...
            'Faces - Scrambled | EEG',          'Faces-Scrambled_EEG'; ...
            'Faces - Scrambled | MEG',          'Faces-Scrambled_MEG'; ...
            'Famous - Unfamiliar | EEG',        'Famous-Unfamiliar_EEG'; ...
            'Famous - Unfamiliar | MEG',        'Famous-Unfamiliar_MEG'};
% Move group by group
for i = 1:size(ListTags,1)
    % Find names of files to move
    iTag = find(~cellfun(@(c)isempty(strfind(c, ListTags{i,1})), {sProjSrc.Comment}));
    if isempty(iTag)
        continue;
    end
    % Process: Move files
    bst_process('CallProcess', 'process_movefile', {sProjSrc(iTag).FileName}, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      ListTags{i,2});
end


% ===== LOOP ON MODALITY =====
AllModalities = {'MEG','EEG'};
for iMod = 1:length(AllModalities)
    % Selected modality for this loop
    Mod = AllModalities{iMod};
    
    % ===== |FACES-SCRAMBLED|: MEAN DIFF ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_', Mod]);
    % Process: Weighted Average: By folder (grand average)
    sAbsDiffFaces = bst_process('CallProcess', 'process_average', sDiffFaces, [], ...
        'avgtype',         4, ...  % By folder (grand average)
        'avg_func',        1, ...  % Arithmetic average:  mean(x)
        'weighted',        0, ...
        'scalenormalized', 0);
    % Process: Set comment
    sAbsDiffFaces = bst_process('CallProcess', 'process_set_comment', sAbsDiffFaces, [], ...
        'tag',     ['mean(|Faces-Scrambled|) | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sAbsDiffFaces = bst_process('CallProcess', 'process_movefile', sAbsDiffFaces, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Colormap: Custom max: [-10,+10] Z
    bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces-Scrambled|)']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sAbsDiffFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces-Scrambled|)']);

    % ===== |FACES-SCRAMBLED|: PARAMETRIC CHI2-TEST ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_' Mod]);
    % Process: Chi2-test [all]          H0:(|Zi| = 0)
    sChiParamFaces = bst_process('CallProcess', 'process_test_parametric1', sDiffFaces, [], ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'chi2_onesample', ...  % One-sample Chi2 test     Zi~N(0,1), i=1..nQ = sum(|Zi|^2)     Q~Chi2(n)
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sChiParamFaces = bst_process('CallProcess', 'process_set_comment', sChiParamFaces, [], ...
        'tag',     ['|Faces-Scrambled|=0: Parametric Chi2 test | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sChiParamFaces = bst_process('CallProcess', 'process_movefile', sChiParamFaces, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces-Scrambled|=0: Parametric Chi2-test']);

    % ===== LOG(|FACES-SCRAMBLED|): PARAMETRIC CHI2-TEST ======
    % Process: Select source files in: Group_analysis/Faces-Scrambled_MOD
    sDiffFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces-Scrambled_' Mod]);
    % Process: Duplicate folders: Add tag "_log"
    sDiffFacesLog = bst_process('CallProcess', 'process_duplicate', sDiffFaces, [], ...
        'target', 2, ...  % Duplicate folders
        'tag',    '_log');
    % Process: Run Matlab command
    sDiffFacesLog = bst_process('CallProcess', 'process_matlab_eval', sDiffFacesLog, [], ...
        'matlab',    'Data = log(Data);', ...
        'overwrite', 1);
    % Process: Chi2-test [all]          H0:(|Zi| = 0)
    sChiParamFacesLog = bst_process('CallProcess', 'process_test_parametric1', sDiffFacesLog, [], ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'chi2_onesample', ...  % One-sample Chi2 test     Zi~N(0,1), i=1..nQ = sum(|Zi|^2)     Q~Chi2(n)
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sChiParamFacesLog = bst_process('CallProcess', 'process_set_comment', sChiParamFacesLog, [], ...
        'tag',     ['log(|Faces-Scrambled|)=0: Parametric Chi2 test | ' Mod], ...
        'isindex', 0);
    % Process: Move files
    sChiParamFacesLog = bst_process('CallProcess', 'process_movefile', sChiParamFacesLog, [], ...
        'subjectname', GroupSubjectName, ...
        'folder',      DirIntra);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sChiParamFacesLog, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': log(|Faces-Scrambled|)=0: Parametric Chi2-test']);


    % ===== |FACES|-|SCRAMBLED|: MEAN DIFF ======
    % Process: Select source files in: Group_analysis/Faces_MOD
    sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces_' Mod]);
    % Process: Select source files in: Group_analysis/Scrambled_MOD
    sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Scrambled_' Mod]);
    % Process: Difference of means [abs(mean)]
    sDiffAbsFaces = bst_process('CallProcess', 'process_diff_mean', sAvgFaces, sAvgScrambled, ...
        'avg_func',   2, ...  % Absolute value of average abs(mean(A)) - abs(mean(B))
        'weighted',   0);
    % Process: Set comment
    sDiffAbsFaces = bst_process('CallProcess', 'process_set_comment', sDiffAbsFaces, [], ...
        'tag',     ['mean(|Faces|)-mean(|Scrambled|) | ' Mod], ...
        'isindex', 0);
    % Colormap: Custom max: [-10,+10] Z
    bst_colormaps('SetMaxCustom', 'stat2', [], -10, 10);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces|)-mean(|Scrambled|)']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sDiffAbsFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'threshold',      30, ...
        'Comment',        [Mod ': mean(|Faces|)-mean(|Scrambled|)']);

    % ===== |FACES|-|SCRAMBLED|: PARAMETRIC T-TEST ======
    % Process: Select source files in: Group_analysis/Faces_MOD
    sAvgFaces = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Faces_' Mod]);
    % Process: Select source files in: Group_analysis/Scrambled_MOD
    sAvgScrambled = bst_process('CallProcess', 'process_select_files_results', [], [], ...
        'subjectname',   GroupSubjectName, ...
        'condition',     ['Scrambled_' Mod]);
    % Process: t-test paired [-202ms,900ms]          H0:(A=B), H1:(A<>B)
    sTtestParamFaces = bst_process('CallProcess', 'process_test_parametric2p', sAvgFaces, sAvgScrambled, ...
        'timewindow',    [], ...
        'scoutsel',      {}, ...
        'scoutfunc',     1, ...  % Mean
        'isnorm',        0, ...
        'avgtime',       0, ...
        'Comment',       '', ...
        'test_type',     'ttest_paired', ...  % Paired Student's t-test        (A-B)~N(m,v)t = mean(A-B) / std(A-B) * sqrt(n)      df=n-1
        'tail',          'two');  % Two-tailed
    % Process: Set comment
    sTtestParamFaces = bst_process('CallProcess', 'process_set_comment', sTtestParamFaces, [], ...
        'tag',     ['|Faces|=|Scrambled|: Parametric t-test | ' Mod], ...
        'isindex', 0);
    % Set colormap: global color scale
    bst_colormaps('SetMaxMode', 'stat2', 'global');
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         4, ...  % bottom
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces|-|Scrambled|: Parametric t-test']);
    % Process: Snapshot: Sources (contact sheet)
    bst_process('CallProcess', 'process_snapshot', sTtestParamFaces, [], ...
        'target',         9, ...  % Sources (contact sheet)
        'orient',         2, ...  % right
        'contact_time',   [0.05, 0.4], ...
        'contact_nimage', 16, ...
        'Comment',        [Mod ': |Faces|-|Scrambled|: Parametric t-test']);    
end



% Save report
ReportFile = bst_report('Save', []);
if ~isempty(reports_dir) && ~isempty(ReportFile)
    bst_report('Export', ReportFile, bst_fullfile(reports_dir, ['report_' ProtocolName '_4sources.html']));
end




end

