First of all write down an explicit expression for the summation over all microstates.
Edit
Since you're treating the system classically this includes an integral over phase-space & a summation over all possible spin-configurations.
$$ \sum_{\mu_s} = \sum_{\{s_i\}}\int\frac{d^Np d^Nq}{h^{3N}N!}$$
The second thing is to realize that your Hamiltonian is non-interacting and the canonical density $e^{-\beta H} $ is just a product of one-particle Hamiltonians
$$ e^{-\beta H} =\prod_{i=1}^Ne^{-\beta h_i} $$
where of course
$$ h_i = {\frac{(p_i + \gamma s_i)^2}{2m} -b s_i} $$
(I renamed $\beta$ to $\gamma$, because you don't mean the inverse temp. here)
So you have to evaluate
$$ \frac{1}{N!}\sum_{\{s_i\}}\prod_{i=1}^N\int\frac{dp_i dq_i}{h^{3}} e^{-\beta h_i} = $$
Because $H$ ist non-interacting, the N-particle phase-space-integral factorizes into N integrations over a 1-particle phase-space. Similarily one may interchange the spin-summation and the product (convince yourself that this is true! e.g that one ends up with the same terms, whether you sum over spin first or not.)
That means, instead of summing over all the many-body microstates, one first sums over the possible configurations of a single particle and accounts for the fact, that there are many afterwards. Additionally all $h_i$ are equivalent. They each just carry a different but redundant index:
$$ Z = \frac{1}{N!}\prod_{i=1}^N\sum_{s_i=\pm1}\int\frac{dp_i dq_i}{h^{3}} e^{-\beta h_i} = \frac{1}{N!}\left(\sum_{s=\pm1}\int\frac{dp dq}{h^{3}} e^{-\beta h}\right)^N$$
where $h$ is $h_i$ but without an index, because the reference is not to a specific particle anymore.
/Edit
You'll have to think about, what to do with the momentum integration. I haven't calculated the result, but it might be you won't end up with a solution in closed form. There might be an approximation needed to do the momentum summation. Edit I think a gaussian integration will do the trick. /Edit
Let us know what you end up with!
You are misunderstanding what is meant by a product over monopartiuclar states. This is not a product over the states of the $N$ particles in the system, it is a product over all possible single particles states. The sum is then over the occupation number of those states, i.e. the number of particles actually in that state (which may be $0$). The advantage of this occupation number approach to over keeping track of individual particles is that it automatically takes account of particle indistinguishably.
The total number of particles in a state $\gamma$ is clearly simply the sum of the number of particles in each single particle state
$$
N_\gamma = \sum_i n_i
$$
The total energy is the energy of each single particle state, times the number of particles in that state
$$
E_\gamma = \sum_i \epsilon_i n_i
$$
Now the formula follows from a simple manipulation
\begin{align}
\mathcal{Z} &= \sum_{\gamma} e^{-\beta(E_\gamma - \mu N_\gamma)}\\
&= \sum_{\gamma} e^{-\beta\sum_i n_i(\epsilon_i-\mu)}\\
&= \sum_{n_0}\sum_{n_1}\ldots\; \left(\prod_i e^{\beta n_i(\mu - \epsilon_i)}\right)\\
&= \prod_i \sum_{n_i} \left(ze^{-\beta\epsilon_i}\right)^{n_i}
\end{align}
Best Answer
They are not actually proportional, I think it's meant something else. The partition function $Z(\beta)$ is the Laplace Transform (LT) of the density $g(E)$ of states of energy $E$, i.e. \begin{equation} Z(\beta)=\int_{0}^{\infty}e^{-\beta E}g(E)dE. \end{equation} The partition function $Z$ plays the role of a characteristic function and $g(E)$ the role of probability density function (p.d.f.), that is the p.d.f. is $p(E)=\frac{g(E)}{\int g(E)dE}$.
To be more precise, the partition function $Z$ is a "probability-generating function" like it is the characteristic function ($\varphi_{X}:=\langle e^{itX} \rangle$ and in this case is the fourier transform of the p.d.f.), the moment-generating function $(M_{X}(t):=\langle e^{tX} \rangle)$ and the cumulant-generating function.
The Probability-generating function is defined as \begin{equation*} G(s)=\sum_{n=0}^{\infty}p_{n}s^{n} \end{equation*} where $p_{n}$ is the probability mass and the series is absolutely convergent for $|s|\leq1$. Then \begin{equation*} Z(\beta)=G(e^{-\beta})=\sum_{E}^{\infty}p_{E}e^{-\beta E}. \end{equation*} Note that this can be generalized to \begin{equation} Z_{H}(\beta)=\int e^{-\beta H(\vec{X})}d\mu(\vec{X}), \end{equation} where the integral is over the space where $\vec{X}$ belongs to (usually phase-space) and $d\mu(\vec{X})$ is a measure on this space (usually $d\mu(\vec{X})=g(\vec{X})d\vec{X}$). It's common to use the assumption that the density is uniform (microcanonical ensemble).
This is a nice question since these probability-generating functions appear in many different areas of physics such statistical mechanics, electromagnetism, quantum field theory and quantum mechanics. All with different names but the same mathematical structure which simplifies calculations of moments and correlations.
References: 1. http://en.wikipedia.org/wiki/Partition_function_(mathematics)