Multiqubit state tomography by performing measurement in the same basis For a $n$-qubit state $\rho$ we perform all projective measurement consisting of one-particle measurements in the same basis, that is,
$$p_{i_1i_2\ldots i_n}(\theta,\varphi) = \text{Tr}\left \{ \rho [P_{i_1}(\theta,\varphi)\otimes P_{i_2}(\theta,\varphi) \otimes \ldots\otimes P_{i_n}(\theta,\varphi)]    \right\},$$
where $P_0 = |\psi(\theta,\varphi)\rangle\langle \psi(\theta,\varphi) |$ and $P_1=\mathbb{I}-P_0$.
Does it suffice to recover the state? 
If so, how can one do it?
 A: Working out an explicit case of Alex's construction, you get the two states
$$ \rho_1 = \frac{1}{2} \left(|0\rangle\langle 0| \otimes |+ \rangle \langle + | + |1\rangle\langle 1| \otimes |- \rangle \langle - |\right) = \frac{1}{4}\left(\begin{array}{cccc}1 & 1 && \\ 1 & 1 & & \\ & & \phantom{-}1 & -1 \\ & & -1 & \phantom{-}1\end{array}\right)$$
and its image under the swap operator
$$ \rho_2 = \frac{1}{2}\left(|+\rangle\langle +| \otimes |0 \rangle \langle 0 | + |-\rangle\langle -| \otimes |1 \rangle \langle 1 |\right) = \frac{1}{4}\left(\begin{array}{cccc}1 &  &1& \\ & \phantom{-}1 & & -1\\ 1 & & 1 &  \\ & -1&  & \phantom{-}1\end{array}\right). $$
It's fairly easy to check that you cannot tell the difference between these. Let $|\psi\rangle = \alpha |0\rangle + \beta |1 \rangle$. Then $\langle \psi | = \alpha^* \langle 0 | + \beta^* \langle 1 |$ and $\langle \psi^\perp | = \beta \langle 0 | - \alpha \langle 1 |$. 
Taking inner products, you find the probability of all the outcomes are equal for both $\rho_1$ and $\rho_2$. The outcomes are given by the following table:
$$
\begin{array}{c  c}
\mathrm{outcome} &\mathrm{probability}\\
\langle \psi | \langle \psi | &\textstyle{\frac{1}{4}}\left(|\alpha|^2 |\alpha+\beta|^2  + |\beta|^2 |\alpha-\beta|^2\right)\\
\langle \psi | \langle \psi^\perp | &\textstyle{\frac{1}{4}}\left(|\alpha|^2 |\alpha-\beta|^2  + |\beta|^2 |\alpha+\beta|^2\right)\\
\langle \psi^\perp | \langle \psi | &\textstyle{\frac{1}{4}}\left(|\alpha|^2 |\alpha-\beta|^2  + |\beta|^2 |\alpha+\beta|^2\right)\\
\langle \psi^\perp | \langle \psi^\perp | &\textstyle{\frac{1}{4}}\left(|\alpha|^2 |\alpha+\beta|^2  + |\beta|^2 |\alpha-\beta|^2\right)
\end{array}
$$
It does seem possible that you can distinguish between any two pure states using the OP's measurements, but I haven't checked this carefully. 
A: It does not suffice for $n>1$.
The idea is that density matrix for qubit may be expressed as 
$\frac{1}{2}\bigl(\hat{1}+ \sum_{k=1}^3 v_k \hat\sigma_k\bigr)$, there $v$ is a 3D real unit vector.
Now we should consider span of tensor products of such terms. 
For any fixed $v$ we have $2^n$ tensor products of 
$\hat{P}_i = \bigl(\hat{1} + (-1)^i\sum_{k=1}^3 v_k \hat\sigma_k\bigr)/2$, $i =0,1$. 
Linear span of given products is equal with span generated by $2^n$ different combinations of tensor 
products of $2 \times 2$ unit matrix and (fixed for given $v$) matrix 
$\hat v = \sum_{k=1}^3 v_k \hat\sigma_k$.
We want to know span of all such products for arbitrary $v$. This span is less than maximal $4^n$ 
due to linear dependence between coefficients (see example below).
Already consideration of the case $n=2$ shows, why it does not work. Four projectors relevant for this case may be written as 
$$\frac{1}{4}(\hat{1} \pm \sum_{k=1}^3 v_k \hat\sigma_k) \otimes (\hat{1} \pm \sum_{k=1}^3 v_k \hat\sigma_k)$$
so coefficients for $\hat\sigma_k \otimes \hat\sigma_j$ and $\hat\sigma_j \otimes \hat\sigma_k$ are equal for $k,j = 1,2,3$ in all four cases ($\pm v_k v_j$). They are equal for any $v_k$ if bases for both qubits are equal, so consideration of many different measurements does not change the property. So decomposition with Pauli matrices based on such terms always has three equal coefficients and may not represent decomposition of general state.
A: If it can be done, then you will need at least $4^n$ angles $(\theta_k, \phi_k)$.  However, the code below suggests that this method will generate at most $O\left(2^{\frac32n}\right)$ independent probabilities.  Thus, it only works for $n=1$.

