Feature #1982

Phase information in event file

Added by Cardenzana Josh almost 7 years ago. Updated almost 7 years ago.

Status:ClosedStart date:04/11/2017
Priority:NormalDue date:
Assigned To:Di Venere Leonardo% Done:

100%

Category:-Estimated time:3.00 hours
Target version:1.3.0
Duration:

Description

Add a new tool 'ctphase’ for adding phase information to event files. Allow both information from an xml file and parameters from a 'ctphase.par’ file.

MJD-59000_P0-0_F0-1_F1-1e-8_F2-0.png - Uniformly distributed phases (210 KB) Cardenzana Josh, 04/06/2017 11:44 AM

MJD-59000_P0-0_F0-1_F1-1e-8_F2-1e-8.png - Precision issues in phase computation (128 KB) Cardenzana Josh, 04/06/2017 11:52 AM

Period-3minutes.png - Phases of a source with a 3 minute period (92.1 KB) Cardenzana Josh, 04/06/2017 12:03 PM

model_phase_info.fits - Phase node information (8.44 KB) Cardenzana Josh, 04/06/2017 05:51 PM

model_phase_info.png - Phase node information (35.2 KB) Cardenzana Josh, 04/06/2017 05:57 PM

PhaseFoldedCrab.png - Phase folded distribution of Crab (25.9 KB) Cardenzana Josh, 04/06/2017 06:02 PM

PhaseFoldedCrab_RandomPhaseInfo.png - Phase folded distribution of Crab (Random) (23.1 KB) Cardenzana Josh, 04/06/2017 06:07 PM

test_model.xml Magnifier - Model defining temporal component for Crab model (1.8 KB) Cardenzana Josh, 04/11/2017 05:15 PM

test.py Magnifier - Python script for testing ctphase (2.29 KB) Cardenzana Josh, 04/11/2017 05:15 PM

test.py Magnifier (2.52 KB) Knödlseder Jürgen, 04/12/2017 09:35 AM

Mjd-59000_p0-0_f0-1_f1-1e-8_f2-0 Mjd-59000_p0-0_f0-1_f1-1e-8_f2-1e-8 Period-3minutes Model_phase_info Phasefoldedcrab Phasefoldedcrab_randomphaseinfo

Recurrence

No recurrence.

History

#1 Updated by Cardenzana Josh almost 7 years ago

The tool has been developed and somewhat tested. Phase information is derived using the GModelTemporalPhaseCurve class and the user supplied reference values for:
  • “MJD” (reference MJD)
  • “P0” (phase at reference MJD)
  • “F0” (frequency at reference MJD in Hz)
  • “F1” (first frequency derivative at reference MJD in sec^(-1))
  • “F2” (second frequency derivative at reference MJD in sec^(-2))

The phase information is successfully appended to the events file, however precision issues can become a problem. These issues appear to only occur when a reference date far from the MJD of the data is chosen (more than a few years difference in the dates) or when a large value for the first/second derivate are chosen. For the time being, a warning is printed into the ctphase log file when either (F1 * dt) or (F2 * dt * dt) is greater than 10^8 (dt is seconds between the reference MJD and the MJD of the events). When this occurs the user should check that the distribution of derived phases makes sense. The method should be stable for most sensible choices of the parameters.

Tests of the above issue have been done simulating events with ctobssim and then passing the resulting events file through ctphase. These results are documented below:

Example plot for 3 minute period:
  • Data MJD: 51544.5
  • Reference MJD: 59000
  • P0: 0.0
  • F0: 0.00555556 (1/180)
  • F1: 0.0
  • F2: 0.0
    Phases of a source with a 3 minute period
Example plot for sensible, uniform distribution of phases as a function of event time:
  • Data MJD: 51544.5
  • Reference MJD: 59000
  • P0: 0.0
  • F0: 1.0
  • F1: 1e-8
  • F2: 0
    Uniformly distributed phases
Example plot showing precision issues (generated phases do not appear uniformly distributed)
  • Data MJD: 51544.5
  • Reference MJD: 59000
  • P0: 0.0
  • F0: 1.0
  • F1: 1.0e-8
  • F2: 1.0e-8
    Precision issues in phase computation

#2 Updated by Cardenzana Josh almost 7 years ago

