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

test

parent bfd10df5
Pipeline #4067 failed with stage
in 0 seconds
Hide whitespace changes
Inline Side-by-side
csst_mci_sim/csst_mci_sim.py
View file @ e6ed2c85
...@@ -1401,8 +1401,8 @@ class MCIsimulator(): ...@@ -1401,8 +1401,8 @@ class MCIsimulator():
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']))
crLengths = np.loadtxt(self.information['cosmicraylengths']) crLengths = np.loadtxt(self.information['dir_path']+self.information['cosmicraylengths'])
crDists = np.loadtxt(self.information['cosmicraydistance']) crDists = np.loadtxt(self.information['dir_path']+self.information['cosmicraydistance'])
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])
......
csst_mci_sim/support/cosmicrays.py
View file @ e6ed2c85
...@@ -50,9 +50,9 @@ class cosmicrays(): ...@@ -50,9 +50,9 @@ class cosmicrays():
self.ysize, self.xsize = self.image.shape self.ysize, self.xsize = self.image.shape
#set up the information dictionary, first with defaults and then overwrite with inputs if given #set up the information dictionary, first with defaults and then overwrite with inputs if given
self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat', # self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat',
cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat', # cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat',
exptime=565)) # exptime=565))
if information is not None: if information is not None:
self.information.update(information) self.information.update(information)
...@@ -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):
...@@ -104,23 +104,23 @@ class cosmicrays(): ...@@ -104,23 +104,23 @@ class cosmicrays():
ilo = np.round(x0.copy() - dx) ilo = np.round(x0.copy() - dx)
msk = ilo < 0. msk = ilo < 0.
ilo[msk] = 0 ilo[msk] = 0
ilo = ilo.astype(np.int) ilo = ilo.astype(int)
ihi = 1 + np.round(x0.copy() + dx) ihi = 1 + np.round(x0.copy() + dx)
msk = ihi > self.xsize msk = ihi > self.xsize
ihi[msk] = self.xsize ihi[msk] = self.xsize
ihi = ihi.astype(np.int) ihi = ihi.astype(int)
#pixels in y-directions #pixels in y-directions
jlo = np.round(y0.copy() - dy) jlo = np.round(y0.copy() - dy)
msk = jlo < 0. msk = jlo < 0.
jlo[msk] = 0 jlo[msk] = 0
jlo = jlo.astype(np.int) jlo = jlo.astype(int)
jhi = 1 + np.round(y0.copy() + dy) jhi = 1 + np.round(y0.copy() + dy)
msk = jhi > self.ysize msk = jhi > self.ysize
jhi[msk] = self.ysize jhi[msk] = self.ysize
jhi = jhi.astype(np.int) jhi = jhi.astype(int)
#loop over the individual events #loop over the individual events
for i, luminosity in enumerate(lum): for i, luminosity in enumerate(lum):
...@@ -195,109 +195,113 @@ class cosmicrays(): ...@@ -195,109 +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()
cr_y = self.ysize * 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_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
# np.random.seed()
#find the intercepts # cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
if limit is None:
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) # #find the intercepts
print ('Number of cosmic ray events:', len(self.cr['cr_e'])) # if limit is None:
else: # self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
#limit to electron levels < limit # print ('Number of cosmic ray events:', len(self.cr['cr_e']))
msk = self.cr['cr_e'] < limit # else:
print ('Number of cosmic ray events: %i / %i' % (len(self.cr['cr_e'][msk]), int(np.floor(cr_n)))) # #limit to electron levels < limit
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk], # msk = self.cr['cr_e'] < limit
self.cr['cr_l'][msk], cr_phi[msk]) # 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],
#count the covering factor # self.cr['cr_l'][msk], cr_phi[msk])
area_cr = np.count_nonzero(self.cosmicrayMap)
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \ # #count the covering factor
% (area_cr, 100.*area_cr / (self.xsize*self.ysize)) # area_cr = np.count_nonzero(self.cosmicrayMap)
self.log.info(text) # text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \
print (text) # % (area_cr, 100.*area_cr / (self.xsize*self.ysize))
# self.log.info(text)
# print (text)
def _drawSingleEvent(self, limit=1000, cr_n=1):
"""
Generate a single cosmic ray event and include it to a cosmic ray map (self.cosmicrayMap). # def _drawSingleEvent(self, limit=1000, cr_n=1):
# """
:param limit: limiting energy for the cosmic ray event # Generate a single cosmic ray event and include it to a cosmic ray map (self.cosmicrayMap).
:type limit: float
:param cr_n: number of cosmic ray events to include # :param limit: limiting energy for the cosmic ray event
:type cr_n: int # :type limit: float
# :param cr_n: number of cosmic ray events to include
:return: None # :type cr_n: int
"""
#pseudo-random numbers taken from a uniform distribution between 0 and 1 # :return: None
np.random.seed() # """
luck = np.random.rand(cr_n) # #pseudo-random numbers taken from a uniform distribution between 0 and 1
# np.random.seed()
#draw the length of the tracks # luck = np.random.rand(cr_n)
ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
self.cr['cr_l'] = ius(luck) # #draw the length of the tracks
# ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u'])
#set the energy directly to the limit # self.cr['cr_l'] = ius(luck)
self.cr['cr_e'] = np.asarray([limit, ]*cr_n)
# #set the energy directly to the limit
#Choose the properties such as positions and an angle from a random Uniform dist # self.cr['cr_e'] = np.asarray([limit, ]*cr_n)
np.random.seed()
cr_x = self.xsize * np.random.rand(int(np.floor(cr_n))) # #Choose the properties such as positions and an angle from a random Uniform dist
np.random.seed() # np.random.seed()
cr_y = self.ysize * np.random.rand(int(np.floor(cr_n))) # cr_x = self.xsize * np.random.rand(int(np.floor(cr_n)))
np.random.seed()
cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) # np.random.seed()
# cr_y = self.ysize * np.random.rand(int(np.floor(cr_n)))
#find the intercepts
self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi) # np.random.seed()
# cr_phi = np.pi * np.random.rand(int(np.floor(cr_n)))
#count the covering factor
area_cr = np.count_nonzero(self.cosmicrayMap) # #find the intercepts
text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \ # self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'], cr_x, cr_y, self.cr['cr_l'], cr_phi)
% (area_cr, 100.*area_cr / (self.xsize*self.ysize))
self.log.info(text) # #count the covering factor
print( text) # area_cr = np.count_nonzero(self.cosmicrayMap)
# text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' \
# % (area_cr, 100.*area_cr / (self.xsize*self.ysize))
# self.log.info(text)
# print( text)
def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False): def _drawEventsToCoveringFactor(self, coveringFraction=3.0, limit=1000, verbose=False):
...@@ -340,8 +344,10 @@ class cosmicrays(): ...@@ -340,8 +344,10 @@ class cosmicrays():
#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)))
...@@ -355,38 +361,38 @@ class cosmicrays(): ...@@ -355,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):
...@@ -411,38 +417,8 @@ class cosmicrays(): ...@@ -411,38 +417,8 @@ class cosmicrays():
return self.image return self.image
if __name__ == "__main__": # if __name__ == "__main__":
import sys
sys.path.append('/home/yan/csst-master/')
from support import logger as lg
from support import files as fileIO # print()
from scipy import ndimage
from astropy.io import fits
#set up logger
log = lg.setUpLogger('VISsim.log')
#test section
crImage = np.zeros((2066, 2048), dtype=np.float64)
#cosmic ray instance
cosmics = cosmicrays(log, crImage)
#add cosmic rays up to the covering fraction
CCD_cr = cosmics.addUpToFraction(1.4, limit=None, verbose=True)
print(CCD_cr)
print(type(CCD_cr))
effected = np.count_nonzero(CCD_cr)
print (effected, effected*100./(CCD_cr.shape[0]*CCD_cr.shape[1]))
#save to FITS
fits.writeto(CCD_cr, 'cosmicrayTest.fits', overwrite=True)
#smooth with a charge diffusion kernel
smooth = ndimage.filters.gaussian_filter(CCD_cr, (0.32, 0.32))
fits.writeto(smooth, 'cosmicrayTestSmoothed.fits')
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