The negative values in un-rectified SED of SALT2 (and SALT2-extended) at

examples/testSN.py

@@ -9,6 +9,10 @@
         implementations of extinction using OD94 model.)
     - Band Flux for extincted SED in r Band
     - Band Mag for extincted SED in r Band
+    - rectification of SED: while the SEDs from the model can go negative
+      resulting in negative model fluxes and adus, this has to be prevented in
+      SNObject, particularly for image generation software. Test that if
+      rectifySEDs is set to True, such situations are avoided
 
 SNIaCatalog_tests:
 A Class containing tests to check crictical functionality for SNIaCatalog 
@@ -24,6 +28,7 @@
 from lsst.sims.photUtils import Bandpass
 from lsst.sims.photUtils import BandpassDict
 from lsst.sims.utils import ObservationMetaData
+from lsst.sims.photUtils.PhotometricParameters import PhotometricParameters
 from lsst.sims.utils import spatiallySample_obsmetadata as sample_obsmetadata
 from lsst.sims.catUtils.utils import ObservationMetaDataGenerator
 from lsst.sims.catalogs.generation.db import CatalogDBObject
@@ -78,6 +83,7 @@ def setUp(self):
         self.SNCosmoModel = self.SN_extincted.equivalentSNCosmoModel()
 
         self.lsstBandPass = BandpassDict.loadTotalBandpassesFromFiles()
