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.

    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 

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