Fit data to parametric distribution I have data with nice bell-shaped histogram PDF. However, the Normal distribution fitting (by calculating mean and variance) does not work as the figure below.

My question is that if there are other distributions I should try according to your experience. In other words, by which distribution should a nice bell-shaped data, which is not well fit by Normal distribution, be fit.
My ultimate goal is to have an approximated analytic form of cummulative distribution function to analyze the according probability. Thus any advices towards solving this goal are appreciated.
I include the data (space as delimiter):
Data to fit.
Update: q-q plot

 A: The histogram, as presented by the OP, gives the impression that the data is symmetrical. Given that the data is noticeably more peaked than Normal, and if the data is roughly symmetrical, then a natural suggestion to try is the Student's t  with location parameter $\mu$, scale parameter $\sigma$, and $v$ degrees of freedom, and pdf $f(x)$:
$$f = \frac{1}{\sigma  \sqrt{v} \; B\left(\frac{v}{2},\frac{1}{2}\right)} \left(\frac{v}{v+\frac{(x-\mu )^2}{\sigma ^2}}\right)^{\frac{v+1}{2}} \quad \text{defined on the real line}$$
Student t fit
The following diagram shows a sample fit using the Student's t, with $\mu = 5.45$, $\sigma = 6.61$ and $v = 2.97$:

In the diagram:


*

*the dashed red curve is the fitted Student's t pdf

*the squiggly blue curve is the empirical pdf (frequency polygon) of the raw data
On the upside, this appears to be a significantly better fit than the Normal, using the same raw data set provided.
On the possible downside, I am not sure I would fully agree with the OP's opening statement:  "I have data with nice bell-shaped histogram PDF". In particular, if one looks more closely at your data set (which contains 100,000 samples), the maximum is 37.45, while the minimum is -910. Moreover, there is not just one large negative value, but a whole bunch of them. This suggests that your data set is not symmetrical, but negatively skewed ... and that there other things going on in the tails, and if so, other distributions may perhaps be better suited. Zooming out, again with the same Student's t fit, we can see this feature of the data, in the right and left tails:

A: In short: your two plots show a big discrepancy, the smallest values shown in the histogram is about $-30$, while the qqplot shows values down to around $-900$.  All those lom-tailed outliers is about 0.7% of the sample, but dominates the qqplot. So you need to ask yourself what produces those outliers! and that should guide you to what to do with your data. If I make a qqplot after eliminating that long tail, it looks much closer to normal, but not perfect. Look at these:
mean(Y)
[1] 3.9657
mean(Y[Y>= -30])
[1] 4.414797

but the effect on standard deviation is larger:
sd(Y)
[1] 10.92237
sd(Y[Y>= -30])
[1] 8.006223

and that explains the strange form of your first plot (histogram): the fitted normal  curve you shows is influenced by that long tail you omitted from the plot. 
A: You might try a Gaussian mixture, which is easy using Mclust in the mclust library of R.
library(mclust)
mc.fit = Mclust(data$V1)
summary(mc.fit,parameters=TRUE)
This gives a three-component Gaussian mixture (8 parameters total), with components
1: N(-69.269908, 6995.71627), p1 = 0.003970506 
2: N( -4.314187,  171.76873), p2 = 0.115329209 
3: N(  5.380137,   46.26587), p3 = 0.880700285
The log likelihood is -352620.4, which you can use to compare other possible fits such as those suggested.  
The long left tail is captured by the first two components, especially the first.
The cumulative distribution estimate at "x" is (in R form)
p1*pnorm(x, -69.269908, sqrt(6995.71627)) + p2*pnorm(x, -4.314187, sqrt(171.76873)) 
  + p3*pnorm(x, 5.380137, sqrt(46.26587)) 
I tried various quantiles (x) from .0001 to .9999 and the accuracy of the estimate seems reasonable to me. 
