Change request #2892

Allow fitting of instrument scaling factors

Added by Knödlseder Jürgen almost 5 years ago. Updated almost 5 years ago.

Status:ClosedStart date:05/22/2019
Priority:NormalDue date:
Assigned To:Knödlseder Jürgen% Done:

100%

Category:-
Target version:1.6.0
Duration:

Description

It should be possible to fit the instrument scaling factors.

This can probably achieved by adding in GOptimizerPars GModels::pars(void) the instrument scaling factors to the list of parameters to optimise. It remains to be seen if more needs to be changed.

Recurrence

No recurrence.

History

#1 Updated by Knödlseder Jürgen almost 5 years ago

The method GResponse::eval_prob needs to take into account the scale factor gradient.

#2 Updated by Knödlseder Jürgen almost 5 years ago

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

I did the following modifications:

The GModel::scales() method was added that returns the number of instrument scales for a given model. In addition, GModel::scale() access operators were added to access the scale with a given index. Finally, scale parameters allocated by GModel now have a scale parameter gradient.

The GModelSky::set_pointers() method now appends any scale factors after the spatial, spectral, and temporal parameters to the model list. In that way, the instrument scale parameters are counted as normal model parameters. The GModelSky::print() method was modified so that it outputs now also the scale parameters for each model, as shown below:

2019-05-22T20:50:05: === GModelSky ===
2019-05-22T20:50:05:  Name ......................: Crab
2019-05-22T20:50:05:  Instruments ...............: all
2019-05-22T20:50:05:  Observation identifiers ...: all
2019-05-22T20:50:05:  Model type ................: PointSource
2019-05-22T20:50:05:  Model components ..........: "PointSource" * "PowerLaw" * "Constant" 
2019-05-22T20:50:05:  Number of parameters ......: 7
2019-05-22T20:50:05:  Number of spatial par's ...: 2
2019-05-22T20:50:05:   RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
2019-05-22T20:50:05:   DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
2019-05-22T20:50:05:  Number of spectral par's ..: 3
2019-05-22T20:50:05:   Prefactor ................: 5.42372069310751e-16 +/- 4.15467688757817e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2019-05-22T20:50:05:   Index ....................: -2.45717293803165 +/- 0.0300356983499498 [-0,-5]  (free,scale=-1,gradient)
2019-05-22T20:50:05:   PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
2019-05-22T20:50:05:  Number of temporal par's ..: 1
2019-05-22T20:50:05:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2019-05-22T20:50:05:  Number of scale par's .....: 1
2019-05-22T20:50:05:   CTAOnOff .................: 1.05046099933027 +/- 0.0886387854346198 [0.1,10]  (free,scale=1,gradient)

The instrument scale factor gradient computation for 3D analysis was implemented in the GResponse::eval_prob() method. And the GCTAOnOffObservation::N_gamma() method was adapted to achieve the same thing for the On/Off analysis.

Tests were done with two Crab observations, where the first observation was labeled for CTA and the second for HESS. The model fit output for the On/Off analysis is shown above, the one for an unbinned 3D analysis below.

2019-05-22T22:06:35: === GModelSky ===
2019-05-22T22:06:35:  Name ......................: Crab
2019-05-22T22:06:35:  Instruments ...............: all
2019-05-22T22:06:35:  Observation identifiers ...: all
2019-05-22T22:06:35:  Model type ................: PointSource
2019-05-22T22:06:35:  Model components ..........: "PointSource" * "PowerLaw" * "Constant" 
2019-05-22T22:06:35:  Number of parameters ......: 7
2019-05-22T22:06:35:  Number of spatial par's ...: 2
2019-05-22T22:06:35:   RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
2019-05-22T22:06:35:   DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
2019-05-22T22:06:35:  Number of spectral par's ..: 3
2019-05-22T22:06:35:   Prefactor ................: 5.61765554070147e-16 +/- 3.71298776013594e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2019-05-22T22:06:35:   Index ....................: -2.46875191786324 +/- 0.0260350540789234 [-0,-5]  (free,scale=-1,gradient)
2019-05-22T22:06:35:   PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
2019-05-22T22:06:35:  Number of temporal par's ..: 1
2019-05-22T22:06:35:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2019-05-22T22:06:35:  Number of scale par's .....: 1
2019-05-22T22:06:35:   CTA ......................: 1.01569396284157 +/- 0.0727807512411943 [0.1,10]  (free,scale=1,gradient)

For comparison we show also the unbinned 3D fit without an instrument scaling factor below. There are now only 6 instead of 7 parameters, and the statistical uncertainties are smaller.

2019-05-22T22:07:55: === GModelSky ===
2019-05-22T22:07:55:  Name ......................: Crab
2019-05-22T22:07:55:  Instruments ...............: all
2019-05-22T22:07:55:  Observation identifiers ...: all
2019-05-22T22:07:55:  Model type ................: PointSource
2019-05-22T22:07:55:  Model components ..........: "PointSource" * "PowerLaw" * "Constant" 
2019-05-22T22:07:55:  Number of parameters ......: 6
2019-05-22T22:07:55:  Number of spatial par's ...: 2
2019-05-22T22:07:55:   RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)
2019-05-22T22:07:55:   DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)
2019-05-22T22:07:55:  Number of spectral par's ..: 3
2019-05-22T22:07:55:   Prefactor ................: 5.69574504609791e-16 +/- 1.01198119103208e-17 [1e-23,1e-13] ph/cm2/s/MeV (free,scale=1e-16,gradient)
2019-05-22T22:07:55:   Index ....................: -2.47322510783642 +/- 0.0154943435176592 [-0,-5]  (free,scale=-1,gradient)
2019-05-22T22:07:55:   PivotEnergy ..............: 300000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)
2019-05-22T22:07:55:  Number of temporal par's ..: 1
2019-05-22T22:07:55:   Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
2019-05-22T22:07:55:  Number of scale par's .....: 0

Finally, the unit test for the model was adapted to take into account the additional parameters.

#3 Updated by Knödlseder Jürgen almost 5 years ago

For the record, here the XML model definition file for the 3D analysis tests:

<?xml version="1.0" standalone="no"?>
<source_library title="source library">
  <source name="Crab" type="PointSource">
    <spectrum type="PowerLaw">
       <parameter name="Prefactor"   scale="1e-16" value="5.7"  min="1e-07" max="1000.0" free="1"/>
       <parameter name="Index"       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"/>
    </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>
    <scaling>
      <instrument name="CTA" scale="1.0" min="0.1" max="10.0" value="1.0" free="1"/>
    </scaling>
  </source>
  <source name="CTA background model" type="CTAIrfBackground" instrument="CTA">
    <spectrum type="PowerLaw">  
      <parameter name="Prefactor"   scale="1.0"  value="1.0"  min="1e-3" max="1e+3"   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 name="HESS background model" type="CTAIrfBackground" instrument="HESS">
    <spectrum type="PowerLaw">
      <parameter name="Prefactor"   scale="1.0"  value="1.0"  min="1e-3" max="1e+3"   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>

#4 Updated by Knödlseder Jürgen almost 5 years ago

  • Status changed from Pull request to Closed

Code was merged into devel.

Also available in: Atom PDF