Extended Source Response Generator

[hiyoneda]

Regarding #129 and #267, I added new functions into FullDetectorResponse to generate an extended source response from a full detector response.

• get_extended_source_response : it generates an extended source response using a full detector response and cosipy.spacecraftfile.SpacecraftFile.
• get_point_source_response_per_image_pixel : it generates a point source response for a specific HEALPix pixel using a full detector response and cosipy.spacecraftfile.SpacecraftFile.
• merge_psr_to_extended_source_response : it creates an extended source response by merging multiple point source responses.
• _setup_extended_source_response_params : it is a supplementary function to setup NSIDE parameters for an image and a scatt map.

Basically, an extended source response can be generated like this.

from cosipy.spacecraftfile import SpacecraftFile
from cosipy.response import FullDetectorResponse, ExtendedSourceResponse

full_detector_response = FullDetectorResponse.open("SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5")

orientation = SpacecraftFile.parse_from_file("20280301_3_month_with_orbital_info.ori")

extended_source_response = full_detector_response.get_extended_source_response(orientation,
                                                                               coordsys='galactic',
                                                                               nside_scatt_map=None,
                                                                               Earth_occ=True)

extended_source_response.write("extended_source_response_continuum.h5", overwrite = True)

As default, the NSIDE of the scatt map is the same as that of the full detector response, but when we want to use a finer binned scatt map, its NSIDE parameter can be modified by setting the input parameternside_scatt_map.

By using get_point_source_response_per_image_pixel and merge_psr_to_extended_source_response, we can generate a response on each pixel and merge the generated point source responses into an extended source response later. It can be used to generate a large extended source response using multiple nodes in a computer cluster.
On each node,

from cosipy.spacecraftfile import SpacecraftFile
from cosipy.response import FullDetectorResponse, ExtendedSourceResponse

full_detector_response = FullDetectorResponse.open("SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5")

orientation = SpacecraftFile.parse_from_file("20280301_3_month_with_orbital_info.ori")

ipix_image_list = [int(_) for _ in sys.argv[1:]]

basename = "psr/psr_"

for ipix_image in ipix_image_list:

    psr = full_detector_response.get_point_source_response_per_image_pixel(ipix_image, orientation, 
                                                                           coordsys='galactic',
                                                                           nside_image=None,
                                                                           nside_scatt_map=None,
                                                                           Earth_occ=True)
    
    psr.write(f"{basename}{ipix_image:08}.h5",overwrite = True)

When we merge them, the code can be like this.

from cosipy.spacecraftfile import SpacecraftFile
from cosipy.response import FullDetectorResponse, ExtendedSourceResponse

full_detector_response = FullDetectorResponse.open("SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5")

basename = "psr/psr_"

extended_source_response = full_detector_response.merge_psr_to_extended_source_response(basename, coordsys = 'galactic', nside_image = None)

extended_source_response.write("extended_source_response_continuum_merged.h5", overwrite = True)

[codecov]

Codecov Report

Attention: Patch coverage is 93.75000% with 3 lines in your changes missing coverage. Please review.

Project coverage is 79.74%. Comparing base (aea580b) to head (f1366ae).
Report is 6 commits behind head on develop.

Files with missing lines	Patch %	Lines
cosipy/response/FullDetectorResponse.py	93.61%	3 Missing ⚠️

Files with missing lines	Coverage Δ
cosipy/response/ExtendedSourceResponse.py	91.17% <100.00%> (+0.26%)	⬆️
cosipy/response/FullDetectorResponse.py	49.71% <93.61%> (+10.42%)	⬆️

[hiyoneda]

I've compared the generated extended source response with that used in DC2. I used "SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5" as a full detector response, and "20280301_3_month_with_orbital_info.ori" as an orientation file. Then, I compared the response generated with the new function with "psr_gal_continuum_DC2.h5". Here, the earth occultation is not considered.

They are almost the same, while there is a slight difference of about 0.1% when I projected them onto the Ei axis. I suppose the slight difference may be caused by a different parameter setup, maybe the NSIDE of the scatt map.

When I generated the response by considering the earth occultation, the effective area gets smaller in a higher energy band. I think that it is a correct behavior because a higher energy gamma ray can penetrate the shield and then the effective area gets larger if the earth occultation is ignored.

Also, I confirmed that the response is identical with the standard method when I used get_point_source_response_per_image_pixel + merge_psr_to_extended_source_response.

[hiyoneda]

@ckarwin can you check the new functions roughly? If there is no critical issue, I'll continue to make unit tests.

[hiyoneda]

With nside_scatt_map=full_detector_response.nside*2, the generated response becomes identical to psr_gal_continuum_DC2.h5. So, we can understand that the slight difference above was caused by the different nside_scatt_map.

[ckarwin]

