Optical Photon Simulation

It finally seems to be working without too much of a problem!  This is great news, as the ARI conference is in a week. I haven't really had time to go through the results yet, but they seem to be mostly aligned with the measurements. For the GS20 we are only losing about 8% of the photons, which seems low to me - but the mounting conditions usually account for around a 5-10% change in the peak channel number, so if we are only losing 8% then it kinda works, if you squint hard enough.

GS20 Simulation with Teflon Reflector, BK7 photocathode and silicone optical grease

The simulation is very slow, so I need to work on regaining the parallel capability that I discarded to avoid the code clutter.

In another news we have a big summer league game coming up tomorrow against my future roommate. And, for a limited time, I'm back on Facebook!

Layered Detector Fabrication

The past few days I have been working on measuring semi-large scale layered detector systems. These are layers of plexiglass, and detectors, with light collection only from one end.  So far the results have been a little bit disappointing, but that is the thug research life.

Looking down the fishtail light guide for the first sample (EJ-426 HD2)

Wrapped EJ-426 HD2

 The following is an example of the somewhat poor optical coupling that was achieved. You can see that there are thin bubbles of air in the optical grease that are trapped.  I couldn't figure out how to get rid of them without applying too much optical grease, which tends to degrade detector performance.

Example of poor optical coupling

Layered carborane detector.  The layers are coupled with optical grease.

Example of coupling to the fishtail light guide (carborane detector)

So we will see how these detectors work - in general I think the performance will be much lower than expected and simulated.

Other projects:

•  trying to figure out what's giving me so many gammas in the thermal neutron subtraction (I don't know why this wasn't done before - perhaps when the irradiator was designed and built by a previous student)
• Thesis Proposal
• Project Website - you should really check it out in case you need a bigger picture of this project

MARS!

Super Awesome Range Simulations!

Yay!  I can simulate ranges accurately with GEANT4!  What I did was compare my GEANT4 simulated ranges to those published by the NIST in ASTAR and ESTAR.  A ratio of 1.0 indicates perfect agreement. For the most part (asides from being lazy and only running a low number of Monte Carlo trials) I have very good agreement to the NIST values.

Alpha

Electrons

Super Cool Energy Deposition

A while ago I did a short post about the positional energy deposition in polymeric films.  I recently ran the simulation for 25 cm thick films, and it is cool that it is possible to see a lot of fundamental physics in the simulation.

Neutron Energy Deposition. (I wish I was good enough such that you could rotate the figure)

First, you can see for very thick films pratically of the energy is deposited in the film (the neutron reaction has a q-value of 4.78 MeV). There is some backscattering lost of the front edge of the film, shown by the number of events in the first thickness bin.  This effect is decreased as one goes into the film. Secondly, you can see the exponential attenuation of neutrons through the film by th linear decrease in the frequency of counts as you go through the film at 4.78 MeV -  just taking the slice at 4.78 MeV.  This decrease is linear on a log scale, and therfore follows an exponential attenutation.  From this it would then be possible to esimate the macroscopic absoribition cross section, which determines the rate of attenuation in the medium. Pretty cool, huh?

Compton Scattering & Secondary Electron Kinetic Energy Distrbution

Every once and a while it is good to do basic physics again. This problem was for the validation (again) of my GEANT4 code that simulates the energy deposition and distribution of secondary electrons in a thin film in order to make sure that I am grabbing the right particles / using the correct tracking.

The physics involved are the energy of the Compton scattering of the photon
 \[E_e= E - \frac{E m_e c^2}{m_e c^2 + E (1-\cos\theta)}\]
