And thinking a little more, this looked suspiciously like a superposition of Gaussians. Does theory justify such an interpretation?

In[31]:= fit2 = 
  NonlinearModelFit[data, 
   a Exp[-(-b + x)^2/c] + 
    a1 Exp[-(-b1 + x)^2/c1], {{a, 700}, {b, 3490}, {c, 20}, {a1, 
     700}, {b1, 3490}, {c1, 20}}, x];

In[32]:= fit2["BestFitParameters"]

Out[32]= {a -> 456.927, b -> 3490.81, c -> 106.369, a1 -> 441.68, 
 b1 -> 3470.54, c1 -> 766.738}

In[33]:= Show[Plot[fit2[x], {x, 3400, 3520}, PlotStyle -> Red], 
 ListPlot@data, PlotRange -> All]

The exponent of the Gaussian needs to be negative. I used NonlinearModelFit below, but FindFit would work.

In[25]:= fit = 
  NonlinearModelFit[data, 
   a Exp[-(-b + x)^2/c], {{a, 700}, {b, 3490}, {c, 20}}, x];

In[26]:= fit["BestFitParameters"]

Out[26]= {a -> 622.26, b -> 3479.8, c -> 728.5}

In[27]:= Show[Plot[fit[x], {x, 3400, 3520}, PlotStyle -> Red], 
 ListPlot@data, PlotRange -> All]

You need "Exp" rather than "exp". (But the plot of the data suggests that a such a curve will not be a very good fit.)

Thanks, Marco -- I was not thinking at all clearly about this problem. We are given a PDF -- to fit a distribution we have to reconstitute the samples themselves, as you did with data2. Following you, I tiried to fit this to a superposition of normal distributions, with the same result.

In[241]:= candidate2 = 
 MixtureDistribution[{a1, a2}, {NormalDistribution[u1, s1], 
   NormalDistribution[u2, s2]}]

Out[241]= MixtureDistribution[{a1, a2}, {NormalDistribution[u1, s1], 
  NormalDistribution[u2, s2]}]

In[242]:= par2 = FindDistributionParameters[
  data2,
  candidate2,
  {{a1, 1000}, {a2, 1000}, {u1, 3460}, {s1, 10}, {u2, 3490}, {s2, 10}}]

During evaluation of In[242]:= FindMaximum::eit: The algorithm does not converge to the tolerance of 4.806217383937354`*^-6 in 500 iterations. The best estimated solution, with feasibility residual, KKT residual, or complementary residual of {1.83244*10^-9,0.00273008,4.85286*10^-10}, is returned. >>

Out[242]= {a1 -> 1187.93, a2 -> 866.089, u1 -> 3462.92, s1 -> 21.0563,
  u2 -> 3488.61, s2 -> 9.17502}

In[243]:= dist = candidate2 /. par2;

In[244]:= Show[
 Plot[Length[data2] PDF[dist][x], {x, 3400, 3520}, PlotStyle -> Red], 
 ListPlot@data, PlotRange -> All]

While FindDistributionParameters might provide a reasonable "fit" that well describes the observed data (a mixture of two or more normals for example), there is no random sampling from a mixture of normal distributions or a mixture of Poisson distributions. The counts likely have a random component (approximately Poisson) but not the particle energy variable (as those values are uniformly spaced and not randomly distributed).

It appears that what you have is a Poisson regression with the particle's energy as the fixed and known predictor and the counts as the dependent variable. As the plot of the data shows, this is not a simple linear or quadratic relationship. I think what you need is a generalized additive model (gam) but I'm unaware of a Mathematica-supplied function for that. However, you could call the appropriate function in R from Mathematica using R's mgcv library.

An update: I'll backtrack a bit (well, a whole lot) if the equal-spaced energy values and the counts represent a histogram with the particle energy values parked in the appropriate bins. Then you are in the mode of fitting a distribution rather than a Poisson regression. But if the particle energy levels are specifically set and increased from the lowest value to the highest, then a Poisson regression is more appropriate.

This is a spectrum of alpha particles. On the horizontal axis are integer values proportional to the particle's energy, on the vertical axis are counts. This is actually an overlap of two gaussians/poisson distributions, because there are two decay processes, an excited and an unexcited one.

The most interesting information from the data are the amplitude, peak position and variance of the peak with the higher energy. Therefore FindDistributionParameters seems to be a good choice.