Feature #1948
Add smooth broken power law spectrum model
| Status: | Closed | Start date: | 03/10/2017 | |
|---|---|---|---|---|
| Priority: | Normal | Due date: | ||
| Assigned To: |  Cardenzana Josh | % Done: | 100% | |
| Category: | - | |||
| Target version: | 1.3.0 | |||
| Duration: |
Description
A smooth broken power law spectrum model should be implemented, following the Fermi-LAT syntax (see https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html#SmoothBrokenPowerLaw and https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#SmoothBrokenPowerLaw).
The XML format should be
<spectrum type="SmoothBrokenPowerLaw"> <parameter free="1" max="1e10" min="0.0" name="Prefactor" scale="1e-06" value="1.0"/> <parameter free="1" max="-1.0" min="-5.0" name="Index1" scale="1.0" value="-2.0"/> <parameter free="0" max="2000.0" min="30.0" name="Scale" scale="1.0" value="100.0"/> <parameter free="1" max="-1.0" min="-5.0" name="Index2" scale="1.0" value="-2.0"/> <parameter free="1" max="5e5" min="20" name="BreakValue" scale="1.0" value="1e3"/> <parameter free="1" max="10" min="0.01" name="Beta" scale="1.0" value="0.2"/> </spectrum>
The following alternative format, being more consistent in the namings, should also be supported:
<spectrum type="SmoothBrokenPowerLaw"> <parameter free="1" max="1e10" min="0.0" name="Prefactor" scale="1e-06" value="1.0"/> <parameter free="1" max="-1.0" min="-5.0" name="Index1" scale="1.0" value="-2.0"/> <parameter free="0" max="2000.0" min="30.0" name="PivotEnergy" scale="1.0" value="100.0"/> <parameter free="1" max="-1.0" min="-5.0" name="Index2" scale="1.0" value="-2.0"/> <parameter free="1" max="5e5" min="20" name="BreakEnergy" scale="1.0" value="1e3"/> <parameter free="1" max="10" min="0.01" name="BreakSmoothness" scale="1.0" value="0.2"/> </spectrum>
To implement the model, the best is to start from the GModelSpectralBrokenPlaw class, rename the .hpp, .cpp and .i files to GModelSpectralSmoothBrokenPlaw (see https://cta-redmine.irap.omp.eu/attachments/download/1919/6th-coding-sprint.pdf) and modify the files as needed.
prefactor.png (32.3 KB)
index1.png (35.5 KB)
index2.png (35.5 KB)
breakenergy.png (38 KB)
smoothness.png (40.3 KB)
prefactor_nosmooth.png (40.4 KB)
index1_nosmooth.png (39.2 KB)
index2_nosmooth.png (39.8 KB)
breakenergy_nosmooth.png (37.9 KB)
prefactor_smooth2.0.png (33.8 KB)
index1_smooth2.0.png (31.2 KB)
index2_smooth2.0.png (31.5 KB)
breakenergy_smooth2.0.png (34.7 KB)
prefactor_smooth0.02.png (40.3 KB)
index1_smooth0.02.png (39.1 KB)
index2_smooth0.02.png (39.4 KB)
breakenergy_smooth0.02.png (41.7 KB)
Recurrence
No recurrence.
History
#1
Updated by Knödlseder Jürgen about 7 years ago
- Status changed from New to In Progress
- % Done changed from 0 to 50
#2
Updated by Knödlseder Jürgen about 7 years ago
- File prefactor.png added
- File index1.png added
- File index2.png added
- File breakenergy.png added
- File smoothness.png added
I created 100 pulls for the following model definition XML file:
Below the resulting pull histograms:
<?xml version="1.0" standalone="no"?>
<source_library title="source library">
<source name="Crab" type="PointSource">
<spectrum type="SmoothBrokenPowerLaw">
<parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
<parameter name="Index1" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
<parameter name="PivotEnergy" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
<parameter name="Index2" scale="-1" value="2.70" min="0.01" max="+5.0" free="1"/>
<parameter name="BreakEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="1"/>
<parameter name="BreakSmoothness" scale="1.0" value="0.2" min="0.01" max="10.0" free="1"/>
</spectrum>
<spatialModel type="PointSource">
<parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
<parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
</spatialModel>
</source>
<source name="Background" type="CTAIrfBackground" instrument="CTA">
<spectrum type="PowerLaw">
<parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e3" free="1"/>
<parameter name="Index" scale="1.0" value="0.0" min="-5.0" max="+5.0" free="1"/>
<parameter name="PivotEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
</spectrum>
</source>
</source_library>
Below the resulting pull histograms:
 |
 |
 |
 |
 |
#3
Updated by Knödlseder Jürgen about 7 years ago
- File prefactor_nosmooth.png added
- File index1_nosmooth.png added
- File index2_nosmooth.png added
- File breakenergy_nosmooth.png added
I now tried fixing the smoothness factor. This did change the pull distributions, but did not lead to satisfactory pull histograms:
 |
 |
 |
 |
#4
Updated by Knödlseder Jürgen about 7 years ago
- File prefactor_smooth2.0.png added
- File index1_smooth2.0.png added
- File index2_smooth2.0.png added
- File breakenergy_smooth2.0.png added
I now made the smoothness larger:
This had a dramatic impact on the pull distributions:
<?xml version="1.0" standalone="no"?>
<source_library title="source library">
<source name="Crab" type="PointSource">
<spectrum type="SmoothBrokenPowerLaw">
<parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
<parameter name="Index1" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
<parameter name="PivotEnergy" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
<parameter name="Index2" scale="-1" value="2.70" min="0.01" max="+5.0" free="1"/>
<parameter name="BreakEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="1"/>
<parameter name="BreakSmoothness" scale="1.0" value="2.0" min="0.01" max="10.0" free="0"/>
</spectrum>
<spatialModel type="PointSource">
<parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
<parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
</spatialModel>
</source>
<source name="Background" type="CTAIrfBackground" instrument="CTA">
<spectrum type="PowerLaw">
<parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e3" free="1"/>
<parameter name="Index" scale="1.0" value="0.0" min="-5.0" max="+5.0" free="1"/>
<parameter name="PivotEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
</spectrum>
</source>
</source_library>
This had a dramatic impact on the pull distributions:
 |
 |
 |
 |
#5
Updated by Knödlseder Jürgen about 7 years ago
- File prefactor_smooth0.02.png added
- File index1_smooth0.02.png added
- File index2_smooth0.02.png added
- File breakenergy_smooth0.02.png added
Setting the smoothness parameters to a small value (0.02) improves things, but the break energy is still not perfect.
 |
 |
 |
 |
#6
Updated by Knödlseder Jürgen about 7 years ago
Note that a simple simulation of 30 min of data using the initial XML model, and fitting the simulated events using ctlike indicates that there are convergence problems. This typical happens for invalid gradients:
2017-03-28T13:53:07: +=================================+ 2017-03-28T13:53:07: | Maximum likelihood optimisation | 2017-03-28T13:53:07: +=================================+ 2017-03-28T13:53:07: >Iteration 0: -logL=157024.509, Lambda=1.0e-03 2017-03-28T13:53:07: Parameter "BreakEnergy" drives optimization step (step=0.497168) 2017-03-28T13:53:07: Iteration 1: -logL=157024.509, Lambda=1.0e-03, delta=-696.181, max(|grad|)=-38674.620105 [BreakEnergy:6] (stalled) 2017-03-28T13:53:07: Parameter "BreakEnergy" does not drive optimization step anymore. 2017-03-28T13:53:07: Iteration 2: -logL=157024.509, Lambda=1.0e-02, delta=-16.823, max(|grad|)=283.058569 [Index2:4] (stalled) 2017-03-28T13:53:07: >Iteration 3: -logL=157022.062, Lambda=1.0e-01, delta=2.446, max(|grad|)=11.343717 [Index:10] 2017-03-28T13:53:08: Iteration 4: -logL=157022.062, Lambda=1.0e-02, delta=-6.516, max(|grad|)=128.891134 [Index2:4] (stalled) 2017-03-28T13:53:08: >Iteration 5: -logL=157021.979, Lambda=1.0e-01, delta=0.083, max(|grad|)=4.062394 [Index1:3] 2017-03-28T13:53:08: Iteration 6: -logL=157021.979, Lambda=1.0e-02, delta=-2.265, max(|grad|)=66.600685 [Index2:4] (stalled) 2017-03-28T13:53:08: >Iteration 7: -logL=157021.925, Lambda=1.0e-01, delta=0.054, max(|grad|)=2.610925 [Index1:3] 2017-03-28T13:53:08: Iteration 8: -logL=157022.951, Lambda=1.0e-02, delta=-1.027, max(|grad|)=43.167592 [Index2:4] (stalled) 2017-03-28T13:53:08: >Iteration 9: -logL=157021.701, Lambda=1.0e-01, delta=1.250, max(|grad|)=3.155219 [Index2:4] 2017-03-28T13:53:08: Iteration 10: -logL=157021.722, Lambda=1.0e-02, delta=-0.020, max(|grad|)=14.470936 [BreakSmoothness:7] (stalled) 2017-03-28T13:53:08: >Iteration 11: -logL=157021.597, Lambda=1.0e-01, delta=0.124, max(|grad|)=0.701246 [Index2:4] 2017-03-28T13:53:08: >Iteration 12: -logL=157021.577, Lambda=1.0e-02, delta=0.020, max(|grad|)=-9.644692 [BreakEnergy:6] 2017-03-28T13:53:08: Iteration 13: -logL=157021.577, Lambda=1.0e-03, delta=-100.522, max(|grad|)=-2958.702580 [BreakEnergy:6] (stalled) 2017-03-28T13:53:08: >Iteration 14: -logL=157021.504, Lambda=1.0e-02, delta=0.073, max(|grad|)=-8.585555 [BreakEnergy:6] 2017-03-28T13:53:08: Iteration 15: -logL=157021.504, Lambda=1.0e-03, delta=-77.307, max(|grad|)=-2626.302597 [BreakEnergy:6] (stalled) 2017-03-28T13:53:08: >Iteration 16: -logL=157021.447, Lambda=1.0e-02, delta=0.057, max(|grad|)=-9.376360 [BreakEnergy:6] 2017-03-28T13:53:08: Iteration 17: -logL=157021.447, Lambda=1.0e-03, delta=-53.416, max(|grad|)=-2163.691984 [BreakEnergy:6] (stalled) 2017-03-28T13:53:08: >Iteration 18: -logL=157021.385, Lambda=1.0e-02, delta=0.061, max(|grad|)=-10.625401 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 19: -logL=157021.385, Lambda=1.0e-03, delta=-38.693, max(|grad|)=-1886.824989 [BreakEnergy:6] (stalled) 2017-03-28T13:53:09: >Iteration 20: -logL=157021.307, Lambda=1.0e-02, delta=0.079, max(|grad|)=-13.005929 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 21: -logL=157021.307, Lambda=1.0e-03, delta=-26.422, max(|grad|)=-1607.011273 [BreakEnergy:6] (stalled) 2017-03-28T13:53:09: >Iteration 22: -logL=157021.169, Lambda=1.0e-02, delta=0.137, max(|grad|)=-17.670655 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 23: -logL=157021.169, Lambda=1.0e-03, delta=-11.325, max(|grad|)=-1036.219873 [BreakEnergy:6] (stalled) 2017-03-28T13:53:09: >Iteration 24: -logL=157020.836, Lambda=1.0e-02, delta=0.334, max(|grad|)=-27.476821 [BreakEnergy:6] 2017-03-28T13:53:09: >Iteration 25: -logL=157020.782, Lambda=1.0e-03, delta=0.054, max(|grad|)=-143.289056 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 26: -logL=157020.782, Lambda=1.0e-04, delta=-12.308, max(|grad|)=-1518.463100 [BreakEnergy:6] (stalled) 2017-03-28T13:53:09: >Iteration 27: -logL=157020.424, Lambda=1.0e-03, delta=0.358, max(|grad|)=-73.785714 [BreakEnergy:6] 2017-03-28T13:53:09: >Iteration 28: -logL=157020.337, Lambda=1.0e-04, delta=0.087, max(|grad|)=-20.992382 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 29: -logL=157021.414, Lambda=1.0e-05, delta=-1.077, max(|grad|)=-422.701433 [BreakEnergy:6] (stalled) 2017-03-28T13:53:09: >Iteration 30: -logL=157021.237, Lambda=1.0e-04, delta=0.177, max(|grad|)=-286.616084 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 31: -logL=157021.237, Lambda=1.0e-05, delta=-1.526, max(|grad|)=-679.670560 [BreakEnergy:6] (stalled) 2017-03-28T13:53:09: >Iteration 32: -logL=157020.652, Lambda=1.0e-04, delta=0.585, max(|grad|)=-212.471407 [BreakEnergy:6] 2017-03-28T13:53:09: Iteration 33: -logL=157020.652, Lambda=1.0e-05, delta=-18.470, max(|grad|)=-804.677134 [BreakEnergy:6] (stalled) 2017-03-28T13:53:10: >Iteration 34: -logL=157020.632, Lambda=1.0e-04, delta=0.021, max(|grad|)=-162.477223 [BreakEnergy:6] 2017-03-28T13:53:10: Iteration 35: -logL=157020.632, Lambda=1.0e-05, delta=-2.698, max(|grad|)=-768.319088 [BreakEnergy:6] (stalled) 2017-03-28T13:53:10: Iteration 36: -logL=157020.632, Lambda=1.0e-04, delta=-0.133, max(|grad|)=-250.945203 [BreakEnergy:6] (stalled) 2017-03-28T13:53:10: >Iteration 37: -logL=157020.332, Lambda=1.0e-03, delta=0.300, max(|grad|)=-11.121790 [BreakEnergy:6] 2017-03-28T13:53:10: >Iteration 38: -logL=157020.330, Lambda=1.0e-04, delta=0.002, max(|grad|)=-4.292704 [BreakEnergy:6]
#7
Updated by Knödlseder Jürgen about 7 years ago
I replaced the analytical gradient computations by numerical gradient computations which did not change the results. This indicates that the problem is possibly inherent to the model.
I therefore change the slope of the second index to a stepper value as follows:
<?xml version="1.0" standalone="no"?>
<source_library title="source library">
<source name="Crab" type="PointSource">
<spectrum type="SmoothBrokenPowerLaw">
<parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
<parameter name="Index1" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
<parameter name="PivotEnergy" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
<parameter name="Index2" scale="-1" value="3.48" min="0.01" max="+5.0" free="1"/>
<parameter name="BreakEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="1"/>
<parameter name="BreakSmoothness" scale="1.0" value="0.2" min="0.01" max="10.0" free="1"/>
</spectrum>
<spatialModel type="PointSource">
<parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
<parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
</spatialModel>
</source>
<source name="Background" type="CTAIrfBackground" instrument="CTA">
<spectrum type="PowerLaw">
<parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e3" free="1"/>
<parameter name="Index" scale="1.0" value="0.0" min="-5.0" max="+5.0" free="1"/>
<parameter name="PivotEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
</spectrum>
</source>
</source_library>
Now the fit seems to converge nicely:
2017-03-28T14:11:14: +=================================+ 2017-03-28T14:11:14: | Maximum likelihood optimisation | 2017-03-28T14:11:14: +=================================+ 2017-03-28T14:11:14: >Iteration 0: -logL=156311.871, Lambda=1.0e-03 2017-03-28T14:11:14: >Iteration 1: -logL=156310.512, Lambda=1.0e-03, delta=1.360, max(|grad|)=2.459834 [Index:10] 2017-03-28T14:11:14: >Iteration 2: -logL=156310.510, Lambda=1.0e-04, delta=0.002, max(|grad|)=-0.124483 [BreakEnergy:6]
#8
Updated by Knödlseder Jürgen about 7 years ago
- Status changed from In Progress to Feedback
- % Done changed from 50 to 100
Added user documentation to GammaLib and ctools.
#9
Updated by Knödlseder Jürgen almost 7 years ago
- Status changed from Feedback to Closed