Physicists never work this way, this is an extreme caricature which would be a joke if it weren't tragic. Veneziano was not fitting experimental data to a mathematical function, he was using general principles to deduce what form a tree-level self-consistent scattering amplitude in a theory of infinitely many particles on straight line Regge trajectories, obeying Dolen-Horn-Schmidt duality, must have.
This story is something that people made up afterwards to disguise origin of string theory in the S-matrix theory of the 1960s. The S-matrix people that made string theory all were political losers, their theory was too advanced for people without internet to understand and disseminate. Thanks to the foresight and political clout of Witten, the situation reversed in 1984, when Witten understood that string theory actually made sense. But most of the people who took up string theory in the 1980s were still dismissive of Regge theory and S-matrix theory, because the 1960s battle between S-matrix and field theory ended in 1974 with total victory for field theory, and all these people were on the winning side.
Most of the original S-matrix folks jumped ship too, out of political necessity (and the undeniable fact that nonabelian gauge field theory was interesting and correct for the strong interactions, and also there was a lot of low-hanging fruit after t'Hooft and Veltman's revolutionary work). Most of the original string theorists that developed proto-string theory in 1968-1974 changed fields: they studied field theory from 1974 to 1984. This includes Veneziano, Gross, and Susskind. All these folks re-embraced string theory again in 1984 when it was politically ok to do so. You can't really blame them, as anybody studying S-matrix was dismissed in the 1970s, as part of the baby-boomer conservative counter-revolution. Anything smelling of logical-positivism was thrown out of academia, including S-matrix theory.
The fake story that "Veneziano looked in a math book, found the beta-function, and noticed it was a scattering amplitude" was made up by people who didn't understand what Veneziano was doing. This is not true, and Veneziano has been bitter about his work (which took hard thinking and some of the most inspired mathematical and physical reasoning in physics history) being misrepresented in this way.
Veneziano derived the open-string tachyon scattering amplitude from principles of Regge theory and S-matrix theory, and was forced to use the beta-function to make the following conditions all work:
- The scattering amplitude (at tree level) should have a pole whenever the energy hits a particle mass.
- There should be infinitely many particles lying along a trajectory function $\alpha(s)$ which relates the spin of the particles to the square-mass by a straight line with universal slope (all the mesons lie on these lines).
- The residue of the pole should be positive (required unitarity).
- The residue of the pole in s as a function of t should be the appropriate polynomial in t, to reproduce the required angular dependence of scattering of particles of a varying spin.
- The scattering amplitude should be analytic away from the poles (to tree order).
- The amplitude should be crossing invariant (this is why Veneziano's amplitude is a sum of three beta functions).
- The scattering amplitude should obey the apprpriate regge scaling in s, so that the width of the scattering at high energies should shrink appropriately.
- The mass of the particle that is coming in to do the scattering should be identifiable with the particles in the theory (you should be able to put any of the exchanged intermediate particles on the outgoing legs).
and the most important property (inducted from experiments on the strong interactions by Dolen Horn and Schmidt):
- The amplitude should not be a sum over t and s exchange of particles, like in field theory, but the t-exchange should be dual to s-exchange, in that summing over all particles exchanged in the t-channel must produce the s-channel poles.
These conditions are so stringent (especially the last), that people were claiming to prove that there is no solution in 1968. Veneziano put an end to these claims, by giving what is essentially the unique solution to these constraints. Starting with the B function (more or less the beta function)
$$ B(s,t) = {\Gamma(-(s-\alpha(s)))\Gamma(-(t-\alpha(t)))\over\Gamma(-(s+t-\alpha(s+t)))} $$
and symmetrize by crossing to make a scattering amplitude
$$ A(s,t) = B(s,t) + B(u,t) + B(s,u) $$
(note that $B(x,y)=B(y,x)$) Where s,t,u are the Mandelstam variables, and $\alpha(x) = ax+b$ is the straight line Regge trajectory, you solve the self-consistency conditions at tree level.
This history is summarized in two excellent review articles, both appearing in "Dual Models" (a reprint collection from 1984). One article is Mandelstam's Dual-Resonance Models from 1974.
Veneziano found the function B(x,y) by his own methods, which took deep profound Regge theory intuition. The first thing is that a Regge trajectory has poles at infinitely many positive positions. The gamma function
$$ \Gamma(s) = \int_0^\infty x^{s-1} e^{-x} dx $$
has the property that it has poles at all negative integer positions, which follows from the functional equation:
$$ \Gamma(s+1)=s\Gamma(s) $$
Which defines $\Gamma$ as the generalization of the factorial function to noninteger values: $\Gamma(n)=(n-1)!$. Using the functional equation for negative values of s, you get a division by 0 as soon as you get to 0, and you get poles with a residue that increases factorially at negative positions.
This was well known in Regge theory world, and the Gamma function was used for phenomenological Regge amplitudes throughout the late 1960s. But these phenomenological fits didn't work to make a complete consistent theory, just to get a few scattering amplitudes more or less right to lowest order in Regge exchange. You can find these types of amplitudes everywhere in the literature in the late 60s. These types of amplitudes work to give the near-beam scattering, which has some paradoxical properties--- it shrinks in angle at high energies according to a superposition of power-laws, each separate power-law corresponds to the exchange of a sequence of regularly spaced particles of higher mass and spin, which are related by $m^2 = a J + b$ (a straight line of mass-squared to spin) with the coefficient a the same for all trajectories (except the pomeron, but ignore this), and b different from trajectory to trajectory.
The integral form of the gamma function was also known to be the Regge theory analog of the Schwinger representation of particle propagators by proper time (but this is not quite right--- the proper generalization is the world-sheet, but this came much later). Nobody before Veneziano had a consistent generalization using Gamma functions of the field theory propagators, since the particles were non scattering fully consistently at all energies.
The beta function can be thought of as (the reciprocal of) a generalization of the binomial coefficient to noninteger values:
$$ {a +b \choose b} = {(a+b)!\over a!b!} $$
Which, in gamma form, and dropping some obvious factors which just shift things by one unit, is more or less
$$ {\Gamma(x+y) \over \Gamma(x)\Gamma(y)} = {1\over \beta(x,y)}$$
The thing that happens here is that the double zeros you expect when x and y are both negative integers are smoothed out to single zeros because they are cancelled by the $\Gamma(x+y)$. This cancellation is required by unitarity, because the reciprocal of this is the amplitude, and you can't have a double pole in a scattering amplitude. You need to cancel the double poles, and this essentially uniquely determines the third gamma factor (since the poles of an analytic function determine the function up to asymptotics).
This form is essentially one of only two simple amplitudes that obey the condition that it can be thought of as tree level exchange linear Regge trajectories (the other is the closed string amplitude of Koba and Nielsen, which involves a similar product of three gamma functions).
In his early papers, Veneziano noticed that Euler had already defined this function (for other reasons), and gave a bunch of interesting identities for it, including
$$ \beta(x,y) = \int_0^1 q^{x-1} (1-q)^{y-1} dq$$
This identity might be what he found in a math book (maybe, maybe not, you can ask him), and this identity was all-important: it allowed Veneziano to make a 1-d picture of the scattering, where there is an internal sigma variable which controls "where" each contribution to the scattering was happening. By inserting the appropriate factor at the right position $q$, Veneziano could generalize his scattering amplitude to find the right amplitude for scattering the other particles in the theory. This amplitude isn't arbitrary, because he already knew where all the particles in the theory were in terms of mass and spin, from the tachyon amplitude. But the q picture gave him the idea to introduce local operators to change the type of particles.
So Venziano and Fubini introduced vertex operators, which are localized in $q$. The variable $q$ is now understood to paramertrize the open-string endpoint position where the particles enter the disk after a conformal map, but this took a while to develop. The vertex operator formalism is essentially equivalent to modern string world-sheet, but it required recognizing that there is an internal string in the theory, something which isn't obvious. This was done by Susskind, Nielsen, and Nambu, and also by Ramond each of which had a slightly different picture of what the internal string picture was like. Susskind and Nambu had the right idea, and then string theory develops logically from there.
It is impossible to learn the proper history of string theory from field theorists, or from popular works. Reading the original articles is essential. There was also a string of review and history articles by the original authors in recent years on arxiv, which helps sort the history out properly.