  We present Next-to-Leading Order (NLO) predictions for
$\gamma\gamma\gamma\gamma$ final states, including the effects of photon
fragmentation. Our results are calculated fully analytically using the
techniques of $D$-dimensional unitarity, and we discuss some refinements to
existing methods, focusing particularly on the role of three-mass triangle
coefficients. Using these techniques we are able to produce a numerically
stable Monto Carlo code that runs fully in double precision. We investigate the
$4\gamma$ signal at colliders, both for Run II of the LHC and at future
machines. Our results are implemented in MCFM.
