Use of the Dirac delta as a sampling operator As we are all taught in signal processing classes, the Dirac delta is not a normal function but a generalized function. It is sometimes defined as either the limit of some normal function whose support goes to zero as some parameter goes to zero, while keeping unit area. This apparently leads to the usual representation of the Dirac delta function as being defined and only having meaning under the so-called sifting integral,
$$
\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-t\right)d\tau=x\left(t\right)
$$

Yet, pick up any book or paper on signal sampling and you will see the sampling process defined as a modulation process with the delta standing alone as a normal function, multiplying the analog signal to be sampled, like this:
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)
$$
Why is this legal? Why isn't the sampling process modeled as
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau
$$
or
$$
x\left(nT\right)=\intop_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau\text{?}
$$

Similarly, in the frequency domain, we know that a linear system's response to a sinusoidal input is a common measurement. Using the common definitions for $X$, $Y$, and $H$, we know that
$$
Y\left(f\right)=X\left(f\right)H\left(f\right).
$$
If 
$$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
is a sinusoid, then
$$
X\left(f\right)=\delta\left(f-f_{0}\right)
$$
and
$$
Y\left(f\right)=H\left(f_{0}\right)
$$
where we have once again used the stand-alone delta as a sampling operator.

Can someone please clearly explain my confusion on this matter? Why do we seem to be so cavalier about how we use the Dirac delta when modeling sampling?

(While we're in the neighborhood, maybe someone could explain why we're even allowed to think about the Fourier transform of a sinusoid since it it is neither absolutely integrable nor square integrable.)

By the way, I find the discussion at Why is dirac delta used in continuous signal sampling? to be unsatisfying. Can someone do better?
The OP's "equation"
$$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\tag{1}
$$
is incorrect; the sum on the right side of $(1)$ is the sampled signal $x_s(t)$, not $x(nT)$ at all.  That is,
$$x_s(t) = \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right).\tag{2}$$
This is not an ordinary function since it is a weighted sum of Dirac deltas, but if we pretend that it is indeed an ordinary function, then the  "value" of $x_s(t)$ at $t=mT$ is
$$x_s(mT) = \sum_{n=-\infty}^{\infty}x\left(nT\right)\delta\left(mT-nT\right) = x(mT)\delta(0),$$
while if $t$ is not an integer multiple of $T$, then the "value" of $x_s(t)$ is $0$. Note the quotation marks around the word value and avoid the usual nonsense of saying that $\delta(0) = \infty$ that is found in so many low-level DSP books.  $x_s(t)$ is not a function and it does not have any value at any real number $t$. What the sampled signal
$x_s(t)$ as defined in $(2)$ is is an impulse train; a collection of impulses spaced $T$ seconds apart with coefficient or amplitude of the impulse at time $t=mT$ being $x(mT)$, that is, the sampled value of $x(t)$ at time $mT$. As stated earlier, $x_s(t)$ is not an ordinary function of time; it is, as all the other answers have pointed out, a functional, that is, a function that operates on functions as a whole, not on individual specific values that the function being operated on takes on as ordinary functions do.   Now, it is a standard notion that
$$x(t)\delta(t-t_0) = x(t_0)\delta(t-t_0)$$ provided that $x(t)$ is continuous at $t=t_0$, and so if we re-write $(2)$  as
$$x_s(t) = \sum_{n=-\infty}^\infty x\left(nT\right)\delta\left(t-nT\right)
\tag{3}$$ we see that $(3)$ formally resembles the commonly used sum
$$x_m(t) = \sum_{n=-\infty}^{\infty}x\left(nT\right)g\left(t-nT\right)\tag{4}$$
which describes the modulation of a stream of pulses (of the form $g(t)$ spaced $T$ seconds apart) by the sample values of $x(t)$ spaced $T$ seconds apart. Small wonder then that textbook writers notice the analogy between $x_s(t)$ as given in $(3)$ and $x_m(t)$ in $(4)$ and impulsively describe the process that leads to $(2)$ or $(3)$ as a modulation process instead of a multiplication process. (At this point, it is worth reflecting that several modulation processes (e.g. AM-DSBSC which can be described as $x(t) \mapsto 
x(t)\cos(2\pi f_ct)$) are indeed just multiplication of a carrier signal by a baseband signal).

So, why is this rigamarole useful to anyone?  Well, let's look at the Fourier transform of $x_s(t)$:

\begin{align}
 \int_{-\infty}^\infty x_s(t) \exp(-j2\pi ft) \,\mathrm dt
