In principle in QFT you want to calculate the interaction picture unitary time evolution operator. Then you can use the operator to evolve quantum states just like you would in normal quantum mechanics. So, for example, suppose you set up your particle accelerator to generate a certain state $|i\rangle$ at a time $t_0$, and then after some time $t_1 - t_0$ has passed you want to know the probability it has evolved into a state $|f\rangle$, you calculate the probability amplitude:
$$
\langle f|U(t_1,t_0)|i \rangle
$$
The expression for $U(t_1,t_0)$ is given by Dyson's formula:
$$
U(t_1,t_0) = \mathcal{T}\left\{\exp\left(-i\int_{t_0}^{t_1}H_{\mathrm{int}}(t) \mathrm{d}t\right)\right\}
$$
where the $\mathcal{T}$ stands for time ordering and the expansion of the exponential is given by a Dyson series, the $n$th term of which has the form:
$$
U_{n}(t_1,t_0) = (-i)^{n}\int_{t_0}^{t_1}dt\int_{t_0}^{t}dt'\int_{t_0}^{t'}dt'' \cdots\int_{t_0}^{t^{(n-1)}}dt^{(n)}\mathcal{T}\left\{H_{\mathrm{int}}(t)H_{\mathrm{int}}(t') \cdots H_{\mathrm{int}}(t^{(n)})\right\}
$$
Typically the Hamiltonian has a coupling constant at the front, and so the Dyson series becomes an expansion in powers of that coupling constant (the power of which is equal to the number of vertices in the corresponding Feynman diagram). I don't really know what he means by powers of $\hbar$.
Typically the above integral is impossible to compute independently like you would in normal QM, but we can use a bit of magic called the LSZ reduction formula to write a probability amplitude like:
$$
\langle f|U(-\infty,\infty)|i \rangle \sim \langle 0|\mathcal{T}\left\{\phi(x)\cdots\phi(x^{(b)})\phi(y)\cdots\phi(y^{(b')})\ U(-\infty,\infty)\right\}|0 \rangle
$$
where $|0\rangle$ is the vacuum state of the free theory, $\phi$ are the fields and you have one for each particle in the initial state (labelled with $x$) and one for each in the final state (labelled with $y$). I've left out some integrals and constants that you'll get to eventually but aren't important for this high-level overview. Notice that this only holds in the $\pm\infty$ time limit, which is called the 'adiabatic limit'.
Once you have this, you can substitute in your $n$th order term for $U$ and get the $n$th order approximation for your probability amplitude. Remember that the $H_{\mathrm{int}}$ are typically products of fields themselves, so substituting just gives a larger product of time ordered fields (I've labelled the fields from the interaction terms with $z$):
$$
\langle 0|\mathcal{T}\left\{\phi(x)\cdots\phi(x^{(b)})\phi(y)\cdots\phi(y^{(b')})\phi(z)\cdots\phi(z^{(b'')})\right\}|0 \rangle
$$
Then you use a bit of technical mathematics called Wick's Theorem to show that this thing reduces to a product of two-point Green's functions (a.k.a. propagators):
$$
\langle 0|\mathcal{T}\left\{\phi\phi\right\}|0 \rangle\langle 0|\mathcal{T}\left\{\phi\phi\right\}|0 \rangle\langle 0|\mathcal{T}\left\{\phi\phi\right\}|0 \rangle\cdots
$$
Now, obviously there's more than one way to pair off the $b+b'+b''$ fields you have into Green's functions, which is why you have multiple Feynman diagrams for each order in the coupling constant. Each one of these Green's functions in the product corresponds to a line in a Feynman diagram, and each way to pair the fields off corresponds to a different diagram (up to symmetry factors etc). Loops are when you get Green's functions where both fields are at the same spacetime point:
$$
\langle 0|\mathcal{T}\left\{\phi(z)\phi(z)\right\}|0 \rangle
$$
because obviously you get a line from a point back to itself -- a loop. This results in unconstrained momentum because a typical Green's function has the form:
$$
\langle 0|\mathcal{T}\left\{\phi(x)\phi(y)\right\}|0 \rangle = \int \mathrm{d}^{3}p\ \frac{i}{p^2 - m^2}e^{-ip(x-y)}
$$
The exponential then reduces to a delta function when you integrate over position, allowing you to easily integrate over momentum to essentially select the momentum you want (this comes from one of the integrals I left out earlier). But if $x=y$ then the exponential reduces to 1, meaning your momentum is unconstrained and generally this also leaves the integral divergent. This is the root of ultraviolet divergences in QFT, and we use regularisation to renormalise certain constants (like the mass) and absorb the infinities into unmeasureable quantities.
Edit: I notice you say you're reading P&S as well as Schwartz. I don't know Schwartz's book at all, but in P&S they derive the LSZ Reduction formula for scalar fields in Chapter 7. It's quite instructive to look through it because it gives you a much better appreciation for the form of the Feynman rules.