Suppose for the moment that it can be done; here is how to reconstruct the state from a set of angles.  Arrange the angles into one index set $\Lambda = \{(\theta_k,\phi_k)| k\in\{0,1,2,\ldots,4^n\}\}$. So now you have an overcomplete (and necessarily non-orthogonal) basis of operators $\{P(\lambda) | \lambda\in \Lambda\}$.  In other branches of mathematics, this thing is called a frame.  To recover the state from the set of probabilities $\{p(\lambda)|\lambda\in\Lambda\}$ you need to compute a dual frame, call it $\{Q(\lambda)|\lambda\in\Lambda\}$.  Then the state is given by
$$
\rho = \sum_\lambda p(\lambda) Q(\lambda) \stackrel{\text{or}}{=} \sum_\lambda \text{Tr}(\rho Q(\lambda)) P(\lambda).
$$
(Just like the $P$ and $Q$ function in quantum optics -- except in reverse notation.)
Unfortunately, you cannot compute a dual frame for an object that is not a frame (i.e. it doesn't span the operator space).  The MATLAB code suggests that this is the case.  It picks a bunch of random angles and computes the measurement operators for them.  Then it checks if they are linearly independent and tries to compute a dual for them.  Running it thousands of times produced no variation in the result that this procedure will produce roughly $O\left(2^{\frac32n}\right)$ linearly independent measurement operators.
% Number of qubits
n = 3;

% ===> Generate random qubit states
num_states = 2^n+2;

% Pick random Bloch sphere directions
phi = 2*pi*rand(1,num_states);
theta = acos(2*rand(1,num_states)-1);
state = [cos(theta/2) ; exp(1i*phi).*sin(theta/2)];
orth_state = [cos((pi-theta)/2) ; exp(1i*(phi+pi)).*sin((pi-theta)/2)];

% Now construct the tensor product states (there should be 2^n for each
% single qubit state since there are this many distinct outcomes) 
big_state = zeros(2^n,2^n*num_states);
for state_idx = 1:num_states
    for perm_idx = 1:2^n
        temp_state = 1;
        binrep = repmat(de2bi(perm_idx-1,n),2,1);
        to_tensor = repmat(state(:,state_idx),1,n).^binrep.*...
                    repmat(orth_state(:,state_idx),1,n).^(1-binrep);
        for sys_idx = 1:n
            temp_state = kron(temp_state,to_tensor(:,sys_idx));
        end
        big_state(:,(state_idx-1)*(2^n)+perm_idx) = temp_state;
    end
end

% ===> Calculate the frame operator

% Vectorize the projectors onto the states
P = zeros(4^n,2^n*num_states);
for state_idx = 1:2^n*num_states
    P(:,state_idx) = kron(big_state(:,state_idx)...
        ,conj(big_state(:,state_idx)));
end

% The frame operator
S = P*P';

% ===> Check for informational completeness and compute dual

% Compute the rank of S
fprintf('Rank of S is %d.\n', rank(S));

% The dual frame
Q = S\P;

% ===> Verify the reconstruction formula

% pick a random state
rho = eye(2^n)/2^n;
rho = rho(:);

% Compute reconstruction error
err = sum(abs(rho - ...
    sum(repmat(rho'*P,4^n,1).*Q,2)...
    ));
fprintf('Error in reconstruction was %d.\n', err);

A: It seems this cannot be done in general for the following reason: We can regard the process as a global rotation $G = R_X^{\otimes n}(\phi) R_Z^{\otimes n}(\theta)$ (i.e. a local rotation through the same angles on each qubit) followed by a $Z$ basis measurement. As $G$ is symmetric, it can be decomposed in terms of symmetric sums of Pauli operators. Since there is exactly one symmetric sum for each choice of $n_I$, $n_X$, $n_Y$ and $n_Z$, the numbers of $I$, $X$, $Y$ and $Z$ local operators in the Pauli terms. Since $n = n_I + n_X + n_Y + n_Z$, there are $\binom{n+3}{3}$ independent terms in including the identity. Since each of these can have a complex coefficient, there are $2\binom{n+3}{3}$ parameters. We have some constraints on these to enforce unitarity, and to avoid counting the global phase. I'll ignore these constraints here, so I will be over counting. Since we make a $Z$ measurement of each qubit and we can consider correlations from arbitrary combinations of these ($2^n - 1$ independent operators), at most we are determining $2(2^n - 1)\binom{n+3}{3}$ parameters.
Note that $$2(2^n - 1)\binom{n+3}{3} = \frac{2(2^n -1) (n+3)!}{n!3!} = \frac{1}{3}(2^n-1)(n+3)(n+2)(n+1)$$
however for an $n$ qubit state there are $4^n - 1$ free parameters, so the number of free parameters in the quantum state rapidly overtakes the maximum number of independent parameters you can measure with such measurements, since if you divide the former by the latter you get $\frac{(n+1)(n+2)(n+3)}{3(2^n+1)}$, which rapidly drops below $1$ (between $n=8$ and $n=9$).
A: I want to give you a general argument for the negative result: A  general $n$-qubit density matrix cannot be contructed from this type of measurements.
The following observation will be needed:
Let $G$ be a compact Lie group acting on a vector space $V$ via a representation $\pi$. Let $\rho$ be a density matrix on $V$.
Consider the quantum characteristic function:
$\phi(g) = \mathrm{Tr}(\rho \pi(g)), \   \ g\in G$, 
When $\pi$ is irreducible, then the density matrix can be reconstructed by virtue of the Peter-Weyl theorem:
$\rho_{ij} = \int \phi(g)\pi_{ij} (g) d\mu(g)\   \ i,j = 1,.,.,.,\mathrm{dim}(V)$
However, in the reducible case, there exists a basis in which all matrix representatives are block diagonal.Thus, in this basis, the matrix elements of the density matrix between basis vectors of different irreducible factors cannot be reconstructed, even given the quantum characteristic function at all points of the group manifold.
To apply that to our case consider the group $G = \bigotimes_{i=1}^{n} U(1) \bigotimes SU(2)$.
Where every Abelian factor acts on a different qubit by phase multiplication and the $SU(2)$ factor acts diagonally on the qubits.
Observing that an SU(2) group element acting on a qubit can be parameterized in terms of the Euler angles:
$g = e^{i\frac{\psi}{2}}P_1(\theta, \phi) +e^{-i\frac{\psi}{2}}P_0(\theta, \phi)$.
Then the quantum characteristic function can be obtained from the series of measurement expectations by a Fourier transform with respect to the Abelian coordinates and $\psi$ and vice versa.
But the group $G$ acts reducibly on the $n$-qubit, since every Abelian factor acts reducibly and $SU(2)$ acts irreducibly only on symmetric powers of qubits (thus reducibly on the full tensor power), which completes the argument.
A: I think that a concise way to understand why one-particle measurements cannot always recover any mixed state is to make use of the following fact: If a POVM does not span the space of operators on the Hilbert space in question, then there always exist two distinct states that give rise to exactly the same measurement probabilities. 
(See http://www.cs.uwaterloo.ca/~watrous/CS766/ProblemSets/solutions1.pdf for a proof). 
For a system of $n$ qubits, the space of operators has dimension $4^n$ but a one-particle measurement has only $2^n$ operators, so it cannot possibly span the space and hence not all mixed states can be distinguished unambiguously.
