Source code for ixpeobssim.irfgen.alphas

#!/usr/bin/env python
#
# Copyright (C) 2020--2021, the ixpeobssim team.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

"""Small module providing simulation tools for the alpha particles emitted by
the Be window.
"""

from __future__ import print_function, division

import numpy

from ixpeobssim.irfgen.gpd import WINDOW_BE_THICKNESS, WINDOW_AL_THICKNESS
from ixpeobssim.irfgen.gpd import ABSORPTION_GAP_THICKNESS
from ixpeobssim.instrument.gpd import GPD_PHYSICAL_HALF_SIDE_X
from ixpeobssim.core.geometry import xPoint, xLine, xRay
from ixpeobssim.irfgen.astar import load_alpha_stopping_power_data
from ixpeobssim.irfgen.astar import load_dme_alpha_stopping_power_data
from ixpeobssim.irfgen.constants import BE_DENSITY, AL_DENSITY
from ixpeobssim.utils.matplotlib_ import plt
from ixpeobssim.core.hist import xHistogram1d
from ixpeobssim.utils.logging_ import logger
from ixpeobssim.utils.math_ import modulo_2pi


WINDOW_HALF_SIDE = 4. #cm
ASIC_HALF_SIDE = GPD_PHYSICAL_HALF_SIDE_X / 10. #cm
# Setup the geometry---all dimensions in cm.
Z_ASIC = 0.
Z_AL_BOT = Z_ASIC + ABSORPTION_GAP_THICKNESS
Z_BE_BOT = Z_AL_BOT + WINDOW_AL_THICKNESS
Z_BE_TOP = Z_BE_BOT + WINDOW_BE_THICKNESS
logger.info('Basic detector geometry...')
logger.info('ASIC at %.6f cm', Z_ASIC)
logger.info('Window Al bottom at %.8f cm', Z_AL_BOT)
logger.info('Window Be bottom at %.8f cm', Z_BE_BOT)
logger.info('Window Be top at %.8f cm', Z_BE_TOP)

BE_ASTAR_DATA = load_alpha_stopping_power_data('Be', BE_DENSITY)
AL_ASTAR_DATA = load_alpha_stopping_power_data('Al', AL_DENSITY)
DME_ASTAR_DATA = load_dme_alpha_stopping_power_data()



def _inside(val, minimum, maximum):
    """Convenience function returning True if a value is within a min and max.
    """
    return val >= minimum and val <= maximum

def _zinside(val):
    """
    """
    return _inside(val, Z_ASIC, Z_AL_BOT)

def _xyinside(val):
    """
    """
    return _inside(val, -ASIC_HALF_SIDE, ASIC_HALF_SIDE)

