Hi, I believe in this case the units of X and Y are [cm] as you show a projection of the DPA in the XY plane of your detector.
The secondary particle tracks you see now are merely a result of not collecting enough statistics by not running enough primary beam particles. In the Monte Carlo approach , if you collect results from more primary particle runs you will eventually arrive at a statistically converged distribution (see the examples in these slides).
What you can see from your plot though is that the beam size you define is very small, at coordinates (0,0) the dpa value is very high compared to anywhere else in the detector. This could be a result of also simply be a result of running a very small amount of primary particles. I suggest you restrict the USRBIN ranges to +/- 1mm in the X and Y directions and within this range vary the number of bins, this you can do by defining multiple USRBIN scoring meshes (setting the amount of bins differently each time: for example 10, 50, 100, 200). Compare the results after running more particles and plot the results for example using the “1D Max” option in the plotting tab in Flair. I think in light of this, this discussion might also be useful for your work.
You can gather results for dpa caused by neutron, proton and pion beams by choosing these particles as particle source using the BEAM card, just as you did for neutrons in the inputfile you attached above. Getting the energy dependence will mean you need to run different simulations at different beam energies by also specifying this in the BEAM card.
Hope this is helpful,
Andreas