Let's first test which distribution is more likely.

    FindDistribution[data[[All, 2]], 5, All]
[![enter image description here][1]][1]
What I did was ignore the binning with the command data[[All, 2]], and doing that doesn't really harm anything other than shift the data by about 23 to the left. This shows us that a normal distribution is likely good enough, and it would take a 3 parameter Weibull distribution to be competitive. Now this may be a three parameter Weibull distribution. In any case, treating it like a normal distribution one can

    f1 = FindFit[data,c0 PDF[NormalDistribution[mu,sigma], x], {{c0, 1}, {mu, 23}, {sigma, 0.11}}, x]
    Show[ListPlot[data, Filling -> Axis, FillingStyle -> Thickness[0.018],PlotStyle -> Directive[Blue, Opacity[.3]]], Plot[Evaluate[c0 PDF[NormalDistribution[mu, sigma], x] /. f1], {x, 18, 30}, PlotStyle -> Red]]

  [![enter image description here][2]][2]

Which shows the data with its fit function. There is likely some skewness, so one could try a three parameter Weibull distribution. Given the magnitude of the $x$-values, without a location parameter, I don't think a two-parameter Weibull would work as well. 

**Edit**
Worrying about ML or not ML is not a first consideration. Before that let's further examine which distribution is in play. So how did Sjoerd C. de Vries, in his answer, know to try a three parameter Weibull distribution other by doing something like I did above? Here is another way of showing a three parameter Weibull.

    fit=NonlinearModelFit[data2,kk PDF[WeibullDistribution[a0,a1,a2],x],{kk,a0,a1,a2},x,WorkingPrecision->30,MaxIterations->1000];
    Show[ListPlot[data2],Plot[fit[x],{x,18,30}],Frame->True]
    fit["BestFitParameters"]



[![enter image description here][3]][3]

{kk->0.964661496866093819527120602321,a0->4.72300143536200858818774078537,a1->6.80075698971527535007544086964,a2->17.0803838666758819512633107980}

Further steps are to model the data **as is**, e.g., with a Poisson loss function as it is count data (see JimB's answer) and note that we haven't proven anything here. We still need to test distributions for goodness of fit to gain more insight.

Also note this post https://mathematica.stackexchange.com/a/104279/42558, which includes

    (* GerneralizedLinearModelFit - logically the same as maximum
    likelihood estimate *)
    glm = GeneralizedLinearModelFit[data, xx^2, xx,
    ExponentialFamily -> "Poisson", LinearOffsetFunction -> 1,
    IncludeConstantBasis -> False, LinkFunction -> "IdentityLink",
    WorkingPrecision -> 100];

That is, if we have a Poisson process, and model it as such, the claim is that it will not differ from ML.

To do ML, we can use @JimB's guess as to what the number of realizations are in each bin, duplicate that many samples in each bin at the bin location, and just use FindDistributionParameters, which is an ML procedure:

    data3 = Rationalize[{{18.6, 1}, {18.84, 1}, {19.32, 1}, {19.56, 
    5}, {19.8, 3}, {20.04, 3}, {20.28, 11}, {20.52, 7}, {20.76, 
    4}, {21, 12}, {21.24, 13}, {21.48, 9}, {21.72, 22}, {21.96, 
    18}, {22.2, 37}, {22.44, 18}, {22.68, 22}, {22.92, 42}, {23.16, 
    20}, {23.4, 43}, {23.64, 33}, {23.88, 30}, {24.12, 36}, {24.36, 
    28}, {24.6, 32}, {24.84, 19}, {25.08, 18}, {25.32, 10}, {25.56, 
    9}, {25.8, 6}, {26.04, 10}, {26.28, 11}, {26.52, 4}, {26.76, 
    3}, {27, 2}, {27.24, 3}, {27.48, 4}, {27.72, 3}, {27.96, 
    1}, {28.92, 1}, {29.16, 1}}, 0]; 
    datamess = Flatten[Table[data3[[i, 1]], {i, Length[data3]}, {j, 1, data3[[i, 2]]}]];

    H1 = FindDistributionParameters[datamess, WeibullDistribution[a0, a1, a2], WorkingPrecision -> 30]
    QQ1 = QuantilePlot[datamess, WeibullDistribution[a0, a1, a2] /. H1, PlotLabel -> "3 parm Weibull", ImageSize -> Medium];
    H2 = FindDistributionParameters[datamess, WeibullDistribution[a0, a1], WorkingPrecision -> 30]
    QQ2 = QuantilePlot[datamess, WeibullDistribution[a0, a1] /. H2, PlotLabel -> "2 parm Weibull", ImageSize -> Medium];
    H3 = FindDistributionParameters[datamess, NormalDistribution[a0, a1], WorkingPrecision -> 30]
    QQ3 = QuantilePlot[datamess, NormalDistribution[a0, a1] /. H3, PlotLabel -> "Normal", ImageSize -> Medium];
    Print[QQ1, QQ2, QQ3]

Which gives us
[![enter image description here][4]][4]
From the QQ plots, we can see the problem with just duplicating samples, we get a dotted step pattern $...\cdots$. Despite this, it is quite clear that the Normal distribution (see left tail) fits slightly better than the 3 parameter Weibull, which in turn fits much better than the 2 parameter Weibull.

  [1]: https://i.sstatic.net/BzGwb.png
  [2]: https://i.sstatic.net/2ueJgm.png
  [3]: https://i.sstatic.net/0XIjnm.png
  [4]: https://i.sstatic.net/DK1NU.png