The apparently random pattern of the distribution of the zeroes of the Riemann zeta function (all on the critical line, according to the Riemann hypothesis) has led to the suggestion that it might be related to the spectrum of a `Hamiltonian'. It has also been known for some time that the statistical properties of the eigenvalue distribution of an ensemble of random matrices resemble those of the zeroes of the zeta function. With the objective to get to this operator, we assume Riemann hypothesis, start with the zeroes, and construct a unitary matrix model (UMM) for the zeta function. Our approach to the problem, however, is `piecemeal'. That is, we consider each factor in the Euler product representation of the zeta function to get a UMM for each prime. This helps to write the partition function as a trace and suggests an operator. We then assemble these to get a matrix model for the full zeta function. Similar construction seems to work for a family of related functions.