Why using the impulse respond of an LTI to characterize its output?
For the two special classes of systems:
$$ \textbf{Linear time and shift invariant systems } $$ An LTI and LSI system is one that obeys superposition and time and shift invariance. The LTI systems are described by linear ordinary differential equations and are your analog systems. The LSI systems are described by your linear constant coefficient difference equations. There are many analogues between them, as is to be expected as the discrete systems can be taken to be the discrete versions of the analog ones. 

The dirac pulse and unit impulse sequence completely characterise the systems in both the time and frequency domain. 

Let us first define what a system is, it is a relation that takes in an input signal $ x(t) $ or $ x[n] $ and produces an output $y(t)$ or $y[n]$. As is the case for analogue and discrete ones. The system is thus a relation, or an operator on the input to produce the output. 

$$
y(t) = \textbf{T} \big \{x(t) \big \} 
$$

$$
y[n] = \textbf{T} \big \{x[n] \big \} 
$$
We define the response to a dirac pulse and a unit impulse respectively for system classes 1 and 2 as: 
$$
h(t) = \textbf{T} \big \{ \delta(t) \big \} 
$$

And the unit impulse response for discrete linear shift invariant systems as: 
$$
h[n] = \textbf{T} \big \{ \delta[n] \big \} 
$$

The impulse response for these two classes of systems is hence, the output of such systems when the input is a dirac pulse or a unit impulse sequence. Note the following properties of the dirac pulse and unit impulse sequence: 


Every function can be expressed as a superposition or sum of scaled and shifted dirac pulses or unit impulses. 
$$
y(t) = \displaystyle \int_{- \infty}^{\infty} y(\tau) \cdot \delta(\tau - t) \,\,\,\, \text{d}t
$$


$$
y[n]= \displaystyle \sum_{ k \to - \infty}^{k \to \infty} x[k] \cdot \delta[k - n] 
$$


The continous time fourier transform and the discrete fourier transform of the dirac pulse and unit impulse are both 1, hence, in the frequency domain, the pulse subjects your system to all possible frequencies, giving a fingerprint of the system. 


$$
\textbf{The transfer function of a system (frequency domain or s domain) }
$$

It is the fourier transform of the impulse response, or the laplace transform of the impulse response (s domain). 

$$
H(j \omega) = \displaystyle \int_{- \infty}^{\infty} h(t) \cdot e^{-j(\omega t) } \,\,\,\, \text{d}t 
$$

And

$$
 H(s) = \displaystyle \int_{- \infty}^{\infty} h(t) \cdot e^{-s (t) } \,\,\,\, \text{d}t 
$$

Now we have to talk about convolution and multiplication time duality. 
Multiplying two functions in the frequency domain is equal to convolving them in the frequency domain. 
Convolving in the time domain is the same as multiplying in the frequency domain. 

Now to put it all togther, since your input function can be expressed as a sum of scaled and shifted dirac pulse or unit impulses, your output, due to linearity and time invariance, would be equally scaled and shifted impulse responses: 
$$ a h(t - t_{0} ) = \textbf{T} \big \{ a \delta(t - t_{0} ) \big \} $$
This is where convolution comes in, you express the input function as a sum of scaled and shifted dirac pulses or unit impulses, you then input these into your system, and to get the output you sum up or integrate all the scaled and shifted impulse responses, to get your output. 

$$
y(t) = \displaystyle \int_{-\infty}^{\infty} f( \tau) \cdot h(t - \tau) \,\,\,\, \text{d} \tau 
$$

$$
y[n] = \displaystyle \sum_{k \to - \infty}^{ k \to \infty} x[k] \cdot h[n - k] 
$$
And in the frequency domain due to the theorems this is equal to
$$
Y(s) = H(s) \cdot X(s)
$$