Added ability to obtain the phase information from a model xml file. User needs to specify two parameters:
  • inmodel (model XML file containing source information)
  • source (Name of the source in 'inmodel’ from which phase information will be taken)
If either of these parameters is not specified, then the phase information is derived from the parameters in the ctphase.par file. An exception is thrown if either of the following occurs:
  • Model XML file does not define a source with name 'source’
  • 'source’ does not define a temporal model with phase information.

This functionality was tested using a temporal model for the Crab of:

  <source name="Crab" type="PointSource">
    <spatialModel type="PointSource">
      <parameter name="RA"  scale="1" value="83.6331" min="-360" max="360" free="0"/>
      <parameter name="DEC" scale="1" value="22.0145" min="-90"  max="90"  free="0"/>
    </spatialModel>
    <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>
    <temporal type="PhaseCurve" file="model_phase_info.fits">
      <parameter name="Normalization" scale="1" value="1.0"     min="0.0" max="1000.0"   free="0"/>
      <parameter name="MJD"           scale="1" value="51544.5" min="0.0" max="100000.0" free="0"/>
      <parameter name="Phase"         scale="1" value="0.0"     min="0.0" max="1.0"      free="0"/>
      <parameter name="F0"            scale="1" value="1.0"     min="0.0" max="1000.0"   free="0"/>
      <parameter name="F1"            scale="1" value="0.1"     min="0.0" max="1000.0"   free="0"/>
      <parameter name="F2"            scale="1" value="0.01"    min="0.0" max="1000.0"   free="0"/>
    </temporal>
  </source>

'model_phase_info.fits’ is attached and defines node information in the following way:
Phase node information

Generating events from only the above Crab model and processing through ctphase produces the following distribution of events (as a function of phase):
Phase folded distribution of Crab

Processing the same data through ctphase, but using random values for the phase information produces the following (seemingly random) distribution of events. This is the expected behavior.
Phase folded distribution of Crab (Random)

I presume that the features requested are now implemented and a pull of this code can be done. Note, a pull of the associated gammaLib branch is also necessary for appropriate functionality, otherwise errors will be produced by GModelTemporalPhaseCurve.

Relevant repositories for pull request:

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

  • Target version set to 1.3.0
  • Start date set to 04/11/2017

I integrated the GammaLib branch 1982-add-phase-info.

I’ve seen that you allowed an empty filename in the GModelTemporalPhaseCurve::read() method, however I’m worried about the possible side effects of that. In particular, a phase object without any nodes does not really make sense. What is the Use Case of this functionality? For the moment I have not merged in that part.

#4 Updated by Cardenzana Josh almost 7 years ago

We originally were going to just set the parameters individually using the dedicated "::phase()”, "::f0()”, etc... methods, but ran into issues where any parameter values passed that were outside the range defined in “init_members()” resulted in an exception (this was an issue for negative values of F1 and F2).

However, we could control the Min/Max values if we pass the parameters inside a “GXmlElement” object. But then there’s an error that gets produced in the “GModelTemporalPhaseCurve::read(const GXmlElement& xml)” method when the value in the “file” attribute doesn’t exist. When computing the phases (to the best of my knowledge) we only need to know the time of an event, thus the node information from the file is never used. An empty file attribute still produces an error in an analysis where the node information is needed. But on second thought the original “invalid file” error is probably more useful for the user than a seg-fault later on during execution.

Alternative approach:
Expand the ranges of the f1/f2 parameters in “GModelTemporalPhaseCurve::init_members()” to be [-1000, 1000] instead of [0,1000] and then set all phase/frequency parameters via their dedicated call methods (doesn’t require an additional constructor and the read method can be reset to requiring a legitimate file attribute). This should probably be done anyway, since the frequency of a pulsar will most likely decrease over time (i.e. f1 will most likely be negative).

Unless you can think of a reason not to, I will make this change, rerun the tests to make sure everything still works as expected, and post once I’ve verified the results haven’t changed.

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

Indeed, I had not recognized that the phase file is not needed ;)

I prefer expanding the default range of F1 and F2 to [-1000,1000] so that the coherence of the XML format is always checked. Otherwise a user may eventually omit the file attribute and does not understand why the model is not working.

#6 Updated by Knödlseder Jürgen almost 7 years ago

  • Status changed from Pull request to Closed

I merged the ctphase tool into devel.

I also added a page to the reference manual.

#7 Updated by Cardenzana Josh almost 7 years ago

I’ve made the above change and verified that the results are consistent with those presented previously. I’ve attached my test python script (with supplemental files) that produces the events, then stores the phase information and plots the distribution of events in “phase”. Note that the 'model_phase_info.fits’ file attached in a previous post is also necessary to get the node information in the pulse profile.

Relevant repositories for pull request:

#8 Updated by Knödlseder Jürgen almost 7 years ago

Thanks. Note that I modified two parameter names so that they are compliant with the names of other tools:

source -> srcname
p0     -> phase

I modified a bit the script to avoid using astropy (which I don’t have installed): test.py

Also available in: Atom PDF