There isn't a general process of selecting an activation function. A lot of it is emperical, once they satisfy certain properties..
Although, we are interested in having the network compute interesting functions, so if you were to use for example a linear function as activation function (i.e. $f(x) = x$) then you're network wouldn't be able to model non-lineariites that might be present in your dataset. As a result, you would want non-linearities in your activation function, which are present in ReLU, sigmoid, tanh, etc.
The reason you see ReLU $f(x) = max(0, x)$ being used by default is because it enabled gradients to flow when the input to the ReLU function is positive, and does not have the saturation problems of sigmoid/tanh. Then some subequent papers they saw that: "Oh! units die on the left half" (where it's flat), and so they introduced modifications like PReLU which is $f(x) = max(-\alpha x, x)$ so now you also have gradients flowing when the input is negative (i.e. learning won't stop).
Now as to which one you should use for your network, you would have to run experiments!