We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
Problem of fit with sherpa
Hi @cdeil and @joleroi , 
I've got some problem to do a spectral fit with sherpa on CTA simulated data. I don't really understand why it's happening. Everything worked fine ~2 weeks ago.

Here is a small example to illustrate the problem: 
```
import astropy.units as u

from gammapy.spectrum.models import PowerLaw
from gammapy.scripts.cta_utils import CTAObservationSimulation, Target, ObservationParameters

from gammapy.scripts import CTAPerf

# observation parameters
alpha = 0.2 * u.Unit('')
livetime = 100 * u.h
emin = 0.05 * u.TeV
emax = 5. * u.TeV
obs_param = ObservationParameters(alpha=alpha,
                                  livetime=livetime,
                                  emin=emin,
                                  emax=emax)

# target
pwl = PowerLaw(amplitude=1e-11 * u.Unit('1/(s cm2 TeV)'),
               index=2.,
               reference=1. * u.TeV)
target = Target(name='bob', model=pwl)

# perf
dir_south = './CTA-Performance-South-20170323/'
perf = CTAPerf.read('$GAMMAPY_EXTRA/datasets/cta/perf_prod2/point_like_non_smoothed/South_5h.fits.gz')

# simu
obs_id = 0
simu = CTAObservationSimulation.simulate_obs(perf=perf,
                                             target=target,
                                             obs_param=obs_param,
                                             obs_id=obs_id)

# simu output
print(simu)
stats = simu.stats_table()
print(stats[['energy_min', 'energy_max', 'background', 'excess', 'sigma']])
simu.write(outdir='./', use_sherpa=True)


from sherpa.astro.datastack import DataStack
import sherpa.astro.datastack as sh
from sherpa_model import Absorption

simu_output_dir = './'
filename = simu_output_dir + 'pha_obs' + str(obs_id) + '.fits'
sh.load_data(id=obs_id, filename=filename)
sherpa_model = 'powlaw1d.pwl'
sh.set_source(obs_id, sherpa_model)
    
pwl.ampl = 1.e-12 * u.Unit('1/(s cm2 TeV)').to('1/(s cm2 keV)')
pwl.ref = 1 * u.TeV.to('keV')
pwl.gamma = 2.

sh.notice(emin.to('keV')*0.99,emax.to('keV')*1.01)
print('JLK: threshold for spectra: ({},{})'.format(emin.to('keV')*0.99,emax.to('keV')*1.01))

sh.set_stat('WStat')
sh.fit()
sh.conf()
```

which results in weird NaN values and implies problem for the fit:
```
Dataset               = 0
Method                = levmar
Statistic             = wstat
Initial fit statistic = nan
Final fit statistic   = nan at function evaluation 33
Data points           = 12
Degrees of freedom    = 10
Reduced statistic     = nan
Change in statistic   = nan
   pwl.gamma      2           
   pwl.ampl       1e-21       
```
 
Do you have any idea what could cause that? OGIP format writting or something like that?

Thanks in advance ++

