Point source sensitivity calculator
This notebook shows how to compute a point source sensitivity for a specific model. In the current status of cosipy, we are limited by the available responses, i.e at 511, 1157, 1173, 1333 and 1809 keV and a coarse continuum. Keep in mind that due to the coarseness of the current responses, the sensitivity might be underestimated.
The collaboration is currently developping a neural network approach for the response to use with an unbinned analysis. This will allow us to fully exploit COSI power.
[11]:
import os
import logging
import sys
logger = logging.getLogger('cosipy')
logger.setLevel(logging.INFO)
logger.addHandler(logging.StreamHandler(sys.stdout))
os.environ["OMP_NUM_THREADS"] = "1"
os.environ["MKL_NUM_THREADS"] = "1"
os.environ["NUMEXPR_NUM_THREADS"] = "1"
from cosipy.spacecraftfile import SpacecraftHistory
from cosipy.response.FullDetectorResponse import FullDetectorResponse
from cosipy.util import fetch_wasabi_file
from cosipy.statistics import PoissonLikelihood
from cosipy.background_estimation import FreeNormBinnedBackground
from cosipy.interfaces import ThreeMLPluginInterface
from cosipy.response import BinnedThreeMLModelFolding, BinnedInstrumentResponse, BinnedThreeMLPointSourceResponse
from cosipy.data_io import EmCDSBinnedData
from cosipy.source_injector import SourceInjector
from histpy import Histogram
from scoords import SpacecraftFrame
from astropy.time import Time
import astropy.units as u
from astropy.coordinates import SkyCoord, Galactic
import numpy as np
import matplotlib.pyplot as plt
from threeML import *
from astromodels import *
from pathlib import Path
%matplotlib inline
Download the data
[3]:
data_path = Path("/localscratch/sgallego/linkToXauron/COSIpyData/DC4_data") # /path/to/files. Current dir by default
[19]:
#Download the Orientation file
fetch_wasabi_file('COSI-SMEX/develop/Data/Orientation/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.fits', checksum = 'e86df2407eb052cf0c1db4a8e7598727')
Downloading cosi-pipeline-public/COSI-SMEX/develop/Data/Orientation/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.fits
[ ]:
#Download the response
#Choose the response depending on your model
#Model for a 511 keV line
#fetch_wasabi_file('COSI-SMEX/develop/Data/Responses/Response511.o4.e509_513.s20881894470591.m2555.filtered.nonsparse.binnedimaging.imagingresponse.h5',checksum = 'd15cbd03b2cd3d50f3411cd6b8ee362e')
#Model for a 1157 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Responses/Response44Ti.o4.e1154_1160.s9607532021290.m1215.filtered.nonsparse.binnedimaging.imagingresponse.h5',checksum = 'e965fdf8a7ee50cd6b300bd83778237b')
#Model for a 1173 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Responses/Response60FeLow.o4.e1170_1176.s9552269354945.m1188.filtered.nonsparse.binnedimaging.imagingresponse.h5',checksum = '76ea57aef5c73d1471e77703935657f1')
#Model for a 1333 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Responses/Response60FeHigh.o4.e1329_1336.s10201526728102.m1287.filtered.nonsparse.binnedimaging.imagingresponse.h5', checksum = 'b6fe847dd29570a991fc5316d39b3195')
#Model for a 1809 keV line
fetch_wasabi_file('COSI-SMEX/DC4/Data/Responses/Response26Al.o4.e1805_1812.s10036231691364.m1045.filtered.nonsparse.binnedimaging.imagingresponse.h5',checksum = '4de7c15f95c62698ff139585f5f5cdd1')
#Model for an other line or continuum source
#fetch_wasabi_file('COSI-SMEX/develop/Data/Responses/ResponseContinuum.o3.e100_10000.b10log.s10396905069491.m2284.filtered.nonsparse.binnedimaging.imagingresponse.h5',checksum = '7121f094be50e7bfe9b31e53015b0e85')
[ ]:
#Download the background model
#Choose the file depending on your model
#Model for a 511 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Backgrounds/gal_DC4_bkg_511_binned_smoothed_fwhm_10.0deg.hdf5',checksum = '1b437be9b4d9e85a230dfe4f4d2ca67e')
#Model for a 1157 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Backgrounds/gal_DC4_bkg_Ti44_binned_smoothed_fwhm_10.0deg.hdf5',checksum = 'bd62bfde46628a6e68434a6c68c7167d')
#Model for a 1173 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Backgrounds/gal_DC4_bkg_Fe60_low_binned_smoothed_fwhm_10.0deg.hdf5',checksum = '57f972b0a1c2ece4626c85ced3cc050f')
#Model for a 1333 keV line
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Backgrounds/gal_DC4_bkg_Fe60_high_binned_smoothed_fwhm_10.0deg.hdf5',checksum = 'db56a9fd9d8fc37cb03fba5e907a80c3')
#Model for a 1809 keV line
fetch_wasabi_file('COSI-SMEX/DC4/Data/Backgrounds/gal_DC4_bkg_Al26_binned_smoothed_fwhm_10.0deg.hdf5',checksum = 'bd6a32f8abee5b989b87149111487faa')
#Model for an other line or continuum source
#fetch_wasabi_file('COSI-SMEX/DC4/Data/Backgrounds/gal_DC4_bkg_continuum_binned_smoothed_fwhm_10.0deg.hdf5',checksum = '1acceb619b0d7e43eea8cb29ffe74279')
Define your source model
For this example we are using a Gaussian but other functions can be used (see https://astromodels.readthedocs.io/en/latest/notebooks/function_list.html#Included-1D-Functions)
Additionaly you can also provide your own spectrum template from a file. See Template (Table) Models in : https://threeml.readthedocs.io/en/stable/notebooks/spectral_models.html
[13]:
#Position in Galactic coord of the source
l = 184.56
b = -5.78
F = 0.45e-5 / u.cm / u.cm / u.s
mu = 1809*u.keV
sigma = 0.3*u.keV
spectrum = Gaussian()
spectrum.F.value = F.value
spectrum.F.unit = F.unit
spectrum.mu.value = mu.value
spectrum.mu.unit = mu.unit
spectrum.sigma.value = sigma.value
spectrum.sigma.unit = sigma.unit
# Set spectral parameters for fitting:
spectrum.F.free = True
spectrum.mu.free = False
spectrum.sigma.free = False
#set limit for the param
spectrum.F.min_value = 0
spectrum.F.max_value = 1
source = PointSource("source", # Name of source (arbitrary, but needs to be unique)
l=l, # Longitude (deg)
b=b, # Latitude (deg)
spectral_shape=spectrum) # Spectral model
model = Model(source) # Model with single source. If we had multiple sources, we would do Model(source1, source2, ...)
Open the Response, Orientation and background files
[5]:
#open the Response file
#replace the background and response file by the one you want to use
responsefile = "Response26Al.o4.e1805_1812.s10036231691364.m1045.filtered.nonsparse.binnedimaging.imagingresponse.h5"
bkgfile = "gal_DC4_bkg_Al26_binned_smoothed_fwhm_10.0deg.hdf5"
sc_orientation_path = "DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth_SAA.fits"
dr = FullDetectorResponse.open(data_path / responsefile)
#Create the point source injector
injector = SourceInjector(response_path=data_path / responsefile, response_frame = "spacecraftframe")
#open the background model
bg_model = Histogram.open(data_path / bkgfile)
#open the Orientation file
sc_orientation = SpacecraftHistory.open(data_path / sc_orientation_path)
WARNING: XDG_CACHE_HOME is set to '/localscratch/.cache/sgallego', but the default location, /uni-mainz.de/homes/sgallego/.astropy/cache, already exists, and takes precedence. This environment variable will be ignored. [astropy.config.paths]
WARNING: XDG_CACHE_HOME is set to '/localscratch/.cache/sgallego', but the default location, /uni-mainz.de/homes/sgallego/.astropy/cache, already exists, and takes precedence. This environment variable will be ignored. [astropy.config.paths]
Compute the sensitivity
The Likelihood (LH) fit is ran hundred times (more would be better but depending on how fast it is you could increase this number). Each time a new signal dataset is created by convolving the model with the response and applying some poisson statistic (also called the source injector). Similarly for the background, a dataset is created by taking a poisson sample of the background model.
For each iteration the LH ratio test is computed :
\[\text{LHtest} = -2 * (L_0 - L_1)\]
\(L_0\) being the LH value for the background only hypothesis and \(L_1\) for background+signal . Assuming the Wilk’s theorem is applying, the LH ratio test should follow a \(\chi^2\) distribution with dof = nb of parameters. The median value of the LH ratio test distribution should tell us if COSI is sensitive or not to this flux value.
For our example we have only one parameter we test : the signal strength . So for a 3 \(\sigma\) sensitivity we want to have the median superior or equal to 9.
Keep in mind we are testing only 3 months of observations but we will scale the flux later to 2 years of observation if it is already sensitive after 3 months.
⚠ !! WARNING !! ⚠ In this example we used a predifined Gaussian function from Threeml so the parameter value is F and is call with spectrum.F.value . If you changed the spectrum shape model, please change in the code below F with the new name of your parameter. To know it, you can use model.display()
[15]:
# Calculate the response once
spectrum.F.value = F.value
signal = injector.inject_model(model = model, orientation = sc_orientation , fluctuate=True)
bkg_poissonsamp = Histogram(bg_model.axes)
bkg_poissonsamp[:] += np.random.poisson(lam = bg_model.to_dense()[:])
data = EmCDSBinnedData(signal.project('Em', 'Phi', 'PsiChi') + bkg_poissonsamp.project('Em', 'Phi', 'PsiChi'))
instrument_response = BinnedInstrumentResponse(dr,data)
psr = BinnedThreeMLPointSourceResponse(data = data,
instrument_response = instrument_response,
sc_history=sc_orientation,
energy_axis = dr.axes['Ei'],
polarization_axis = dr.axes['Pol'] if 'Pol' in dr.axes.labels else None,
nside = 2*data.axes['PsiChi'].nside)
response = BinnedThreeMLModelFolding(data = data, point_source_response = psr)
[19]:
# We start with 100 but higher would be better
# To adjust depending on how fast the LH fit is
# You can also run a script in parallel and save the likelihood ratio test
nbiter = 100 #put this number to 100
LHratiotest = []
for i in range(nbiter) :
print(f"ITERATION {i}")
#create a signal dataset with the source injector
spectrum.F.value = F.value
spectrum.F.free = True
signal = injector.inject_model(model = model, orientation = sc_orientation ,make_spectrum_plot = False ,data_save_path = None,fluctuate=True)
#create a bkg dataset from the bkg model
bkg_poissonsamp = Histogram(bg_model.axes)
bkg_poissonsamp[:] += np.random.poisson(lam = bg_model.to_dense()[:])
# Set background parameter, which is used to fit the amplitude of the background, and instantiate the COSI 3ML plugin
# here we define only one bkg model but we could add more
bkg_dist = {"Background":bg_model.project('Em', 'Phi', 'PsiChi')}
# Workaround to avoid inf values. Out bkg should be smooth, but currently it's not.
for bckfile in bkg_dist.keys() :
bkg_dist[bckfile] += sys.float_info.min
#combine the data + the bck like we would get for real data
data_i = EmCDSBinnedData(signal.project('Em', 'Phi', 'PsiChi')
+ bkg_poissonsamp.project('Em', 'Phi', 'PsiChi')
)
# check psr caching
# print('previous psr._data:', np.sum( psr._data._data ) )
psr._data = data_i
# print('current psr._data:', np.sum( psr._data._data ) )
# print('data_i:', np.sum( data_i.data ))
bkg = FreeNormBinnedBackground(bkg_dist,
sc_history=sc_orientation,
copy = True)
like_fun = PoissonLikelihood(data_i, response, bkg)
cosi = ThreeMLPluginInterface('cosi',
like_fun,
response,
bkg)
# Nuisance parameter guess, bounds, etc.
cosi.bkg_parameter['Background'] = Parameter("Background", # background parameter
1, # initial value of parameter
min_value=0, # minimum value of parameter
max_value=10, # maximum value of parameter
delta=0.05, # initial step used by fitting engine
unit = u.Hz
)
# Set the source Model
cosi.set_model(model)
plugins = DataList(cosi) # If we had multiple instruments, we would do e.g. DataList(cosi, lat, hawc, ...)
like = JointLikelihood(model, plugins, verbose = False)
#do the fit
like.fit(compute_covariance = False,quiet=True)
# like.fit()
#get L1
L1 = cosi.get_log_like()
#do the fit again for null hypothesis (no signal)
spectrum.F.value = 0
spectrum.F.free = False
like.fit(compute_covariance = False,quiet=True)
# like.fit()
#Get L0
L0 = cosi.get_log_like()
#save LH ratio test
LHratiotest.append( -2 * (L0-L1) )
ITERATION 0
11:20:17 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 1
11:20:23 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 2
11:20:28 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 3
11:20:34 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 4
11:20:41 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 5
11:20:47 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 6
11:20:53 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 7
11:20:59 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 8
11:21:05 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 9
11:21:11 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 10
11:21:17 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 11
11:21:23 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 12
11:21:29 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 13
11:21:35 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 14
11:21:41 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 15
11:21:47 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 16
11:21:53 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 17
11:22:00 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 18
11:22:06 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 19
11:22:12 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 20
11:22:18 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 21
11:22:24 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 22
11:22:31 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 23
11:22:37 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 24
11:22:43 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 25
11:22:49 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 26
11:22:55 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 27
11:23:01 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 28
11:23:07 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 29
11:23:13 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 30
11:23:19 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 31
11:23:25 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 32
11:23:31 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 33
11:23:37 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 34
11:23:43 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 35
11:23:50 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 36
11:23:56 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 37
11:24:02 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 38
11:24:08 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 39
11:24:15 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 40
11:24:21 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 41
11:24:27 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 42
11:24:33 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 43
11:24:39 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 44
11:24:45 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 45
11:24:52 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 46
11:24:58 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 47
11:25:04 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 48
11:25:10 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 49
11:25:16 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 50
11:25:22 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 51
11:25:29 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 52
11:25:35 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 53
11:25:41 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 54
11:25:47 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 55
11:25:53 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 56
11:25:59 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 57
11:26:06 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 58
11:26:12 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 59
11:26:18 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 60
11:26:25 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 61
11:26:31 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 62
11:26:37 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 63
11:26:43 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 64
11:26:50 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 65
11:26:56 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 66
11:27:02 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 67
11:27:08 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 68
11:27:14 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 69
11:27:21 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 70
11:27:27 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 71
11:27:33 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 72
11:27:39 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 73
11:27:46 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 74
11:27:52 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 75
11:27:58 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 76
11:28:04 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 77
11:28:10 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 78
11:28:17 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 79
11:28:23 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 80
11:28:31 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 81
11:28:38 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 82
11:28:44 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 83
11:28:50 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 84
11:28:56 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 85
11:29:03 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 86
11:29:10 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 87
11:29:16 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 88
11:29:22 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 89
11:29:27 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 90
11:29:34 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 91
11:29:40 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 92
11:29:46 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 93
11:29:52 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 94
11:29:58 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 95
11:30:05 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 96
11:30:11 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 97
11:30:17 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 98
11:30:23 INFO set the minimizer to minuit joint_likelihood.py:1017
ITERATION 99
11:30:29 INFO set the minimizer to minuit joint_likelihood.py:1017
Plot the LH ratio test
When testing high flux value for your model, some “islands” in the TS distribution could appears (peaks with very high TS). This happens when a signal event lands in a bin where the estimated background was originally 0, and then filled with an arbitrary floor value. This is not a problem with the sensitivity estimation algorithm as this doesn’t affect significantly the median of a TS distribution of a flux near threshold, so the estimated sensitivity is still sound. See this issue for more informations : https://github.com/cositools/cosipy/issues/598
[18]:
if np.median(LHratiotest) <= 9 :
print(f"COSI is not sensitive to this model with a flux of {F.value} ph/cm2/s. Please try a higher flux")
else :
print(f"COSI is 3sigma sensitive to {F.value} ph/cm2/s. You can continue to the next cell")
fig,ax = plt.subplots()
ax.hist(LHratiotest,bins=100)
ax.axvline(np.median(LHratiotest),label=f"median {np.median(LHratiotest):.3f}",color="r",linestyle="--")
ax.legend()
ax.set_xlabel("LH ratio test")
ax.set_ylabel("counts")
ax.set_title(f"Model - flux {F.value} ph/cm2/s")
COSI is not sensitive to this model with a flux of 4.5e-06 ph/cm2/s. Please try a higher flux
[18]:
Text(0.5, 1.0, 'Model - flux 4.5e-06 ph/cm2/s')
Extrapolate the sensitivity to 2 years
DC4 simulations only includes 3 months of orbit. To get the sensitivity after 2 years (nominal life time of COSI) as a first approximation we can scale the flux by \(\frac{1}{\sqrt{T}}\)
[39]:
# 3 * 8 = 24 months = 2 years
sensi_2y = F/np.sqrt(8)
print(f"3sigma sensitivity after 2 years of COSI observation : {sensi_2y}")
3sigma sensitivity after 2 years of COSI observation : 1.0606601717798212e-05 1 / (s cm2)