+        self.photParams = PhotometricParameters()
         self.SNCosmoBP = sncosmo.Bandpass(wave=self.lsstBandPass['r'].wavelen,
                                           trans=self.lsstBandPass['r'].sb,
                                           wave_unit=astropy.units.Unit('nm'),
@@ -95,6 +101,53 @@ def test_SNstatenotEmpty(self):
         for key in myDict.keys():
             assert myDict[key] is not None
 
+    def test_attributeDefaults(self):
+        """
+        Check the defaults and the setter properties for rectifySED and modelOutSideRange
+        """
+        snobj = SNObject(ra=30., dec=-60., source='salt2')
+        assert snobj.rectifySED == True
+        assert snobj.modelOutSideRange == 'zero'
+
+        snobj.rectifySED = False
+        assert snobj.rectifySED == False
+        assert snobj.modelOutSideRange == 'zero'
+        # Should run but do nothing
+        snobj.modelOutSideRange = 'False'
+        assert snobj.modelOutSideRange == 'zero'
+
+
+
+
+
+    def test_rectifiedSED(self):
+        """
+        Check for an extreme case that the SN seds are being rectified. This is
+        done by setting up an extreme case where there will be negative seds, and
+        checking that this is indeed the case, and checking that they are not
+        negative if rectified.
+        """
+
+        snobj = SNObject(ra=30., dec=-60., source='salt2')
+        snobj.set(z=0.96, t0=self.mjdobs, x1=-3., x0=1.8e-6)
+        snobj.rectifySED  = False
+        times = np.arange(self.mjdobs - 50., self.mjdobs + 150., 1.)
+        badTimes = []
+        for time in times:
+            sed = snobj.SNObjectSED(time=time,
+                                    bandpass=self.lsstBandPass['r'])
+            if any(sed.flambda < 0.):
+                badTimes.append(time)
+        # Check that there are negative SEDs
+        assert(len(badTimes) > 0)
+        snobj.rectifySED = True
+        for time in badTimes:
+            sed = snobj.SNObjectSED(time=time,
+                                    bandpass=self.lsstBandPass['r'])
+            assert not sed.calcADU(bandpass=self.lsstBandPass['r'],
+                                   photParams=self.photParams) < 0.
+            assert not any(sed.flambda < 0.)
+
     def test_ComparebandFluxes2photUtils(self):
         """
         The SNObject.catsimBandFlux computation uses the sims.photUtils.sed
@@ -386,9 +439,9 @@ def test_redrawingCatalog(self):
 
         for key in oldlcs.groups.keys():
             df_old = oldlcs.get_group(key)
-            df_old.sort(['time', 'band'], inplace=True)
+            df_old.sort_values(['time', 'band'], inplace=True)
             df_new = newlcs.get_group(key)
-            df_new.sort(['time', 'band'], inplace=True)
+            df_new.sort_values(['time', 'band'], inplace=True)
             s = "Testing equality for SNID {0:8d} with {1:2d} datapoints" 
             print(s.format(df_new.snid.iloc[0], len(df_old)))
             assert_frame_equal(df_new, df_old)

python/lsst/sims/catUtils/exampleCatalogDefinitions/phoSimCatalogExamples.py

@@ -115,7 +115,7 @@ class PhoSimCatalogPoint(PhosimInputBase, PhoSimAstrometryStars, EBVmixin):
 
     default_columns = [('redshift', 0., float),('shear1', 0., float), ('shear2', 0., float),
                        ('kappa', 0., float), ('raOffset', 0., float), ('decOffset', 0., float),
-                       ('galacticExtinctionModel', 'CCM', (str,3)),
+                       ('galacticExtinctionModel', 'CCM', (str,3)), ('galacticRv', 3.1, float),
                        ('internalExtinctionModel', 'none', (str,4))]
 
     default_formats = {'S':'%s', 'f':'%.9g', 'i':'%i'}
@@ -137,7 +137,7 @@ class PhoSimCatalogZPoint(PhosimInputBase, PhoSimAstrometryGalaxies, EBVmixin):
     default_columns = [('shear1', 0., float), ('shear2', 0., float), ('kappa', 0., float),
                        ('raOffset', 0., float), ('decOffset', 0., float), ('spatialmodel', 'ZPOINT', (str, 6)),
                        ('galacticExtinctionModel', 'CCM', (str,3)),
-                       ('galacticAv', 0.1, float),
+                       ('galacticAv', 0.1, float), ('galacticRv', 3.1, float),
                        ('internalExtinctionModel', 'none', (str,4))]
 
     default_formats = {'S':'%s', 'f':'%.9g', 'i':'%i'}
@@ -159,7 +159,8 @@ class PhoSimCatalogSersic2D(PhoSimCatalogZPoint):
                       'internalExtinctionModel','internalAv','internalRv']
 
     default_columns = [('shear1', 0., float), ('shear2', 0., float), ('kappa', 0., float),
-                       ('raOffset', 0., float), ('decOffset', 0., float), ('galacticAv', 0.1, float),
+                       ('raOffset', 0., float), ('decOffset', 0., float),
+                       ('galacticAv', 0.1, float), ('galacticRv', 3.1, float),
                        ('galacticExtinctionModel', 'CCM', (str,3)),
                        ('internalExtinctionModel', 'CCM', (str,3)), ('internalAv', 0., float),
                        ('internalRv', 3.1, float) ]

python/lsst/sims/catUtils/mixins/EBVmixin.py

@@ -14,27 +14,25 @@ class EBVmixin(EBVbase):
     This mixin class contains the getters which a catalog object will use to call
     calculateEbv in the EBVbase class
     """
-
-
+
     #and finally, here is the getter
     @cached
     def get_EBV(self):
         """
         Getter for the InstanceCatalog framework
         """
-        
+
         galacticCoordinates=numpy.array([self.column_by_name('glon'),self.column_by_name('glat')])
-  
+
         EBV_out=numpy.array(self.calculateEbv(galacticCoordinates=galacticCoordinates,interp=True))
         return EBV_out
-    
-    @cached    
-    def get_galacticRv(self):
+
+    @cached
+    def get_galacticAv(self):
         """
-        Returns galactic RV by getting galacticAv and EBV and assuming Rv = Av/EBV"
+        Getter to return galactic(Milky Way) Av based on EBV values
+        from getter (reading dustmaps) and the RV value from galacticRv
         """
-        Av=self.column_by_name('galacticAv')
-        EBV=self.column_by_name('EBV')
-        return numpy.array(Av/EBV)
-
+        return self.column_by_name('EBV')*self.column_by_name('galacticRv')
+
 

python/lsst/sims/catUtils/supernovae/snObject.py

@@ -58,6 +58,7 @@ class SNObject(sncosmo.Model):
         Therefore, the value must be set explicitly to 0. to get unextincted
         quantities.
 
+    rectifySED : Bool, True by Default
     Methods
     -------
 
@@ -121,8 +122,32 @@ def __init__(self, ra=None, dec=None, source='salt2-extended'):
         self.ebvofMW = None
         if self._hascoords:
             self.mwEBVfromMaps()
+
+
+        # Behavior of model outside range :
+        # if 'zero' then all fluxes outside the temporal range of the model
+        # are set to 0.
+        self._modelOutSideRange = 'zero'
+
+        # SED will be rectified to 0. for negative values of SED if this
+        # attribute is set to True
+        self.rectifySED = True
         return
 
+    @property
+    def modelOutSideRange(self):
+        """
+        Defines the behavior of the model when sampled at times beyond the
+        model definition.
+        """
+        return self._modelOutSideRange
+
+    @modelOutSideRange.setter
+    def modelOutSideRange(self, value):
+        if value != 'zero':
+            print ('Model not implemented, defaulting to zero method\n')
+        return self._modelOutSideRange
+
     @property
     def SNstate(self):
         """
@@ -210,6 +235,7 @@ def equivalentSNCosmoModel(self):
         sncosmoParams['mwebv'] = snState['MWE(B-V)']
         sncosmoModel.set(**sncosmoParams)
         return sncosmoModel
+
     @staticmethod
     def equivsncosmoParamDict(SNstate, SNCosmoModel):
         """
@@ -261,7 +287,7 @@ def sncosmoParamDict(SNstate, SNCosmoModel):
 
 
     def summary(self):
-        '''
+        """
         summarizes the current state of the SNObject class in a returned
         string.
 
@@ -277,7 +303,8 @@ def summary(self):
         --------
         >>> t = SNObject()
         >>> print t.summary()
-        '''
+        """
+
         state = '  SNObject Summary      \n'
 
         state += 'Model = ' + '\n'
@@ -390,7 +417,7 @@ def mwEBVfromMaps(self):
 
     def SNObjectSED(self, time, wavelen=None, bandpass=None,
                     applyExtinction=True):
-        '''
+        """
         return a `lsst.sims.photUtils.sed` object from the SN model at the
         requested time and wavelengths with or without extinction from MW
         according to the SED extinction methods. The wavelengths may be
@@ -431,7 +458,7 @@ def SNObjectSED(self, time, wavelen=None, bandpass=None,
         Examples
         --------
         >>> sed = SN.SNObjectSED(time=0., wavelen=wavenm)
-        '''
+        """
 
         if wavelen is None and bandpass is None:
             raise ValueError('A non None input to either wavelen or\
@@ -476,6 +503,13 @@ def SNObjectSED(self, time, wavelen=None, bandpass=None,
             flambda[mask] = self.flux(time=time, wave=wave)
             flambda[mask] = flambda[mask] * 10.0
 
+        if not (self.modelOutSideRange=='zero'):
+            raise ValueError('Other values of modelOutSideRange not\
+                              implemented')
+        # Rectify
+        if self.rectifySED:
+            flambda = np.where(flambda > 0., flambda, 0.)
+
         SEDfromSNcosmo = Sed(wavelen=wavelen, flambda=flambda)
 
         if not applyExtinction:

python/lsst/sims/catUtils/utils/CatalogTestUtils.py

@@ -402,7 +402,7 @@ class cartoonStars(InstanceCatalog,AstrometryStars,EBVmixin,VariabilityStars,car
     #I need to give it the name of an actual SED file that spans the expected wavelength range
     defSedName = 'km30_5250.fits_g00_5370'
     default_columns = [('sedFilename', defSedName, (str,len(defSedName))), ('glon', 180., float),
-                       ('glat', 30., float), ('galacticAv', 0.1, float)]
+                       ('glat', 30., float), ('galacticAv', 0.1, float), ('galacticRv', 3.1, float)]
 
 class cartoonStarsOnlyI(InstanceCatalog, AstrometryStars ,EBVmixin, VariabilityStars, PhotometryStars):
     catalog_type = 'cartoonStarsOnlyI'
@@ -411,7 +411,7 @@ class cartoonStarsOnlyI(InstanceCatalog, AstrometryStars ,EBVmixin, VariabilityS
     #I need to give it the name of an actual SED file that spans the expected wavelength range
     defSedName = 'km30_5250.fits_g00_5370'
     default_columns = [('sedFilename', defSedName, (str,len(defSedName))), ('glon', 180., float),
-                       ('glat', 30., float), ('galacticAv', 0.1, float)]
+                       ('glat', 30., float), ('galacticAv', 0.1, float), ('galacticRv', 3.1, float)]
 
     @compound('cartoon_u','cartoon_g','cartoon_r','cartoon_i','cartoon_z')
     def get_magnitudes(self):
@@ -563,7 +563,7 @@ class testStars(InstanceCatalog, EBVmixin, VariabilityStars, TestVariabilityMixi
     'EBV','varParamStr']
     defSedName = 'sed_flat.txt'
     default_columns = [('sedFilename', defSedName, (str,len(defSedName))), ('glon', 180., float),
-                       ('glat', 30., float), ('galacticAv', 0.1, float)]
+                       ('glat', 30., float), ('galacticAv', 0.1, float), ('galacticRv', 3.1, float)]
 
 class testGalaxies(InstanceCatalog,EBVmixin,VariabilityGalaxies,TestVariabilityMixin,PhotometryGalaxies,testDefaults):
     """

tests/testPhotometryMixins.py

@@ -368,7 +368,7 @@ def testAlternateBandpassesStars(self):
         self.assertEqual(len(mags), len(test_cat.cartoonBandpassDict))
         self.assertGreater(len(mags), 0)
         for j in range(len(mags)):
-            self.assertAlmostEqual(mags[j],test_cat.magnitudeMasterList[i][j],5)
+            self.assertAlmostEqual(mags[j],test_cat.magnitudeMasterList[i][j],4)
 
         with open(catName, 'r') as input_file:
             lines = input_file.readlines()