New non-monte-carlo simulation for TUDA

The idea: Quickly produce energy spectra for all strips of a segmented detector, for reactions that populate many excited states in the final nucleus.

Background: During the TUDA 12C+12C experiment, I wrote a monte carlo code to do just this, but to get reasonable statistics in each strip, you need to run 100000 events (about 3 minutes run time on my laptop) per excited state, which means that it's quite labor-intensive to produce spectra with many excited states and several reactions...and of course making mistakes is both inevitable and a pain. Götz Ruprecht asked why it was necessary to do a full event-by-event monte carlo simulation, instead of just finding peak locations and widths from energy loss/straggling codes; I couldn't think of a good reason why not, but I didn't have time to try out that way of doing things--until now.


Approach:

• Calculate the beam's energy loss up to half way through the target. This gives the average beam energy when the reaction takes place. For the sigma of the beam energy distribution, use the energy loss itself (because (a) this is easy and (b) it should be approximately right: the beam will lose somewhere between no energy at all and twice the average energy, and the straggling will be something like 10% of the energy loss).
  
• For a detector at a certain position (I started with a LEDA at 10 cm; I've tried moving it farther back, and am planning to add a forward S2 and a backward LEDA), find the angles corresponding to the inner and outer edges of each strip.
• For each excited energy, cycle through all strips: do the reaction using the inner and outer angles for the strip as the angle of the ejectile; find the initial energy of the ejectile; find the energy lost by the ejectile as it passes through the remaining half of the target (and again use sigma = energy loss); find its energy loss in the dead layer of the detector (this is probably negligible but it is included for completeness; sigma = 10% energy loss since there is no variability in the thickness in this case).
• Now for each strip we have θ_inner, θ_outer, E_{final inner}, E_{final outer}, and σ_E (which is the sum in quadrature of σ for the energy loss of beam in target, ejectile in target, and ejectile in dead layer) . E_{final inner} is always bigger than E_{final outer} (for forward angles), so make the energy spectra as follows: assign a channel number to E_{final inner} (so far I'm using 1 keV per channel), and for higher channel numbers calculate a gaussian distribution with amplitude 1000, mean E_{final inner}, and sigma; for channels below E_{final outer}, calculate a gaussian with mean E_{final outer}; between E_{final inner} and E_{final outer}, write the peak value of the gaussians. (Since both gaussians have the same amplitude and sigma, their peaks will have to be the same size; just in case, there are a couple of lines to calculate both peaks and choose the greater.) The result is a rectangle with softened edges. The shape represents the number of particles of a given energy hitting the strip. The flat top is not quite accurate, since (particularly for a close detector for which each strip spans a broader range of angles) there will be a variation in intensity across the strip; also, the cross section will actually change from strip to strip, which is not reflected in the code. The results, then, do not reflect the peak shapes or relative heights that can be expected, but should reflect the locations and maximum widths.
  
• For additional exitation energies, the spectra are incremented instead of being overwritten, so that a composite of all excitation energies is formed. It would be possible to do multiple reactions in a single run, but I think that it's more useful to do each reaction separately so that you can (for example) plot the alpha particle spectra in blue and the proton spectra in red (or the other way around); having separate spectra for separate reactions makes particle identification easier.

Below I compare the results of this simulation with the event-by-event simulation. The peak locations are the same; the widths are similar (although the new simulation necessarily gives an upper bound on the width and in fact is a slight overestimate because currently the strip width of LEDA is taken to be 5 mm instead of ...4.9 mm? It doesn't account for the inactive gaps between the strips, anyhow.); but the run time for the new simulation was practically instant whereas the event-by-event simulation took ~ 6 minutes. The spectra for the innermost strip of a LEDA at 10 cm is plotted in both cases. The reaction is 80 MeV 12C on a 50 ug/cm2 12C target --> proton ejectile.



And here's a comparison of two detector locations: 10 cm vs 50 cm. The farther position has narrower peaks because of the smaller angular range spanned by each strip, and of course the energies change because the angles are smaller for the farther detector.