&amp;= \int_{-\infty}^\infty \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt\\
&amp;= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) \exp(-j2\pi ft) \,\mathrm dt
&amp;{\scriptstyle{\text{sifting property of Dirac delta}}}\\
&amp;= \sum_{n=-\infty}^{\infty} x(nT)\exp(j2\pi f(nT)) \tag{5}
\end{align}
which is the DTFT of the discrete-time sequence $\big\{x(nT)\big\}_{n=-\infty}^\infty$ of sample values of $x(t)$ spaced $T$ seconds apart. More generally, if we have a (linear time-invariant BIBO) filter whose impulse response (dirty word!) is $g(t)$, then applying the input $x_s(t)$ to this filter, we get that the output is
\begin{align}
\int_{-\infty}^\infty \left( \sum_{n=-\infty}^{\infty}x\left(t\right)\delta\left(t-nT\right)\right) g(\tau-t)
\,\mathrm dt
&amp;= \sum_{n=-\infty}^{\infty}\int_{-\infty}^\infty x\left(t\right)\delta\left(t-nT\right) g(\tau-t)
\,\mathrm dt\\
&amp;= \sum_{n=-\infty}^{\infty} x(nT) g(\tau-nT)\\
&amp;= x_m(\tau)&amp; \text{as defined in }(4)
\end{align}
So, Equation $(2)$ does have some utility whereas the jury is still out on the OP's $(1)$.

The OP asks:


  Why isn't the sampling process modeled as
  $$
x\left(nT\right)=\sum_{n=-\infty}^{\infty}\intop_{-\infty}^{\infty}x\left(\tau\right)\delta\left(\tau-nT\right)d\tau ??
$$


Well, quite apart from the royal confusion that the OP is creating by using $n$, a parameter on the left side of the OP's alleged definition, as an index of summation, note that the sifting integral property says that the value of that integral is $x(nT)$.  Consequently, the OP's purported definition reduces to
$$x(nT) = \sum_{n=-\infty}^\infty x(nT),$$
which is nonsensical, and matters don't improve very much if a different symbol is used for the parameter on the left side because all one gets is
$$x(mT) = \sum_{n=-\infty}^\infty x(nT),$$
which suggests that $x(nT) = 0$ for all integers $n$ that are not equal to the chosen $m$.  That's no way to run a railroad!



Turning to the frequency domain, the OP asserts that


  If 
  $$x\left(t\right)=\exp\left(j2\pi f_{0}t\right)$$
  is a sinusoid, then
  $$
X\left(f\right)=\delta\left(f-f_{0}\right) \tag{4}
$$
  and
  $$
Y\left(f\right)=H\left(f_{0}\right).\tag{5}
$$


While $(4)$ is correct, $(5)$ is not. Actually,
$$Y(f) = H(f)X(f) = H(f)\delta(f-f_0) = H(f_0)\delta(f-f_0)\tag{6}$$ which does not equal
$H(f_0)$ as the OP claims.  Indeed, in $(5)$ the left side is a function of $f$ while the right side is a constant so that the standard interpretation of $(5)$ is that $Y(f)$ is a constant, which in turn implies that
$y(t)$ is an impulse or Dirac delta!  Since $H(f)$ is arbitrarily chosen, what $(5)$ seems to be implying is that the response of any filter to a complex sinusoid is an impulse, which is manifestly nonsense.  What is true is that if we assume the system with transfer function $H(f)$ (a.k.a frequency response) has input $\exp(j2\pi f_0 t)$, then the output is 
$$y(t) = \mathcal F^{-1}(Y(f)) = \int_{\infty}^\infty H(f)\delta(f-f_0)\exp(j2\pi tf)\, \mathrm df = H(f_0)\exp(j2\pi f_0t)$$
which is the epitome of the definition of the frequency response function: 


  If for each real number $f_0$, the complex sinusoidal input $\exp(j2\pi f_0t)$ at frequency $f_0$ produces complex sinusoidal output $H(f_0)\exp(j2\pi f_0t) = |H(f_0)|\exp\big(j(2\pi f_0t+\angle H(f_0))\big)$ (which is also at frequency $f_0$ but has amplitude
  $|H(f)|$ and phase $\angle H(f_0)$), then $H(f)$ is called the frequency response function or transfer function of the linear time-invariant system.


In short, the OP's whole question is based on massaging "alternative facts" (that do not withstand scrutiny) to arrive at incorrect conclusions.