An event generator for diphoton ($gammagamma$) production in hadron collisions that includes associated jet production up to two jets has been developed using a subtraction method based on the LLL subtraction. The parton shower (PS) simulation to restore the subtracted divergent components involves both QED and QCD radiation, and QED radiation at very small $Q^{2}$ are simulated by referring to a fragmentation function (FF). The PS/FF simulation has the ability to enforce the radiation of a given number of energetic photons. The generated events can be fed to PYTHIA to obtain particle (hadron)-level event information, which enables us to perform realistic simulations of photon isolation and hadron-jet reconstruction. The simulated events, in which the loop-mediated $gg rightarrow gammagamma$ process is involved, reasonably reproduce the diphoton kinematics measured at the LHC. Using the developed simulation, we found that the 2-jet processes significantly contribute to diphoton production. A large 2-jet contribution can be considered as a common feature in electroweak-boson production in hadron collisions although the reason is yet to be understood. Discussion concerning the treatment of the underlying events in photon isolation is necessary for future higher precision measurements.