Add unpolarized unbinned source injector capability [Yong2Sheng's]

[israelmcmc]

I'm bringing @Yong2Sheng from my fork to develop, now that the interfaces branch is merged.

Summary

This PR introduces an initial implementation for unpolarized unbinned continuum point source injector in the spacecraft frame. It follows the existing line injector workflow, but with additional steps to (1) compute the integrated detected rate across an energy band and (2) sample detected photon energies across that band.

Compute incident photons and simulate events

1) Define a continuum spectrum (differential photon flux)

A power-law spectrum is used as the working example in the tutorial:

$$S(E)=K\left(\frac{E}{\mathrm{piv}}\right)^{\mathrm{index}}, \qquad [S(E)] = \frac{\mathrm{ph}}{\mathrm{cm^2}\times\mathrm{s}\times\mathrm{keV}} $$

where $K$ is the normalization (differential flux at the pivot energy $\mathrm{piv}$).

2) Convolve the spectrum with effective area

The detected count-rate density is

$$ w(E) \equiv S(E)\times A_{\text{eff}}(E), \qquad [w(E)] = \frac{\mathrm{ph}}{\mathrm{s}\times\mathrm{keV}}. $$

3) Integrate $w(E)$ to get the total detected rate and total counts

The integrated detected photon rate is

$$ I_0 = \int_{E_\mathrm{min}}^{E_\mathrm{max}} w(E) dE, \qquad [I_0] = \frac{\mathrm{ph}}{\mathrm{s}}. $$

In practice, the spectrum might not be analytical. And the effective area is also usually not analytical as well. So we do not perform the integral directly. Instead, we will use a very fine energy grid to obtain the integral by trapezoids.

The expected detected counts is

$$ \mu = T\times I_0, $$

and the realized detected counts are sampled as

$$ N \sim \mathrm{Poisson}(\mu). $$

4) Build an inverse-CDF sampler to draw detected photon energies

• Construct a discrete CDF from $w(E)$ on the energy grid via cumulative trapezoidal integration, and normalized it by $I_0$ to $[0,1]$.
• Draw a uniform distribution from [0,1): $u \sim \mathrm{Unif}(0,1)$.
• Apply inverse transform sampling: $E_i = CDF^{-1}(u)$, where $CDF^{-1}$ is implemented using scipy.interpolate.interp1d with linear interpolation.

5) Simulate each sampled energy

Feed the sampled photon energies (and their multiplicities, although they should show up only once) into UnpolarizedIdealComptonIRF.random_events(...) to produce the event list.

Implementation notes

1. A new class RandomEventDataFromContinuumInSCFrame is responsible for continuum injection.
2. For performance and clarity, I cache the following intermediates with @cached_property:
   • energy_grid
   • unpolarized_aeff
   • unpol_counts
   • incident_photons
3. A structured return container for detected-count (unpol_counts) bookkeeping:
   • UnpolCountsSummary(rate_density, rate, total_expected_counts)

What is not included yet

• Polarization handling for the continuum branch (this PR is unpolarized only).

[codecov]

Codecov Report

❌ Patch coverage is 85.81925% with 193 lines in your changes missing coverage. Please review.
✅ Project coverage is 71.35%. Comparing base (6765236) to head (3e7dd63).
⚠️ Report is 196 commits behind head on develop.

Files with missing lines	Patch %	Lines
cosipy/response/ideal_response.py	0.00%	93 Missing ⚠️
...background_estimation/ContinuumEstimationInterp.py	88.18%	26 Missing ⚠️
.../background_estimation/ml/ContinuumEstimationNN.py	86.63%	25 Missing ⚠️
...e_deconvolution/data_interfaces/dataIF_Parallel.py	25.00%	15 Missing ⚠️
...ution/data_interfaces/data_interface_collection.py	85.96%	8 Missing ⚠️
...deconvolution/algorithms/RichardsonLucyAdvanced.py	94.95%	6 Missing ⚠️
...ge_deconvolution/algorithms/RichardsonLucyBasic.py	94.36%	4 Missing ⚠️
...y/image_deconvolution/algorithms/RichardsonLucy.py	96.72%	2 Missing ⚠️
...lution/algorithms/accelerators/accelerator_base.py	95.45%	2 Missing ⚠️
...ution/algorithms/accelerators/build_accelerator.py	87.50%	2 Missing ⚠️
... and 8 more

Files with missing lines	Coverage Δ
cosipy/__init__.py	92.59% <100.00%> (+3.70%)	⬆️
cosipy/background_estimation/__init__.py	100.00% <100.00%> (ø)
cosipy/event_selection/good_time_interval.py	93.38% <100.00%> (+1.15%)	⬆️
cosipy/image_deconvolution/__init__.py	100.00% <100.00%> (ø)
cosipy/image_deconvolution/algorithms/__init__.py	100.00% <100.00%> (ø)
..._deconvolution/algorithms/accelerators/__init__.py	100.00% <100.00%> (ø)
...on/algorithms/accelerators/max_step_accelerator.py	100.00% <100.00%> (ø)
...olution/algorithms/deconvolution_algorithm_base.py	94.11% <100.00%> (ø)
...osipy/image_deconvolution/algorithms/prior_base.py	100.00% <ø> (ø)
...py/image_deconvolution/algorithms/prior_entropy.py	100.00% <100.00%> (ø)
... and 39 more

... and 10 files with indirect coverage changes

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Yong2Sheng]

I will add the unit tests after the changes have been reviewed and approved. I’d like to avoid adding them too early, since they may need to be updated if the code changes during review.

[israelmcmc]

See my comment in #528

This is also related to #499: if the user passes a UnbinnedThreeMLModelFoldingInterface, we can inject event data right from the same class.

It's fine for now, but eventually this capability should be moved to the SourceInjector class.

[read-the-docs-community]

Documentation build overview

📚 cositools-cosipy | 🛠️ Build #32212080 | 📁 Comparing 3e7dd63 against latest (92bb057)

  🔍 Preview build  

Show files changed (10 files in total): 📝 8 modified | ➕ 2 added | ➖ 0 deleted

File	Status
index.html	📝 modified
api/background_estimation.html	➕ added
api/index.html	📝 modified
tutorials/index.html	📝 modified
tutorials/image_deconvolution/511keV-Galactic-ImageDeconvolution.html	📝 modified
tutorials/source_injector/Unbinned_Source_Injector.html	➕ added
tutorials/background_estimation/continuum_estimation/BG_estimationNN_example.html	📝 modified
tutorials/image_deconvolution/optional/scatt_binned_data_w_local_coordinate/511keV-ScAtt-DataReduction.html	📝 modified
tutorials/image_deconvolution/optional/scatt_binned_data_w_local_coordinate/511keV-ScAtt-ImageDeconvolution.html	📝 modified
tutorials/image_deconvolution/optional/time_binned_data_w_local_coordinate/GRB-ScAtt-ImageDeconvolution.html	📝 modified

[Yong2Sheng]

Hi @ldigesu , I have added the saving function. Could you please take a look? Thank you!

[ldigesu]

Hello Yong,
The tutorial runs ok for me; it simulates and saves the output as a fits file correctly. My main comments are about explanations in the tutorial and sanity checks.

I suggest adding an introductory paragraph at the beginning of the notebook, describing what this simulator does and what it is expected to produce as outputs. Something like the summary written in this PR should be ok.

Secondly, I thinks the tutorial is missing some sanity checks of the outputs. All the plots run correctly to show the simulated data, but there is no comparison with what was input, thus no checks that what was produced by the simulator is indeed what we expected. The saved fits output cannot be used in a binned analysis: the binning function does not work, because it can't find the Time Tags of the photons. Thus, my suggestion as a sanity check is the following: instead of showing as an example a random point source, you could take one of the point sources of DC3 e.g., a bright GRB or the Crab nebula, and simulate exactly the same position, spectrum, source duration, ecc.. Then you can produce plots from your simulator and from the fits file produced by Megalib and show that they are consistent. Let me know if this could work.

Thirdly, I remind you Israel's comment to PR #528 (here above), that is relevant to this PR.

Finally, a last minor comment: I find the last plot showing the distribution of the photons in the CDS not particularly explanatory. I understand that there are no Galactic coordinates in these data, thus no way to visualize a ring around the expected source position. However, maybe you could try to you apply some cuts in Energy or Phi angles before plotting and see if point clusters. I suggest in any case to add a note for inexperienced COSI users explaining that in local/spacecraft coordinates, no ring/cone can be visualized.

Thanks! Let me know If everything is clear,
Laura

[israelmcmc]

Thanks for the comments @ldigesu

About this:

Finally, a last minor comment: I find the last plot showing the distribution of the photons in the CDS not particularly explanatory. I understand that there are no Galactic coordinates in these data, thus no way to visualize a ring around the expected source position. However, maybe you could try to you apply some cuts in Energy or Phi angles before plotting and see if point clusters. I suggest in any case to add a note for inexperienced COSI users explaining that in local/spacecraft coordinates, no ring/cone can be visualized.

Maybe we can do a better job for that plot and the explanation, but the point of the plot is that if you take all the event with a given phi (color coded) then they should make a circle with radius=phi around the location of the source (the cross). It's a visualization of the Compton cone, using a colorbar for Phi instead ofa "z-axis".