[docs] def gas_cell_intersection_points(ray): """Calculate the intersecting points of a given ray with the gas cell. This is a horrible list of conditions, and should be """ w = ASIC_HALF_SIDE east = ray.yintersect(ASIC_HALF_SIDE) west = ray.yintersect(-ASIC_HALF_SIDE) north = ray.xintersect(ASIC_HALF_SIDE) south = ray.xintersect(-ASIC_HALF_SIDE) top = ray.zintersect(Z_AL_BOT) bottom = ray.zintersect(Z_ASIC) points = [] for point in [east, west]: if _zinside(point.z) and _xyinside(point.x): points.append(point) for point in [north, south]: if _zinside(point.z) and _xyinside(point.y): points.append(point) for point in [top, bottom]: if _xyinside(point.x) and _xyinside(point.y): points.append(point) num_points = len(points) # If there are no valid intersection points, return None and handle this # downstream. if num_points == 0: return None # At this pooint the intersection points must be exactly 2. assert num_points == 2 # Since the alpha is travelling downward, make sure that the z of the first # point is higher than the z of the second one. p1, p2 = points if p1.z < p2.z: points = [p2, p1] return points
[docs] def generate_decays(alpha_energy=5., num_alphas=1000000, plot=True, energy_cut=0.1): """Generate the initial alpha particles in the Be window volume. This is returning a list of xRay object and a numpy array of energies describing the alpha particles that emerge from the window in a direction that is intersecting the fiducial cylinder enclosing the GPD active gas volume. This first part of the program is vectorized, as only a small fraction of the particles survive the initial cut, and the whole thing would be excruciatingly slow, if stuffed into a plain Python loop. """ logger.info('Generating %d initial alpha decays...', num_alphas) # Generate the z decay poition within the Be window. z = numpy.random.uniform(Z_BE_BOT, Z_BE_TOP, size=num_alphas) # Generate the z cosine director uniformly in the half-sphere. zdir = numpy.random.random(size=num_alphas) # Calculate the path length in the Berillium be_path_length = (z - Z_BE_BOT) / zdir logger.info('Average decay z: %.8f cm', z.mean()) # Calculate the range of an alpha particle in the Be. be_range = BE_ASTAR_DATA.csda_range(alpha_energy) logger.info('Alpha range in Be: %.8f um', be_range * 1.e4) # Initial debug plots. if plot: plt.figure('Conversion point') binning = numpy.linspace(Z_BE_BOT, Z_BE_TOP, 100) h = xHistogram1d(binning, xlabel='z conversion point [cm]').fill(z) h.plot() plt.figure('Path length') binning = numpy.linspace(0., 0.025, 100) h = xHistogram1d(binning, xlabel='Path length in Be [cm]').fill(be_path_length) h.plot() plt.vlines(be_range, 0., num_alphas, ls='dashed', color='gray') # Throw away all the events that would not make it out of the Be window. mask = be_path_length < be_range n = mask.sum() # Calculate the fraction of alpha escaping the Be Window # (note the factor of 0.5 is due to the fact that we are generating over # 2pi steradians.) frac = 100. * 0.5 * n / num_alphas logger.info('%d alpha(s) escaping the Be (%.2f%%)...', n, frac) # Select the particles that survive. be_path_length = be_path_length[mask] z = z[mask] zdir = zdir[mask] # Calculate the energy of the surviving alpha particles at the Be exit. be_energy_profile = BE_ASTAR_DATA._energy_profile(alpha_energy) energy = be_energy_profile(be_path_length) if plot: plt.figure('Be energy profile') be_energy_profile.plot(grids=True) # Now the energy loss in the Al layer. al_path_length = WINDOW_AL_THICKNESS / zdir al_energy_loss = AL_ASTAR_DATA.stopping_power(energy) * al_path_length logger.info('Average energy loss in the Al layer: %.2e MeV', al_energy_loss.mean()) mask = energy > al_energy_loss n = mask.sum() frac = 100. * 0.5 * n / num_alphas logger.info('%d alpha(s) escaping the Al (%.2f%%)...', n, frac) z = z[mask] zdir = zdir[mask] energy = energy[mask] # Apply the minimum energy cut. mask = energy > energy_cut n = mask.sum() frac = 100. * 0.5 * n / num_alphas logger.info('%d alpha(s) above %.2f MeV (%.2f%%)...', n, energy_cut, frac) z = z[mask] zdir = zdir[mask] energy = energy[mask] # More debug plots. if plot: plt.figure('Conversion point (gas volume)') binning = numpy.linspace(Z_BE_BOT, Z_BE_TOP, 100) h = xHistogram1d(binning, xlabel='z conversion point [cm]').fill(z) h.plot() plt.figure('Alpha energy (gas volume)') binning = numpy.linspace(0., alpha_energy, 100) h = xHistogram1d(binning, xlabel='Alpha energy [MeV]').fill(energy) h.plot() # Generate the starting decay position within the Be window. x = numpy.random.uniform(-WINDOW_HALF_SIDE, WINDOW_HALF_SIDE, size=n) y = numpy.random.uniform(-WINDOW_HALF_SIDE, WINDOW_HALF_SIDE, size=n) phi = 2. * numpy.pi * numpy.random.random(size=n) # Additional geometrical cut to throw away all the rays that do not intercept # the circle enclosing the active volume of the gas. # See https://math.stackexchange.com/questions/543496 for the basic math. r = ASIC_HALF_SIDE * numpy.sqrt(2.) d = numpy.sqrt(x**2. + y**2) rho = r / d b = d * rho * numpy.sqrt(1. - rho**2.) phi0 = numpy.pi + numpy.arctan2(y, x) delta_phi = numpy.arccos(b / r) mask = numpy.logical_or((rho > 1.), (abs(modulo_2pi(phi - phi0)) < delta_phi)) n = mask.sum() frac = 100. * 0.5 * n / num_alphas logger.info('%d alpha(s) within fiducial cylinder (%.2f%%)...', n, frac) # Prepare the output rays. x = x[mask] y = y[mask] z = numpy.full(x.shape, Z_AL_BOT) theta = numpy.degrees(numpy.arccos(-zdir[mask])) phi = numpy.degrees(phi[mask]) energy = energy[mask] # Even more debug plots. if plot: plt.figure('Alpha energy (fiducial cylinder)') binning = numpy.linspace(0., alpha_energy, 100) h = xHistogram1d(binning, xlabel='Alpha energy [MeV]').fill(energy) h.plot() if plot: plt.show() rays = [] for _x, _y, _z, _theta, _phi in zip(x, y, z, theta, phi): p = xPoint(_x, _y, _z) ray = xRay(p, _theta, _phi) rays.append(ray) return tuple(rays), energy
[docs] def plot_rectangle(width, height, center=(0., 0.), **kwargs): """Facility to plot a generic rectangle. """ x0, y0 = center x1 = x0 - 0.5 * width x2 = x0 + 0.5 * width y1 = y0 - 0.5 * height y2 = y0 + 0.5 * height plt.plot([x1, x2, x2, x1, x1], [y1, y1, y2, y2, y1], **kwargs)
[docs] def plot_event(ray, enery, track_start=None, track_end=None): """Rudimentary single-event display. """ p0 = ray.origin p1 = ray.zintersect(Z_ASIC) ray_fmt = dict(ls='dashed', color='orange', lw=1.) window_fmt = dict(ls='solid', color='gray', lw=1.) gas_fmt = dict(ls='solid', color='blue', lw=1.) fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, num='Event display', figsize=(9., 9.)) plt.sca(ax3) plot_rectangle(2. * WINDOW_HALF_SIDE, 2. * WINDOW_HALF_SIDE, **window_fmt) plot_rectangle(2. * ASIC_HALF_SIDE, 2. * ASIC_HALF_SIDE, **gas_fmt) plt.axis([-5., 5., -5., 5.]) plt.gca().set_aspect('equal') plt.xlabel('x [cm]') plt.ylabel('y [cm]') plt.plot([p0.x, p1.x], [p0.y, p1.y], **ray_fmt) for point in [track_start, track_end]: plt.plot([track_start.x, track_end.x], [track_start.y, track_end.y], 'o') plt.sca(ax1) plot_rectangle(2. * WINDOW_HALF_SIDE, ABSORPTION_GAP_THICKNESS, center=(0, 0.5 * ABSORPTION_GAP_THICKNESS), **window_fmt) plot_rectangle(2. * ASIC_HALF_SIDE, ABSORPTION_GAP_THICKNESS, center=(0, 0.5 * ABSORPTION_GAP_THICKNESS), **gas_fmt) plt.axis([-5., 5., -5., 5.]) plt.gca().set_aspect('equal') plt.xlabel('x [cm]') plt.ylabel('z [cm]') plt.plot([p0.x, p1.x], [p0.z, p1.z], **ray_fmt) for point in [track_start, track_end]: plt.plot([track_start.x, track_end.x], [track_start.z, track_end.z], 'o') plt.sca(ax4) plot_rectangle(2. * WINDOW_HALF_SIDE, ABSORPTION_GAP_THICKNESS, center=(0, 0.5 * ABSORPTION_GAP_THICKNESS), **window_fmt) plot_rectangle(2. * ASIC_HALF_SIDE, ABSORPTION_GAP_THICKNESS, center=(0, 0.5 * ABSORPTION_GAP_THICKNESS), **gas_fmt) plt.axis([-5., 5., -5., 5.]) plt.gca().set_aspect('equal') plt.xlabel('y [cm]') plt.ylabel('z [cm]') plt.plot([p0.y, p1.y], [p0.z, p1.z], **ray_fmt) for point in [track_start, track_end]: plt.plot([track_start.y, track_end.y], [track_start.z, track_end.z], 'o')
class AlphaTrack: def __init__(self, start, end, deposited_energy): """Constructor. """ self.start = start self.end = end self.deposited_energy = deposited_energy self.length = self.start.distance_to(self.end) def projected_length(self): """Return the length of the track, projected on the xy plane. """ dx = self.end.x - self.start.x dy = self.end.y - self.start.y return numpy.sqrt(dx**2. + dy**2.) def roi_size(self, padding=0.025, pixel_area=2.17e-5): """ """ dx = abs(self.end.x - self.start.x) + 2. * padding dy = abs(self.end.y - self.start.y) + 2. * padding return int(dx * dy / pixel_area) def track_size(self, width=0.05, pixel_area=2.17e-5): """ """ return int(self.length * width / pixel_area)
[docs] def simulate(alpha_energy=5.0, num_alphas=100000, plot=False, energy_cut=0.1): """Run the actual simulation. """ track_list = [] # Setup the energy loss for the DME. dme_energy_loss = DME_ASTAR_DATA.energy_loss(alpha_energy) # Create the list of particle direction and energy for the alpha particles # that emerge from the window with a direction intersecting the fiducial # cylinder enclosing the active gas cell. logger.info('Propagating tracks...') for ray, energy in zip(*generate_decays(alpha_energy, num_alphas, plot)): points = gas_cell_intersection_points(ray) # If the event does not hit the active gas volume, there is nothing to do. if points is None: continue # Unpack the enter and exit pointo into/from the gas active volume. track_start, track_end = points # Calculate the energy at the entry point in the active gas volume. energy -= dme_energy_loss(energy, ray.origin.distance_to(track_start)) # Does the particle die before reaching the active gas volume? if energy <= energy_cut: continue # Build the actual track and energy loss. track_length = track_start.distance_to(track_end) deposited_energy = dme_energy_loss(energy, track_length) # If the particle looses all the remaining energy, it dies within the # active volume, and we have to go back and re-calculate the track # endpoint. if energy - deposited_energy < 1e-6: track_length = DME_ASTAR_DATA.csda_range(energy) deposited_energy = energy track_end = track_start.move(track_length, ray) track_list.append(AlphaTrack(track_start, track_end, deposited_energy)) if plot: plot_event(ray, energy, track_start, track_end) plt.show() n = len(track_list) frac = 100. * 0.5 * n / num_alphas logger.info('Done, %d good tracks generated (%.2f%%).', n, frac) return track_list
[docs] def generate_distributions(alpha_energy=5.0, num_alphas=10000000): """ """ track_list = simulate(alpha_energy, num_alphas) deposited_energy = [track.deposited_energy for track in track_list] length = [track.length for track in track_list] projected_length = [track.projected_length() for track in track_list] roi_size = [track.roi_size() for track in track_list] track_size = [track.track_size() for track in track_list] binning = numpy.linspace(0., alpha_energy, 100) hde = xHistogram1d(binning, xlabel='Deposited energy [MeV]').fill(deposited_energy) binning = numpy.linspace(0., 2., 100) htl = xHistogram1d(binning, xlabel='Track length [cm]').fill(length) binning = numpy.linspace(0., 2., 100) hptl = xHistogram1d(binning, xlabel='Projected track length [cm]').fill(projected_length) binning = numpy.linspace(0., 40000, 100) hrs = xHistogram1d(binning, xlabel='ROI size [pixels]').fill(roi_size) binning = numpy.linspace(0., 10000, 100) hts = xHistogram1d(binning, xlabel='Track size [pixels]').fill(track_size) return hde, htl, hptl, hrs, hts
[docs] def plot_distributions(alpha_energy=5.0, num_alphas=1000000): """ """ hde, htl, hptl, hrs, hts = generate_distributions(alpha_energy, num_alphas) plt.figure('Deposited energy') hde.plot() plt.figure('Track length') htl.plot() plt.figure('Projected track length') hptl.plot() plt.figure('ROI size') hrs.plot() plt.figure('Track size') hts.plot()
[docs] def plot_raindrops(alpha_energy=5.0, num_alphas=2000, max_roi_size=5000): """ """ plot_rectangle(2. * ASIC_HALF_SIDE, 2. * ASIC_HALF_SIDE) plt.axis([-1., 1., -1., 1.]) plt.gca().set_aspect('equal') plt.xlabel('x [cm]') plt.ylabel('y [cm]') for track in simulate(alpha_energy, num_alphas): if track.roi_size() > max_roi_size: fmt = dict(lw=1., ls='dashed', color='lightgray') else: fmt = dict(ls='solid', color='gray') plt.plot([track.start.x, track.end.x], [track.start.y, track.end.y], **fmt)
if __name__ == '__main__': plot_distributions() plt.figure('Event raindrops') plot_raindrops() plt.show()