Hi Rick,

Thank you for your help. Using 3dcalc was a bit cumbersome as the number of reference waveforms can vary and is often a large number
I also needed to extract the associated beta value.
I found another way that I think works using 3dTstat. Also, would you use the full R^2 or Stim#1 R^2 ?
Here is what i wrote:

foreach n (`seq -f %02g 0 $n_max`)
3dDeconvolve -mask ${mask}+orig -input ${data}+orig -polort 1 -num_stimts 1 \
-stim_file 1 ${ref}_multi.1D\[$n\] -rout -bout \
-bucket tmp \
-jobs 16
3dbucket -aglueto R+orig tmp+orig\[0\]
3dbucket -aglueto B+orig tmp+orig\[4\]
end

3dTstat -prefix R_best -max R+orig
3dTstat -prefix I_best -argmax R+orig

3dcalc -a B+orig -b I_best+orig -expr "a*equals(l,b)" -prefix tmp_B_best
3dTstat -prefix B_best -max tmp_B_best+orig