Bug #2121: ctbutterfly does not produce comparable result for GAUSSIAN and ENVELOPE method - ctools - CTA IRAP Project Gateway

Bug #2121

ctbutterfly does not produce comparable result for GAUSSIAN and ENVELOPE method

Added by Knödlseder Jürgen over 7 years ago. Updated over 7 years ago.

Status:

Closed

Start date:

05/31/2017

Priority:

Urgent

Due date:

Assigned To:

Knödlseder Jürgen

% Done:

100%

Category:

Target version:

1.3.0

Duration:

Description

While writing the tutorial for the first Data Challenge I recognised that ctbutterfly does not produce comparable results for the GAUSSIAN and ENVELOPE methods. The model I used included two point sources and a background model component:

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<source_library title="source library">
  <source name="Src001" type="PointSource">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor" value="7.1484636662325" error="0.0952385139246196" scale="5.7e-18" min="0" free="1" />
      <parameter name="Index" value="0.963447086926589" error="0.00391412506678486" scale="-2.48" min="-4.03225806451613" max="4.03225806451613" free="1" />
      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
    </spectrum>
    <spatialModel type="PointSource">
      <parameter name="RA" value="266.421783219514" error="0.000537203925055351" scale="1" free="1" />
      <parameter name="DEC" value="-29.0041386363916" error="0.000470522584227518" scale="1" free="1" />
    </spatialModel>
  </source>
  <source name="Src002" type="PointSource">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor" value="4.12717772057087" error="0.0803499133359769" scale="5.7e-18" min="0" free="1" />
      <parameter name="Index" value="1.01696465985064" error="0.00607963505525929" scale="-2.48" min="-4.03225806451613" max="4.03225806451613" free="1" />
      <parameter name="PivotEnergy" value="1" scale="300000" free="0" />
    </spectrum>
    <spatialModel type="PointSource">
      <parameter name="RA" value="266.836181272963" error="0.000856619928781573" scale="1" free="1" />
      <parameter name="DEC" value="-28.156245095923" error="0.000751283643100217" scale="1" free="1" />
    </spatialModel>
  </source>
  <source name="Background" type="CTAIrfBackground">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor" value="1.17737659577972" error="0.00105486061616639" scale="1" min="0" free="1" />
      <parameter name="Index" value="0.047504495184372" error="0.000542458489367568" scale="1" min="-10" max="10" free="1" />
      <parameter name="PivotEnergy" value="1" scale="1000000" free="0" />
    </spectrum>
  </source>
</source_library>

When using the ENVELOPE method I obtain the following result for Src001:

103514.216667934 5.17912536878112e-16 5.16850866431869e-16 5.18974207324354e-16
110917.48152624 4.39112662303548e-16 4.38212523954033e-16 4.40012800653062e-16
118850.222743702 3.72302109845796e-16 3.71538926646024e-16 3.73065293045569e-16
127350.308101666 3.15656716134079e-16 3.15009650494971e-16 3.16303781773188e-16
136458.313658893 2.67629862430327e-16 2.67081247181144e-16 2.68178477679509e-16
146217.717445672 2.26910246490849e-16 2.26445102503215e-16 2.27375390478482e-16
156675.107010815 1.92386079397032e-16 1.91991706601974e-16 1.9278045219209e-16
...

Using the GAUSSIAN method I obtained the following result for Src001:

103514.216667934 5.17912536878112e-16 -1.20383747637499 1.203837476375
110917.48152624 4.39112662303548e-16 -1.02067480817947 1.02067480817947
118850.222743702 3.72302109845796e-16 -0.86538015678761 0.865380156787611
127350.308101666 3.15656716134079e-16 -0.733713431310699 0.733713431310699
136458.313658893 2.67629862430327e-16 -0.622079666448653 0.622079666448654
146217.717445672 2.26910246490849e-16 -0.527430867276841 0.527430867276841
156675.107010815 1.92386079397032e-16 -0.447182788250421 0.447182788250421
...

Obviously, the minimum and maximum intensities (columns 3 and 4) are way to large.

