  We compute multi-gluon production in the Color Glass Condensate approach in
dilute-dense collisions, pA, extending previous calculations up to four gluons.
We include the contributions that are leading in the overlap area of the
collision but keep all orders in the expansion in the number of colours. We
develop a diagrammatic technique to write the numerous colour contractions and
exploit the symmetries to group the diagrams and simplify the expressions. To
proceed further, we use the McLerran-Venugopalan and Golec-Biernat-W\"usthoff
models for the projectile and target averages, respectively. We use a form of
the Lipatov vertices that leads to the Wigner function approach for the
projectile previously employed, that we generalise to take into account quantum
correlations in the projectile wave function. We provide analytic expressions
for integrated and differential two gluon cumulants and show a smooth
dependence on the parameters defining the projectile and target Wigner function
and dipole, respectively. For four gluon correlations we find that the second
order four particle cumulant is negative, so a sensible second Fourier
azimuthal coefficient can be defined. The effect of correlations in the
projectile on this result results qualitatively and quantitatively large.