and the Klein-Nishma scattering cross section
\[ d\sigma = r_0^2 f(\theta)\sin\theta d\theta d\phi \]
where
\[f(\theta) = \frac{1}{2}\left(\frac{E'}{E}\right)^2 \left(\frac{E'}{E} + \frac{E}{E'}-\sin^2\theta\right)\].
After integrating out the phi component it is possible to arrive at
\[\frac{d\sigma}{d \theta} = 2\pi r_e^2 f(\theta) \sin \theta \]

Then taking the derivative of the electron kinetic energy allows for the relationship between the scattering probability and the electron kinetic energy.
 \[\frac{d\sigma}{dE_e} = 2\pi r_e^2 \sin \theta f(\theta)\left [ \frac{1+\frac{E}{m_e c^2}\left(1-\cos\theta \right)^2}{E^2 \sin \theta} \right ]\]

The difficult part was in correctly sampling the distribution, for this I just used a simple rejection method in which the probability of an energy is proportional the the scattering angle of that energy as the PDF (shown above) would be difficult to invert. For this a python script was written. The first part was to make sure that I could reproduce the Compton Scattering cross section; i.e. make sure that my math was correct (Hey, I had to do some fancy chain rule stuff and it has been awhile).
 

Scattering Cross Section as a function of Photon Energy (Co-60)

Finally, I sample from the above distribution and binned it according to a logarithmic energy binning, as shown below. The comparison to the GEANT4 data is also shown as the final figure.

Klein-Nishma Electron Kinetic Energy Spectra

Pretty Good, but I don't have the peak slightly below 1 MeV. 

After all of this it is kinda depressing - I spent so much time just making sure I was right that the simulation code wasn't really worth the effort if in two days you can do a quick analytical check - but hey, it is pretty sweet that I can reproduce things.

Positional Dependance of Energy Deposition

Often times I find myself making vague arguments about the influence of an neutron interaction position and it's energy deposition.  For example, if an interaction takes place "close" (say 10 microns) very little of the energy from the neutron interaction will be deposited in the film as the range of the alpha and triton are both in the 30-50 micron range and will thus leave the film. However, these are very hand-wavy arguments, so I quantified them with a GEANT4 simulation, bining the postion of first interaction and the energy deposition resulting from that interaction.  The detector is a cylinder with the impingent neutrons along the axis of the cylinder; thus it is axial symmetric and the only dimension of intrest is the z.

The following plots (however badly made; I'm still trying to figure out how to nicely plot things in ROOT) are the results of these simulations.

25 micron neutron - most of the energy depostion comes from the alpha, which has a shorter range than the triton

250 micron film neutron.  Can observe that interactions near the edges deposit less energy, are are far less likely.

2,500 micron (2.5mm) neutron. Essentially all of the energy is deposited in the film.

The gamma probablities of energy depsotion are pretty flat even for very thick films, with most of the interactions depositing very little energy, as expected.

10,000 micron (10 mm) film Co60

Incidently, if anyone knows how to space the axis label from the axis in ROOT, please let me know.

Particle Reflectors for Photons

My radio-chemistry professor is an old bloke, and like most old professors is special. In our section on detection he advocated using shielding to increase the count rate of a weak source.  The idea is that the shielding will serve to reflect some of the radiation back towards the detector.  I argued against this, as you will be losing spectra information for each of the collided count as energy will be lost in the quanta of radiation incident.

Four separate geometries were simulated in MCNPX.  The first was a bare HPGe crystal (with the aluminum housing) 10 cm away from a mixed Cs-137 and Co-60 point source. The second (side scattering) was with 1 cm thick lead cylinder encasing the air between the source and the detector. The third (all scatter) was with a 1 cm thick lead  cylinder on the sides, and a 1 cm thick slab below the source.  The final was with a 2" (5.08 cm) thick (all scatter thick bottom) slab below the source and a 1 cm thick side scattering cylinder.  1 cm of lead was chosen because 1.07 cm is the half-thickness for lead (Co-60 source).

X-Z profile of side scattering

X-Z profile of back scattering and side scatter (all scattering)

This also provides cute little examples of solid angle calculations. The simplest fraction of the solid angle the detector subtends is the ratio of the area of the detector to the area of the sphere of radius the distance from the source to the detector is.  For example since the simulated detector is a cylinder along the z-axis with a radius of 2.54 cm, and is 10 cm away from the point source, as is shown below.
\begin{align}
\eta &= \frac{A_{detector}}{A_{Sphere}} \newline
&= \frac{\pi(2.54cm)^2}{4\pi(10cm)^2} \newline
&= 0.016
\end{align}
A more accurate solid angle fraction can be calculated as follows:
\begin{align}
\eta &= \frac{2\pi\left(1-\cos(\theta)\right)}{4\pi}\newline
&=\frac{2\pi\left(1-\cos \left(\tan^{-1}\frac{2.54}{10}\right)\right)}{4\pi} \newline
&= 0.015
\end{align}

Okay, so we do get around a 30% increase in the incident flux with reflection.  This was higher than I expected. In addition, we don't take that much of a hit on our detector response, as shown below. I used  a gaussian energy broadening on the pulse height tally in order to accurate simulate the HPGe response, which is something that I haven't done before and was is really nifty.  What could be super nifty would be to compare this to a measurement, but I don't think that is gonna happen.

Optimal Detector Thickness

As you may recall from the previous post I was investigating the making a layered detector, and I plotted the contribution of a particular film as a function of how far the film was in the detector.  This lead to the idea to ask what would be the best use of material - by making the detector thinner we could avoid some of the low contributions of reactions on the back of the detector.

I held the volume of the detector constant, and adjusted the length and width.  If we make our detector half of what the current He-3 footprint is (which means we need to make them twice as wide) our efficiency goes from 3.1 to 4.6.

Interaction Rate's in a DHS-DNDO Detector

My Master's Thesis is "due" tomorrow, and today I was just finishing up my MCNPX simulations.  And while I should have done this about a month ago, the results are still pretty nifty.  The following figure shows the effect of moderating a Cf-252 source with a 0.5 cm of lead and 2.5 cm of HDPE.  What is neat is that the 2.5 cm of HDPE (lead won't moderate neutrons much because of a neutron cannot impart much recoil energy to the a lead nucleus) is that the spectra is essentially thermalized - the peak is around 0.025 eV, which if you recall is the kinetic energy of room temperature particles predicted by the Boltzmann Distribution.

The detector assembly is composed of a bunch of layers of thin (less than 100 microns) films each separated by acrylic.  The following two plots show the effects of the different layers.  The first layer see's the most neutrons, but doesn't have the largest contribution to the reactions because it only sees them from the front - there isn't much effect of the reflector.  As you go farther into the detecor the number of neutrons crossing a film decrease (it is not because the films absorbs them, only about 20 neutrons are absorbed in all 120 of the films) but rather neutron's escape and scatter out of the detector.

This is probably my favorite plot of my thesis

I looks like I will graduate with a Masters either in May or Summer, and then hopefully (finger's crossed) I can get out of here by May 2013.  Have to do work.

Battleship - Game Play

I found out today that someone else did something similar to what I did: http://thevirtuosi.blogspot.com/2011/10/linear-theory-of-battleship.html, but we are a little bit different.

I calculate what happens with a miss by passing a board where that is an imaginary ship - place all of the ships, and then subtract off the imaginary ship board before returning.  I also changed the order in which I place the ships; placing the larger ships first makes it not take as many searches to find a space where a large ship would fit.

Effects of misses (right) on the (left) PDF.  A miss means that a ship cannot be there so the PDF at that location must be zero.

Anyway, I have the rest of the code written but haven't made any visuals of it.  It is essentially just computing the PDF, and then choosing from that PDF for where to place the ships.  It gets a little bit hung up on deciding what ships are remaining; currently I keep on assuming that I don't eliminate any ship until they all are sunk.

In the event of a nuclear disaster don't "Eat Fresh, Eat Local".   And try not to breathe the air.  You get about 500 times more dose* if you eat contaminated food.  Especially don't drink the milk 2 days after the event.

*depending on what sort of isotope, etc.

Battleship (Ship Placement)

The idea is to find the optimal game play of battle ship; one that ensures the minimum number of loses - because of the random nature of the game (not knowing where the opponents ship's are) it is impossible to ensure a win.

Naive Ship Placement
I then wanted to find a probability density function (PDF) from which to pick my positions; want to sample (take more shots) more where there are more shots. This was first approached numerically and then analytically.  The numerical PDF's were calculated by creating 1 million boards and summing them, and then averaging for the number of boards run (1 million).  Ships are placed anywhere on the board under the following conditions:

• The ship is continuous
• Another ship is not present
•  Entire ship fits on the board

The key point in implementation is that it is only necessary to implement ships placement in one orientation, and just transpose the board to for the other orientation.  The board is simply an (n x n) logical array, where the true (1s) mark the locations of ships.

PDF of ships of length [2,3,3,4,5] being placed on a 10x10 board

 We observe that as less ships are added (ignoring the aircraft carrier of size 5) that the PDF flattens more.

PDF of ships of length [2,3,3,4] place on a 10x10 board

These are not true PDF because they are not normalized to 1, but rather to the amount of ships; i.e. the board of ships [2,3,3,4,5] would sum to 17.

Analytically this was approached by counting the number of ways to arrange a ship of a given size on the board.  All possible combinations are worked out for ships of size 2 and 3.

Possible placements of ships of size two, 9 total.

Possible combinations of ships of size three. 17 possible combinations.

 Given a position (i,j) on the board, we can then calculate  (analytically) the probability that a ship will be placed there, but that will have to wait for tomorrow.

Code for generating Boards:

 function [board] = generateBoard(obj)
            % Allocating space for the board
            board = zeros(obj.n,obj.n);

            i = 1;
            while i <= numel(obj.ships)
                ship = obj.ships(i);

                % Determining Orientation
                ort = (rand>0.5);

                % Determining Position
                x = randi(obj.n);
                y = randi(obj.n);
                shipEndPoints = [x x+ship-1];

                % Placing Ship, if it is a valid position
                if canPlaceShip(obj,shipEndPoints,y,board,ort)
                   board = placeShip(shipEndPoints,y,board,ort);
                   i = i +1;
                end
            end
        end 
end 
function validPosition = canPlaceShip(obj,shipEndpoints,col,board,ort)
if ort;
    board = board';
end
% Logical Values for Board Placement
if (shipEndpoints(1) >= 1 && shipEndpoints(2) <= obj.n) && (0 == sum(board(shipEndpoints,col)))
    validPosition = true;
else
    validPosition = false;
end
end
function board = placeShip(shipEndpoints,col,board,ort)
if ort;
    board = board';
end
% Spots are now filled
board(shipEndpoints(1):shipEndpoints(2),col) = ones(abs(diff(shipEndpoints)+1),1);
% Need to Flip Back
if ort;
    board = board';
end
end

Gamma Interaction Simulations

One of the drawbacks of using MCNPX is a lack of control for the physics options, and its refusal to implement certain features.  So I have taken to playing around with Geant4 again, which means I get to make awesome pictures.

50 Events in a Plastic Scintillator

The results even seem correct, so I guess I'm on the correct track! (Get it, track . . .).   Next step is to add a photo cathode behind it and implement hits collections and sensitive detectors. I've also started to write my thesis; I think I am starting to look forward to getting away from Rocky Top.

In other news our IM Ultimate suffered our first defeat of the season after a very sluggish, hammer throwing, who needs defense first half.  We should still have the head to head and point diff, so I am still on track to win a T-Shirt.

Optimal Detector Thickness

6Li has a higher chance of detecting a neutron if the neutron is lower in energy; this is accomplished by moderating the neutron, good moderators are water, plastic, or anything with a bunch of hydrogen (a neutron can lose all of it's kinetic energy to hydrogen).  We can then increase our detection effiency by moderating the neutron's that come into our detector.  This is the front thickness.  You also want a rear moderator / reflector to reflect neutrons back into the detector; this is the rear thickness.

The data sets are grouped by rear thickness.  The interaction rate will  asymptotically approach some limit for the rear shield thickness, but for practical purposes it needs to be a few (5 cm or so) to have light collection.

Anyway, thought it was very cool.

Intrisinic Efficiency of GS20

In order to prove that I can simulate our spectra (using the same input file as Martian, but that is another story) I had to demonstrate why  our GS20 had an higher efficiency.  Our spectra is not really a beam; it is more isotropic because of the moderation of the HDPE.  Using the attenuation coefficients published by Saint Gobain, the following figure shows the intrinsic efficiency of a 0.2 cm GS20 detector, with the GS20 perpendicular to the beam. 

Intrinsic Efficiency of GS20

The above is the reason for the stacking of the following curves in the below figures.  What I did was simulate the detector for neutron  beams incident on different angles.  The important part to note is that for high energies (greater than a few hundred eV) the effect of orientation is seen.  This is because the low energies have a high intrinsic efficiency, so small values don't make a large contribution.

Intrinsic Efficiency of GS20 for different detector angles

Gamma Count Rates

Been spending a bunch of time in Gamma Town.  It's nice, you don't have to interact with anybody*.  Shown below is comparing results to simulated values from the total interactions. 

Source of all errors

Alright, so not perfect. But hey, I'm gonna take it.  It explains a lot of the problems we have been having; doubling the distance about doubles the counts.   I think at the father distances photons bouncing around off the table and the lab might also be a little bit an issue. It is also unsettling how the simulated and calculated cross each other.  Now the problem is going to make sure that we can have a repeatable geometry; so maybe we want to place it around 6" away, where 1/2" in each direction only changes the count rate by 12%.  At 1/4", moving the sample by 1/4" changes the counts by 30%.

*Groan worthy, but I'm proud of it.  Gamma's are photons, and are characterized by a probability of interaction, not like electrons which are always experiencing E&M fields.

Gamma Spectra Simulations

I've tried to simulate the spectra that our detectors see in MCNPX from a Cs137 source (0.6617 MeV), and Co60 (1.17 and 1.33 MeV).  You can observe the Compton Edge (maximum energy of a Compton scattered electron) in both spectra; notice that the Co60 spectra has two compton edges because of the two photons it emits.  You can see the photopeaks, which is caused by when all of the incident energy is deposited in the detector clearly in the Cs137 spectra.

EJ-200 Gamma Response.  Black is Cs137, blue is Co60

GS20 Gamma Response.  Cs137 is black, Cobalt 60 is blue.

The simulated spectra where then compared. The EJ200 is 1/4" thick 2" diameter, while the GS20 is 1" diameter and about 1/8" thick, that is why the EJ200 has a higher tally than the GS20.  You can notice the extra resolution that using a glass scintillator has a much higher resolution (notice the photo peaks that occur at 1.17 and 1.33 MeV in the black, and the absence of them in the the blue).

Co60 Comparison.  The black is the GS20, the blue is EJ200.

The problem is that the count rate's do not match the observed spectra.  However, the shapes are the same.  Gives hope for tomorrow, anyway.

Observed Spectra from GS20.  Green is the Cobalt 60.

Gonna have to trust me that the shapes are the same, since I don't have a spectra readily available on a log ordinate axis.  The channel numbers are analogous to energy, since the channel number corresponds to a voltage, which in turn corresponds to the amount of light collected by the PMT, which in turn corresponds to the energy deposited by the photon.

Simulated

I have been spending a very long time learning how to use MCNPX to simulate detectors, and an even longer time fiddling with my geometry to ensure what I modeled is actually what is there.  It is an iterative procedure; simulate, look at results, learning something about MCNPX, simulate, change a composition / geometry (mostly adding stuff, Martin's input was very back of the envelope), and continue.  Each simulation takes around 40 hours to run on a single processor, so I run it on the NE computing cluster with anywhere from 30 - 50 cores.  (The cluster only haves 150 cores, and I like to leave space for others to run their jobs.)

Using the (n,triton) reaction card (MT 105) I used an F4 (cell flux) tally modified by an FM card to calculate the interaction rate.  For the GS20 (a standard detector with published efficiency values) I also used an F1 (surface particle) tally to find the number of neutrons (of certian energy ranges) that cross the detector.

Results

For GS20 the simulated interaction rate CPS is  431, the simulated reaction rate from efficiency is 423 cps, and the observed count rate is 428 cps.  Less than a 1% difference.  For the PS based LiF loaded film, I calculated a simulated interaction rates of 108 cps, the measured is 96 cps, a 13% error.  However, this is to be expected since the PS film is opaque, and light generated may be quenched or pushed below the LLD setting, so our measurement may not record all of the light generated.

MCNPX simulations

I recently got cluster access, so now I can simulate much larger systems!  (I also have managed to bring the cluster down twice in about a week . . .).  So what I am doing is simulating the neutron flux and photon flux in our irriadiator and comparing it to how our detector would perform if it was in the DHS configuration.  For the most part, we aren't doing too badly.  Sure that looks like a huge difference between the spectra, but our source is stronger than the DHS.  I was really worried that the DHS spectra would have a different energy structure than ours, but that doesn't look like it is the case.

What is true is that we have a lot more photons crossing our detector than the DHS.  This is a problem - currently we have been subtracting two types of metals with (mostly similar) attenuation to get eliminate the photon response, but since we are using cadmium as the second type of metal we also cut down on our low energy neutrons, which from the above is a large majority of them.

The first energy bin highlights the difference between the neutron spectra

Very different gamma cross sections (and poor statistics)

It does appear, however, that since the cross section for the (n,t) is larger a low energies that our gamma correction by tube subtraction actually will prove to be bothersome.

Tomorrow I go North!

Cold Water Glass

I am a pro at become distracted by side projects.  In any case it is very similar to the coffee cup model, expect that the temperatures will switch.  Each simulation just brings more questions, namely:

• What happens if their is condensation of the outside?  How much energy does it take to for the condensation?  Is the convection more efficient for a film of water compared to silica glass?
• Can I add ice cubes?
• How does this compare to experimental data?
• What are the effects of different geometries?  Would the cup that we (Zilla, Nate, plus other) actually be more effective?

Check out the results for a rather simple cup:

Linear Heat Loss - Convection is proportional to the temperature difference

All of the heat flux is through the top of the glass

I see an ear.