Choosing center of histogram bins for fitting I have a bimodal distribution, and if plotted with Mathematica it looks like this:

Now, the lowest value from the actual data is 8196 and 690720, but as seen in the plot, Mathematica lets the data range go from 0 to 744572. Is Mathematica choosing a bad histogram data range?
What is in general a good choice for defining the center of the bins and data range, so that I can fit a distribution through the histogram?
My approach would be:
(1) Calculate (bin width) = (Max-Min)/(number of bins)
    [I'm aware that there are different rules how to choose the optimal number of
    bins depending on the underlying distribution, let's just assume this is 12]

(2) then I have 12 equal bins, starting from 8196 and ending at 690720, each having a width of 56877

(3) The first bin goes then from 8196 to 65073=8196+56877 and so on

(4) As the center of the bin I define the middle between 8196 and 65073 which is 36624.5 and I position my first bin there.

(5) Then I get 12 data pairs of bin center position and number of observations and I can fit a bimodal distribution through it

Am I making a mistake if I do that, or what is the reasoning behind Mathematica's choice of the histogram range exceeding the actual data range?
Edit: I've uploaded the raw data here: raw data
Edit2: For clarity explaining the mysterious frequency of 17 which was pointed out by Stephan: The data is a confocal photoluminescence map where a laser scans an emitter and it looks like this: 

The laser scans row by row so in the middle of each row the emitter lights up which explains the frequency of 17 when the raw data is plotted as it originates from a single list.
 A: As @Nick Cox says, fit your distribution directly to the data. Do not first bin the data into a histogram. Why would you want to do so?
Instead, fit a standard kernel density. I'll use R, because I know it better, but I assume Mathematica has similar functionalities. (If it doesn't, I recommend you learn R.) Below is code that will fit such a density to your data and extract the $x$ value for the second peak.
For added enlightenment, we can assess how sure we are of this coordinate by bootstrapping it. I'm also plotting a bootstrapped 95% quantile. Notice how this is slightly asymmetrical.

dataset <- unlist(read.table("https://files.fm/down.php?i=qucxqxgw"))

foo <- density(dataset)
max.index <- which(foo$x>4e5)[which.max(foo$y[foo$x>4e5])]

plot(foo)
points(foo$x[max.index],foo$y[max.index],pch=19,col="red")
text(foo$x[max.index],foo$y[max.index],round(foo$x[max.index]),pos=3,col="red")

library(boot)
bootstrap <- boot(dataset,statistic=function(dataset,index){
    foo <- density(dataset[index])
    max.index <- which(foo$x>4e5)[which.max(foo$y[foo$x>4e5])]
 foo$x[max.index]
}, R=1e3)

lines(quantile(bootstrap$t,c(0.025,0.975)),rep(foo$y[max.index],2),col="red",lwd=2)

If you want the width of the second peak, you can extract it from the density (and bootstrap it) after you have decided how you define a peak (anything more than 95% of the peak value, or a fixed offset, or something else).
(Yes, in principle we could correct for the fact that your data seem to be all nonnegative, whereas the density estimate goes negative. In practice, since you are only interested in the second peak, I don't really see the point.)
However...
Here is a plot of your original data:

plot(dataset,type="o")

This looks strangely regular. So, after playing around a bit with the frequency parameter, we find the following seasonplot:

library(forecast)
seasonplot(ts(dataset,frequency=17))

Unless you did some very strange sorting to your raw data, your data is actually seasonal with a period of 17. Thus, I'd question whether finding the location of your second mode in such data is really what you want to be doing at all.
A: You can use Mathematica to draw a SmoothHistogram. I use here the data you provided in a flattened list.

data = Flatten[Import["raw.txt", "TSV"]]
Histogram[data]



SmoothHistogram[data]


Now assume you want to fit to a mixture of two distributions, and obtain the parameters and also the mixture, first define a distribution mixture:

distMix = 
   MixtureDistribution[{p, 1 - p}, {NormalDistribution[a, b], 
     NormalDistribution[c, d]}]

And then, obtain the parameters:

params = FindDistributionParameters[data, distMix]

{p -> 0.382944, a -> 516841., b -> 101764., c -> 124503., 
 d -> 81260.4}