--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `gammapy/scripts/cta_utils.py`
Content:
```
1 from __future__ import absolute_import, division, print_function, unicode_literals
2 import numpy as np
3 import astropy.units as u
4 from ..spectrum import SpectrumObservation
5 from ..spectrum.utils import CountsPredictor
6 from ..spectrum.core import PHACountsSpectrum
7 from ..utils.random import get_random_state
8 
9 __all__ = [
10     'Target',
11     'ObservationParameters',
12     'CTAObservationSimulation',
13 ]
14 
15 
16 class Target(object):
17     """Observation target information.
18 
19     Parameters
20     ----------
21     name : `str`
22         Name of the source
23     model : `~gammapy.spectrum.models.SpectralModel`
24         Model of the source
25     """
26 
27     def __init__(self, name=None,
28                  model=None):
29         self.name = name
30         self.model = model
31 
32     def __str__(self):
33         """Target report (`str`)."""
34         ss = '*** Target parameters ***\n'
35         ss += 'Name={}\n'.format(self.name)
36         for par in self.model.parameters.parameters:
37             ss += '{}={} {}\n'.format(par.name, str(par.value), par.unit)
38         return ss
39 
40     def from_fermi_lat_catalogue(name):
41         raise NotImplementedError
42 
43 
44 class ObservationParameters(object):
45     """Container for observation parameters.
46 
47     Parameters
48     ----------
49     alpha : `~astropy.units.Quantity`
50         Normalisation between ON and OFF regions
51     livetime :  `~astropy.units.Quantity`
52         Observation time
53     emin :  `~astropy.units.Quantity`
54         Minimal energy for simulation
55     emax : `~astropy.units.Quantity`
56         Maximal energy for simulation
57     """
58 
59     def __init__(self, alpha=None, livetime=None,
60                  emin=None, emax=None):
61         self.alpha = alpha
62         self.livetime = livetime
63         self.emin = emin
64         self.emax = emax
65 
66     def __str__(self):
67         """Observation summary report (`str`)."""
68         ss = '*** Observation parameters summary ***\n'
69         ss += 'alpha={} [{}]\n'.format(self.alpha.value, self.alpha.unit)
70         ss += 'livetime={} [{}]\n'.format(self.livetime.value,
71                                           self.livetime.unit)
72         ss += 'emin={} [{}]\n'.format(self.emin.value, self.emin.unit)
73         ss += 'emax={} [{}]\n'.format(self.emax.value, self.emax.unit)
74         return ss
75 
76 
77 class CTAObservationSimulation(object):
78     """Simulate observation for one IRF and target.
79 
80     TODO : Should be merge with `~gammapy.spectrum.SpectrumSimulation`
81 
82     Parameters
83     ----------
84     perf : `~gammapy.scripts.CTAPerf`
85         CTA performance
86     target : `~gammapy.scripts.Target`
87         Source
88     """
89 
90     @staticmethod
91     def simulate_obs(perf, target, obs_param, obs_id=0):
92         """
93         Simulate observation with given parameters
94 
95         Parameters
96         ----------
97         perf : `~gammapy.scripts.CTAPerf`
98             CTA performance
99         target : `~gammapy.scripts.Target`
100             Source
101         obs_param : `~gammapy.scripts.ObservationParameters`
102             Observation parameters
103         obs_id : `int`, optional
104             Observation Id
105         """
106         livetime = obs_param.livetime
107         alpha = obs_param.alpha.value
108         emin = obs_param.emin
109         emax = obs_param.emax
110 
111         model = target.model
112 
113         # Compute expected counts
114         reco_energy = perf.bkg.energy
115         bkg_rate_values = perf.bkg.data.data * livetime.to('s')
116         predicted_counts = CountsPredictor(model=model,
117                                            aeff=perf.aeff,
118                                            livetime=livetime,
119                                            edisp=perf.rmf)
120         predicted_counts.run()
121         npred = predicted_counts.npred
122         # set negative values to zero (interpolation issue)
123         idx = np.where(npred.data.data < 0.)
124         npred.data.data[idx] = 0
125 
126         # Randomise counts
127         rand = get_random_state('random-seed')
128         on_counts = rand.poisson(npred.data.data.value)  # excess
129         bkg_counts = rand.poisson(bkg_rate_values.value)  # bkg in ON region
130         off_counts = rand.poisson(
131             bkg_rate_values.value / alpha)  # bkg in OFF region
132 
133         on_counts += bkg_counts  # evts in ON region
134 
135         meta = dict(EXPOSURE=livetime.to('s').value)
136 
137         on_vector = PHACountsSpectrum(
138             data=on_counts,
139             backscal=1,
140             energy_lo=reco_energy.lo,
141             energy_hi=reco_energy.hi,
142             meta=meta,
143         )
144 
145         off_vector = PHACountsSpectrum(energy_lo=reco_energy.lo,
146                                        energy_hi=reco_energy.hi,
147                                        data=off_counts,
148                                        backscal=1. / alpha,
149                                        is_bkg=True,
150                                        )
151 
152         obs = SpectrumObservation(on_vector=on_vector,
153                                   off_vector=off_vector,
154                                   aeff=perf.aeff,
155                                   edisp=perf.rmf)
156         obs.obs_id = obs_id
157 
158         # Set threshold according to the closest energy reco from bkg bins
159         idx_min = np.abs(reco_energy.lo - emin).argmin()
160         idx_max = np.abs(reco_energy.lo - emax).argmin()
161         obs.lo_threshold = reco_energy.lo[idx_min]
162         obs.hi_threshold = reco_energy.lo[idx_max]
163 
164         return obs
165 
166     @staticmethod
167     def plot_simu(simu, target):
168         import matplotlib.pyplot as plt
169         fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,
170                                        figsize=(10, 5))
171 
172         # Spectrum plot
173         energy_range = [0.01 * u.TeV, 100 * u.TeV]
174         target.model.plot(ax=ax1, energy_range=energy_range,
175                           label='Model')
176         plt.text(0.55, 0.65, target.__str__(),
177                  style='italic', transform=ax1.transAxes, fontsize=7,
178                  bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10})
179         ax1.set_xlim([energy_range[0].value, energy_range[1].value])
180         ax1.set_ylim(1.e-17, 1.e-5)
181         ax1.grid(which='both')
182         ax1.legend(loc=0)
183 
184         # Counts plot
185         on_off = simu.on_vector.data.data.value
186         off = 1. / simu.off_vector.backscal * simu.off_vector.data.data.value
187         excess = on_off - off
188         bins = simu.on_vector.energy.lo.value
189         x = simu.on_vector.energy.nodes.value
190         ax2.hist(x, bins=bins, weights=on_off,
191                  facecolor='blue', alpha=1, label='ON')
192         ax2.hist(x, bins=bins, weights=off,
193                  facecolor='green', alpha=1, label='OFF')
194         ax2.hist(x, bins=bins, weights=excess,
195                  facecolor='red', alpha=1, label='EXCESS')
196         ax2.legend(loc='best')
197         ax2.set_xscale('log')
198         ax2.set_xlabel('Energy [TeV]')
199         ax2.set_ylabel('Expected counts')
200         ax2.set_xlim([energy_range[0].value, energy_range[1].value])
201         ax2.set_ylim([0.0001, on_off.max() * (1 + 0.05)])
202         ax2.vlines(simu.lo_threshold.value, 0, 1.1 * on_off.max(),
203                    linestyles='dashed')
204         ax2.grid(which='both')
205         plt.text(0.55, 0.05, simu.__str__(),
206                  style='italic', transform=ax2.transAxes, fontsize=7,
207                  bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10})
208         plt.tight_layout()
209
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/gammapy/scripts/cta_utils.py b/gammapy/scripts/cta_utils.py
--- a/gammapy/scripts/cta_utils.py
+++ b/gammapy/scripts/cta_utils.py
@@ -1,4 +1,4 @@
-from __future__ import absolute_import, division, print_function, unicode_literals
+
 import numpy as np
 import astropy.units as u
 from ..spectrum import SpectrumObservation
@@ -132,22 +132,21 @@
 
         on_counts += bkg_counts  # evts in ON region
 
-        meta = dict(EXPOSURE=livetime.to('s').value)
-
         on_vector = PHACountsSpectrum(
             data=on_counts,
             backscal=1,
             energy_lo=reco_energy.lo,
             energy_hi=reco_energy.hi,
-            meta=meta,
         )
 
+        on_vector.livetime = livetime
         off_vector = PHACountsSpectrum(energy_lo=reco_energy.lo,
                                        energy_hi=reco_energy.hi,
                                        data=off_counts,
                                        backscal=1. / alpha,
                                        is_bkg=True,
                                        )
+        off_vector.livetime = livetime
 
         obs = SpectrumObservation(on_vector=on_vector,
                                   off_vector=off_vector,

{"golden_diff": "diff --git a/gammapy/scripts/cta_utils.py b/gammapy/scripts/cta_utils.py\n--- a/gammapy/scripts/cta_utils.py\n+++ b/gammapy/scripts/cta_utils.py\n@@ -1,4 +1,4 @@\n-from __future__ import absolute_import, division, print_function, unicode_literals\n+\n import numpy as np\n import astropy.units as u\n from ..spectrum import SpectrumObservation\n@@ -132,22 +132,21 @@\n \n         on_counts += bkg_counts  # evts in ON region\n \n-        meta = dict(EXPOSURE=livetime.to('s').value)\n-\n         on_vector = PHACountsSpectrum(\n             data=on_counts,\n             backscal=1,\n             energy_lo=reco_energy.lo,\n             energy_hi=reco_energy.hi,\n-            meta=meta,\n         )\n \n+        on_vector.livetime = livetime\n         off_vector = PHACountsSpectrum(energy_lo=reco_energy.lo,\n                                        energy_hi=reco_energy.hi,\n                                        data=off_counts,\n                                        backscal=1. / alpha,\n                                        is_bkg=True,\n                                        )\n+        off_vector.livetime = livetime\n \n         obs = SpectrumObservation(on_vector=on_vector,\n                                   off_vector=off_vector,\n", "issue": "Problem of fit with sherpa\nHi @cdeil and @joleroi , \r\nI've got some problem to do a spectral fit with sherpa on CTA simulated data. I don't really understand why it's happening. Everything worked fine ~2 weeks ago.\r\n\r\nHere is a small example to illustrate the problem: \r\n```\r\nimport astropy.units as u\r\n\r\nfrom gammapy.spectrum.models import PowerLaw\r\nfrom gammapy.scripts.cta_utils import CTAObservationSimulation, Target, ObservationParameters\r\n\r\nfrom gammapy.scripts import CTAPerf\r\n\r\n# observation parameters\r\nalpha = 0.2 * u.Unit('')\r\nlivetime = 100 * u.h\r\nemin = 0.05 * u.TeV\r\nemax = 5. * u.TeV\r\nobs_param = ObservationParameters(alpha=alpha,\r\n                                  livetime=livetime,\r\n                                  emin=emin,\r\n                                  emax=emax)\r\n\r\n# target\r\npwl = PowerLaw(amplitude=1e-11 * u.Unit('1/(s cm2 TeV)'),\r\n               index=2.,\r\n               reference=1. * u.TeV)\r\ntarget = Target(name='bob', model=pwl)\r\n\r\n# perf\r\ndir_south = './CTA-Performance-South-20170323/'\r\nperf = CTAPerf.read('$GAMMAPY_EXTRA/datasets/cta/perf_prod2/point_like_non_smoothed/South_5h.fits.gz')\r\n\r\n# simu\r\nobs_id = 0\r\nsimu = CTAObservationSimulation.simulate_obs(perf=perf,\r\n                                             target=target,\r\n                                             obs_param=obs_param,\r\n                                             obs_id=obs_id)\r\n\r\n# simu output\r\nprint(simu)\r\nstats = simu.stats_table()\r\nprint(stats[['energy_min', 'energy_max', 'background', 'excess', 'sigma']])\r\nsimu.write(outdir='./', use_sherpa=True)\r\n\r\n\r\nfrom sherpa.astro.datastack import DataStack\r\nimport sherpa.astro.datastack as sh\r\nfrom sherpa_model import Absorption\r\n\r\nsimu_output_dir = './'\r\nfilename = simu_output_dir + 'pha_obs' + str(obs_id) + '.fits'\r\nsh.load_data(id=obs_id, filename=filename)\r\nsherpa_model = 'powlaw1d.pwl'\r\nsh.set_source(obs_id, sherpa_model)\r\n    \r\npwl.ampl = 1.e-12 * u.Unit('1/(s cm2 TeV)').to('1/(s cm2 keV)')\r\npwl.ref = 1 * u.TeV.to('keV')\r\npwl.gamma = 2.\r\n\r\nsh.notice(emin.to('keV')*0.99,emax.to('keV')*1.01)\r\nprint('JLK: threshold for spectra: ({},{})'.format(emin.to('keV')*0.99,emax.to('keV')*1.01))\r\n\r\nsh.set_stat('WStat')\r\nsh.fit()\r\nsh.conf()\r\n```\r\n\r\nwhich results in weird NaN values and implies problem for the fit:\r\n```\r\nDataset               = 0\r\nMethod                = levmar\r\nStatistic             = wstat\r\nInitial fit statistic = nan\r\nFinal fit statistic   = nan at function evaluation 33\r\nData points           = 12\r\nDegrees of freedom    = 10\r\nReduced statistic     = nan\r\nChange in statistic   = nan\r\n   pwl.gamma      2           \r\n   pwl.ampl       1e-21       \r\n```\r\n \r\nDo you have any idea what could cause that? OGIP format writting or something like that?\r\n\r\nThanks in advance ++\n", "before_files": [{"content": "from __future__ import absolute_import, division, print_function, unicode_literals\nimport numpy as np\nimport astropy.units as u\nfrom ..spectrum import SpectrumObservation\nfrom ..spectrum.utils import CountsPredictor\nfrom ..spectrum.core import PHACountsSpectrum\nfrom ..utils.random import get_random_state\n\n__all__ = [\n    'Target',\n    'ObservationParameters',\n    'CTAObservationSimulation',\n]\n\n\nclass Target(object):\n    \"\"\"Observation target information.\n\n    Parameters\n    ----------\n    name : `str`\n        Name of the source\n    model : `~gammapy.spectrum.models.SpectralModel`\n        Model of the source\n    \"\"\"\n\n    def __init__(self, name=None,\n                 model=None):\n        self.name = name\n        self.model = model\n\n    def __str__(self):\n        \"\"\"Target report (`str`).\"\"\"\n        ss = '*** Target parameters ***\\n'\n        ss += 'Name={}\\n'.format(self.name)\n        for par in self.model.parameters.parameters:\n            ss += '{}={} {}\\n'.format(par.name, str(par.value), par.unit)\n        return ss\n\n    def from_fermi_lat_catalogue(name):\n        raise NotImplementedError\n\n\nclass ObservationParameters(object):\n    \"\"\"Container for observation parameters.\n\n    Parameters\n    ----------\n    alpha : `~astropy.units.Quantity`\n        Normalisation between ON and OFF regions\n    livetime :  `~astropy.units.Quantity`\n        Observation time\n    emin :  `~astropy.units.Quantity`\n        Minimal energy for simulation\n    emax : `~astropy.units.Quantity`\n        Maximal energy for simulation\n    \"\"\"\n\n    def __init__(self, alpha=None, livetime=None,\n                 emin=None, emax=None):\n        self.alpha = alpha\n        self.livetime = livetime\n        self.emin = emin\n        self.emax = emax\n\n    def __str__(self):\n        \"\"\"Observation summary report (`str`).\"\"\"\n        ss = '*** Observation parameters summary ***\\n'\n        ss += 'alpha={} [{}]\\n'.format(self.alpha.value, self.alpha.unit)\n        ss += 'livetime={} [{}]\\n'.format(self.livetime.value,\n                                          self.livetime.unit)\n        ss += 'emin={} [{}]\\n'.format(self.emin.value, self.emin.unit)\n        ss += 'emax={} [{}]\\n'.format(self.emax.value, self.emax.unit)\n        return ss\n\n\nclass CTAObservationSimulation(object):\n    \"\"\"Simulate observation for one IRF and target.\n\n    TODO : Should be merge with `~gammapy.spectrum.SpectrumSimulation`\n\n    Parameters\n    ----------\n    perf : `~gammapy.scripts.CTAPerf`\n        CTA performance\n    target : `~gammapy.scripts.Target`\n        Source\n    \"\"\"\n\n    @staticmethod\n    def simulate_obs(perf, target, obs_param, obs_id=0):\n        \"\"\"\n        Simulate observation with given parameters\n\n        Parameters\n        ----------\n        perf : `~gammapy.scripts.CTAPerf`\n            CTA performance\n        target : `~gammapy.scripts.Target`\n            Source\n        obs_param : `~gammapy.scripts.ObservationParameters`\n            Observation parameters\n        obs_id : `int`, optional\n            Observation Id\n        \"\"\"\n        livetime = obs_param.livetime\n        alpha = obs_param.alpha.value\n        emin = obs_param.emin\n        emax = obs_param.emax\n\n        model = target.model\n\n        # Compute expected counts\n        reco_energy = perf.bkg.energy\n        bkg_rate_values = perf.bkg.data.data * livetime.to('s')\n        predicted_counts = CountsPredictor(model=model,\n                                           aeff=perf.aeff,\n                                           livetime=livetime,\n                                           edisp=perf.rmf)\n        predicted_counts.run()\n        npred = predicted_counts.npred\n        # set negative values to zero (interpolation issue)\n        idx = np.where(npred.data.data < 0.)\n        npred.data.data[idx] = 0\n\n        # Randomise counts\n        rand = get_random_state('random-seed')\n        on_counts = rand.poisson(npred.data.data.value)  # excess\n        bkg_counts = rand.poisson(bkg_rate_values.value)  # bkg in ON region\n        off_counts = rand.poisson(\n            bkg_rate_values.value / alpha)  # bkg in OFF region\n\n        on_counts += bkg_counts  # evts in ON region\n\n        meta = dict(EXPOSURE=livetime.to('s').value)\n\n        on_vector = PHACountsSpectrum(\n            data=on_counts,\n            backscal=1,\n            energy_lo=reco_energy.lo,\n            energy_hi=reco_energy.hi,\n            meta=meta,\n        )\n\n        off_vector = PHACountsSpectrum(energy_lo=reco_energy.lo,\n                                       energy_hi=reco_energy.hi,\n                                       data=off_counts,\n                                       backscal=1. / alpha,\n                                       is_bkg=True,\n                                       )\n\n        obs = SpectrumObservation(on_vector=on_vector,\n                                  off_vector=off_vector,\n                                  aeff=perf.aeff,\n                                  edisp=perf.rmf)\n        obs.obs_id = obs_id\n\n        # Set threshold according to the closest energy reco from bkg bins\n        idx_min = np.abs(reco_energy.lo - emin).argmin()\n        idx_max = np.abs(reco_energy.lo - emax).argmin()\n        obs.lo_threshold = reco_energy.lo[idx_min]\n        obs.hi_threshold = reco_energy.lo[idx_max]\n\n        return obs\n\n    @staticmethod\n    def plot_simu(simu, target):\n        import matplotlib.pyplot as plt\n        fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,\n                                       figsize=(10, 5))\n\n        # Spectrum plot\n        energy_range = [0.01 * u.TeV, 100 * u.TeV]\n        target.model.plot(ax=ax1, energy_range=energy_range,\n                          label='Model')\n        plt.text(0.55, 0.65, target.__str__(),\n                 style='italic', transform=ax1.transAxes, fontsize=7,\n                 bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10})\n        ax1.set_xlim([energy_range[0].value, energy_range[1].value])\n        ax1.set_ylim(1.e-17, 1.e-5)\n        ax1.grid(which='both')\n        ax1.legend(loc=0)\n\n        # Counts plot\n        on_off = simu.on_vector.data.data.value\n        off = 1. / simu.off_vector.backscal * simu.off_vector.data.data.value\n        excess = on_off - off\n        bins = simu.on_vector.energy.lo.value\n        x = simu.on_vector.energy.nodes.value\n        ax2.hist(x, bins=bins, weights=on_off,\n                 facecolor='blue', alpha=1, label='ON')\n        ax2.hist(x, bins=bins, weights=off,\n                 facecolor='green', alpha=1, label='OFF')\n        ax2.hist(x, bins=bins, weights=excess,\n                 facecolor='red', alpha=1, label='EXCESS')\n        ax2.legend(loc='best')\n        ax2.set_xscale('log')\n        ax2.set_xlabel('Energy [TeV]')\n        ax2.set_ylabel('Expected counts')\n        ax2.set_xlim([energy_range[0].value, energy_range[1].value])\n        ax2.set_ylim([0.0001, on_off.max() * (1 + 0.05)])\n        ax2.vlines(simu.lo_threshold.value, 0, 1.1 * on_off.max(),\n                   linestyles='dashed')\n        ax2.grid(which='both')\n        plt.text(0.55, 0.05, simu.__str__(),\n                 style='italic', transform=ax2.transAxes, fontsize=7,\n                 bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10})\n        plt.tight_layout()\n", "path": "gammapy/scripts/cta_utils.py"}], "after_files": [{"content": "\nimport numpy as np\nimport astropy.units as u\nfrom ..spectrum import SpectrumObservation\nfrom ..spectrum.utils import CountsPredictor\nfrom ..spectrum.core import PHACountsSpectrum\nfrom ..utils.random import get_random_state\n\n__all__ = [\n    'Target',\n    'ObservationParameters',\n    'CTAObservationSimulation',\n]\n\n\nclass Target(object):\n    \"\"\"Observation target information.\n\n    Parameters\n    ----------\n    name : `str`\n        Name of the source\n    model : `~gammapy.spectrum.models.SpectralModel`\n        Model of the source\n    \"\"\"\n\n    def __init__(self, name=None,\n                 model=None):\n        self.name = name\n        self.model = model\n\n    def __str__(self):\n        \"\"\"Target report (`str`).\"\"\"\n        ss = '*** Target parameters ***\\n'\n        ss += 'Name={}\\n'.format(self.name)\n        for par in self.model.parameters.parameters:\n            ss += '{}={} {}\\n'.format(par.name, str(par.value), par.unit)\n        return ss\n\n    def from_fermi_lat_catalogue(name):\n        raise NotImplementedError\n\n\nclass ObservationParameters(object):\n    \"\"\"Container for observation parameters.\n\n    Parameters\n    ----------\n    alpha : `~astropy.units.Quantity`\n        Normalisation between ON and OFF regions\n    livetime :  `~astropy.units.Quantity`\n        Observation time\n    emin :  `~astropy.units.Quantity`\n        Minimal energy for simulation\n    emax : `~astropy.units.Quantity`\n        Maximal energy for simulation\n    \"\"\"\n\n    def __init__(self, alpha=None, livetime=None,\n                 emin=None, emax=None):\n        self.alpha = alpha\n        self.livetime = livetime\n        self.emin = emin\n        self.emax = emax\n\n    def __str__(self):\n        \"\"\"Observation summary report (`str`).\"\"\"\n        ss = '*** Observation parameters summary ***\\n'\n        ss += 'alpha={} [{}]\\n'.format(self.alpha.value, self.alpha.unit)\n        ss += 'livetime={} [{}]\\n'.format(self.livetime.value,\n                                          self.livetime.unit)\n        ss += 'emin={} [{}]\\n'.format(self.emin.value, self.emin.unit)\n        ss += 'emax={} [{}]\\n'.format(self.emax.value, self.emax.unit)\n        return ss\n\n\nclass CTAObservationSimulation(object):\n    \"\"\"Simulate observation for one IRF and target.\n\n    TODO : Should be merge with `~gammapy.spectrum.SpectrumSimulation`\n\n    Parameters\n    ----------\n    perf : `~gammapy.scripts.CTAPerf`\n        CTA performance\n    target : `~gammapy.scripts.Target`\n        Source\n    \"\"\"\n\n    @staticmethod\n    def simulate_obs(perf, target, obs_param, obs_id=0):\n        \"\"\"\n        Simulate observation with given parameters\n\n        Parameters\n        ----------\n        perf : `~gammapy.scripts.CTAPerf`\n            CTA performance\n        target : `~gammapy.scripts.Target`\n            Source\n        obs_param : `~gammapy.scripts.ObservationParameters`\n            Observation parameters\n        obs_id : `int`, optional\n            Observation Id\n        \"\"\"\n        livetime = obs_param.livetime\n        alpha = obs_param.alpha.value\n        emin = obs_param.emin\n        emax = obs_param.emax\n\n        model = target.model\n\n        # Compute expected counts\n        reco_energy = perf.bkg.energy\n        bkg_rate_values = perf.bkg.data.data * livetime.to('s')\n        predicted_counts = CountsPredictor(model=model,\n                                           aeff=perf.aeff,\n                                           livetime=livetime,\n                                           edisp=perf.rmf)\n        predicted_counts.run()\n        npred = predicted_counts.npred\n        # set negative values to zero (interpolation issue)\n        idx = np.where(npred.data.data < 0.)\n        npred.data.data[idx] = 0\n\n        # Randomise counts\n        rand = get_random_state('random-seed')\n        on_counts = rand.poisson(npred.data.data.value)  # excess\n        bkg_counts = rand.poisson(bkg_rate_values.value)  # bkg in ON region\n        off_counts = rand.poisson(\n            bkg_rate_values.value / alpha)  # bkg in OFF region\n\n        on_counts += bkg_counts  # evts in ON region\n\n        on_vector = PHACountsSpectrum(\n            data=on_counts,\n            backscal=1,\n            energy_lo=reco_energy.lo,\n            energy_hi=reco_energy.hi,\n        )\n\n        on_vector.livetime = livetime\n        off_vector = PHACountsSpectrum(energy_lo=reco_energy.lo,\n                                       energy_hi=reco_energy.hi,\n                                       data=off_counts,\n                                       backscal=1. / alpha,\n                                       is_bkg=True,\n                                       )\n        off_vector.livetime = livetime\n\n        obs = SpectrumObservation(on_vector=on_vector,\n                                  off_vector=off_vector,\n                                  aeff=perf.aeff,\n                                  edisp=perf.rmf)\n        obs.obs_id = obs_id\n\n        # Set threshold according to the closest energy reco from bkg bins\n        idx_min = np.abs(reco_energy.lo - emin).argmin()\n        idx_max = np.abs(reco_energy.lo - emax).argmin()\n        obs.lo_threshold = reco_energy.lo[idx_min]\n        obs.hi_threshold = reco_energy.lo[idx_max]\n\n        return obs\n\n    @staticmethod\n    def plot_simu(simu, target):\n        import matplotlib.pyplot as plt\n        fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,\n                                       figsize=(10, 5))\n\n        # Spectrum plot\n        energy_range = [0.01 * u.TeV, 100 * u.TeV]\n        target.model.plot(ax=ax1, energy_range=energy_range,\n                          label='Model')\n        plt.text(0.55, 0.65, target.__str__(),\n                 style='italic', transform=ax1.transAxes, fontsize=7,\n                 bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10})\n        ax1.set_xlim([energy_range[0].value, energy_range[1].value])\n        ax1.set_ylim(1.e-17, 1.e-5)\n        ax1.grid(which='both')\n        ax1.legend(loc=0)\n\n        # Counts plot\n        on_off = simu.on_vector.data.data.value\n        off = 1. / simu.off_vector.backscal * simu.off_vector.data.data.value\n        excess = on_off - off\n        bins = simu.on_vector.energy.lo.value\n        x = simu.on_vector.energy.nodes.value\n        ax2.hist(x, bins=bins, weights=on_off,\n                 facecolor='blue', alpha=1, label='ON')\n        ax2.hist(x, bins=bins, weights=off,\n                 facecolor='green', alpha=1, label='OFF')\n        ax2.hist(x, bins=bins, weights=excess,\n                 facecolor='red', alpha=1, label='EXCESS')\n        ax2.legend(loc='best')\n        ax2.set_xscale('log')\n        ax2.set_xlabel('Energy [TeV]')\n        ax2.set_ylabel('Expected counts')\n        ax2.set_xlim([energy_range[0].value, energy_range[1].value])\n        ax2.set_ylim([0.0001, on_off.max() * (1 + 0.05)])\n        ax2.vlines(simu.lo_threshold.value, 0, 1.1 * on_off.max(),\n                   linestyles='dashed')\n        ax2.grid(which='both')\n        plt.text(0.55, 0.05, simu.__str__(),\n                 style='italic', transform=ax2.transAxes, fontsize=7,\n                 bbox={'facecolor': 'white', 'alpha': 1, 'pad': 10})\n        plt.tight_layout()\n", "path": "gammapy/scripts/cta_utils.py"}]}