Thanks for testing that @hiyoneda. I'm going to test your code by producing the 511 PSR for DC3 on the NASA cluster. I should be able to do this early next week. If everything works ok then I can merge your PR.

[hiyoneda]

@israelmcmc Can I ask to make a full detector response in a dense matrix for the unit test? Currently, the full detector response for the unit test is in the old format; its axes include 'SigmaTau' and 'Dist', and the response matrix is sparse. Actually, with it, get_point_source_response, which is used for the extended response generation, raises an error because Coord + scatt_map is currently only supported for dense responses.

[ckarwin]

@hiyoneda I have been testing your PR and it works well. Very nice job! I have some suggestions before merging:

1. Please add unit tests. I guess your waiting for the detector response from @israelmcmc.
2. The examples that you provided above work well and are very helpful. I think you should add this to the examples directory. It might not work so well as a Jupyter notebook, unless you use the same response from the unit tests. Or maybe you can just add python scripts. What do you thin @israelmcmc, is there a good place in the docs directory for example scripts? Or is it better to do Jupyter notebooks?
3. Your example for get_point_source_response_per_image_pixel needs import sys.

Currently, I'm calculating the 511 keV PSR for DC3 on the NASA cluster. If it looks ok, I think we can merge this PR after addressing the above comments.

[israelmcmc]

Sorry for the delay in getting the test response. I'll work on it to not delay further the unit test.

Yes, I think it's time to start adding some scripts. The scripts and Jupyter notebooks serve a different purpose though, and they are not mutually exclusive. Jupyter notebook are meant to explain to the user how to use the library code, include more verbose than a script, and don't need to be very efficient. In general, we could have both a script and a Jupyter doing similar things.

In this particular case, I think this should be a script for production, not an example --i.e. it doesn't require the extra work to hold the users' hand and doesn't need to be in the docs folders.

We'll have a discussion next Thursday (Feb 6) about scripts. It was prompted by the SDOC group, but it also applies to this thread.

I think the way to go is to create scripts using setuptools' entry points. I wrote an example a while ago. See the entry_points item in setup.py and the corresponding function. I think it's out of date now, but enough to get an idea of how to proceed.

[augustus-thomas]

As a note on creating scripts that serve similar purposes to Jupyter Notebooks. If you want to save time, I have found that to convert a notebook to script quickly you can use jupyter nbconvert --to script "path/to/notebook.ipynb" . As @israelmcmc mentioned though, scripts should be more efficient and not exactly a copy-paste. Also this would not follow the entry-points guidance. But it could save development time and I wasn't sure if people were aware of it.

[israelmcmc]

Thanks, @augustus-thomas! I was not aware of it. It does seem useful as a starting point for a script. I'm tagging @ldigesu who's working on converting some of our jupyter nb to executables in order the build the automatic pipelines. Please see above, Laura.

[israelmcmc]

@hirokiyoneda Sorry for the delay. You can find the test file in PR #287

I changed the tests to use this file instead. I also committed a small change to the code that allowed me to create this file in an easier way.

If you can review it while you are at it I'd appreciate it 😁

[hiyoneda]

@israelmcmc Thanks for making the new test response file. I've checked it and PR #287; they look great. I want to ask others who worked on the full detector response before; from my side, I think it can be merged once open_rsp is covered by the unit test. For this PR, I’ll make the unit test using the new dense test response file.

[israelmcmc]

All sounds good, thanks @hiyoneda. I'll come back to the open_rsp test, although it's not a high priority at the moment.

[ckarwin]

@hiyoneda I tested your PR by generating the 511 keV PSR on the NASA cluster. Everything worked fine, and I can generate the PSR very quickly (currently using 1000 parallel cores). As a sanity check, below is an example image taking a slice and projecting onto PsiChi:
NuLambda: 1440 (Galactic center)
Em: 0
Phi: 5

I believe that we still need to update COSILike in order to seamlessly use the precomputed PSR for analyzing point sources. I didn't see an issue for this, but it's documented in the code comments. I don't think it needs to be done for this PR, but maybe we can at least open an issue (I didn't see one). Please let me know what you think.

Once you finish adding the unit tests I think this PR can be merged.

[israelmcmc]

Pretty cool, thank you both!

I'll open the issue about pre-computed PSR for point sources.

[israelmcmc]

Here it is #288

[hiyoneda]

@ckarwin Thanks for testing the code with the NASA cluster. I am happy to see the attached figure. I have added the unit tests with the new test response which @israelmcmc generated. I have also added the example python scripts in docs/tutorials/response/ (hope that I added them in the correct directory). I think that we've done things mostly.

[ckarwin]

Looks good, and all checks have passed. @israelmcmc Are you ok with me merging this? Maybe we can implement the scripts as entry points later on?

[israelmcmc]

Yeah, that's all right. We can move the scripts around later. Thanks!

[ckarwin]

All looks good now. Merging.