Two and three point function of primary fields of arbitrary fields I was looking at this paper hep-th/0011040 and I found the following equation:
$$ \langle C_{\mu_1 \dots \mu_l} \mathcal{O}^{\mu_1 \dots \mu_l}(x_1) D_{\nu_1 \dots \nu_l} \mathcal{O}^{\nu_1 \dots \nu_l}(x_2)\rangle $$$$= \dfrac{1}{|x_1 - x_2|^{2\eta}} C_{\mu_1 \dots \mu_l}I^{\mu_1 \nu_1}(x_{12}) \ldots I^{\mu_l \nu_l}(x_{12}) D_{\nu_1 \dots \nu_l}, \tag{2.6}$$ 
where the $\mathcal{O}$ are spin $\ell$ symmetric traceless primary operators of dimension $\eta$ and 
$$ I^{\mu \nu}(x) = \delta^{\mu \nu} - 2\dfrac{x^\mu x^\nu}{x^2}. \tag{2.7}$$
Could anyone tell me how to arrive at this result, or any reference where this has been given? I'm also interested in understanding the analogous derivation for $3$-point functions. 
 A: Note: see comments for generalization from scalar fields to higher spins.
I shall do the $2$-point case and leave the $3$-point one to you as an exercise. The canonical reference which I'm using is Di Francesco, Mathieu and Senechal. This is pretty much obligatory reading if you're interested in conformal field theory.
Let $\phi_1$ and $\phi_2$ be two primary fields. Then by definition their conformal transformations are
$$\phi_1(x_1)=\left|\frac{\partial x'}{\partial x}\right|_{x=x_1}^{\Delta_1/d}\phi_1(x_1')$$
where $d$ is the spacetime dimension and $\Delta_1$ the scaling dimension of the field. There's a similar formula for $\phi_2$. Now since the measure and action in the functional integral are conformally invariant we can promote the transformation above to one of the correlation function, viz.
$$\langle\phi_1(x_1)\phi_2(x_2)\rangle=\left|\frac{\partial x'}{\partial x}\right|_{x=x_1}^{\Delta_1/d}\left|\frac{\partial x'}{\partial x}\right|_{x=x_2}^{\Delta_2/d}\langle\phi_1(x_1')\phi_2(x_2')\rangle$$
Now we start specialising to specific symmetries. Already from rotational and translational symmetry we know
$$\langle\phi_1(x_1)\phi_2(x_2)\rangle=f(|x_1-x_2|)$$
Now using a scale transformation $x\to \lambda x$ in our formula above produces 
$$f(|x_1-x_2|)=\lambda^{\Delta_1+\Delta_2}f(\lambda|x_1-x_2|)$$
but this fixes the correlation function up to a constant, that is
$$\langle\phi_1(x_1)\phi_2(x_2)\rangle=\frac{C}{|x_1-x_2|^{\Delta_1+\Delta_2}}$$
The final ingredient is invariance under special conformal transformations. These have 
$$\left|\frac{\partial x'}{\partial x}\right|=\frac{1}{(1-2b\cdot x+b^2x^2)^d}$$
Substituting this into our formulae above fixes $\Delta_1=\Delta_2$. Therefore we have exactly proved your first result above
$$\langle\phi_1(x_1)\phi_2(x_2)\rangle=\frac{C}{|x_1-x_2|^{2\Delta_1}}$$
if your operators have the same scaling dimension $\Delta_1$. 
Now try the $3$-point case using the same symmetries as above! Also, have a think about why this method fails at $4$-point. Hint: you can make conformally invariant cross-ratios when you have $4$ positions.
Endnote on Scaling Dimensions
Remember that the scaling dimension is not necessarily the naive engineering dimension of your field, even in a CFT. This is because loop corrections can introduce divergences even though there's global conformal symmetry. Various theorems guarantee that these don't renormalize masses or couplings but they do cause field strength renormalization $\phi \to \sqrt{Z}\phi$. Typically $Z$ carries some dimensionality, hence we get anomalous dimensions. This argument only falls down if the conformal symmetry is local (as in string theory) or if the operator is protected by supersymmetry.
