diff --git a/csst_mci_sim/csst_mci_sim.py b/csst_mci_sim/csst_mci_sim.py index 39d26d88a54d1971cccb3fa5cfd81b4aff5e4476..e3a3c03fc55e3b5aa7e7a5cdd4c6e6ace3e8dce8 100644 --- a/csst_mci_sim/csst_mci_sim.py +++ b/csst_mci_sim/csst_mci_sim.py @@ -1401,8 +1401,8 @@ class MCIsimulator(): self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], self.information['cosmicraydistance'])) - crLengths = np.loadtxt(self.information['cosmicraylengths']) - crDists = np.loadtxt(self.information['cosmicraydistance']) + crLengths = np.loadtxt(self.information['dir_path']+self.information['cosmicraylengths']) + 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], cr_v=crDists[:, 0], cr_cde=crDists[:, 1], cr_cden=np.shape(crDists)[0]) diff --git a/csst_mci_sim/support/cosmicrays.py b/csst_mci_sim/support/cosmicrays.py index 3b8e4bafe679feaf540d6dd8afd24ce67797af9f..2a36ec4c64c2c72fd5f138f3cd5e9e6045c7d9c0 100644 --- a/csst_mci_sim/support/cosmicrays.py +++ b/csst_mci_sim/support/cosmicrays.py @@ -50,9 +50,9 @@ class cosmicrays(): self.ysize, self.xsize = self.image.shape #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', - cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat', - exptime=565)) + # self.information = (dict(cosmicraylengths='/home/yan/csst-master/data/cdf_cr_length.dat', + # cosmicraydistance='/home/yan/csst-master/data/cdf_cr_total.dat', + # exptime=565)) if information is not None: self.information.update(information) @@ -62,18 +62,18 @@ class cosmicrays(): self._readCosmicrayInformation() - def _readCosmicrayInformation(self): - self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], - self.information['cosmicraydistance'])) - #read in the information from the files - crLengths = np.loadtxt(self.information['cosmicraylengths']) - crDists = np.loadtxt(self.information['cosmicraydistance']) + # def _readCosmicrayInformation(self): + # self.log.info('Reading in cosmic ray information from %s and %s' % (self.information['cosmicraylengths'], + # self.information['cosmicraydistance'])) + # #read in the information from the files + # crLengths = np.loadtxt(self.information['cosmicraylengths']) + # crDists = np.loadtxt(self.information['cosmicraydistance']) - #set up the cosmic ray information dictionary - 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]) + # #set up the cosmic ray information dictionary + # 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]) - return self.cr + # return self.cr def _cosmicRayIntercepts(self, lum, x0, y0, l, phi): @@ -104,23 +104,23 @@ class cosmicrays(): ilo = np.round(x0.copy() - dx) msk = ilo < 0. ilo[msk] = 0 - ilo = ilo.astype(np.int) + ilo = ilo.astype(int) ihi = 1 + np.round(x0.copy() + dx) msk = ihi > self.xsize ihi[msk] = self.xsize - ihi = ihi.astype(np.int) + ihi = ihi.astype(int) #pixels in y-directions jlo = np.round(y0.copy() - dy) msk = jlo < 0. jlo[msk] = 0 - jlo = jlo.astype(np.int) + jlo = jlo.astype(int) jhi = 1 + np.round(y0.copy() + dy) msk = jhi > self.ysize jhi[msk] = self.ysize - jhi = jhi.astype(np.int) + jhi = jhi.astype(int) #loop over the individual events for i, luminosity in enumerate(lum): @@ -195,109 +195,113 @@ class cosmicrays(): return crImage - def _drawCosmicRays(self, limit=None): - """ - 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. - """ - #estimate the number of cosmics - cr_n = self.xsize * self.ysize * 0.014 / 43.263316 * 2. - #scale with exposure time, the above numbers are for the nominal 565s exposure - cr_n *= (self.information['exptime'] / 565.0) - - #assume a power-law intensity distribution for tracks - fit = dict(cr_lo=1.0e3, cr_hi=1.0e5, cr_q=2.0e0) - fit['q1'] = 1.0e0 - fit['cr_q'] - fit['en1'] = fit['cr_lo'] ** fit['q1'] - fit['en2'] = fit['cr_hi'] ** fit['q1'] - - #pseudo-random numbers taken from a uniform distribution between 0 and 1 - np.random.seed() - luck = np.random.rand(int(np.floor(cr_n))) - - #draw the length of the tracks - if self.cr['cr_cdfn'] > 1: - ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) - self.cr['cr_l'] = ius(luck) - else: - self.cr['cr_l'] = np.sqrt(1.0 - luck ** 2) / luck - - #draw the energy of the tracks - if self.cr['cr_cden'] > 1: - ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) - self.cr['cr_e'] = ius(luck) - else: - np.random.seed() - self.cr['cr_e'] = (fit['en1'] + (fit['en2'] - fit['en1']) * - 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 - np.random.seed() - 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() - cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) - - #find the intercepts - if limit is None: - 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'])) - else: - #limit to electron levels < 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)))) - self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk], - self.cr['cr_l'][msk], cr_phi[msk]) - - #count the covering factor - 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 _drawSingleEvent(self, limit=1000, cr_n=1): - """ - 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 - :type limit: float - :param cr_n: number of cosmic ray events to include - :type cr_n: int - - :return: None - """ - #pseudo-random numbers taken from a uniform distribution between 0 and 1 - np.random.seed() - luck = np.random.rand(cr_n) - - #draw the length of the tracks - ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) - self.cr['cr_l'] = ius(luck) - - #set the energy directly to the limit - self.cr['cr_e'] = np.asarray([limit, ]*cr_n) - - #Choose the properties such as positions and an angle from a random Uniform dist - np.random.seed() - 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() - cr_phi = np.pi * 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) - - #count the covering factor - 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 _drawCosmicRays(self, limit=None): + # """ + # 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. + # """ + # #estimate the number of cosmics + # cr_n = self.xsize * self.ysize * 0.014 / 43.263316 * 2. + # #scale with exposure time, the above numbers are for the nominal 565s exposure + # cr_n *= (self.information['exptime'] / 565.0) + + # #assume a power-law intensity distribution for tracks + # fit = dict(cr_lo=1.0e3, cr_hi=1.0e5, cr_q=2.0e0) + # fit['q1'] = 1.0e0 - fit['cr_q'] + # fit['en1'] = fit['cr_lo'] ** fit['q1'] + # fit['en2'] = fit['cr_hi'] ** fit['q1'] + + # #pseudo-random numbers taken from a uniform distribution between 0 and 1 + # np.random.seed() + # luck = np.random.rand(int(np.floor(cr_n))) + + # #draw the length of the tracks + # if self.cr['cr_cdfn'] > 1: + # ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) + # self.cr['cr_l'] = ius(luck) + # else: + # self.cr['cr_l'] = np.sqrt(1.0 - luck ** 2) / luck + + # #draw the energy of the tracks + # if self.cr['cr_cden'] > 1: + # ius = InterpolatedUnivariateSpline(self.cr['cr_cde'], self.cr['cr_v']) + # self.cr['cr_e'] = ius(luck) + # else: + # np.random.seed() + # self.cr['cr_e'] = (fit['en1'] + (fit['en2'] - fit['en1']) * + # 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 + # np.random.seed() + # 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() + # cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) + + # #find the intercepts + # if limit is None: + # 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'])) + # else: + # #limit to electron levels < 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)))) + # self.cosmicrayMap = self._cosmicRayIntercepts(self.cr['cr_e'][msk], cr_x[msk], cr_y[msk], + # self.cr['cr_l'][msk], cr_phi[msk]) + + # #count the covering factor + # 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 _drawSingleEvent(self, limit=1000, cr_n=1): + # """ + # 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 + # :type limit: float + # :param cr_n: number of cosmic ray events to include + # :type cr_n: int + + # :return: None + # """ + # #pseudo-random numbers taken from a uniform distribution between 0 and 1 + # np.random.seed() + # luck = np.random.rand(cr_n) + + # #draw the length of the tracks + # ius = InterpolatedUnivariateSpline(self.cr['cr_cdf'], self.cr['cr_u']) + # self.cr['cr_l'] = ius(luck) + + # #set the energy directly to the limit + # self.cr['cr_e'] = np.asarray([limit, ]*cr_n) + + # #Choose the properties such as positions and an angle from a random Uniform dist + # np.random.seed() + # 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() + # cr_phi = np.pi * 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) + + # #count the covering factor + # 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): @@ -340,8 +344,10 @@ class cosmicrays(): #Choose the properties such as positions and an angle from a random Uniform dist np.random.seed() 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() cr_phi = np.pi * np.random.rand(int(np.floor(cr_n))) @@ -355,38 +361,38 @@ class cosmicrays(): text = 'The cosmic ray covering factor is %i pixels i.e. %.3f per cent' % (area_cr, covering) self.log.info(text) - if verbose: - print( text) + # if verbose: + # print( text) - def addCosmicRays(self, limit=None): - """ - Include cosmic rays to the image given. + # def addCosmicRays(self, limit=None): + # """ + # Include cosmic rays to the image given. - :return: image with cosmic rays - :rtype: ndarray - """ - self._drawCosmicRays(limit=limit) + # :return: image with cosmic rays + # :rtype: ndarray + # """ + # self._drawCosmicRays(limit=limit) - #paste cosmic rays - self.image += self.cosmicrayMap + # #paste cosmic rays + # self.image += self.cosmicrayMap - return self.image + # return self.image - def addSingleEvent(self, limit=None): - """ - Include a single cosmic ray event to the image given. + # def addSingleEvent(self, limit=None): + # """ + # Include a single cosmic ray event to the image given. - :return: image with cosmic rays - :rtype: ndarray - """ - self._drawSingleEvent(limit=limit) + # :return: image with cosmic rays + # :rtype: ndarray + # """ + # self._drawSingleEvent(limit=limit) - #paste cosmic rays - self.image += self.cosmicrayMap + # #paste cosmic rays + # self.image += self.cosmicrayMap - return self.image + # return self.image def addUpToFraction(self, coveringFraction, limit=None, verbose=False): @@ -411,38 +417,8 @@ class cosmicrays(): return self.image -if __name__ == "__main__": - import sys - sys.path.append('/home/yan/csst-master/') - +# if __name__ == "__main__": - from support import logger as lg - from support import files as fileIO - 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)) + +# print() - fits.writeto(smooth, 'cosmicrayTestSmoothed.fits')