CMRO2_computation.tcsh

#!/bin/tcsh -xef

# to execute via bash: 
#   tcsh -xef XXX.tcsh 2>&1 | tee output.XXX
echo "auto-generated by afni_proc.py, Tue Nov 29 09:27:30 2022"
echo "(version 7.12, April 14, 2020)"
echo "execution started: `date`"

# =========================== block: init ============================

# assign output & input directory names
set suffixe = "classic"
set input_directory = "regression_$suffixe"
set input_directory = "regression_factorization"   # erase me
set output_directory = "CMRO2calc_$suffixe"
set input_dir = common_processing
set temp_dir = $input_dir
set dataFolder =  'XXXX/data'
set HOME = 'Github folder'
set SNRsavingdirectory = "$HOME/SNR_results_$suffixe"
set imageDirectory = "Images/$suffixe"

set subjects =  ('test')

# regression parametres
set pValthreshold = 1/100
set REML_flag = 0 

set beginning_offsetTR_BOLD = 1 # how many TRs must be eliminated from the beginning of BOLD processed data (mostly because of stable state in fMRI)
set ending_offsetTR_BOLD    = 0 # how many TRs must be eliminated from the ending of BOLD processed data (mostly due to filtering)
set beginning_offsetTR_PW   = 4 # how many TRs must be eliminated from the beginning of ASL processed data
set ending_offsetTR_PW      = 2 # how many TRs must be eliminated from the ending of ASL processed data


# Signal to noise analysis parametres
set methods = ( "PW" "SurSub" "MEICA")

# Multi-echos parametres
set echo_list = (`count -digits 2 1 4`)
set echo_times = ( 9.1 25.0 39.6 54.3 )
set fave_echo = '01'



mkdir -p $SNRsavingdirectory
mkdir -p $imageDirectory

# initialise SNR files
# ------ Hick ------
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tSNR_Hick_tlrc.txt
# CNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/CNR_Hick_tlrc.txt
# tCNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tCNR_Hick_tlrc.txt
# calculating Percent active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/percentVoxAct_Hick_tlrc.txt
# calculating mean t-scores in GM for active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/GM_tscores_Hick_tlrc.txt
# PVE
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/PVE_Hick_tlrc.txt
# mean CMRO2 in Left MC
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/CMRO2_TvsR_L_MC_tlrc.txt
# mean CMRO2 in right MC
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/CMRO2_TvsR_R_MC_tlrc.txt

# ------ tapping ------
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tSNR_tapping_tlrc.txt
# CNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/CNR_tapping_tlrc.txt
# tCNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tCNR_tapping_tlrc.txt
# calculating Percent active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/percentVoxAct_tapping_tlrc.txt
# calculating mean t-scores in GM for active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/GM_tscores_tapping_tlrc.txt
# PVE
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/PVE_tapping_tlrc.txt
# mean CMRO2 in Left MC
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/CMRO2_tapping_L_MC_tlrc.txt
# mean CMRO2 in right MC
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/CMRO2_tapping_R_MC_tlrc.txt
# Percent signal change due to Hypercapnia in left MC
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/PerctChg_tapping_L_MC_tlrc.txt
# Percent signal change due to Hypercapnia in right MC
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/PerctChg_tapping_R_MC_tlrc.txt

# HRF
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tSNR_HRF_tlrc.txt
# CNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/CNR_HRF_tlrc.txt
# tCNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tCNR_HRF_tlrc.txt
# calculating Percent active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/percentVoxAct_HRF_tlrc.txt
# calculating mean t-scores in GM for active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/GM_tscores_HRF_tlrc.txt
# PVE
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/PVE_HRF_tlrc.txt
# Percent signal change due to Hypercapnia
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/HC_HRF_tlrc.txt
# Percent signal change due to Hypercapnia in GM
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/HC_GM_HRF_tlrc.txt


# RRF
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tSNR_RRF_tlrc.txt
# CNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/CNR_RRF_tlrc.txt
# tCNR
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/tCNR_RRF_tlrc.txt
# calculating Percent active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/percentVoxAct_RRF_tlrc.txt
# calculating mean t-scores in GM for active voxels
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/GM_tscores_RRF_tlrc.txt
# PVE
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/PVE_RRF_tlrc.txt
# Percent signal change due to Hypercapnia
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/HC_RRF_tlrc.txt
# Percent signal change due to Hypercapnia in GM
echo "Subject MEC E2 PW SurSub MEICA" > $SNRsavingdirectory/HC_GM_RRF_tlrc.txt



# Mean M in GM
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/mean_M_GM_tlrc.txt
# Mean M in L_MC
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/mean_M_L_MC_tlrc.txt
# Mean M in R_MC
echo "Subject HRF RRF E2_HRF E2_RRF" > $SNRsavingdirectory/mean_M_R_MC_tlrc.txt

# =========================== auto block: setup ============================
# script setup

# take note of the AFNI version
afni -ver

# check that the current AFNI version is recent enough
afni_history -check_date 27 Jun 2019
if ( $status ) then
    echo "** this script requires newer AFNI binaries (than 27 Jun 2019)"
    echo "   (consider: @update.afni.binaries -defaults)"
    exit
endif

# =================== SUBJECT loop =====================
# the user may specify a single subject to run with
foreach subj ( $subjects )

    set output_dir = $dataFolder/${subj}/Analysis/


    cd $output_dir

    mkdir -p $output_directory
    
    # if (! -d ${input_directory}) then
    #     unzip -q ${input_directory}.zip
    # endif
    set brainmask = "$input_dir/mask_epi_anat.${subj}+tlrc.HEAD"
    set GMmask = "$input_dir/mask_GM_resam+tlrc.HEAD"
    # cp ${input_dir}/mask_epi_anat.$subj+tlrc* $dataFolder/$subj/Analysis/$output_directory
    # cp $input_dir/mask_GM_resam+tlrc.* ../$output_dir/

    set TR01 = "`3dinfo -nt $input_dir/pb05.$subj.r01.blur+tlrc`"
    set TR02 = "`3dinfo -nt $input_dir/pb05.$subj.r02.blur+tlrc`"

    # suppressing 4 TRs at beginning and 2TRs at end (eliminates artifacts (low values), from filtering)
    @ TR01end_BOLD = ($TR01 - 1 - $ending_offsetTR_BOLD )
    @ TR02end_BOLD = ($TR02 - 1 - $ending_offsetTR_BOLD )
    @ TR01end_PW = ($TR01 - 1 - $ending_offsetTR_PW )
    @ TR02end_PW = ($TR02 - 1 - $ending_offsetTR_PW )

    @ common_begin_offset = ( $beginning_offsetTR_PW - $beginning_offsetTR_BOLD )
    @ common_TR02end = ( $TR02end_PW - $beginning_offsetTR_BOLD )


    whereami -prefix "$output_directory/MNI_Glasser_HCP_v1.0.L_Primary_Motor_Cortex" -mask_atlas_region MNI_Glasser_HCP_v1.0:L_Primary_Motor_Cortex
    3dresample -master "$brainmask" -prefix "$output_directory/LMC_mask" -input "$output_directory/MNI_Glasser_HCP_v1.0.L_Primary_Motor_Cortex+tlrc"
    whereami -prefix "$output_directory/MNI_Glasser_HCP_v1.0.R_Primary_Motor_Cortex" -mask_atlas_region MNI_Glasser_HCP_v1.0:R_Primary_Motor_Cortex
    3dresample -master "$brainmask" -prefix "$output_directory/RMC_mask" -input "$output_directory/MNI_Glasser_HCP_v1.0.R_Primary_Motor_Cortex+tlrc"


    set RMC_mask = "${output_directory}/RMC_mask+tlrc"
    set LMC_mask = "${output_directory}/LMC_mask+tlrc"

  # ====================  BREATHHOLD CMRO2 calulation  ========================
   # ------------------  MEC tSNR / CNR / PVE calculations --------------------
    # calculating std of input signal (both RRF and HRF)
    echo $PWD
    3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.MEC1.all" "$input_directory/dataset.MEC.r01+tlrc"

    # HRF
    # calculating std of residuals
    3dTstat -nzstdev -prefix "$output_directory/rm.noise.MEC1.hrf" "$input_directory/errts_select_MECh+tlrc"
    # calculating -- tSNR --
    3dcalc -a "$input_directory/meanMECh+tlrc"                                             \
           -b "$output_directory/rm.noise.MEC1.hrf+tlrc"                                             \
           -c $GMmask                                \
           -expr 'c*a/b' -prefix "$output_directory/TSNR.MEC.hrf.${subj}"
    # calculating -- CNR --
    3dcalc -a "$output_directory/rm.signalSTD.MEC1.all+tlrc"                                            \
           -b "$output_directory/rm.noise.MEC1.hrf+tlrc"                                             \
           -c $GMmask                           \
           -expr 'c*a/b' -prefix "$output_directory/CNR.MEC.hrf.${subj}"
    # calculating tCNR for HRF
        3dcalc -a "$input_directory/best_coefMECh+tlrc.HEAD[0]"                                            \
               -b "$output_directory/rm.noise.MEC1.hrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/tCNR.MEC.hrf.${subj}"
        
    # calculating -- Percent active voxels --
    set tempVariable = (`p2dsetstat -quiet -inset "$input_directory/stats.MEC.0TR.HRF.${subj}_REML+tlrc.HEAD[7]" -pval $pValthreshold -bisided`)
    3dcalc -a "$input_directory/best_TtestsMECh+tlrc[0]" \
           -b $GMmask                               \
           -expr "astep(a,${tempVariable})*b" -prefix "$output_directory/ActVoxGMMask.MECh"
    # mean t-score in active voxels of GM      
    set GMtscoreMECh = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.MECh+tlrc.HEAD" "$input_directory/best_TtestsMECh+tlrc[0]"`
    # Percent of signal increase
    3dcalc -a "$input_directory/meanMECh+tlrc" \
           -b "$input_directory/best_coefMECh+tlrc.HEAD[0]" \
           -expr '100*b/a' -prefixe "$output_directory/MEC_HC_HRF"

    # RRF
    # calculating std of residuals
    3dTstat -nzstdev -prefix "$output_directory/rm.noise.MEC1.rrf" "$input_directory/errts_select_MECr+tlrc"
    # calculating tSNR for HRF
    3dcalc -a "$input_directory/meanMECr+tlrc"                                             \
           -b "$output_directory/rm.noise.MEC1.rrf+tlrc"                                             \
           -c $GMmask                                \
           -expr 'c*a/b' -prefix "$output_directory/TSNR.MEC.rrf.${subj}"
    # calculating CNR for HRF
    3dcalc -a "$output_directory/rm.signalSTD.MEC1.all+tlrc"                                            \
           -b "$output_directory/rm.noise.MEC1.rrf+tlrc"                                             \
           -c $GMmask                           \
           -expr 'c*a/b' -prefix "$output_directory/CNR.MEC.rrf.${subj}"
    # calculating tCNR for HRF
        3dcalc -a "$input_directory/best_coefMECr+tlrc.HEAD[0]"                                            \
               -b "$output_directory/rm.noise.MEC1.rrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/tCNR.MEC.rrf.${subj}"
        
    
    # calculating Percent active voxels for HRF
    set tempVariable = (`p2dsetstat -quiet -inset "$input_directory/stats.MEC.0TR.RRF.${subj}_REML+tlrc.HEAD[7]" -pval $pValthreshold -bisided`)
    3dcalc -a "$input_directory/best_TtestsMECr+tlrc[0]" \
           -b $GMmask                                  \
           -expr "astep(a,${tempVariable})*b" -prefix "$output_directory/ActVoxGMMask.MECr"
    # mean t-score in active voxels of GM      
    set GMtscoreMECr = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.MECr+tlrc.HEAD" "$input_directory/best_coefMECr+tlrc.HEAD[1]"`
              
    # Percent of signal increase
    3dcalc -a "$input_directory/meanMECr+tlrc" \
           -b "$input_directory/best_coefMECr+tlrc.HEAD[0]" \
           -expr '100*b/a' -prefixe "$output_directory/MEC_HC_RRF"

   # ------------------  E2 tSNR / CNR / PVE calculations --------------------
    # calculating std of input signal (both RRF and HRF)
    3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.E21.all" "$input_directory/dataset.E2.r01+tlrc"

    # HRF
    # calculating std of residuals
    3dTstat -nzstdev -prefix "$output_directory/rm.noise.E21.hrf" "$input_directory/errts_select_E2h+tlrc"
    # calculating -- tSNR --
    3dcalc -a "$input_directory/meanE2h+tlrc"                                             \
           -b "$output_directory/rm.noise.E21.hrf+tlrc"                                             \
           -c $GMmask                                \
           -expr 'c*a/b' -prefix "$output_directory/TSNR.E2.hrf.${subj}"
    # calculating -- CNR --
    3dcalc -a "$output_directory/rm.signalSTD.E21.all+tlrc"                                            \
           -b "$output_directory/rm.noise.E21.hrf+tlrc"                                             \
           -c $GMmask                           \
           -expr 'c*a/b' -prefix "$output_directory/CNR.E2.hrf.${subj}"
    # calculating tCNR for HRF
        3dcalc -a "$input_directory/best_coefE2h+tlrc.HEAD[0]"                                            \
               -b "$output_directory/rm.noise.E21.hrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/tCNR.E2.hrf.${subj}"
        
    
    # calculating -- Percent active voxels --
    set tempVariable = (`p2dsetstat -quiet -inset "$input_directory/stats.E2.0TR.HRF.${subj}_REML+tlrc.HEAD[7]" -pval $pValthreshold -bisided`)
    3dcalc -a "$input_directory/best_TtestsE2h+tlrc[0]" \
           -b $GMmask                               \
           -expr "astep(a,${tempVariable})*b" -prefix "$output_directory/ActVoxGMMask.E2h"
    # mean t-score in active voxels of GM      
    set GMtscoreE2h = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.E2h+tlrc.HEAD" "$input_directory/best_coefE2h+tlrc.HEAD[1]"`

    # Percent of signal increase
    3dcalc -a "$input_directory/meanE2h+tlrc" \
           -b "$input_directory/best_coefE2h+tlrc.HEAD[0]" \
           -expr '100*b/a' -prefixe "$output_directory/E2_HC_HRF"

    # RRF
    # calculating std of residuals
    3dTstat -nzstdev -prefix "$output_directory/rm.noise.E21.rrf" "$input_directory/errts_select_E2r+tlrc"
    # calculating tSNR for HRF
    3dcalc -a "$input_directory/meanE2r+tlrc"                                             \
           -b "$output_directory/rm.noise.E21.rrf+tlrc"                                             \
           -c $GMmask                                \
           -expr 'c*a/b' -prefix "$output_directory/TSNR.E2.rrf.${subj}"
    # calculating CNR for HRF
    3dcalc -a "$output_directory/rm.signalSTD.E21.all+tlrc"                                            \
           -b "$output_directory/rm.noise.E21.rrf+tlrc"                                             \
           -c $GMmask                           \
           -expr 'c*a/b' -prefix "$output_directory/CNR.E2.rrf.${subj}"
    # calculating tCNR for HRF
        3dcalc -a "$input_directory/best_coefE2r+tlrc.HEAD[0]"                                            \
               -b "$output_directory/rm.noise.E21.rrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/tCNR.E2.rrf.${subj}"
        
    # calculating Percent active voxels for HRF
    set tempVariable = (`p2dsetstat -quiet -inset "$input_directory/stats.E2.0TR.RRF.${subj}_REML+tlrc.HEAD[7]" -pval $pValthreshold -bisided`)
    3dcalc -a "$input_directory/best_TtestsE2r+tlrc[0]" \
           -b $GMmask                                  \
           -expr "astep(a,${tempVariable})*b" -prefix "$output_directory/ActVoxGMMask.E2r"
    # mean t-score in active voxels of GM      
    set GMtscoreE2r = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.E2r+tlrc.HEAD" "$input_directory/best_coefE2r+tlrc.HEAD[1]"`

    # Percent of signal increase
    3dcalc -a "$input_directory/meanE2r+tlrc" \
           -b "$input_directory/best_coefE2r+tlrc.HEAD[0]" \
           -expr '100*b/a' -prefixe "$output_directory/E2_HC_RRF"
   # ------------------  ASL tSNR / CNR / PVE calculations --------------------  
    # calculating means and std for signal
    foreach method ( $methods )
        # calculating std of input signal (both RRF and HRF)
        3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.${method}1.all" "$input_directory/dataset.${method}.r01+tlrc"

        # HRF
        # calculating std of residuals        
        3dTstat -nzstdev -prefix "$output_directory/rm.noise.${method}1.hrf" "$input_directory/errts_select_${method}h+tlrc"
        # calculating tSNR 
        3dcalc -a "$input_directory/mean${method}h+tlrc"                                             \
               -b "$output_directory/rm.noise.${method}1.hrf+tlrc"                                             \
               -c $GMmask                                \
               -expr 'c*a/b' -prefix "$output_directory/TSNR.${method}1.hrf.${subj}"
        # calculating CNR 
        3dcalc -a "$output_directory/rm.signalSTD.${method}1.all+tlrc"                                            \
               -b "$output_directory/rm.noise.${method}1.hrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/CNR.${method}1.hrf.${subj}"
        # calculating tCNR for HRF
        3dcalc -a "$input_directory/best_coef${method}h+tlrc.HEAD[0]"                                            \
               -b "$output_directory/rm.noise.${method}1.hrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/tCNR.${method}.hrf.${subj}"
        
        # calculating Percent active voxels
        set tempVariable = (`p2dsetstat -quiet -inset "$input_directory/stats.${method}.0TR.HRF.${subj}+tlrc.HEAD[7]" -pval $pValthreshold -bisided`)
        3dcalc -a "$input_directory/best_Ttests${method}h+tlrc[0]" \
               -b $GMmask                                       \
               -expr "astep(a,${tempVariable})*b" -prefix "$output_directory/ActVoxGMMask.${method}h"
        # mean t-score in active voxels of GM      
        set GMtscore${method}h = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.${method}h+tlrc.HEAD" "$input_directory/best_coef${method}h+tlrc.HEAD[1]"`

        # Percent of signal increase
        3dcalc -a "$input_directory/mean${method}h+tlrc" \
               -b "$input_directory/best_coef${method}h+tlrc.HEAD[0]" \
               -expr '100*b/a' -prefixe "$output_directory/${method}_HC_HRF"

        # RRF
        # calculating std and sos for residuals for RRF
        3dTstat -nzstdev -prefix "$output_directory/rm.noise.${method}1.rrf" "$input_directory/errts_select_${method}r+tlrc.HEAD"
        # calculating tSNR for HRF
        3dcalc -a "$input_directory/mean${method}r+tlrc"                                             \
               -b "$output_directory/rm.noise.${method}1.rrf+tlrc"                                             \
               -c $GMmask                                \
               -expr 'c*a/b' -prefix "$output_directory/TSNR.${method}.rrf.${subj}"
        # calculating CNR for HRF
        3dcalc -a "$output_directory/rm.signalSTD.${method}1.all+tlrc"                                            \
               -b "$output_directory/rm.noise.${method}1.rrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/CNR.${method}.rrf.${subj}"
        # calculating tCNR for HRF
        3dcalc -a "$input_directory/best_coef${method}r+tlrc.HEAD[0]"                                            \
               -b "$output_directory/rm.noise.${method}1.rrf+tlrc"                                             \
               -c $GMmask                           \
               -expr 'c*a/b' -prefix "$output_directory/tCNR.${method}.rrf.${subj}"
        
        # calculating Percent active voxels for RRF
        set tempVariable = (`p2dsetstat -quiet -inset "$input_directory/stats.${method}.0TR.RRF.${subj}+tlrc.HEAD[7]" -pval $pValthreshold -bisided`)
        3dcalc -a "$input_directory/best_Ttests${method}r+tlrc[0]" \
               -b $GMmask                               \
               -expr "astep(a,${tempVariable})*b" -prefix "$output_directory/ActVoxGMMask.${method}r"
        # mean t-score in active voxels of GM      
        set GMtscore${method}r = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.${method}r+tlrc.HEAD" "$input_directory/best_coef${method}r+tlrc.HEAD[1]"`

        # Percent of signal increase
        3dcalc -a "$input_directory/mean${method}r+tlrc" \
               -b "$input_directory/best_coef${method}r+tlrc.HEAD[0]" \
               -expr '100*b/a' -prefixe "$output_directory/${method}_HC_RRF"
    end

   # --------- compiling results in one txt file for copy paste in csv ---------
    # RRF
    3dTcat -output $output_directory/tSNRall.RRF \
        "$output_directory/TSNR.MEC.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.E2.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.PW.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.SurSub.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.MEICA.rrf.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tSNRall.RRF+tlrc.HEAD"`                          >> $SNRsavingdirectory/tSNR_RRF_tlrc.txt

    3dTcat -output $output_directory/CNRall.RRF \
        "$output_directory/CNR.MEC.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.E2.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.PW.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.SurSub.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.MEICA.rrf.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/CNRall.RRF+tlrc.HEAD"`                           >> $SNRsavingdirectory/CNR_RRF_tlrc.txt

    3dTcat -output $output_directory/tCNRall.RRF \
        "$output_directory/tCNR.MEC.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.E2.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.PW.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.SurSub.rrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.MEICA.rrf.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tCNRall.RRF+tlrc.HEAD"`                           >> $SNRsavingdirectory/tCNR_RRF_tlrc.txt

    3dTcat -output $output_directory/ActVoxAll.RRF \
        "$output_directory/ActVoxGMMask.MECr+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.E2r+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.PWr+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.SurSubr+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.MEICAr+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/ActVoxAll.RRF+tlrc.HEAD"`                        >> $SNRsavingdirectory/percentVoxAct_RRF_tlrc.txt

    echo $subj $GMtscoreMECr $GMtscoreE2r $GMtscorePWr $GMtscoreSurSubr $GMtscoreMEICAr                             >> $SNRsavingdirectory/GM_tscores_RRF_tlrc.txt

    3dTcat -output $output_directory/PVEall.RRF \
        "$input_directory/best_R2MECr+tlrc.HEAD" \
        "$input_directory/best_R2E2r+tlrc.HEAD" \
        "$input_directory/best_R2PWr+tlrc.HEAD" \
        "$input_directory/best_R2SurSubr+tlrc.HEAD" \
        "$input_directory/best_R2MEICAr+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/PVEall.RRF+tlrc.HEAD"`                           >> $SNRsavingdirectory/PVE_RRF_tlrc.txt



    3dTcat -output $output_directory/HC.RRF \
        "$output_directory/MEC_HC_RRF+tlrc.HEAD" \
        "$output_directory/E2_HC_RRF+tlrc.HEAD" \
        "$output_directory/PW_HC_RRF+tlrc.HEAD" \
        "$output_directory/SurSub_HC_RRF+tlrc.HEAD" \
        "$output_directory/MEICA_HC_RRF+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $brainmask "$output_directory/HC.RRF+tlrc.HEAD"`                            >> $SNRsavingdirectory/HC_RRF_tlrc.txt

    3dTcat -output $output_directory/HC_GM.RRF \
        "$output_directory/MEC_HC_RRF+tlrc.HEAD" \
        "$output_directory/E2_HC_RRF+tlrc.HEAD" \
        "$output_directory/PW_HC_RRF+tlrc.HEAD" \
        "$output_directory/SurSub_HC_RRF+tlrc.HEAD" \
        "$output_directory/MEICA_HC_RRF+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/HC.RRF+tlrc.HEAD"`                               >> $SNRsavingdirectory/HC_GM_RRF_tlrc.txt
    
# motor cortex doesn't concern BH tasks
    # 3dTcat -output $output_directory/HC_L_MC.RRF \
    #     "$output_directory/MEC_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/E2_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/PW_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/SurSub_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/MEICA_HC_RRF+tlrc.HEAD"
    # echo $subj `3dmaskave -q -mask $LMC_mask "$output_directory/HC_L_MC.RRF+tlrc.HEAD"`                        >> $SNRsavingdirectory/PerctChg_tapping_L_MC_RRF_tlrc.txt
    
    # 3dTcat -output $output_directory/HC_R_MC.RRF \
    #     "$output_directory/MEC_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/E2_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/PW_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/SurSub_HC_RRF+tlrc.HEAD" \
    #     "$output_directory/MEICA_HC_RRF+tlrc.HEAD"
    # echo $subj `3dmaskave -q -mask $RMC_mask "$output_directory/HC_R_MC.RRF+tlrc.HEAD"`                         >> $SNRsavingdirectory/PerctChg_tapping_R_MC_RRF_tlrc.txt    


    # HRF
    3dTcat -output $output_directory/tSNRall.HRF \
        "$output_directory/TSNR.MEC.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.E2.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.PW1.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.SurSub1.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.MEICA1.hrf.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tSNRall.HRF+tlrc.HEAD"`                           >> $SNRsavingdirectory/tSNR_HRF_tlrc.txt

    3dTcat -output $output_directory/CNRall.HRF \
        "$output_directory/CNR.MEC.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.E2.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.PW1.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.SurSub1.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.MEICA1.hrf.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/CNRall.HRF+tlrc.HEAD"`                            >> $SNRsavingdirectory/CNR_HRF_tlrc.txt

    3dTcat -output $output_directory/tCNRall.HRF \
        "$output_directory/tCNR.MEC.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.E2.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.PW1.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.SurSub1.hrf.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.MEICA1.hrf.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tCNRall.HRF+tlrc.HEAD"`                            >> $SNRsavingdirectory/tCNR_HRF_tlrc.txt

    3dTcat -output $output_directory/ActVoxAll.HRF \
        "$output_directory/ActVoxGMMask.MECh+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.E2h+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.PWh+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.SurSubh+tlrc.HEAD" \
        "$output_directory/ActVoxGMMask.MEICAh+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/ActVoxAll.HRF+tlrc.HEAD"`                         >> $SNRsavingdirectory/percentVoxAct_HRF_tlrc.txt

    echo $subj $GMtscoreMECh $GMtscoreE2h $GMtscorePWh $GMtscoreSurSubh $GMtscoreMEICAh                              >> $SNRsavingdirectory/GM_tscores_HRF_tlrc.txt

    3dTcat -output $output_directory/PVEall.HRF \
        "$input_directory/best_R2MECh+tlrc.HEAD" \
        "$input_directory/best_R2E2h+tlrc.HEAD" \
        "$input_directory/best_R2PWh+tlrc.HEAD" \
        "$input_directory/best_R2SurSubh+tlrc.HEAD" \
        "$input_directory/best_R2MEICAh+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/PVEall.HRF+tlrc.HEAD"`                            >> $SNRsavingdirectory/PVE_HRF_tlrc.txt




    3dTcat -output $output_directory/HC.HRF \
        "$output_directory/MEC_HC_HRF+tlrc.HEAD" \
        "$output_directory/E2_HC_HRF+tlrc.HEAD" \
        "$output_directory/PW_HC_HRF+tlrc.HEAD" \
        "$output_directory/SurSub_HC_HRF+tlrc.HEAD" \
        "$output_directory/MEICA_HC_HRF+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $brainmask "$output_directory/HC.HRF+tlrc.HEAD"`                            >> $SNRsavingdirectory/HC_HRF_tlrc.txt

    3dTcat -output $output_directory/HC_GM.HRF \
        "$output_directory/MEC_HC_HRF+tlrc.HEAD" \
        "$output_directory/E2_HC_HRF+tlrc.HEAD" \
        "$output_directory/PW_HC_HRF+tlrc.HEAD" \
        "$output_directory/SurSub_HC_HRF+tlrc.HEAD" \
        "$output_directory/MEICA_HC_HRF+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/HC.HRF+tlrc.HEAD"`                               >> $SNRsavingdirectory/HC_GM_HRF_tlrc.txt
    
# Motor cortex doesn't concern BH task, but rather the two others
   #  3dTcat -output $output_directory/HC_L_MC.HRF \
   #      "$output_directory/MEC_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/E2_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/PW_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/SurSub_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/MEICA_HC_HRF+tlrc.HEAD"
   #  echo $subj `3dmaskave -q -mask $LMC_mask "$output_directory/HC_L_MC.HRF+tlrc.HEAD"`                        >> $SNRsavingdirectory/PerctChg_L_MC_HRF_tlrc.txt

   # # this part should be for tapping, not BH    
   #  3dTcat -output $output_directory/HC_R_MC.HRF \
   #      "$output_directory/MEC_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/E2_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/PW_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/SurSub_HC_HRF+tlrc.HEAD" \
   #      "$output_directory/MEICA_HC_HRF+tlrc.HEAD"
   #  echo $subj `3dmaskave -q -mask $RMC_mask "$output_directory/HC_R_MC.HRF+tlrc.HEAD"`                         >> $SNRsavingdirectory/PerctChg_R_MC_HRF_tlrc.txt

  # ============================ M computation ==================================
   # calculating M for MEC and PW using HRF
    3dcalc -prefix "$output_directory/M_hrf" \
        -a "$input_directory/best_coefPWh+tlrc.HEAD[0]" \
        -b "$input_directory/meanPWh+tlrc" \
        -c "$input_directory/best_coefMECh+tlrc.HEAD[0]"\
        -d "$input_directory/meanMECh+tlrc"\
        -expr '(c/d) / (1 - ((a+b)/b)^(-8/10))' 
    3danisosmooth -prefix "$output_directory/blurredM_hrf" -mask $brainmask -3D -iters 2 "$output_directory/M_hrf+tlrc"

   # calculating M for MEC and PW using RRF
    3dcalc -prefix "$output_directory/M_rrf" \
        -a "$input_directory/best_coefPWr+tlrc.HEAD[0]" \
        -b "$input_directory/meanPWr+tlrc" \
        -c "$input_directory/best_coefMECr+tlrc.HEAD[0]"\
        -d "$input_directory/meanMECr+tlrc"\
        -expr '(c/d) / (1 - ((a+b)/b)^(-8/10))' 
    3danisosmooth -prefix "$output_directory/blurredM_rrf" -mask $brainmask -3D -iters 2 "$output_directory/M_rrf+tlrc"

   # calculating M for E2 and PW using HRF
    3dcalc -prefix "$output_directory/M_E2_hrf" \
        -a "$input_directory/best_coefPWh+tlrc.HEAD[0]" \
        -b "$input_directory/meanPWh+tlrc" \
        -c "$input_directory/best_coefE2h+tlrc.HEAD[0]"\
        -d "$input_directory/meanE2h+tlrc"\
        -expr '(c/d) / (1 - ((a+b)/b)^(-8/10))' 
    3danisosmooth -prefix "$output_directory/blurredM_E2_hrf" -mask $brainmask -3D -iters 2 "$output_directory/M_E2_hrf+tlrc"

   # calculating M for E2 and PW using RRF
    3dcalc -prefix "$output_directory/M_E2_rrf" \
        -a "$input_directory/best_coefPWr+tlrc.HEAD[0]" \
        -b "$input_directory/meanPWr+tlrc" \
        -c "$input_directory/best_coefE2r+tlrc.HEAD[0]"\
        -d "$input_directory/meanE2r+tlrc"\
        -expr '(c/d) / (1 - ((a+b)/b)^(-8/10))' 
    3danisosmooth -prefix "$output_directory/blurredM_E2_rrf" -mask $brainmask -3D -iters 2 "$output_directory/M_E2_rrf+tlrc"

    # Mean M in GM
    3dTcat -output $output_directory/M_GM \
        "$output_directory/M_hrf+tlrc.HEAD" \
        "$output_directory/M_rrf+tlrc.HEAD" \
        "$output_directory/M_E2_hrf+tlrc.HEAD" \
        "$output_directory/M_E2_rrf+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/M_GM+tlrc.HEAD"`                            >> $SNRsavingdirectory/mean_M_GM_tlrc.txt

    # Mean M in L_MC
    3dTcat -output $output_directory/M_LMC \
        "$output_directory/M_hrf+tlrc.HEAD" \
        "$output_directory/M_rrf+tlrc.HEAD" \
        "$output_directory/M_E2_hrf+tlrc.HEAD" \
        "$output_directory/M_E2_rrf+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $LMC_mask "$output_directory/M_LMC+tlrc.HEAD"`                         >> $SNRsavingdirectory/mean_M_L_MC_tlrc.txt

    # Mean M in R_MC
    3dTcat -output $output_directory/M_RMC \
        "$output_directory/M_hrf+tlrc.HEAD" \
        "$output_directory/M_rrf+tlrc.HEAD" \
        "$output_directory/M_E2_hrf+tlrc.HEAD" \
        "$output_directory/M_E2_rrf+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $RMC_mask "$output_directory/M_RMC+tlrc.HEAD"`                         >> $SNRsavingdirectory/mean_M_R_MC_tlrc.txt

  # ============================= Hick CMRO2 ====================================
   # ------------------------------ MEC ----------------------------------
    # calculating means and std
    3dTcat -prefix "$output_directory/rm.mean.signal.MEC2.all" "$input_directory/stats.MEC2.${subj}_REML+tlrc[2]"
    3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.MEC2.all" "$input_directory/dataset.MEC.r02+tlrc"
    3dTstat -nzstdev -prefix "$output_directory/rm.noise.MEC2.all" "$input_directory/errts.MEC2.${subj}_REML+tlrc.HEAD"
   
    # calculating tSNR
    3dcalc -a "$output_directory/rm.mean.signal.MEC2.all+tlrc"                                             \
           -b "$output_directory/rm.noise.MEC2.all+tlrc"                                             \
           -c $GMmask                                \
           -expr 'c*a/b' -prefix "$output_directory/TSNR.MEC2.all.${subj}"
   
    # calculating CNR
    3dcalc -a "$output_directory/rm.signalSTD.MEC2.all+tlrc"                                            \
           -b "$output_directory/rm.noise.MEC2.all+tlrc"                                             \
           -c $GMmask                           \
           -expr 'c*a/b' -prefix "$output_directory/CNR.MEC2.all.${subj}"

    # calculating tCNR
    3dcalc -a "$input_directory/stats.MEC2.${subj}_REML+tlrc[6]" \
           -c "$input_directory/stats.MEC2.${subj}_REML+tlrc[10]" \
           -d "$input_directory/stats.MEC2.${subj}_REML+tlrc[14]" \
           -b $GMmask                                                                  \
           -e "$output_directory/rm.noise.MEC2.all+tlrc"                               \
           -expr "b*(a + c + d)/(3*e)" -prefix"$output_directory/tCNR.MEC2.all.${subj}"

    # calculating Percent of variance explained
    3dcalc -a "$input_directory/stats.MEC2.${subj}_REML+tlrc[8]" \
           -c "$input_directory/stats.MEC2.${subj}_REML+tlrc[12]" \
           -d "$input_directory/stats.MEC2.${subj}_REML+tlrc[16]" \
           -b $GMmask                               \
           -expr "b*(a + c + d)/3" -prefix "$output_directory/mean_cond_R2.MEC2"

    # calculating Percent active voxels
    # We choose to calculate the percent of pVal on all three conditions in one mean mask
    3dcalc -a "$input_directory/pVal_stats.MEC2.${subj}_REML+tlrc[7]" \
           -b $GMmask                               \
           -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH0.MEC2"
    set H0ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH0.MEC2+tlrc"`
    3dcalc -a "$input_directory/pVal_stats.MEC2.${subj}_REML+tlrc[11]" \
           -b $GMmask                               \
           -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH1.MEC2"
    set H1ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH1.MEC2+tlrc"`
    3dcalc -a "$input_directory/pVal_stats.MEC2.${subj}_REML+tlrc[15]" \
           -b $GMmask                               \
           -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH2.MEC2"
    set H2ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH2.MEC2+tlrc"`
    set MEC2ActVox = `ccalc -eval "($H0ActPrect + $H1ActPrect + $H2ActPrect)/3"`

    # mean t-score in active voxels of GM      
    set H0GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH0.MEC2+tlrc" "$input_directory/stats.MEC2.${subj}_REML+tlrc[7]"`
    set H1GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH1.MEC2+tlrc" "$input_directory/stats.MEC2.${subj}_REML+tlrc[11]"`
    set H2GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH2.MEC2+tlrc" "$input_directory/stats.MEC2.${subj}_REML+tlrc[15]"`
    set GMtscoreMEC2 = `ccalc -eval "($H0GMtscore + $H1GMtscore + $H2GMtscore)/3"`

   # ------------------------------ E2 ----------------------------------
    # calculating means and std
    3dTcat -prefix "$output_directory/rm.mean.signal.E22.all" "$input_directory/stats.E22.${subj}_REML+tlrc[2]"
    3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.E22.all" "$input_directory/dataset.E2.r02+tlrc"
    3dTstat -nzstdev -prefix "$output_directory/rm.noise.E22.all" "$input_directory/errts.E22.${subj}_REML+tlrc.HEAD"
   
    # calculating tSNR
    3dcalc -a "$output_directory/rm.mean.signal.E22.all+tlrc"                                             \
           -b "$output_directory/rm.noise.E22.all+tlrc"                                             \
           -c $GMmask                                \
           -expr 'c*a/b' -prefix "$output_directory/TSNR.E22.all.${subj}"
   
    # calculating CNR
    3dcalc -a "$output_directory/rm.signalSTD.E22.all+tlrc"                                            \
           -b "$output_directory/rm.noise.E22.all+tlrc"                                             \
           -c $GMmask                           \
           -expr 'c*a/b' -prefix "$output_directory/CNR.E22.all.${subj}"

    # calculating tCNR
    3dcalc -a "$input_directory/stats.E22.${subj}_REML+tlrc[6]" \
           -c "$input_directory/stats.E22.${subj}_REML+tlrc[10]" \
           -d "$input_directory/stats.E22.${subj}_REML+tlrc[14]" \
           -b $GMmask                                                                  \
           -e "$output_directory/rm.noise.E22.all+tlrc"                               \
           -expr "b*(a + c + d)/(3*e)" -prefix"$output_directory/tCNR.E2.all.${subj}"


    # calculating Percent of variance explained
    3dcalc -a "$input_directory/stats.E22.${subj}_REML+tlrc[8]" \
           -c "$input_directory/stats.E22.${subj}_REML+tlrc[12]" \
           -d "$input_directory/stats.E22.${subj}_REML+tlrc[16]" \
           -b $GMmask                               \
           -expr "b*(a + c + d)/3" -prefix "$output_directory/mean_cond_R2.E22"

    # calculating Percent active voxels
    # We choose to calculate the percent of pVal on all three conditions in one mean mask
    3dcalc -a "$input_directory/pVal_stats.E22.${subj}_REML+tlrc[7]" \
           -b $GMmask                               \
           -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH0.E22"
    set H0ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH0.E22+tlrc"`
    3dcalc -a "$input_directory/pVal_stats.E22.${subj}_REML+tlrc[11]" \
           -b $GMmask                               \
           -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH1.E22"
    set H1ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH1.E22+tlrc"`
    3dcalc -a "$input_directory/pVal_stats.E22.${subj}_REML+tlrc[15]" \
           -b $GMmask                               \
           -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH2.E22"
    set H2ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH2.E22+tlrc"`
    set E22ActVox = `ccalc -eval "($H0ActPrect + $H1ActPrect + $H2ActPrect)/3"`

    # mean t-score in active voxels of GM      
    set H0GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH0.E22+tlrc" "$input_directory/stats.E22.${subj}_REML+tlrc[7]"`
    set H1GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH1.E22+tlrc" "$input_directory/stats.E22.${subj}_REML+tlrc[11]"`
    set H2GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH2.E22+tlrc" "$input_directory/stats.E22.${subj}_REML+tlrc[15]"`
    set GMtscoreE22 = `ccalc -eval "($H0GMtscore + $H1GMtscore + $H2GMtscore)/3"`

   # ------------------------------ PW ----------------------------------
    foreach method ( $methods )
        3dTcat -prefix "$output_directory/rm.mean.signal.${method}2.all" "$input_directory/stats.${method}2.$subj+tlrc[2]"
        3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.${method}2.all" "$input_directory/dataset.${method}.r02+tlrc"
        3dTstat -nzstdev -prefix "$output_directory/rm.noise.${method}2.all" "$input_directory/errts.${method}2.${subj}+tlrc"
        
        3dcalc -a "$output_directory/rm.mean.signal.${method}2.all+tlrc"                                             \
               -b  "$output_directory/rm.noise.${method}2.all+tlrc"                                              \
               -c $GMmask                               \
               -expr 'c*a/b' -prefix "$output_directory/TSNR.${method}2.all.${subj}"

        3dcalc -a "$output_directory/rm.signalSTD.${method}2.all+tlrc"                                             \
               -b "$output_directory/rm.noise.${method}2.all+tlrc"                                              \
               -c $GMmask                               \
               -expr 'c*a/b' -prefix "$output_directory/CNR.${method}2.all.${subj}"

        3dcalc -a "$input_directory/stats.${method}2.${subj}+tlrc[6]" \
               -c "$input_directory/stats.${method}2.${subj}+tlrc[10]" \
               -d "$input_directory/stats.${method}2.${subj}+tlrc[14]" \
               -b $GMmask                                                \
               -e "$output_directory/rm.noise.${method}.all+tlrc"                        \
               -expr "b*(a + c + d)/(3*e)" -prefix "$output_directory/tCNR.${method}2.all.${subj}"

       # PVE
        3dcalc -a "$input_directory/stats.${method}2.${subj}+tlrc[8]" \
               -c "$input_directory/stats.${method}2.${subj}+tlrc[12]" \
               -d "$input_directory/stats.${method}2.${subj}+tlrc[16]" \
               -b $GMmask                               \
               -expr "b*(a + c + d)/3" -prefix "$output_directory/mean_cond_R2.${method}2"

       # calculating Percent active voxels
        # We choose to calculate the percent of pVal on all three conditions in one mean mask       
        3dcalc -a "$input_directory/pVal_stats.${method}2.${subj}+tlrc[7]" \
               -b $GMmask                               \
               -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH0.${method}2"
        set H0ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH0.${method}2+tlrc"`
        3dcalc -a "$input_directory/pVal_stats.${method}2.${subj}+tlrc[11]" \
               -b $GMmask                               \
               -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH1.${method}2"
        set H1ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH1.${method}2+tlrc"`
        3dcalc -a "$input_directory/pVal_stats.${method}2.${subj}+tlrc[15]" \
               -b $GMmask                               \
               -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMaskH2.${method}2"
        set H2ActPrect = `3dmaskave -q -mask $GMmask "$output_directory/ActVoxGMMaskH2.${method}2+tlrc"`
        set ${method}2ActVox = `ccalc -eval "($H0ActPrect + $H1ActPrect + $H2ActPrect)/3"`

        # mean t-score in active voxels of GM      
        set H0GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH0.${method}2+tlrc" "$input_directory/stats.${method}2.${subj}+tlrc[7]"`
        set H1GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH1.${method}2+tlrc" "$input_directory/stats.${method}2.${subj}+tlrc[11]"`
        set H2GMtscore = `3dmaskave -q -mask "$output_directory/ActVoxGMMaskH2.${method}2+tlrc" "$input_directory/stats.${method}2.${subj}+tlrc[15]"`
        set GMtscore${method}2 = `ccalc -eval "($H0GMtscore + $H1GMtscore + $H2GMtscore)/3"`
    end

   # --------- compiling results in one txt file for copy paste in csv ---------
    3dTcat -output $output_directory/tSNRall.Hick \
        "$output_directory/TSNR.MEC2.all.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.E22.all.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.PW2.all.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.SurSub2.all.${subj}+tlrc.HEAD" \
        "$output_directory/TSNR.MEICA2.all.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tSNRall.Hick+tlrc.HEAD"`                           >> $SNRsavingdirectory/tSNR_Hick_tlrc.txt

    3dTcat -output $output_directory/CNRall.Hick \
        "$output_directory/CNR.MEC2.all.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.E22.all.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.PW2.all.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.SurSub2.all.${subj}+tlrc.HEAD" \
        "$output_directory/CNR.MEICA2.all.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/CNRall.Hick+tlrc.HEAD"`                            >> $SNRsavingdirectory/CNR_Hick_tlrc.txt

    3dTcat -output $output_directory/tCNRall.Hick \
        "$output_directory/tCNR.MEC2.all.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.E22.all.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.PW2.all.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.SurSub2.all.${subj}+tlrc.HEAD" \
        "$output_directory/tCNR.MEICA2.all.${subj}+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tCNRall.Hick+tlrc.HEAD"`                           >> $SNRsavingdirectory/tCNR_Hick_tlrc.txt

    echo $subj $MEC2ActVox $E22ActVox $PW2ActVox $SurSub2ActVox $MEICA2ActVox                                    >> $SNRsavingdirectory/percentVoxAct_Hick_tlrc.txt

    echo $subj $GMtscoreMEC2 $GMtscoreE22 $GMtscorePW2 $GMtscoreSurSub2 $GMtscoreMEICA2                          >> $SNRsavingdirectory/GM_tscores_Hick_tlrc.txt

    3dTcat -output $output_directory/PVEall.Hick \
        "$output_directory/mean_cond_R2.MEC2+tlrc.HEAD" \
        "$output_directory/mean_cond_R2.E22+tlrc.HEAD" \
        "$output_directory/mean_cond_R2.PW2+tlrc.HEAD" \
        "$output_directory/mean_cond_R2.SurSub2+tlrc.HEAD" \
        "$output_directory/mean_cond_R2.MEICA2+tlrc.HEAD"
    echo $subj `3dmaskave -q -mask $GMmask "$output_directory/PVEall.Hick+tlrc.HEAD"`                            >> $SNRsavingdirectory/PVE_Hick_tlrc.txt


  # ------------------------------ CMRO2 calculation ----------------------------------
   # CMRO2 calculation using MEC 
    # cond1 HRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond1_hrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[6]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.MEC2.${subj}_REML+tlrc.HEAD[6]"\
               -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
               -e "$output_directory/M_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond1_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H0' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_H0_HRF.bmp" -sDEVICE=bmp256
       
    # cond1 RRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond1_rrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[6]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.MEC2.${subj}_REML+tlrc.HEAD[6]"\
               -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
               -e "$output_directory/M_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond1_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H0' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_H0_RRF.bmp" -sDEVICE=bmp256
         
    # cond2 HRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond2_hrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[10]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.MEC2.${subj}_REML+tlrc.HEAD[10]"\
               -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
               -e "$output_directory/M_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond2_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H1' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_H1_HRF.bmp" -sDEVICE=bmp256
        
    # cond2 RRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond2_rrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[10]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.MEC2.${subj}_REML+tlrc.HEAD[10]"\
               -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
               -e "$output_directory/M_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond2_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H1' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_H1_RRF.bmp" -sDEVICE=bmp256
        
    # cond4 HRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond4_hrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[14]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.MEC2.${subj}_REML+tlrc.HEAD[14]"\
               -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
               -e "$output_directory/M_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond4_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H2' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_H2_HRF.bmp" -sDEVICE=bmp256

    # cond4 RRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond4_rrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[14]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.MEC2.${subj}_REML+tlrc.HEAD[14]"\
               -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
               -e "$output_directory/M_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond4_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H2' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_H2_RRF.bmp" -sDEVICE=bmp256
        
    # Task vs rest HRF
        3dcalc -prefix "$output_directory/CMRO2_TvsR_hrf" \
               -a "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[6]" \
               -b "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[2]" \
               -c "$input_directory/stats.MEC2_TvsR.${subj}_REML+tlrc.HEAD[6]"\
               -d "$input_directory/stats.MEC2_TvsR.${subj}_REML+tlrc.HEAD[2]"\
               -e "$output_directory/M_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_TvsR_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 Task versus rest' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_TvsR_HRF.bmp" -sDEVICE=bmp256
    # Task vs rest RRF
        3dcalc -prefix "$output_directory/CMRO2_TvsR_rrf" \
               -a "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[6]" \
               -b "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[2]" \
               -c "$input_directory/stats.MEC2_TvsR.${subj}_REML+tlrc.HEAD[6]"\
               -d "$input_directory/stats.MEC2_TvsR.${subj}_REML+tlrc.HEAD[2]"\
               -e "$output_directory/M_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_TvsR_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H2' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_TvsR_RRF.bmp" -sDEVICE=bmp256

   # CMRO2 calculation using E2 
    # cond1 HRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond1_E2_hrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[6]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.E22.${subj}_REML+tlrc.HEAD[6]"\
               -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
               -e "$output_directory/M_E2_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond1_E2_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H0' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_H0_HRF.bmp" -sDEVICE=bmp256
       
    # cond1 RRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond1_E2_rrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[6]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.E22.${subj}_REML+tlrc.HEAD[6]"\
               -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
               -e "$output_directory/M_E2_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond1_E2_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H0' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_H0_RRF.bmp" -sDEVICE=bmp256
         
    # cond2 HRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond2_E2_hrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[10]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.E22.${subj}_REML+tlrc.HEAD[10]"\
               -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
               -e "$output_directory/M_E2_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond2_E2_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H1' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_H1_HRF.bmp" -sDEVICE=bmp256
        
    # cond2 RRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond2_E2_rrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[10]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.E22.${subj}_REML+tlrc.HEAD[10]"\
               -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
               -e "$output_directory/M_E2_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond2_E2_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H1' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_H1_RRF.bmp" -sDEVICE=bmp256
        
    # cond4 HRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond4_E2_hrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[14]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.E22.${subj}_REML+tlrc.HEAD[14]"\
               -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
               -e "$output_directory/M_E2_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond4_E2_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H2' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_H2_HRF.bmp" -sDEVICE=bmp256

    # cond4 RRF
        3dcalc -prefix "$output_directory/CMRO2_hick_cond4_E2_rrf" \
               -a "$input_directory/stats.PW2.${subj}+tlrc.HEAD[14]" \
               -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
               -c "$input_directory/stats.E22.${subj}_REML+tlrc.HEAD[14]"\
               -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
               -e "$output_directory/M_E2_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_hick_cond4_E2_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H2' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_H2_RRF.bmp" -sDEVICE=bmp256
        
    # Task vs rest HRF
        3dcalc -prefix "$output_directory/CMRO2_TvsR_E2_hrf" \
               -a "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[6]" \
               -b "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[2]" \
               -c "$input_directory/stats.E22_TvsR.${subj}_REML+tlrc.HEAD[6]"\
               -d "$input_directory/stats.E22_TvsR.${subj}_REML+tlrc.HEAD[2]"\
               -e "$output_directory/M_E2_hrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_TvsR_E2_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 Task versus rest' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_TvsR_HRF.bmp" -sDEVICE=bmp256
    # Task vs rest RRF
        3dcalc -prefix "$output_directory/CMRO2_TvsR_E2_rrf" \
               -a "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[6]" \
               -b "$input_directory/stats.PW2_TvsR.${subj}+tlrc.HEAD[2]" \
               -c "$input_directory/stats.E22_TvsR.${subj}_REML+tlrc.HEAD[6]"\
               -d "$input_directory/stats.E22_TvsR.${subj}_REML+tlrc.HEAD[2]"\
               -e "$output_directory/M_E2_rrf+tlrc" \
               -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -mask $brainmask -nbin 500 -notitle -min -2 -max 3 "$output_directory/CMRO2_TvsR_E2_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 H2' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_E2_TvsR_RRF.bmp" -sDEVICE=bmp256

# CMRO2 times series calculation
        # CMRO2(t) HRF
        3dcalc -prefix "$output_directory/CMRO2_time_serie_Hick_TvsR_hrf" \
            -a "$input_directory/dataset.PW.r02+tlrc" \
            -b "$input_directory/stats.PW2.${subj}+tlrc[2]" \
            -c "$input_directory/dataset.MEC.r02+tlrc[$common_begin_offset..$common_TR02end]"\
            -d "$input_directory/stats.MEC2.${subj}_REML+tlrc[2]"\
            -e "$output_directory/M_hrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_time_serie_Hick_TvsR_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_hrf_R_MC.Hick_TvsR.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_time_serie_Hick_TvsR_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_hrf_L_MC.Hick_TvsR.1D"
        1dplot -ps -one -xlabel 'CMRO2 Hick task versus Rest' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_TaskVsRest.E2.1D" "$output_directory/CMRO2_t_hrf_L_MC.Hick_TvsR.1D" "$output_directory/CMRO2_t_hrf_R_MC.Hick_TvsR.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_TvsR_hrf_MC.bmp" -sDEVICE=bmp256
        # CMRO2(t) E2 HRF
        3dcalc -prefix "$output_directory/CMRO2_E2_time_serie_Hick_TvsR_hrf" \
            -a "$input_directory/dataset.PW.r02+tlrc" \
            -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
            -c "$input_directory/dataset.E2.r02+tlrc[$common_begin_offset..$common_TR02end]"\
            -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
            -e "$output_directory/M_E2_hrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_E2_time_serie_Hick_TvsR_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_hrf_R_MC.Hick_TvsR.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_E2_time_serie_Hick_TvsR_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_hrf_L_MC.Hick_TvsR.1D"
        1dplot -ps -one -xlabel 'CMRO2 Hick task versus Rest' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_TaskVsRest.E2.1D" "$output_directory/CMRO2_t_E2_hrf_L_MC.Hick_TvsR.1D" "$output_directory/CMRO2_t_E2_hrf_R_MC.Hick_TvsR.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_TvsR_E2_hrf_MC.bmp" -sDEVICE=bmp256

        # CMRO2(t) RRF
        3dcalc -prefix "$output_directory/CMRO2_time_serie_Hick_TvsR_rrf" \
            -a "$input_directory/dataset.PW.r02+tlrc" \
            -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
            -c "$input_directory/dataset.MEC.r02+tlrc[$common_begin_offset..$common_TR02end]"\
            -d "$output_directory/rm.mean.signal.MEC2.all+tlrc"\
            -e "$output_directory/M_rrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_time_serie_Hick_TvsR_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_rrf_R_MC.Hick_TvsR.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_time_serie_Hick_TvsR_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_rrf_L_MC.Hick_TvsR.1D"
        1dplot -ps -one -xlabel 'CMRO2 Hick task versus Rest' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_TaskVsRest.E2.1D" "$output_directory/CMRO2_t_rrf_L_MC.Hick_TvsR.1D" "$output_directory/CMRO2_t_rrf_R_MC.Hick_TvsR.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_TvsR_rrf_MC.bmp" -sDEVICE=bmp256

        # CMRO2(t) E2 RRF
        3dcalc -prefix "$output_directory/CMRO2_E2_time_serie_Hick_TvsR_rrf" \
            -a "$input_directory/dataset.PW.r02+tlrc" \
            -b "$output_directory/rm.mean.signal.PW2.all+tlrc" \
            -c "$input_directory/dataset.E2.r02+tlrc[$common_begin_offset..$common_TR02end]"\
            -d "$output_directory/rm.mean.signal.E22.all+tlrc"\
            -e "$output_directory/M_E2_rrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_E2_time_serie_Hick_TvsR_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_rrf_R_MC.Hick_TvsR.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_E2_time_serie_Hick_TvsR_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_rrf_L_MC.Hick_TvsR.1D"
        1dplot -ps -one -xlabel 'CMRO2 Hick task versus Rest' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_TaskVsRest.E2.1D" "$output_directory/CMRO2_t_E2_rrf_L_MC.Hick_TvsR.1D" "$output_directory/CMRO2_t_E2_rrf_R_MC.Hick_TvsR.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_E2_TvsR_rrf_MC.bmp" -sDEVICE=bmp256

        # mean CMRO2 in Left MC
        3dTcat -output $output_directory/M_LMC.TvsR \
            "$output_directory/CMRO2_TvsR_hrf+tlrc" \
            "$output_directory/CMRO2_TvsR_rrf+tlrc" \
            "$output_directory/CMRO2_TvsR_E2_hrf+tlrc" \
            "$output_directory/CMRO2_TvsR_E2_rrf+tlrc"
        echo $subj `3dmaskave -q -mask $LMC_mask "$output_directory/M_LMC.TvsR+tlrc.HEAD"`                           >> $SNRsavingdirectory/CMRO2_TvsR_L_MC_tlrc.txt

        # mean CMRO2 in right MC
        3dTcat -output $output_directory/M_RMC.TvsR \
            "$output_directory/CMRO2_TvsR_hrf+tlrc" \
            "$output_directory/CMRO2_TvsR_rrf+tlrc" \
            "$output_directory/CMRO2_TvsR_E2_hrf+tlrc" \
            "$output_directory/CMRO2_TvsR_E2_rrf+tlrc"
        echo $subj `3dmaskave -q -mask $RMC_mask "$output_directory/M_RMC.TvsR+tlrc.HEAD"`                           >> $SNRsavingdirectory/CMRO2_TvsR_R_MC_tlrc.txt

  # =========================== tapping CMRO2 ===============================
   if ( -d ../tapping ) then 
      # MEC analysis

        set TR03 = "`3dinfo -nt $input_dir/pb05.$subj.r03.blur+tlrc`"
        @ TR03end_BOLD = ($TR03 - 1 - $ending_offsetTR_BOLD )
        @ TR03end_PW = ($TR03 - 1 - $ending_offsetTR_PW )
        set tr03 = "`3dinfo -tr dataset.MEC.r03+tlrc`"
        @ common_TR03end = ( $TR03end_PW - $beginning_offsetTR_BOLD  )


       # chose to compute math on non-zero voxels
        3dTcat -prefix "$output_directory/rm.mean.signal.MEC3.all" "$input_directory/stats.MEC3.${subj}_REML+tlrc[2]"
        3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.MEC3.all" "$input_directory/dataset.MEC.r03+tlrc"
        3dTstat -nzstdev -prefix "$output_directory/rm.noise.MEC3.all" "$input_directory/errts.MEC3.${subj}_REML+tlrc.HEAD"
        
       # tSNR
        3dcalc -a "$output_directory/rm.mean.signal.MEC3.all+tlrc"                         \
               -b "$output_directory/rm.noise.MEC3.all+tlrc"                               \
               -c $GMmask                               \
               -expr 'c*a/b' -prefix "$output_directory/TSNR.MEC3.all.${subj}"

       # CNR
        3dcalc -a "$output_directory/rm.signalSTD.MEC3.all+tlrc"                                             \
               -b "$output_directory/rm.noise.MEC3.all+tlrc"                                              \
               -c $GMmask                               \
               -expr 'c*a/b' -prefix "$output_directory/CNR.MEC3.all.${subj}"
         
       # Percent of active voxels
        3dcalc -a "$input_directory/pVal_stats.MEC3.${subj}_REML+tlrc[7]" \
               -b $GMmask                               \
               -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMask.MEC3"
        # mean t-score in active voxels of GM      
        set GMtscoreMEC3 = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.MEC3+tlrc.HEAD" "$input_directory/stats.MEC3.${subj}_REML+tlrc[7]"`
        # Percent of signal increase
        3dcalc -a "$output_directory/rm.mean.signal.MEC3.all+tlrc" \
               -b "$input_directory/stats.MEC3.${subj}_REML+tlrc[6]" \
               -expr '100*b/a' -prefixe "$output_directory/MEC3_PerctChg"

      # E2 analysis
       # chose to compute math on non-zero voxels
        3dTcat -prefix "$output_directory/rm.mean.signal.E23.all" "$input_directory/stats.E23.${subj}_REML+tlrc[2]"
        3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.E23.all" "$input_directory/dataset.E2.r03+tlrc"
        3dTstat -nzstdev -prefix "$output_directory/rm.noise.E23.all" "$input_directory/errts.E23.${subj}_REML+tlrc.HEAD"
        
       # tSNR
        3dcalc -a "$output_directory/rm.mean.signal.E23.all+tlrc"                         \
               -b "$output_directory/rm.noise.E23.all+tlrc"                               \
               -c $GMmask                               \
               -expr 'c*a/b' -prefix "$output_directory/TSNR.E23.all.${subj}"

       # CNR
        3dcalc -a "$output_directory/rm.signalSTD.E23.all+tlrc"                                             \
               -b "$output_directory/rm.noise.E23.all+tlrc"                                              \
               -c $GMmask                               \
               -expr 'c*a/b' -prefix "$output_directory/CNR.E23.all.${subj}"
         
       # Percent of active voxels
        3dcalc -a "$input_directory/pVal_stats.E23.${subj}_REML+tlrc[7]" \
               -b $GMmask                               \
               -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMask.E23"

        # mean t-score in active voxels of GM      
        set GMtscoreE23 = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.E23+tlrc.HEAD" "$input_directory/stats.E23.${subj}_REML+tlrc[7]"`

        # Percent of signal increase
        3dcalc -a "$output_directory/rm.mean.signal.E23.all+tlrc" \
               -b "$input_directory/stats.E23.${subj}_REML+tlrc[6]" \
               -expr '100*b/a' -prefixe "$output_directory/E23_PerctChg"            

      # ASL analysis
        foreach method ( $methods )
            3dTcat -prefix "$output_directory/rm.mean.signal.${method}3.all" "$input_directory/stats.${method}3.$subj+tlrc[2]"
            3dTstat -nzstdev -prefix "$output_directory/rm.signalSTD.${method}3.all" "$input_directory/dataset.${method}.r03+tlrc"
            3dTstat -nzstdev -prefix "$output_directory/rm.noise.${method}3.all" "$input_directory/errts.${method}3.${subj}+tlrc"
            3dcalc -a "$output_directory/rm.mean.signal.${method}3.all+tlrc"                                             \
                   -b "$output_directory/rm.noise.${method}3.all+tlrc"                                              \
                   -c $GMmask                                \
                   -expr 'c*a/b' -prefix "$output_directory/TSNR.${method}3.all.${subj}"

            3dcalc -a "$output_directory/rm.signalSTD.${method}3.all+tlrc"                                             \
                   -b "$output_directory/rm.noise.${method}3.all+tlrc"                                              \
                   -c $GMmask                               \
                   -expr 'c*a/b' -prefix "$output_directory/CNR.${method}3.all.${subj}"           

            # Percent of active voxels
            3dcalc -a "$input_directory/pVal_stats.${method}3.${subj}+tlrc[7]" \
                   -b $GMmask  \
                   -expr "ispositive($pValthreshold-a)*b" -prefix "$output_directory/ActVoxGMMask.${method}3"

            # mean t-score in active voxels of GM      
            set GMtscore${method}3 = `3dmaskave -q -mask "$output_directory/ActVoxGMMask.${method}3+tlrc.HEAD" "$input_directory/stats.${method}3.${subj}+tlrc[7]"`
            # Percent of signal increase
            3dcalc -a "$output_directory/rm.mean.signal.${method}3.all+tlrc" \
                   -b "$input_directory/stats.${method}3.${subj}+tlrc[6]" \
                   -expr '100*b/a' -prefixe "$output_directory/${method}3_PerctChg"
        end

      # compiling results in one txt file for copy paste in csv
        3dTcat -output $output_directory/tSNRall.tapping \
            "$output_directory/TSNR.MEC3.all.${subj}+tlrc.HEAD" \
            "$output_directory/TSNR.E23.all.${subj}+tlrc.HEAD" \
            "$output_directory/TSNR.PW3.all.${subj}+tlrc.HEAD" \
            "$output_directory/TSNR.SurSub3.all.${subj}+tlrc.HEAD" \
            "$output_directory/TSNR.MEICA3.all.${subj}+tlrc.HEAD"
        echo $subj `3dmaskave -q -mask $GMmask "$output_directory/tSNRall.tapping+tlrc.HEAD"`                          >> $SNRsavingdirectory/tSNR_tapping_tlrc.txt

        3dTcat -output $output_directory/CNRall.tapping \
            "$output_directory/CNR.MEC3.all.${subj}+tlrc.HEAD" \
            "$output_directory/CNR.E23.all.${subj}+tlrc.HEAD" \
            "$output_directory/CNR.PW3.all.${subj}+tlrc.HEAD" \
            "$output_directory/CNR.SurSub3.all.${subj}+tlrc.HEAD" \
            "$output_directory/CNR.MEICA3.all.${subj}+tlrc.HEAD"
        echo $subj `3dmaskave -q -mask $GMmask "$output_directory/CNRall.tapping+tlrc.HEAD"`                           >> $SNRsavingdirectory/CNR_tapping_tlrc.txt

        3dTcat -output $output_directory/ActVoxAll.tapping \
            "$output_directory/ActVoxGMMask.MEC3+tlrc.HEAD" \
            "$output_directory/ActVoxGMMask.E23+tlrc.HEAD" \
            "$output_directory/ActVoxGMMask.PW3+tlrc.HEAD" \
            "$output_directory/ActVoxGMMask.SurSub3+tlrc.HEAD" \
            "$output_directory/ActVoxGMMask.MEICA3+tlrc.HEAD"
        echo $subj `3dmaskave -q -mask $GMmask "$output_directory/ActVoxAll.tapping+tlrc.HEAD"`                        >> $SNRsavingdirectory/percentVoxAct_tapping_tlrc.txt

        echo $subj $GMtscoreMEC3 $GMtscoreE23 $GMtscorePW3 $GMtscoreSurSub3 $GMtscoreMEICA3                                 >> $SNRsavingdirectory/GM_tscores_tapping_tlrc.txt

        3dTcat -output $output_directory/PVEall.tapping \
            "$input_directory/stats.MEC3.${subj}_REML+tlrc.HEAD[8]" \
            "$input_directory/stats.E23.${subj}_REML+tlrc.HEAD[8]" \
            "$input_directory/stats.PW3.${subj}+tlrc.HEAD[8]" \
            "$input_directory/stats.SurSub3.${subj}+tlrc.HEAD[8]" \
            "$input_directory/stats.MEICA3.${subj}+tlrc.HEAD[8]"
        echo $subj `3dmaskave -q -mask $GMmask "$output_directory/PVEall.tapping+tlrc.HEAD"`                           >> $SNRsavingdirectory/PVE_tapping_tlrc.txt

       # this part should be for tapping, not BH
        3dTcat -output $output_directory/HC_L_MC.HRF \
            "$output_directory/MEC3_PerctChg+tlrc.HEAD" \
            "$output_directory/E23_PerctChg+tlrc.HEAD" \
            "$output_directory/PW3_PerctChg+tlrc.HEAD" \
            "$output_directory/SurSub3_PerctChg+tlrc.HEAD" \
            "$output_directory/MEICA3_PerctChg+tlrc.HEAD"
        echo $subj `3dmaskave -q -mask $LMC_mask "$output_directory/HC_L_MC.HRF+tlrc.HEAD"`                           >> $SNRsavingdirectory/PerctChg_L_MC_tlrc.txt

       # this part should be for tapping, not BH    
        3dTcat -output $output_directory/HC_R_MC.HRF \
            "$output_directory/MEC3_PerctChg+tlrc.HEAD" \
            "$output_directory/E23_PerctChg+tlrc.HEAD" \
            "$output_directory/PW3_PerctChg+tlrc.HEAD" \
            "$output_directory/SurSub3_PerctChg+tlrc.HEAD" \
            "$output_directory/MEICA3_PerctChg+tlrc.HEAD"
        echo $subj `3dmaskave -q -mask $RMC_mask "$output_directory/HC_R_MC.HRF+tlrc.HEAD"`                           >> $SNRsavingdirectory/PerctChg_R_MC_tlrc.txt

       # CMRO2 calculations
        # CMRO2 calculation HRF using MEC
        3dcalc -prefix "$output_directory/CMRO2_tapping_hrf" \
            -a "$input_directory/stats.PW3.${subj}+tlrc.HEAD[6]" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/stats.MEC3.${subj}_REML+tlrc.HEAD[6]"\
            -d "$output_directory/rm.mean.signal.MEC3.all+tlrc"\
            -e "$output_directory/M_hrf+tlrc" \
            -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -omit 0 -mask $brainmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_hrf+tlrc" > "hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_HRF.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $GMmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_hrf+tlrc" > "$output_directory/hist_GM.1D"
        1dplot -ps -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" "$output_directory/hist_GM.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_HRF_GM_wholebrain.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $LMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_hrf+tlrc" > "$output_directory/hist_LMC.1D"
        3dhistog -omit 0 -mask $RMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_hrf+tlrc" > "$output_directory/hist_RMC.1D"
        1dplot -ps -dashed 1:2:2 -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist_GM.1D[1]" "$output_directory/hist_LMC.1D[1]" "$output_directory/hist_RMC.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_HRF_GM_MC.bmp" -sDEVICE=bmp256


       
        # CMRO2 calculation RRF using MEC
        3dcalc -prefix "$output_directory/CMRO2_tapping_rrf" \
            -a "$input_directory/stats.PW3.${subj}+tlrc.HEAD[6]" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/stats.MEC3.${subj}_REML+tlrc.HEAD[6]"\
            -d "$output_directory/rm.mean.signal.MEC3.all+tlrc"\
            -e "$output_directory/M_rrf+tlrc" \
            -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -omit 0 -mask $brainmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_RRF.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $GMmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_rrf+tlrc" > "hist_GM.1D"
        1dplot -ps -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" "$output_directory/hist_GM.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_RRF_GM_wholebrain.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $LMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_rrf+tlrc" > "$output_directory/hist_LMC.1D"
        3dhistog -omit 0 -mask $RMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_rrf+tlrc" > "$output_directory/hist_RMC.1D"
        1dplot -ps -dashed 1:2:2 -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist_GM.1D[1]" "$output_directory/hist_LMC.1D[1]" "$output_directory/hist_RMC.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_RRF_GM_MC.bmp" -sDEVICE=bmp256

        

        # CMRO2 calculation HRF using E2
        3dcalc -prefix "$output_directory/CMRO2_tapping_E2_hrf" \
            -a "$input_directory/stats.PW3.${subj}+tlrc.HEAD[6]" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/stats.E23.${subj}_REML+tlrc.HEAD[6]"\
            -d "$output_directory/rm.mean.signal.E23.all+tlrc"\
            -e "$output_directory/M_E2_hrf+tlrc" \
            -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -omit 0 -mask $brainmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_hrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_E2_HRF.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $GMmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_hrf+tlrc" > "$output_directory/hist_GM.1D"
        1dplot -ps -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" "$output_directory/hist_GM.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_E2_HRF_GM_wholebrain.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $LMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_hrf+tlrc" > "$output_directory/hist_LMC.1D"
        3dhistog -omit 0 -mask $RMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_hrf+tlrc" > "$output_directory/hist_RMC.1D"
        1dplot -ps -dashed 1:2:2 -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist_GM.1D[1]" "$output_directory/hist_LMC.1D[1]" "$output_directory/hist_RMC.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_E2_HRF_GM_MC.bmp" -sDEVICE=bmp256
        
        # CMRO2 calculation RRF using E2
        3dcalc -prefix "$output_directory/CMRO2_tapping_E2_rrf" \
            -a "$input_directory/stats.PW3.${subj}+tlrc.HEAD[6]" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/stats.E23.${subj}_REML+tlrc.HEAD[6]"\
            -d "$output_directory/rm.mean.signal.E23.all+tlrc"\
            -e "$output_directory/M_E2_rrf+tlrc" \
            -expr '(1 - (c/d)/e)*(((a+b)/b)^(8/10))' 
        3dhistog -omit 0 -mask $brainmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_rrf+tlrc" > "$output_directory/hist.1D"
        1dplot -ps -hist -x "$output_directory/hist.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_E2_RRF.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $GMmask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_rrf+tlrc" > "$output_directory/hist_GM.1D"
        1dplot -ps -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist.1D[1]" "$output_directory/hist_GM.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_E2_RRF_GM_wholebrain.bmp" -sDEVICE=bmp256
        3dhistog -omit 0 -mask $LMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_rrf+tlrc" > "$output_directory/hist_LMC.1D"
        3dhistog -omit 0 -mask $RMC_mask -nbin 100 -notitle -min -2 -max 3 "$output_directory/CMRO2_tapping_E2_rrf+tlrc" > "$output_directory/hist_RMC.1D"
        1dplot -ps -dashed 1:2:2 -one -norm2 -hist -x "$output_directory/hist_GM.1D[0]" -xlabel 'CMRO2 tapping' -ylabel 'histo' "$output_directory/hist_GM.1D[1]" "$output_directory/hist_LMC.1D[1]" "$output_directory/hist_RMC.1D[1]" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_tap_E2_HRF_GM_MC.bmp" -sDEVICE=bmp256
        
       # CMRO2 times series calculation
        # CMRO2(t) HRF
        3dcalc -prefix "$output_directory/CMRO2_time_serie_tapping_hrf" \
            -a "$input_directory/dataset.PW.r03+tlrc" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/dataset.MEC.r03+tlrc[$common_begin_offset..$common_TR03end]"\
            -d "$output_directory/rm.mean.signal.MEC3.all+tlrc"\
            -e "$output_directory/M_hrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_time_serie_tapping_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_hrf_R_MC.tapping.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_time_serie_tapping_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_hrf_L_MC.tapping.1D"
        1dplot -ps -one -xlabel 'CMRO2 tapping' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_tapping.MEC.1D" "$output_directory/CMRO2_t_hrf_L_MC.tapping.1D" "$output_directory/CMRO2_t_hrf_R_MC.tapping.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_tapping_hrf_MC.bmp" -sDEVICE=bmp256

        # CMRO2(t) E2 HRF
        3dcalc -prefix "$output_directory/CMRO2_E2_time_serie_tapping_hrf" \
            -a "$input_directory/dataset.PW.r03+tlrc" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/dataset.E2.r03+tlrc[$common_begin_offset..$common_TR03end]"\
            -d "$output_directory/rm.mean.signal.E23.all+tlrc"\
            -e "$output_directory/M_E2_hrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_E2_time_serie_tapping_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_hrf_R_MC.tapping.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_E2_time_serie_tapping_hrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_hrf_L_MC.tapping.1D"
        1dplot -ps -one -xlabel 'CMRO2 tapping' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_tapping.MEC.1D" "$output_directory/CMRO2_t_E2_hrf_L_MC.tapping.1D" "$output_directory/CMRO2_t_E2_hrf_R_MC.tapping.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_E2_tapping_hrf_MC.bmp" -sDEVICE=bmp256

        # CMRO2(t) RRF
        3dcalc -prefix "$output_directory/CMRO2_time_serie_tapping_rrf" \
            -a "$input_directory/dataset.PW.r03+tlrc" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/dataset.MEC.r03+tlrc[$common_begin_offset..$common_TR03end]"\
            -d "$output_directory/rm.mean.signal.MEC3.all+tlrc"\
            -e "$output_directory/M_rrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_time_serie_tapping_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_rrf_R_MC.tapping.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_time_serie_tapping_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_rrf_L_MC.tapping.1D"
        1dplot -ps -one -xlabel 'CMRO2 tapping' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_tapping.MEC.1D" "$output_directory/CMRO2_t_rrf_L_MC.tapping.1D" "$output_directory/CMRO2_t_rrf_R_MC.tapping.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_tapping_rrf_MC.bmp" -sDEVICE=bmp256

        # CMRO2(t) E2 RRF
        3dcalc -prefix "$output_directory/CMRO2_E2_time_serie_tapping_rrf" \
            -a "$input_directory/dataset.PW.r03+tlrc" \
            -b "$output_directory/rm.mean.signal.PW3.all+tlrc" \
            -c "$input_directory/dataset.E2.r03+tlrc[$common_begin_offset..$common_TR03end]"\
            -d "$output_directory/rm.mean.signal.E23.all+tlrc"\
            -e "$output_directory/M_E2_rrf+tlrc" \
            -expr '(1 - ((c-d)/d)/e)*(((a)/b)^(8/10))'
        echo `3dmaskave -q -mask $RMC_mask "$output_directory/CMRO2_E2_time_serie_tapping_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_rrf_R_MC.tapping.1D"
        echo `3dmaskave -q -mask $LMC_mask "$output_directory/CMRO2_E2_time_serie_tapping_rrf+tlrc.HEAD"`   >  "$output_directory/CMRO2_t_E2_rrf_L_MC.tapping.1D"
        1dplot -ps -one -xlabel 'CMRO2 tapping' -ylabel 'CMRO2 in motor cortex' "$input_directory/ideal_tapping.MEC.1D" "$output_directory/CMRO2_t_E2_rrf_L_MC.tapping.1D" "$output_directory/CMRO2_t_E2_rrf_R_MC.tapping.1D" | gs -r100 -sOutputFile="$imageDirectory/CMRO2_t_E2_tapping_rrf_MC.bmp" -sDEVICE=bmp256

        # mean CMRO2 in Left MC
        3dTcat -output $output_directory/M_LMC.tapping \
            "$output_directory/CMRO2_tapping_hrf+tlrc" \
            "$output_directory/CMRO2_tapping_rrf+tlrc" \
            "$output_directory/CMRO2_tapping_E2_hrf+tlrc" \
            "$output_directory/CMRO2_tapping_E2_rrf+tlrc"
        echo $subj `3dmaskave -q -mask $LMC_mask "$output_directory/M_LMC.tapping+tlrc.HEAD"`                           >> $SNRsavingdirectory/CMRO2_tapping_L_MC_tlrc.txt

        # mean CMRO2 in right MC
        3dTcat -output $output_directory/M_RMC.tapping \
            "$output_directory/CMRO2_tapping_hrf+tlrc" \
            "$output_directory/CMRO2_tapping_rrf+tlrc" \
            "$output_directory/CMRO2_tapping_E2_hrf+tlrc" \
            "$output_directory/CMRO2_tapping_E2_rrf+tlrc"
        echo $subj `3dmaskave -q -mask $RMC_mask "$output_directory/M_RMC.tapping+tlrc.HEAD"`                           >> $SNRsavingdirectory/CMRO2_tapping_R_MC_tlrc.txt


      # saving usefull data in output directory
        cp $input_directory/*stats*PW3*             $output_directory/
        cp $input_directory/*stats*MEC3*_REML+tlrc* $output_directory/
        cp $input_directory/*stats*E23*_REML+tlrc*  $output_directory/
        
    else
        # tSNR
        echo $subj NaN NaN NaN NaN NaN >> $SNRsavingdirectory/tSNR_tapping_tlrc.txt
        # CNR
        echo $subj NaN NaN NaN NaN NaN >> $SNRsavingdirectory/CNR_tapping_tlrc.txt
        # calculating Percent active voxels
        echo $subj NaN NaN NaN NaN NaN >> $SNRsavingdirectory/percentVoxAct_tapping_tlrc.txt
        # PVE
        echo $subj NaN NaN NaN NaN NaN >> $SNRsavingdirectory/PVE_tapping_tlrc.txt 
        # mean M in GM
        echo $subj NaN NaN NaN NaN NaN >> $SNRsavingdirectory/percentVoxAct_tapping_tlrc.txt
        # mean CMRO2 in LRC
        echo $subj NaN NaN NaN NaN     >> $SNRsavingdirectory/CMRO2_tapping_L_MC_tlrc.txt 
        # mean CMRO2 in RRC
        echo $subj NaN NaN NaN NaN     >> $SNRsavingdirectory/CMRO2_tapping_R_MC_tlrc.txt
        
        echo $subj NaN NaN NaN NaN NaN >> $SNRsavingdirectory/GM_tscores_tapping_tlrc.txt

    endif

    cp $temp_dir/mask_GM_resam+tlrc*               $output_directory/
    cp $temp_dir/mask_epi_anat.${subj}+tlrc*       $output_directory/
    cp $input_directory/*stats.PW.0TR*             $output_directory/
    cp $input_directory/*stats*MEC.0TR*_REML+tlrc* $output_directory/
    cp $input_directory/*stats*E2.0TR*_REML+tlrc*  $output_directory/
    cp $input_directory/*stats*PW2*                $output_directory/
    cp $input_directory/*stats*MEC2*_REML+tlrc*    $output_directory/
    cp $input_directory/*stats*E22*_REML+tlrc*     $output_directory/

 # ======================== end of subject loop =============================
	# rm -R $input_directory
	# remove temporary files
	\rm $output_directory/rm.*

end