#!/bin/tcsh

if ( $#argv > 0 ) then
    set subjects = ( $argv )
else
    set subjects = ED
endif

#===========================================================================
# Above command will run script for all our subjects - ED, EE, EF - one after
# the other if, when we execute the script, we type: ./@analyze_ht05 ED EE EF.
# If we type ./@analyze_ht05 or tcsh @analyze_ht05, it'll run the script only
# for subject ED.  The user will then have to go back and edit the script so
# that 'set subjects' = EE and then EF, and then run the script for each subj.
#===========================================================================
    
foreach subj ($subjects)
    
    cd $subj
    
    #======================================================================
    # volume register and time shift our datasets, and remove the first
    # two time points
    #======================================================================    
 
    foreach run ( `count -digits 1 1 10` )
	3dvolreg -verbose			\
	    -base {$subj}_r1+orig'[2]'	        \
	    -tshift 0				\
	    -prefix  {$subj}_r{$run}_vr		\
	    {$subj}_r{$run}+orig'[2..137]'

	# store run data in runs_orig directory

    #======================================================================
    # smooth data with 3dmerge.
    #======================================================================    

	3dmerge -1blur_rms 4			\
		-doall 				\
		-prefix {$subj}_r{$run}_vr_bl	\
		{$subj}_r{$run}_vr+orig
    
    end


    #======================================================================
    # create masks for each run using 3dAutomask
    #====================================================================== 

    foreach run ( `count -digits 1 1 10` )
        3dAutomask -prefix mask_r{$run} {$subj}_r{$run}_vr_bl+orig
    end

    #======================================================================
    # create a mask enveloping masks of the individual runs
    #====================================================================== 
    
    3dcalc -a mask_r1+orig -b mask_r2+orig -c mask_r3+orig	\
	   -d mask_r4+orig -e mask_r5+orig -f mask_r6+orig	\
	   -g mask_r7+orig -h mask_r8+orig -i mask_r9+orig	\
	   -j mask_r10+orig					\
	   -expr 'step(a+b+c+d+e+f+g+h+i+j)'			\
	   -prefix full_mask

    #======================================================================
    # re-scale each run's baseline to 100.
    # E.g., If baseline is 100, and result of 3dcalc on one voxel is 106, then
    # we can say that at that voxel shows a 6% increase is signal activity
    # relative to baseline.
    # Use full_mask to remove non-brain 
    #====================================================================== 

    foreach run ( `count -digits 1 1 10` )
        3dTstat -prefix mean_r{$run} {$subj}_r{$run}_vr_bl+orig
    
        3dcalc -a {$subj}_r{$run}_vr_bl+orig   \
               -b mean_r{$run}+orig            \
               -c full_mask+orig               \
	       -expr "(a/b * 100) * c"         \
               -prefix scaled_r{$run}
    
        \rm mean_r{$run}+orig*
    
    end

    #======================================================================
    # Now we can concatenate our 10 normalized runs with 3dTcat.
    #======================================================================

    3dTcat -prefix {$subj}_all_runs  \
        scaled_r1+orig               \
        scaled_r2+orig               \
        scaled_r3+orig               \
        scaled_r4+orig               \
        scaled_r5+orig               \
        scaled_r6+orig               \
        scaled_r7+orig               \
        scaled_r8+orig               \
        scaled_r9+orig               \
        scaled_r10+orig

    #======================================================================   
    # mv unneeded run data into separate directories
    #======================================================================    

    mkdir runs_orig runs_temp

    mv {$subj}_r*_vr* scaled* runs_temp
    mv {$subj}_r* runs_orig

    #======================================================================   
    # run deconvolution analysis
    #======================================================================    

    3dDeconvolve -input {$subj}_all_runs+orig -num_stimts 4                   \
        -stim_file 1  ../misc_files/all_stims.1D'[0]' -stim_label 1 ToolMovie \
    		 -stim_minlag 1 0 -stim_maxlag 1 14 -stim_nptr 1 2            \
        -stim_file 2  ../misc_files/all_stims.1D'[1]' -stim_label 2 HumanMovie\
    	         -stim_minlag 2 0 -stim_maxlag 2 14 -stim_nptr 2 2            \
        -stim_file 3  ../misc_files/all_stims.1D'[2]' -stim_label 3 ToolPoint \
    	         -stim_minlag 3 0 -stim_maxlag 3 14 -stim_nptr 3 2            \
        -stim_file 4  ../misc_files/all_stims.1D'[3]' -stim_label 4 HumanPoint\
    	         -stim_minlag 4 0 -stim_maxlag 4 14 -stim_nptr 4 2            \
        -glt 4 ../misc_files/contrast1.1D  -glt_label 1 FullF                 \
        -glt 1 ../misc_files/contrast2.1D  -glt_label 2 HvsT                  \
        -glt 1 ../misc_files/contrast3.1D  -glt_label 3 MvsP                  \
        -glt 1 ../misc_files/contrast4.1D  -glt_label 4 HMvsHP                \
        -glt 1 ../misc_files/contrast5.1D  -glt_label 5 TMvsTP                \
        -glt 1 ../misc_files/contrast6.1D  -glt_label 6 HPvsTP                \
        -glt 1 ../misc_files/contrast7.1D  -glt_label 7 HMvsTM                \
        -iresp 1 TMirf -iresp 2 HMirf -iresp 3 TPirf -iresp 4 HPirf           \
        -full_first -fout -tout -nobout -polort 2			      \
        -concat ../misc_files/runs.1D                                         \
        -progress 1000							      \
        -bucket {$subj}_func

    #======================================================================   
    # make slim dataset.  Too many sub-bricks in our bucket dataset.
    # Use 3dbucket to slim it down and include sub-bricks of interest only.
    #======================================================================    

    3dbucket -prefix {$subj}_func_slim -fbuc {$subj}_func+orig'[125..151]'
    
    #======================================================================    
    # Remember IRF datasets created by 3dDeconvolve? 
    # There are 15 time lags in each voxel.  Remove lags 0-3 and 10-15 b/c not 
    # interesting.  Then find mean percent signal change for lags 4-9 in each 
    # voxel with '3dTstat'.  
    # Then transform to Talairach coordinates with 'adwarp'.
    #======================================================================    
    
    foreach cond (TM HM TP HP)
         3dTstat -prefix {$subj}_{$cond}_irf_mean  {$cond}irf+orig'[4..9]'
    
         adwarp  -apar {$subj}spgr+tlrc  -dpar {$subj}_{$cond}_irf_mean+orig
    end
    
    cd ..
end    

    #======================================================================
    # End of script!  
    # Take the {$subj}_{$cond}_irf_mean+tlrc datasets and input into 3dANOVA2.
    #======================================================================