After some discussions with @jingli and @schwab, we realized that the peak of the wavelength distribution \phi(\lambda) was not at the same position as the energy distribution \phi(E). This was because, after the change of variables, the corresponding jacobian was missing. This was properly implemented for the angular distribution.
Consider that |\phi(E) dE| = |\phi(\lambda) d\lambda|, and \lambda [\textup{~\AA}] = \sqrt{81.82/E [\mathrm{meV}]}, it is possible to see that dE = |81.82 \cdot (-2) \cdot \lambda^{-3}| d\lambda.
So, for each bin between \lambda_{\mathrm{min}} and \lambda_{\mathrm{max}}, it would be necessary to divide by the integral of the jacobian, i.e. |81.82\cdot(\lambda^{-2}_{\mathrm{min}}-\lambda^{-2}_{\mathrm{max}})|.
For doing this, a new WavelengthAxis class was implemented to plot the distributions with the proper normalization (same idea as AngleAxis). Also, the tracks-vs-tallies.inp example was changed to see the new shape of the wavelength distributions (see attached).
And, with these implementations, now the peak of both distributions is at the same position