butterfly-envelope.png (58.4 KB) Knödlseder Jürgen, 05/31/2017 02:33 AM

butterfly-gaussian.png (55.9 KB) Knödlseder Jürgen, 05/31/2017 02:34 AM

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen over 7 years ago

Not sure why the code was as follows:

        // Get covariance matrix.
        GObservations::likelihood likelihood = m_obs.function();
        if (m_method == "GAUSSIAN") {
            m_covariance = likelihood.curvature()->invert();
        }
        else {   
            m_covariance = likelihood.covariance();
        }

Changing to

        // Get covariance matrix.
        GObservations::likelihood likelihood = m_obs.function();
        if (m_method == "GAUSSIAN") {
            //m_covariance = likelihood.curvature()->invert();
            m_covariance = likelihood.covariance();
        }
        else {   
            m_covariance = likelihood.covariance();
        }

resulted in the following butterfly, much closer to the ENVELOPE method:

103514.216667934 5.17912536878112e-16 5.07371907944126e-16 5.28453165812097e-16
110917.48152624 4.39112662303548e-16 4.30415229157282e-16 4.47810095449814e-16
118850.222743702 3.72302109845796e-16 3.65128153593063e-16 3.7947606609853e-16
127350.308101666 3.15656716134079e-16 3.09741364334773e-16 3.21572067933386e-16
136458.313658893 2.67629862430327e-16 2.62753780439808e-16 2.72505944420846e-16
146217.717445672 2.26910246490849e-16 2.22891893953578e-16 2.30928599028119e-16
156675.107010815 1.92386079397032e-16 1.89075270886635e-16 1.9569688790743e-16

The difference is that in the likelihood.covariance() method the parameter scale factors are multiplied in the covariance matrix while in the likelihood.curvature()->invert() method they are not multiplied in.

#2 Updated by Knödlseder Jürgen over 7 years ago

File butterfly-envelope.png added
File butterfly-gaussian.png added

It turned out that the lines in the code were inverted, the correct code should be

       GObservations::likelihood likelihood = m_obs.function();
        if (m_method == "GAUSSIAN") {
            m_covariance = likelihood.covariance();
        }
        else {   
            m_covariance = likelihood.curvature()->invert();
        }

The ENVELOPE method gives

103514.216667934 5.17912536878112e-16 5.0207116627315e-16 5.34078772239613e-16
110917.48152624 4.39112662303548e-16 4.26032532260246e-16 4.52443686789782e-16
118850.222743702 3.72302109845796e-16 3.61505849317682e-16 3.83291100225765e-16
127350.308101666 3.15656716134079e-16 3.06748353256624e-16 3.24711933595808e-16
136458.313658893 2.67629862430327e-16 2.60281303845876e-16 2.75089109048519e-16
146217.717445672 2.26910246490849e-16 2.2085011527933e-16 2.33053258032221e-16
156675.107010815 1.92386079397032e-16 1.87389509204483e-16 1.97444053530927e-16

while the GAUSSIAN method gives

103514.216667934 5.17912536878112e-16 5.07371907944126e-16 5.28453165812097e-16
110917.48152624 4.39112662303548e-16 4.30415229157282e-16 4.47810095449814e-16
118850.222743702 3.72302109845796e-16 3.65128153593063e-16 3.7947606609853e-16
127350.308101666 3.15656716134079e-16 3.09741364334773e-16 3.21572067933386e-16
136458.313658893 2.67629862430327e-16 2.62753780439808e-16 2.72505944420846e-16
146217.717445672 2.26910246490849e-16 2.22891893953578e-16 2.30928599028119e-16
156675.107010815 1.92386079397032e-16 1.89075270886635e-16 1.9569688790743e-16

which are very close.

Graphically:

ENVELOPE	GAUSSIAN

#3 Updated by Knödlseder Jürgen over 7 years ago

Status changed from New to Pull request
Assigned To set to Knödlseder Jürgen
% Done changed from 0 to 100

#4 Updated by Knödlseder Jürgen over 7 years ago

Status changed from Pull request to Closed

Also available in: Atom PDF

ctools

Issues

Tags