Skip to content
GitLab
Commit 7134f230 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

redate ifs_sim

parent b2028a91
Pipeline #3989 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_ifs_sim/CTI/CTI.py
View file @ 7134f230
...@@ -14,14 +14,6 @@ parameters in parallel and serial direction. ...@@ -14,14 +14,6 @@ parameters in parallel and serial direction.
""" """
import numpy as np import numpy as np
# try:
# import cdm03bidir
# #import cdm03bidirTest as cdm03bidir #for testing purposes only
# except ImportError:
# print('import CTI module')
# #print ('No CDM03bidir module available, please compile it: f2py -c -m cdm03bidir cdm03bidir.f90')
#CDM03bidir #CDM03bidir
class CDM03bidir(): class CDM03bidir():
...@@ -59,7 +51,7 @@ class CDM03bidir(): ...@@ -59,7 +51,7 @@ class CDM03bidir():
self.params.update(self.values) self.params.update(self.values)
#read in trap information #read in trap information
trapdata = np.loadtxt(self.values['parallelTrapfile']) trapdata = np.loadtxt(self.values['dir_path']+self.values['paralleltrapfile'])
if trapdata.ndim > 1: if trapdata.ndim > 1:
self.nt_p = trapdata[:, 0] self.nt_p = trapdata[:, 0]
self.sigma_p = trapdata[:, 1] self.sigma_p = trapdata[:, 1]
...@@ -70,7 +62,7 @@ class CDM03bidir(): ...@@ -70,7 +62,7 @@ class CDM03bidir():
self.sigma_p = [trapdata[1],] self.sigma_p = [trapdata[1],]
self.taur_p = [trapdata[2],] self.taur_p = [trapdata[2],]
trapdata = np.loadtxt(self.values['serialTrapfile']) trapdata = np.loadtxt(self.values['dir_path']+self.values['serialtrapfile'])
if trapdata.ndim > 1: if trapdata.ndim > 1:
self.nt_s = trapdata[:, 0] self.nt_s = trapdata[:, 0]
self.sigma_s = trapdata[:, 1] self.sigma_s = trapdata[:, 1]
...@@ -99,98 +91,6 @@ class CDM03bidir(): ...@@ -99,98 +91,6 @@ class CDM03bidir():
self.logger = False self.logger = False
def radiateFullCCD(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.zeros((xdim, ydim))
#transpose the data, because Python has different convention than Fortran
data = self.data.transpose().copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[0:self.values['xsize'], 0:self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[0:self.values['xsize'], 0:self.values['ysize']] = tmp
elif quad == 1:
d = data[self.values['xsize']:, :self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, :self.values['ysize']] = tmp
elif quad == 2:
d = data[:self.values['xsize'], self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['xsize'], self.values['ysize']:] = tmp
elif quad == 3:
d = data[self.values['xsize']:, self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, self.values['ysize']:] = tmp
else:
print( 'ERROR -- too many quadrants!!' )
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out.transpose()
def radiateFullCCD2(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.empty((ydim, xdim))
#transpose the data, because Python has different convention than Fortran
data = self.data.copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[:self.values['ysize'], :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], :self.values['xsize']] = tmp
elif quad == 1:
d = data[:self.values['ysize'], self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], self.values['xsize']:] = tmp
elif quad == 2:
d = data[self.values['ysize']:, :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, :self.values['xsize']] = tmp
elif quad == 3:
d = data[self.values['ysize']:, self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, self.values['xsize']:] = tmp
else:
print( 'ERROR -- too many quadrants!!')
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out
def applyRadiationDamage(self, data, iquadrant=0): def applyRadiationDamage(self, data, iquadrant=0):
""" """
Apply radian damage based on FORTRAN CDM03 model. The method assumes that Apply radian damage based on FORTRAN CDM03 model. The method assumes that
...@@ -255,7 +155,6 @@ class CDM03bidir(): ...@@ -255,7 +155,6 @@ class CDM03bidir():
################################################################################# #################################################################################
###modify ###modify
import sys
#sys.path.append('../so') #sys.path.append('../so')
from ifs_so import cdm03bidir from ifs_so import cdm03bidir
# from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir # from ifs_so.cdm03.cpython-38-x86_64-linux-gnu import cdm03bidir
...@@ -275,220 +174,6 @@ class CDM03bidir(): ...@@ -275,220 +174,6 @@ class CDM03bidir():
class CDM03():
"""
Class to run CDM03 CTI model, class Fortran routine to perform the actual CDM03 calculations.
:param data: input data to be radiated
:type data: ndarray
:param input: input parameters
:type input: dictionary
:param log: instance to Python logging
:type log: logging instance
"""
def __init__(self, input, data, log=None):
"""
Class constructor.
:param data: input data to be radiated
:type data: ndarray
:param input: input parameters
:type input: dictionary
:param log: instance to Python logging
:type log: logging instance
"""
try:
import cdm03
except ImportError:
print( 'No CDM03 module available, please compile it: f2py -c -m cdm03 cdm03.f90')
self.data = data
self.values = dict(quads=(0,1,2,3), xsize=2048, ysize=2066, dob=0.0, rdose=8.0e9)
self.values.update(input)
self.log = log
self._setupLogger()
def _setupLogger(self):
"""
Set up the logger.
"""
self.logger = True
if self.log is None:
self.logger = False
def radiateFullCCD(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.zeros((xdim, ydim))
#transpose the data, because Python has different convention than Fortran
data = self.data.transpose().copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[0:self.values['xsize'], 0:self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[0:self.values['xsize'], 0:self.values['ysize']] = tmp
elif quad == 1:
d = data[self.values['xsize']:, :self.values['ysize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, :self.values['ysize']] = tmp
elif quad == 2:
d = data[:self.values['xsize'], self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['xsize'], self.values['ysize']:] = tmp
elif quad == 3:
d = data[self.values['xsize']:, self.values['ysize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['xsize']:, self.values['ysize']:] = tmp
else:
print ('ERROR -- too many quadrants!!')
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out.transpose()
def radiateFullCCD2(self):
"""
This routine allows the whole CCD to be run through a radiation damage mode.
The routine takes into account the fact that the amplifiers are in the corners
of the CCD. The routine assumes that the CCD is using four amplifiers.
There is an excess of .copy() calls, which should probably be cleaned up. However,
given that I had problem with the Fortran code, I have kept the calls. If memory
becomes an issue then this should be cleaned.
:return: radiation damaged image
:rtype: ndarray
"""
ydim, xdim = self.data.shape
out = np.empty((ydim, xdim))
#transpose the data, because Python has different convention than Fortran
data = self.data.copy()
for quad in self.values['quads']:
if self.logger:
self.log.info('Adding CTI to Q%i' % quad)
if quad == 0:
d = data[:self.values['ysize'], :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], :self.values['xsize']] = tmp
elif quad == 1:
d = data[:self.values['ysize'], self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[:self.values['ysize'], self.values['xsize']:] = tmp
elif quad == 2:
d = data[self.values['ysize']:, :self.values['xsize']].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, :self.values['xsize']] = tmp
elif quad == 3:
d = data[self.values['ysize']:, self.values['xsize']:].copy()
tmp = self.applyRadiationDamage(d, iquadrant=quad).copy()
out[self.values['ysize']:, self.values['xsize']:] = tmp
else:
print ('ERROR -- too many quadrants!!')
self.log.error('Too many quadrants! This method allows only four quadrants.')
return out
def applyRadiationDamage(self, data, iquadrant=0):
"""
Apply radian damage based on FORTRAN CDM03 model. The method assumes that
input data covers only a single quadrant defined by the iquadrant integer.
:param data: imaging data to which the CDM03 model will be applied to.
:type data: ndarray
:param iquandrant: number of the quadrant to process
:type iquandrant: int
cdm03 - Function signature::
sout = cdm03(sinp,iflip,jflip,dob,rdose,in_nt,in_sigma,in_tr,[xdim,ydim,zdim])
Required arguments:
sinp : input rank-2 array('d') with bounds (xdim,ydim)
iflip : input int
jflip : input int
dob : input float
rdose : input float
in_nt : input rank-1 array('d') with bounds (zdim)
in_sigma : input rank-1 array('d') with bounds (zdim)
in_tr : input rank-1 array('d') with bounds (zdim)
Optional arguments:
xdim := shape(sinp,0) input int
ydim := shape(sinp,1) input int
zdim := len(in_nt) input int
Return objects:
sout : rank-2 array('d') with bounds (xdim,ydim)
.. Note:: Because Python/NumPy arrays are different row/column based, one needs
to be extra careful here. NumPy.asfortranarray will be called to get
an array laid out in Fortran order in memory. Before returning the
array will be laid out in memory in C-style (row-major order).
:return: image that has been run through the CDM03 model
:rtype: ndarray
"""
#read in trap information
trapdata = np.loadtxt(self.values['trapfile'])
nt = trapdata[:, 0]
sigma = trapdata[:, 1]
taur = trapdata[:, 2]
iflip = iquadrant / 2
jflip = iquadrant % 2
if self.logger:
self.log.info('nt=' + str(nt))
self.log.info('sigma= ' + str(sigma))
self.log.info('taur= ' + str(taur))
self.log.info('dob=%f' % self.values['dob'])
self.log.info('rdose=%e' % self.values['rdose'])
self.log.info('xsize=%i' % data.shape[1])
self.log.info('ysize=%i' % data.shape[0])
self.log.info('quadrant=%i' % iquadrant)
self.log.info('iflip=%i' % iflip)
self.log.info('jflip=%i' % jflip)
# #call Fortran routine
# CTIed = cdm03.cdm03(np.asfortranarray(data),
# iflip, jflip,
# self.values['dob'], self.values['rdose'],
# nt, sigma, taur)
###modify
import sys
sys.path.append('../CTI')
import cdm03
#################################################################################
CTIed = cdm03.cdm03(np.asfortranarray(data),
jflip, iflip,
self.values['dob'], self.values['rdose'],
nt,sigma,taur )
return np.asanyarray(CTIed)
################################################################################################################# #################################################################################################################
......
csst_ifs_sim/csst_ifs_sim.py
View file @ 7134f230
...@@ -5,10 +5,6 @@ Created on Thu Apr 11 15:18:57 2024 ...@@ -5,10 +5,6 @@ Created on Thu Apr 11 15:18:57 2024
@author: yan @author: yan
""" """
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# import sys
import os import os
from scipy.interpolate import interp1d from scipy.interpolate import interp1d
import astropy.coordinates as coord import astropy.coordinates as coord
...@@ -18,17 +14,11 @@ import sys ...@@ -18,17 +14,11 @@ import sys
# from csst_ifs_sim.support import cosmicrays # from csst_ifs_sim.support import cosmicrays
# from csst_ifs_sim.support import logger as lg # from csst_ifs_sim.support import logger as lg
# from csst_ifs_sim.CTI import CTI # from csst_ifs_sim.CTI import CTI
sys.path.append('./csst_ifs_sim') sys.path.append('./csst_ifs_sim')
from CTI import CTI from CTI import CTI
from support import logger as lg from support import logger as lg
from support import cosmicrays from support import cosmicrays
from support import IFSinstrumentModel from support import IFSinstrumentModel
# from optparse import OptionParser
import configparser as ConfigParser import configparser as ConfigParser
import cmath import cmath
from scipy import ndimage from scipy import ndimage
...@@ -62,50 +52,31 @@ The approximate sequence of events in the simulator is as follows: ...@@ -62,50 +52,31 @@ The approximate sequence of events in the simulator is as follows:
#. Read in a configuration file, which defines for example, #. Read in a configuration file, which defines for example,
detector characteristics (bias, dark and readout noise, gain, detector characteristics (bias, dark and readout noise, gain,
plate scale and pixel scale, oversampling factor, exposure time etc.). pixel scale,exposure time etc.).
#. Read in another file containing charge trap definitions (for CTI modelling). #. Read in another file containing charge trap definitions (for CTI modelling).
#. Read in a file defining the cosmic rays (trail lengths and #. Read in a file defining the cosmic rays (trail lengths and
cumulative distributions). cumulative distributions).
#. Read in CCD offset information, displace the image, and modify
the output file name to contain the CCD and quadrant information
#. Load the wavefront aberration data used to calculate PSF with defined #. Load the wavefront aberration data used to calculate PSF with defined
wavelength and field of view. wavelength and field of view.
#. Loop over the number of exposures to co-add and for each object in the
object catalog:
* determine the number of electrons an object should have by scaling the
object's magnitude
with the given zeropoint and exposure time.
* determine whether the object lands on to the detector or not and if it is
a star or an extended source (i.e. a galaxy).
* if object is extended determine the size (using a size-magnitude relation)
and scale counts, convolve with the PSF, and finally overlay onto the
detector according to its position.
* if object is a star, scale counts according to the derived scaling
(first step), and finally overlay onto the detector according to its position.
#. Apply calibration unit flux to mimic flat field exposures [optional]. #. Apply calibration unit flux to mimic flat field exposures [optional].
#. Apply a multiplicative flat-field map to emulate pixel-to-pixel #. Apply a multiplicative flat-field map to emulate pixel-to-pixel
non-uniformity [optional]. non-uniformity.
#. Add a charge injection line (horizontal and/or vertical) [optional].
#. Add cosmic ray tracks onto the CCD with random positions but known #. Add cosmic ray tracks onto the CCD with random positions but known
distribution [optional]. distribution.
#. Apply detector charge bleeding in column direction [optional]. #. Apply detector charge bleeding in column direction.
#. Add constant dark current and background light from Zodiacal #. Add constant dark current and background light from Zodiacal
light [optional]. light.
#. Include spatially uniform scattered light to the pixel grid [optional]. #. Include spatially uniform scattered light to the pixel grid.
#. Add photon (Poisson) noise [optional] #. Add photon (Poisson) noise
#. Add cosmetic defects from an input file [optional]. #. Add cosmetic defects from an input file.
#. Add pre- and overscan regions in the serial direction [optional]. #. Add pre- and overscan regions in the serial direction.
#. Apply the CDM03 radiation damage model [optional]. #. Apply the CDM03 radiation damage model [optional].
#. Apply CCD273 non-linearity model to the pixel data [optional]. #. Apply CCD273 non-linearity model to the pixel data.
#. Add readout noise selected from a Gaussian distribution [optional]. #. Add readout noise selected from a Gaussian distribution.
#. Convert from electrons to ADUs using a given gain factor. #. Convert from electrons to ADUs using a given gain factor.
#. Add a given bias level and discretise the counts (the output is going to be #. Add a given bias level and discretise the counts
in 16bit unsigned integers). #. Finally the simulated image is converted to a FITS file.
#. Finally the simulated image is converted to a FITS file, a WCS is assigned
and the output is saved to the current working directory.
Warning:: The code is still work in progress and new features are being added. Warning:: The code is still work in progress and new features are being added.
The code has been tested, but nevertheless bugs may be lurking in corners, so The code has been tested, but nevertheless bugs may be lurking in corners, so
...@@ -120,16 +91,7 @@ Note:: This class is Python 3 compatible. ...@@ -120,16 +91,7 @@ Note:: This class is Python 3 compatible.
2024.3.20 --- update the fits header as defined order 0 data format 2024.3.20 --- update the fits header as defined order 0 data format
""" """
#### functions definition #####
# from astropy.table import Table
###############################################################################
# filterPivotWave = {'nuv': 2875.5, 'u': 3629.6, 'g': 4808.4,
# 'r': 6178.2, 'i': 7609.0, 'z': 9012.9, 'y': 9627.9}
# filterIndex = {'nuv': 0, 'u': 1, 'g': 2, 'r': 3, 'i': 4, 'z': 5, 'y': 6}
def transRaDec2D(ra, dec): def transRaDec2D(ra, dec):
""" """
...@@ -5596,7 +5558,7 @@ class IFSsimulator(): ...@@ -5596,7 +5558,7 @@ class IFSsimulator():
self.applyPoissonNoise() self.applyPoissonNoise()
print('Applying Possion noise finished.......') print('Applying Possion noise finished.......')
################################## #################################
if self.nonlinearity: if self.nonlinearity:
self.applyNonlinearity() self.applyNonlinearity()
......
csst_ifs_sim/support/IFSinstrumentModel.py
View file @ 7134f230
...@@ -51,5 +51,5 @@ def CCDnonLinearityModel(data, beta=6e-7): ...@@ -51,5 +51,5 @@ def CCDnonLinearityModel(data, beta=6e-7):
return out return out
if __name__ == '__main__': # if __name__ == '__main__':
print('IFSinstrument') # print('IFSinstrument')
csst_ifs_sim/support/cosmicrays.py
View file @ 7134f230
...@@ -62,18 +62,18 @@ class cosmicrays(): ...@@ -62,18 +62,18 @@ class cosmicrays():
self._readCosmicrayInformation() self._readCosmicrayInformation()
def _readCosmicrayInformation(self): # def _readCosmicrayInformation(self):
self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], # self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'],
self.information['cosmicraydistance'])) # self.information['cosmicraydistance']))
#read in the information from the files # #read in the information from the files
crLengths = np.loadtxt(self.information['cosmicraylengths']) # crLengths = np.loadtxt(self.information['cosmicraylengths'])
crDists = np.loadtxt(self.information['cosmicraydistance']) # crDists = np.loadtxt(self.information['cosmicraydistance'])
#set up the cosmic ray information dictionary # #set up the cosmic ray information dictionary
self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0], # self.cr = dict(cr_u=crLengths[:, 0], cr_cdf=crLengths[:, 1], cr_cdfn=np.shape(crLengths)[0],
cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) # cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0])
return self.cr # return self.cr
def _cosmicRayIntercepts(self, lum, x0, y0, l, phi): def _cosmicRayIntercepts(self, lum, x0, y0, l, phi):
...@@ -195,113 +195,113 @@ class cosmicrays(): ...@@ -195,113 +195,113 @@ class cosmicrays():
return crImage return crImage
def _drawCosmicRays(self, limit=None): # def _drawCosmicRays(self, limit=None):
""" # """
Add cosmic rays to the arrays based on a power-law intensity distribution for tracks. # Add cosmic rays to the arrays based on a power-law intensity distribution for tracks.
Cosmic ray properties (such as location and angle) are chosen from random Uniform distribution. # Cosmic ray properties (such as location and angle) are chosen from random Uniform distribution.
""" # """
#estimate the number of cosmics # #estimate the number of cosmics
cr_n = self.xsize * self.ysize * 0.014 / 43.263316 * 2. # cr_n = self.xsize * self.ysize * 0.014 / 43.263316 * 2.
#scale with exposure time, the above numbers are for the nominal 565s exposure # #scale with exposure time, the above numbers are for the nominal 565s exposure
cr_n *= (self.information['exptime'] / 565.0) # cr_n *= (self.information['exptime'] / 565.0)
#assume a power-law intensity distribution for tracks # #assume a power-law intensity distribution for tracks
fit = dict(cr_lo=1.0e3, cr_hi=1.0e5, cr_q=2.0e0) # fit = dict(cr_lo=1.0e3, cr_hi=1.0e5, cr_q=2.0e0)
fit['q1'] = 1.0e0 - fit['cr_q'] # fit['q1'] = 1.0e0 - fit['cr_q']
fit['en1'] = fit['cr_lo'] ** fit['q1'] # fit['en1'] = fit['cr_lo'] ** fit['q1']
fit['en2'] = fit['cr_hi'] ** fit['q1'] # fit['en2'] = fit['cr_hi'] ** fit['q1']
#pseudo-random numbers taken from a uniform distribution between 0 and 1 # #pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed() # np.random.seed()
luck = np.random.rand(int(np.floor(cr_n))) # luck = np.random.rand(int(np.floor(cr_n)))
#draw the length of the tracks # #draw the length of the tracks
if self.cr['cr_cdfn'] > 1: # if self.cr['cr_cdfn'] > 1:
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) # ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck) # self.cr['cr_l'] = ius(luck)
else: # else:
self.cr['cr_l'] = np.sqrt(1.0 - luck ** 2) / luck # self.cr['cr_l'] = np.sqrt(1.0 - luck ** 2) / luck
#draw the energy of the tracks # #draw the energy of the tracks
if self.cr['cr_cden'] > 1: # if self.cr['cr_cden'] > 1:
ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) # ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v'])
self.cr['cr_e'] = ius(luck) # self.cr['cr_e'] = ius(luck)
else: # else:
np.random.seed() # np.random.seed()
self.cr['cr_e'] = (fit['en1'] + (fit['en2'] - fit['en1']) * # self.cr['cr_e'] = (fit['en1'] + (fit['en2'] - fit['en1']) *
np.random.rand(int(np.floor(cr_n)))) ** (1.0 / fit['q1']) # np.random.rand(int(np.floor(cr_n)))) ** (1.0 / fit['q1'])
#Choose the properties such as positions and an angle from a random Uniform dist # #Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed() # np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) # cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed() # np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) # cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed() # np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) # cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts # #find the intercepts
if limit is None: # if limit is None:
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) # self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
print ('Number of cosmic ray events:', len(self.cr['cr_e'])) # print ('Number of cosmic ray events:', len(self.cr['cr_e']))
else: # else:
#limit to electron levels < limit # #limit to electron levels < limit
msk = self.cr['cr_e'] < limit # msk = self.cr['cr_e'] < limit
print ('Number of cosmic ray events: %i / %i' % (len(self.cr['cr_e'][msk]), int(np.floor(cr_n)))) # print ('Number of cosmic ray events: %i / %i' % (len(self.cr['cr_e'][msk]), int(np.floor(cr_n))))
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk], # self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk],
self.cr['cr_l'][msk], cr_phi[msk]) # self.cr['cr_l'][msk], cr_phi[msk])
#count the covering factor # #count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap) # area_cr = np.count_nonzero(self.cosmicrayMap)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \ # text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \
% (area_cr, 100.*area_cr / (self.xsize*self.ysize)) # % (area_cr, 100.*area_cr / (self.xsize*self.ysize))
self.log.info(text) # self.log.info(text)
print (text) # print (text)
def _drawSingleEvent(self, limit=1000, cr_n=1): # def _drawSingleEvent(self, limit=1000, cr_n=1):
""" # """
Generate a single cosmic ray event and include it to a cosmic ray map (self.cosmicrayMap). # Generate a single cosmic ray event and include it to a cosmic ray map (self.cosmicrayMap).
:param limit: limiting energy for the cosmic ray event # :param limit: limiting energy for the cosmic ray event
:type limit: float # :type limit: float
:param cr_n: number of cosmic ray events to include # :param cr_n: number of cosmic ray events to include
:type cr_n: int # :type cr_n: int
:return: None # :return: None
""" # """
#pseudo-random numbers taken from a uniform distribution between 0 and 1 # #pseudo-random numbers taken from a uniform distribution between 0 and 1
np.random.seed() # np.random.seed()
luck = np.random.rand(cr_n) # luck = np.random.rand(cr_n)
#draw the length of the tracks # #draw the length of the tracks
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) # ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck) # self.cr['cr_l'] = ius(luck)
#set the energy directly to the limit # #set the energy directly to the limit
self.cr['cr_e'] = np.asarray([limit, ]*cr_n) # self.cr['cr_e'] = np.asarray([limit, ]*cr_n)
#Choose the properties such as positions and an angle from a random Uniform dist # #Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed() # np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) # cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed() # np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) # cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
np.random.seed() # np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) # cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#find the intercepts # #find the intercepts
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) # self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
#count the covering factor # #count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap) # area_cr = np.count_nonzero(self.cosmicrayMap)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \ # text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \
% (area_cr, 100.*area_cr / (self.xsize*self.ysize)) # % (area_cr, 100.*area_cr / (self.xsize*self.ysize))
self.log.info(text) # self.log.info(text)
print( text) # print( text)
def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False): def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False):
...@@ -361,38 +361,38 @@ class cosmicrays(): ...@@ -361,38 +361,38 @@ class cosmicrays():
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering) text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering)
self.log.info(text) self.log.info(text)
if verbose: # if verbose:
print( text) # print( text)
def addCosmicRays(self, limit=None): # def addCosmicRays(self, limit=None):
""" # """
Include cosmic rays to the image given. # Include cosmic rays to the image given.
:return: image with cosmic rays # :return: image with cosmic rays
:rtype: ndarray # :rtype: ndarray
""" # """
self._drawCosmicRays(limit=limit) # self._drawCosmicRays(limit=limit)
#paste cosmic rays # #paste cosmic rays
self.image += self.cosmicrayMap # self.image += self.cosmicrayMap
return self.image # return self.image
def addSingleEvent(self, limit=None): # def addSingleEvent(self, limit=None):
""" # """
Include a single cosmic ray event to the image given. # Include a single cosmic ray event to the image given.
:return: image with cosmic rays # :return: image with cosmic rays
:rtype: ndarray # :rtype: ndarray
""" # """
self._drawSingleEvent(limit=limit) # self._drawSingleEvent(limit=limit)
#paste cosmic rays # #paste cosmic rays
self.image += self.cosmicrayMap # self.image += self.cosmicrayMap
return self.image # return self.image
def addUpToFraction(self, coveringFraction, limit=None, verbose=False): def addUpToFraction(self, coveringFraction, limit=None, verbose=False):
...@@ -417,8 +417,8 @@ class cosmicrays(): ...@@ -417,8 +417,8 @@ class cosmicrays():
return self.image return self.image
if __name__ == "__main__": # if __name__ == "__main__":
print() # print()
ifs_data/IFS_sim_C90.config
View file @ 7134f230
...@@ -124,7 +124,7 @@ sky_noise=yes ...@@ -124,7 +124,7 @@ sky_noise=yes
cosmetics = yes cosmetics = yes
#apply radiation damage model? #apply radiation damage model?
radiationDamage = no radiationDamage = yes
#add cosmic rays? #add cosmic rays?
cosmicRays = yes cosmicRays = yes
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment