New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add exposure image computation #511

Merged
merged 20 commits into from Apr 18, 2016

Conversation

Projects
None yet
2 participants
@JouvinLea
Contributor

JouvinLea commented Apr 14, 2016

I added the computation of the exposure for the images.
For the moment I define it in two ways make_exposure_one_obs() and make_exposure_one_obs2(). I wanted to discuss with you which is the best.
I also added the excess images.

  1. for make_exposure_one_obs(), excess/exposure will give the integrated flux over the energy range used to compute the images.

  2. for make_exposure_one_obs2(), excess/exposure will give the differential flux at the center of the energy_range used to compute the images.

We think the best is to use the integrated flux but maybe you have some ideas about that?

counts_map.fill(events=events)
counts_map.data = counts_map.data.value
self.maps["counts"] = counts_map
if (len(events) != 0):

This comment has been minimized.

@JouvinLea

JouvinLea Apr 14, 2016

Contributor

For low energy range, there are runs without any events at these energies. I save the id number of these observations to not compute the bkg or the exposure for these runs...

This comment has been minimized.

@cdeil

cdeil Apr 14, 2016

Member

I think keeping a list of runs to skip and using it across method calls is unnecessarily complex.

My suggestion would be to just remove this, i.e. process obs with zero counts.
They should have zero background and everything should work just fine when stacking.

Also, speed shouldn't be big issue here ... this is premature optimisation, no?
If you want to avoid computing the background model for runs with zero counts,
I'd suggest putting the check for zero counts at the start of that function (or in a method _has_counts(obs_id))
instead of pre-computing and caching these flags in a self.pass_run member.

Sometimes such caches are needed for performance, but they are often error-prone, especially if implemented in this way, and in this case I don't think it's needed at all. Just checking for zero counts is much faster than background model computation, so it's OK if the zero-counts check is re-computed a few times (if that is even the case).

OK or am I missing something?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 16, 2016

Contributor

I also want to change that in this PR if you think it is necessary!!
I agree with you that this is quite complex. Ok to add a _has_counts(obs_id) method but in this case I will have also to check that for the exposure so I was afraid that it could be a computation time problem...
Le 14 avr. 2016 à 15:07, Christoph Deil notifications@github.com a écrit :

In gammapy/scripts/image_pipe.py:

@@ -64,9 +67,12 @@ def make_counts(self, obs_id):
obs = self.data_store.obs(obs_id=obs_id)
events = obs.events.select_energy(self.energy_band)
events = events.select_offset(self.offset_band)

  •    counts_map.fill(events=events)
    
  •    counts_map.data = counts_map.data.value
    
  •    self.maps["counts"] = counts_map
    
  •    if (len(events) != 0):
    
    I think keeping a list of runs to skip and using it across method calls is unnecessarily complex.

My suggestion would be to just remove this, i.e. process obs with zero counts.
They should have zero background and everything should work just fine when stacking.

Also, speed shouldn't be big issue here ... this is premature optimisation, no?
If you want to avoid computing the background model for runs with zero counts,
I'd suggest putting the check for zero counts at the start of that function (or in a method _has_counts(obs_id))
instead of pre-computing and caching these flags in a self.pass_run member.

Sometimes such caches are needed for performance, but they are often error-prone, especially if implemented in this way, and in this case I don't think it's needed at all. Just checking for zero counts is much faster than background model computation, so it's OK if the zero-counts check is re-computed a few times (if that is even the case).

OK or am I missing something?


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

I think _has_counts(obs_id) is OK.

Checking counts is just a sum and should be blazing fast ... I think even if you do it multiple times per obs you won't even notice the slowdown really.

exposure_map.data += self.maps["exposure"].data
exposure_map2.data += self.maps["exposure2"].data
self.maps["total_exposure"] = exposure_map
self.maps["total_exposure2"] = exposure_map2

This comment has been minimized.

@JouvinLea

JouvinLea Apr 14, 2016

Contributor

This is just for my test that I added the total_exposure2 but once we choose which way we define the exposure we will only have a total_exposure in the skymap!

@cdeil cdeil changed the title from Exposure to Exposure image computation Apr 14, 2016

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 14, 2016

@JouvinLea - Thanks for working on this!

I discussed this with @adonath ... we want the exposure to be computed so that

differential flux = excess / exposure

(because it's slightly simpler and also that's how it's done in HAP / the HESS survey).

An alternative option to implement exposure images would be to always first compute an exposure cube
http://docs.gammapy.org/en/latest/api/gammapy.cube.exposure_cube.html

and then use this method to compute the integral exposure:
http://docs.gammapy.org/en/latest/api/gammapy.cube.SpectralCube.html?highlight=spectralcube#gammapy.cube.SpectralCube.integral_flux_image
It should be generalised a bit and renames energy_integral_power_law_weighted or something.

If you want to do this I can explain how in detail.
If you want to continue with what you have here, that's OK as well for now, but please leave a TODO comment that this could / should be re-implemented using the exposure_cube function and the spectral_cube.integral_flux_image method.

@cdeil cdeil added the feature label Apr 14, 2016

@cdeil cdeil added this to the 0.4 milestone Apr 14, 2016

@cdeil cdeil self-assigned this Apr 14, 2016

@@ -189,6 +189,13 @@ def nbins(self):
return self.size - 1
@property
def lin_centers(self):

This comment has been minimized.

@cdeil

cdeil Apr 14, 2016

Member

Why introduce lin_center. For the application you use (and any use case I'm aware of), the existing log_center is better, no?

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

You're still calling lin_centers.
I've argued before that that's never useful, you'll always get better results by using log_center.
Did you want to stick with lin_centers? Did I miss your argument why?

@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 15, 2016

@cdeil
Hi,
thanks for your suggestion! I will not have the time to re-implement it with the spectral cube for now so this is ok for you I prefer leave it as it is and add a todo comment.
I'm not sure about the exposure but maybe I can name one of the method make_exposure() and the other one make_exposure_power_law_weighted() ? Like this the user can choose what he wants?

I though that it was easier to cut the energy range of the map in linear energy bins since this is only there to compute the integral. You are not agree?

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 15, 2016

thanks for your suggestion! I will not have the time to re-implement it with the spectral cube for now so this is ok for you I prefer leave it as it is and add a todo comment.

OK

I'm not sure about the exposure but maybe I can name one of the method make_exposure() and the other one make_exposure_power_law_weighted() ? Like this the user can choose what he wants?

You pick names!
(the ones I was suggesting were for the more general SpectralCube class, here just make_exposure is OK, no?)

I though that it was easier to cut the energy range of the map in linear energy bins since this is only there to compute the integral. You are not agree?

Both linear and log equal-width energy bins give the same result for very fine binning.
For a given number of bins, equal-log binning is usually always better, because then quantities are more constant in gamma-ray astronomy within the bins.

But there's several places where equal-log integrals are implemented already.
If this is roughly working for now, OK to leave as-is.
(will be changed to use the equal-log binning integration in SpectralCube.integral_flux_image at some point anyways)

Let me know if you have any further questions or when this is ready to be merged.

@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 16, 2016

ok let’s go this because I don’t really have the time to work a lot on this!!
ok for make exposure that will give the differential flux when you do excess/exposure.
But may be we can live the one that give the integrated flux if you do excess/exposure and by giving it another name? but I don’t have an idea for the name? and maybe let the user choose what he prefers ?
Le 15 avr. 2016 à 20:57, Christoph Deil notifications@github.com a écrit :

thanks for your suggestion! I will not have the time to re-implement it with the spectral cube for now so this is ok for you I prefer leave it as it is and add a todo comment.

OK

I'm not sure about the exposure but maybe I can name one of the method make_exposure() and the other one make_exposure_power_law_weighted() ? Like this the user can choose what he wants?

You pick names!
(the ones I was suggesting were for the more general SpectralCube class, here just make_exposure is OK, no?)

I though that it was easier to cut the energy range of the map in linear energy bins since this is only there to compute the integral. You are not agree?

Both linear and log equal-width energy bins give the same result for very fine binning.
For a given number of bins, equal-log binning is usually always better, because then quantities are more constant in gamma-ray astronomy within the bins.

But there's several places where equal-log integrals are implemented already.
If this is roughly working for now, OK to leave as-is.
(will be changed to use the equal-log binning integration in SpectralCube.integral_flux_image at some point anyways)

Let me know if you have any further questions or when this is ready to be merged.


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 16, 2016

Maybe we can call it make_exposure_weighted_spectrum the one that will give the integrated flux ?

Le 16 avr. 2016 à 18:53, Lea Jouvin ljouvin@apc.in2p3.fr a écrit :

ok let’s go this because I don’t really have the time to work a lot on this!!
ok for make exposure that will give the differential flux when you do excess/exposure.
But may be we can live the one that give the integrated flux if you do excess/exposure and by giving it another name? but I don’t have an idea for the name? and maybe let the user choose what he prefers ?
Le 15 avr. 2016 à 20:57, Christoph Deil notifications@github.com a écrit :

thanks for your suggestion! I will not have the time to re-implement it with the spectral cube for now so this is ok for you I prefer leave it as it is and add a todo comment.

OK

I'm not sure about the exposure but maybe I can name one of the method make_exposure() and the other one make_exposure_power_law_weighted() ? Like this the user can choose what he wants?

You pick names!
(the ones I was suggesting were for the more general SpectralCube class, here just make_exposure is OK, no?)

I though that it was easier to cut the energy range of the map in linear energy bins since this is only there to compute the integral. You are not agree?

Both linear and log equal-width energy bins give the same result for very fine binning.
For a given number of bins, equal-log binning is usually always better, because then quantities are more constant in gamma-ray astronomy within the bins.

But there's several places where equal-log integrals are implemented already.
If this is roughly working for now, OK to leave as-is.
(will be changed to use the equal-log binning integration in SpectralCube.integral_flux_image at some point anyways)

Let me know if you have any further questions or when this is ready to be merged.


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 16, 2016

ok let’s go this

You mean I should merge now.
Or did you still want to make edits here?

The name isn't so important for now.

@JouvinLea

This comment has been minimized.

Contributor

JouvinLea commented Apr 16, 2016

no I will make some edit before you merge!
Le 16 avr. 2016 à 19:14, Christoph Deil notifications@github.com a écrit :

ok let’s go this

You mean I should merge now.
Or did you still want to make edits here?

The name isn't so important for now.


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

Lea Jouvin added some commits Apr 16, 2016

else:
self.maps["counts"] = counts_map
def _has_counts(self, obs_id):

This comment has been minimized.

@JouvinLea

JouvinLea Apr 16, 2016

Contributor

Is it ok now?
And I fill the bkg and exposure of zeros for these obs?

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

Return True or False, i.e. bool.

Implementation looks OK.

If you have a small energy band, there could be no counts, but still some exposure.
So at least for exposure you should always fill it, not look at the counts.
For background you can skip filling it as a speed optimisation, if you like.

This comment has been minimized.

@JouvinLea

JouvinLea Apr 16, 2016

Contributor

But if we don't have any counts for one obs this is false to fill the exposure no ?

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

No. Imagine you have perfect background suppression. And a very week source where you only get one photon per 10 hours. Then to get the correct flux, you need to stack the exposure from the runs also where you had zero counts.

For HESS this can only occur if you pick a very small high-energy band.

But exposure has to be accumulated for all runs.

Makes sense?

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

I take it back.

You have a background modeling technique that doesn't work at all if you have no counts.
So you have to skip runs with zero counts completely, i.e. not accumulate background and exposure.
Maybe emit a warning?

We'll have to think more about how to handle different background modeling methods with the same image / cube analysis pipeline...

This comment has been minimized.

@JouvinLea

JouvinLea Apr 16, 2016

Contributor

yes but If this is fill with zero it's ok no?
you prefer that I add a continue and pass the run instead of filling with zero?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

@cdeil
Just tell me if this is ok to fill with zero or if I have to skip it? after it will be ready to be merge I think!

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

How about this?

Change the no counts statement to a counts < min_counts statement where the user can set min_counts (maybe use 1 or 5 by default), change _has_counts to _low_counts and emit a log warning there:

log.warning('Skipping obs OBS_ID because it only has N_COUNTS counts.')

Then we can think more about how to handle very low counts data later.

OK?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

ok I made the changes.
what do you think?

if self._has_counts(obs_id) == 0:
self.maps["exposure"].data[:] = 0
else:
if exposure_weighted_spectrum:

This comment has been minimized.

@JouvinLea

JouvinLea Apr 16, 2016

Contributor

is it ok for you to have this parameters to choose which king of exposure computation you want?

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

I'd have to think about what a good API is (like this or separate methods).
No time. If this works for you, let's just leave as-is for now?

else:
self.make_exposure(obs_id, spectral_index)
exposure_map.data += self.maps["exposure"].data
self.maps["total_exposure"] = exposure_map

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

Do you document the keys in self.maps somewhere?
Maybe list them in the class-level docstring?

Lea Jouvin
@@ -16,6 +20,16 @@
class ImageAnalysis(object):
"""Gammapy 2D image based analysis.
Fill a `~gammapy.image.SkyMapCollection`:

This comment has been minimized.

@JouvinLea

JouvinLea Apr 16, 2016

Contributor

is it what you want?
How to have new line between each maps description in the docstring?

This comment has been minimized.

@cdeil

cdeil Apr 16, 2016

Member

Yes, this is what I wanted.

You have to put bullets and an empty line before the first bullet for Sphinx to render it correctly.

Like this:

Fill a `~gammapy.image.SkyMapCollection`:

* counts : counts for one obs
* bkg : bkg model for one obs

Lea Jouvin added some commits Apr 16, 2016

Lea Jouvin
Lea Jouvin
Lea Jouvin
exposure.data[offset >= self.offset_band[1]] = 0
self.maps["exposure"] = exposure
def make_exposure_weighted_spectrum(self, obs_id, spectral_index=2.3):

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

So the difference between make_exposure and make_exposure_weighted_spectrum is just that the first returns diff flux at given energy, and the other integral flux in energy band, and the difference is one constant factor for the whole image?

If yes, the second version make_exposure_weighted_spectrum and all the duplicated code there should be removed.
The caller can simply compute and apply that factor herself.
Or you provide a convenience method on this class to do it, but calling the other one or assuming it's called after, not by duplicating the code there.

self.maps["total_excess"].data = excess
def make_exposure(self, obs_id, spectral_index=2.3):
"""Compute the exposure map for one observation. Excess/exposure will give the differential flux

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Put the sentence how flux and exposure are related on a separate line after a blank line.

Also, if you can put a formula how exposure is defined, I think in this case it would be good.
(it's the integral over flux x area x time assuming a PL)
See here as an example how to insert a formula in a docstring.

----------
obs_id : int
Number of the observation
spectral_index : int

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Spectral index is type float and add description, i.e.

spectral_index : float
    Assumed power-law spectral index
else:
return True
def make_total_counts(self, ncounts_min):

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Here you forgot to put a default for ncounts_min.
Generally it's annying for the user to have to pass ncounts_min many times, and it's hard to keep the different defaults consistent.

How about making it a parameter that is set once in the constructor and stored as a data member, that is accessed from _low_counts?

Lea Jouvin added some commits Apr 17, 2016

aeff2d = obs.aeff
offset_tab = Angle(np.linspace(self.offset_band[0].value, self.offset_band[1].value, 10), self.offset_band.unit)
exposure_tab = np.sum(aeff2d.evaluate(offset_tab, energy_bin).to("cm2") * spectrum * energy_band, axis=1)
if exposure_weighted_spectrum:

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

@cdeil one method and I add the norm if the parameters is set at true by the users, is it ok for you?
Sorry for the energy, I changed to log_spacing now!

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Yes it's OK for me.

I only wanted to get rid of the duplicated code, so that we have one version to use and debug.

Merge?

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

just one minute! I wanted to fix the doctring
Le 17 avr. 2016 à 21:25, Christoph Deil notifications@github.com a écrit :

In gammapy/scripts/image_pipe.py:

  •    offset = coord.separation(center)
    
  •    obs = self.data_store.obs(obs_id=obs_id)
    
  •    livetime = obs.observation_live_time_duration
    
  •    # 2D Exposure computation on the self.energy_range and on an offset_tab
    
  •    energy = EnergyBounds.equal_log_spacing(self.energy_band[0].value, self.energy_band[1].value, 100,
    
  •                                            self.energy_band.unit)
    
  •    energy_band = energy.bands
    
  •    energy_bin = energy.log_centers
    
  •    eref = EnergyBounds(self.energy_band).log_centers
    
  •    spectrum = (energy_bin / eref) *\* (-spectral_index)
    
  •    aeff2d = obs.aeff
    
  •    offset_tab = Angle(np.linspace(self.offset_band[0].value, self.offset_band[1].value, 10), self.offset_band.unit)
    
  •    exposure_tab = np.sum(aeff2d.evaluate(offset_tab, energy_bin).to("cm2") \* spectrum \* energy_band, axis=1)
    
  •    if exposure_weighted_spectrum:
    
    Yes it's OK for me.

I only wanted to get rid of the duplicated code, so that we have one version to use and debug.

Merge?


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

Lea Jouvin
EXPOSURE = \int A(E) \phi(E) * time_{obs} \, dE
assuming a power law for the flux :math:`\phi(E) = \phi_{Eref} \times \frac{E}{E_ref}^{spectral_index}`

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

@cdeil
I don't know what is the problem with this line for the docstring? it doesn't appear as a formula?

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

The problem with the docstring is that the backslash has both a special meaning for Sphinx and for LaTeX.
So at the moment Sphinx does something with the backslashes \ in your formulae and then LaTeX sees an incorrect expression and errors out.

The solution is to use a "raw" string for the docstring, i.e. add an "r" character at the front of the docstring.
See here

To debug your LaTeX formulae, running Sphinx repeatedly is annoyingly slow.
For that I use http://www.chachatelier.fr/latexit/ and only once it's correct there, I copy & paste the formula over into the Python docstring.

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

I currently get this:

screen shot 2016-04-17 at 22 27 18

I think it's better to choose symbols T and \Gamma in the formula and define them in the line after the formula.
Also -- please add E_1 and E_2 to the integral, i.e. mention that the integration range is the energy band passed by the caller.

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

Let’s call it for_integral_flux and discuss later what is the best!!
Le 17 avr. 2016 à 22:29, Christoph Deil notifications@github.com a écrit :

In gammapy/scripts/image_pipe.py:

  •    self.maps["total_excess"].data = excess
    
  • def make_exposure(self, obs_id, spectral_index=2.3, exposure_weighted_spectrum=False):
  •    """Compute the exposure map for one observation.
    
  •    Excess/exposure will give the differential flux at the energy Eref at the middle of the self.energy_band
    
  •    If exposure_weighted_spectrum is true, it will give the integrated flux over the self.energy_band
    
  •    Exposure is define as fallow:
    
  •    .. math ::
    
  •        EXPOSURE = \int A(E) \phi(E) \* time_{obs} \, dE
    
  •    assuming a power law for the flux :math:`\phi(E) = \phi_{Eref} \times \frac{E}{E_ref}^{spectral_index}`
    
    I currently get this:

I think it's better to choose symbols T and \Gamma in the formula and define them in the line after the formula.
Also -- please add E_1 and E_2 to the integral, i.e. mention that the integration range is the energy band passed by the caller.


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

OK.

PS: If I find the time, I will use this to compute a HGPS exposure and flux image and compare to the one we computed with HAP and will release with the HGPS paper. This is why I'm excited about the additions you are making like this one ... very soon we will be able to reproduce the HGPS maps and catalog with Gammapy, starting just with event lists and IRFs.

def make_exposure(self, obs_id, spectral_index=2.3, exposure_weighted_spectrum=False):
"""Compute the exposure map for one observation.
Excess/exposure will give the differential flux at the energy Eref at the middle of the self.energy_band

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

In HAP-HD and HGPS, we always use energy_ref = 1 TeV and flux_ref = 1 cm^-2 s^-1 TeV^-1.
Can you please do the same here?
I think that will be simpler / less error-prone for users than to compute energy_ref inside this method.

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

No because otherwise whatever the energy_band you are calculating the same flux when you do total_excess/Exposure whereas you want the differential flux at the center of the self.energy range of the map not at 1 TeV.
Le 17 avr. 2016 à 21:48, Christoph Deil notifications@github.com a écrit :

In gammapy/scripts/image_pipe.py:

 def make_psf(self):
     log.info('Making PSF ...')
  • def make_exposure(self):
  •    log.info('Making EXPOSURE ...')
    
  • def make_total_excess(self):
  •    """Compute excess between counts and bkg map."""
    
  •    self.maps["total_excess"] = SkyMap.empty_like(self.empty_image)
    
  •    excess = self.maps["total_counts"].data - self.maps["total_bkg"].data
    
  •    self.maps["total_excess"].data = excess
    
  • def make_exposure(self, obs_id, spectral_index=2.3, exposure_weighted_spectrum=False):
  •    """Compute the exposure map for one observation.
    
  •    Excess/exposure will give the differential flux at the energy Eref at the middle of the self.energy_band
    
    In HAP-HD and HGPS, we always use energy_ref = 1 TeV and flux_ref = 1 cm^-2 s^-1 TeV^-1.
    Can you please do the same here?
    I think that will be simpler / less error-prone for users than to compute energy_ref inside this method.


You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

I think it's simpler and easier to understand to always use energy_ref = 1 TeV and flux_ref = 1 cm^-2 s^-1 TeV^-1 here. So the user gets back exposure, uses flux = excess / exposure to compute the differential flux at 1 TeV, and then they can compute a differential flux at any energy they like, or an integral flux in any energy band they like.
There are cases where users don't want the integral flux in the measurement energy band, but e.g. above 1 TeV (that's what we do for the flux maps we release for HGPS).

But really it's a matter of taste. If the caller wants the differential flux and the energy band mid-point, then yes, returning exactly that is convenient.

So ... I'll let you make the call here (for today ... will re-discuss with Axel and others).

One more point: the exposure_weighted_spectrum option name isn't very good, no?
Also if it's false, the exposure computed is a power-law spectrum-weighted value, right?

So either rename to for_integral_flux or something with a better name.
Or -- if you follow my suggestion from above to always return diff flux at 1 TeV and let the caller choose a different reference energy or integral energy band, then this option can be removed.

Lea Jouvin added some commits Apr 17, 2016

Lea Jouvin
Lea Jouvin
ok
Lea Jouvin
events = obs.events.select_energy(self.energy_band)
events = events.select_offset(self.offset_band)
if len(events) <= self.ncounts_min:
log.warning('Skipping obs ' + str(obs_id) + ' because it only has ' + str(self.ncounts_min) + ' counts.')

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

self.ncounts_min -> len(events)

@@ -68,12 +85,42 @@ def make_counts(self, obs_id):
counts_map.data = counts_map.data.value
self.maps["counts"] = counts_map
def _low_counts(self, obs_id):

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Based on the name, I would have expected _low_counts to return True if the counts are low.
It's currently the opposite.

Can you change True <-> False ?
Or, if you don't want that, make the name clearer?

def make_total_counts(self):
"""Stack the total count from the observation in the 'DataStore'
Parameters

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Remove this Parameter section (that parameter doesn't exist any more).

If for_integral_flux is true, it will give the integrated flux over the self.energy_band
Exposure is define as fallow:

This comment has been minimized.

@cdeil

cdeil Apr 17, 2016

Member

Typos. Should be: Exposure is defined as follows:

@cdeil

This comment has been minimized.

Member

cdeil commented Apr 17, 2016

I'm offline now.
Please leave another comment when this is ready and I'll merge tomorrow.

Lea Jouvin
self.maps["total_excess"].data = excess
def make_exposure(self, obs_id, spectral_index=2.3, for_integral_flux=False):
r"""Compute the exposure map for one observation.

This comment has been minimized.

@JouvinLea

JouvinLea Apr 17, 2016

Contributor

the r doesn't solve the problem for the docstring...

@@ -3,7 +3,11 @@
unicode_literals)
import logging
import numpy as np
from scipy import interpolate

This comment has been minimized.

@cdeil

cdeil Apr 18, 2016

Member

Imports of optional dependencies like scipy need to be put in the function or method where it's used.
Build fail: https://travis-ci.org/gammapy/gammapy/jobs/123775918#L1428

@cdeil cdeil merged commit bf9efe4 into gammapy:master Apr 18, 2016

1 of 2 checks passed

continuous-integration/travis-ci/pr The Travis CI build failed
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
@cdeil

This comment has been minimized.

Member

cdeil commented Apr 18, 2016

I fixed the scipy import and the exposure docstring in 5637902 and merged this PR.
The issue was extra spaces between :math: and the following backtick.

Thanks!

@cdeil cdeil changed the title from Exposure image computation to Add exposure image computation Apr 20, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment