Pierrick Martin / ctools

Commit 8ef5130628f2c1df7d2ecdc717a91993012fd93f

Authored by L. Tibaldo
1 parent d9c5626a

working

Inline Side-by-side
Showing 2 changed files with 122 additions and 26 deletions
cscripts/csresspec.par
  View file @8ef5130
... ... @@ -31,13 +31,26 @@ outfile, f, a, resspec.fits,,, "Output residual spectrum file"
31 31 #
32 32 # Script parameters
33 33 #==================
  34 +statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic"
34 35 components b, h, no,,, "Calculate model for individual components"
35 36 ebinalg, s, h, LOG,LIN|LOG|FILE,, "Binning algorithm"
36 37 emin, r, a, 0.1,,, "Lower energy limit (TeV)"
37 38 emax, r, a, 100.0,,, "Upper energy limit (TeV)"
38 39 enumbins, i, h, 20,,, "Number of model cube energy bins"
39 40 ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition"
40   -stack, b, a, yes,,, "Stack observations"
  41 +stack, b, a, no,,, "Stack observations"
  42 +coordsys, s, a, CEL,CEL|GAL,, "Coordinate System"
  43 +proj, s, a, CAR,AIT|AZP|CAR|GLS|MER|MOL|SFL|SIN|STG|TAN,, "Projection method"
  44 +xref, r, a, 83.63,0,360, "First coordinate of image/source region center in degrees (RA or galactic l)"
  45 +yref, r, a, 22.01,-90,90, "Second coordinate of image/source region center in degrees (DEC or galactic b)"
  46 +nxpix, i, a, 200,,, "Size of the X axis in pixels"
  47 +nypix, i, a, 200,,, "Size of the Y axis in pixels"
  48 +binsz, r, a, 0.02,,, "Pixel size (deg/pixel)"
  49 +mask b, a, no,,, "Mask data to calculate residuals in ROI"
  50 +ra, r, a, 83.63,0,360, "RA for ROI centre (degrees)"
  51 +dec, r, a, 22.01,-90,90, "Dec for ROI centre (degrees)"
  52 +rad, r, a, 3.0,0,180, "Radius of ROI (degrees)"
  53 +regfile, f, a, "NONE",,, "Input exclusion region file in ds9 format"
41 54  
42 55 #
43 56 # Standard parameters
... ...
cscripts/csresspec.py
  View file @8ef5130
... ... @@ -20,6 +20,7 @@
20 20 # ==========================================================================
21 21 import gammalib
22 22 import ctools
  23 +from cscripts import obsutils
23 24  
24 25  
25 26 # =============== #
... ... @@ -40,9 +41,9 @@ class csresspec(ctools.csobservation):
40 41 argv)
41 42  
42 43 # Initialise class members
43   - self._components = False
44 44 self._use_maps = False
45 45 self._stack = False
  46 + self._mask = False
46 47  
47 48 # Return
48 49 return
... ... @@ -52,10 +53,9 @@ class csresspec(ctools.csobservation):
52 53 """
53 54 Get parameters from parfile and setup the observation
54 55 """
55   - # Set observation if not done before
56   - if self.obs().is_empty():
57   - self._require_inobs('csresspec._get_parameters()')
58   - self.obs(self._get_observations())
  56 + # Setup observations (require response and allow event list as well as
  57 + # counts cube)
  58 + self._setup_observations(self.obs(), True, True, True)
59 59  
60 60 # Set observation statistic
61 61 self._set_obs_statistic(gammalib.toupper(self['statistic'].string()))
... ... @@ -76,26 +76,26 @@ class csresspec(ctools.csobservation):
76 76 n_cta = n_unbinned + n_binned + n_onoff
77 77 n_other = self.obs().size() - n_cta
78 78  
  79 + # Log census of input observations
  80 + self._log_value(gammalib.TERSE, 'Unbinned CTA observations', n_unbinned)
  81 + self._log_value(gammalib.TERSE, 'Binned CTA observations', n_binned)
  82 + self._log_value(gammalib.TERSE, 'On/off CTA observations', n_onoff)
  83 + self._log_value(gammalib.TERSE, 'Other observations', n_other)
  84 + self._log_string(gammalib.TERSE, '\nWARNING: Only CTA observation '
  85 + 'can be handled, all '
  86 + 'non-CTA observation will be ignored.\n')
  87 +
79 88 # Query wether to compute model for individual components
80   - self._components = self['components'].boolean()
  89 + components = self['components'].boolean()
81 90  
82 91 # If there is only one binned observation and no model for individual
83 92 # components is required, query for precomputed model file and set
84 93 # use_maps to True
85   - if self.obs().size() == 1 and n_binned == 1 and not self._components:
  94 + if self.obs().size() == 1 and n_binned == 1 and not components:
86 95 modcube = self['modcube'].filename()
87 96 if modcube != 'NONE':
88 97 self._use_maps = True
89 98  
90   - # If we are not using a precomputed model query input XML model file
91   - if not self._use_maps:
92   - modxml = self['inmodel'].filename()
93   - # if None check whether models are provided in observation
94   - # container, otherwise throw exception and stop here
95   - if modxml == 'NONE' and self.obs().models().is_empty():
96   - msg = 'No model provided. Please specify an input XML model file'
97   - raise RuntimeError(msg)
98   -
99 99 # If there are unbinned observations query the energy binning parameters
100 100 if n_unbinned != 0:
101 101 self['ebinalg'].string()
... ... @@ -105,12 +105,50 @@ class csresspec(ctools.csobservation):
105 105 self['emin'].real()
106 106 self['emax'].real()
107 107 self['enumbins'].integer()
108   -
109   - # If there is more than one observation, and observations are unbinned
110   - # or onoff query user to know if they wish stacked results
  108 + msg = 'User defined energy binning will be used for %s unbinned observations.' % n_unbinned
  109 + self._log_value(gammalib.TERSE, msg)
  110 + if n_cta > n_unbinned:
  111 + n_notunbin = n_cta - n_unbinned
  112 + msg = 'The intrinsic binning will be used for the remaining %s observations.' % (
  113 + n_notunbin)
  114 + self._log_value(gammalib.TERSE, msg)
  115 +
  116 + # If there is more than one observation, and observations are all
  117 + # unbinned or all onoff query user to know if they wish stacked results
111 118 if self.obs().size() > 1 and (
112 119 n_unbinned == self.obs().size() or n_onoff == self.obs().size()):
113   - self['stack'].boolean()
  120 + self._stack = self['stack'].boolean()
  121 + # If we are to stack event lists query parameters for cube creation
  122 + if n_unbinned == self.obs().size():
  123 + self['coordsys'].string()
  124 + self['proj'].string()
  125 + self['xref'].real()
  126 + self['yref'].real()
  127 + self['nxpix'].real()
  128 + self['nypix'].real()
  129 + self['binsz'].real()
  130 +
  131 + # If we are not using a precomputed model query input XML model file
  132 + if not self._use_maps:
  133 + modxml = self['inmodel'].filename()
  134 + # If None check whether models are provided in observation
  135 + # container, otherwise throw exception and stop here
  136 + if modxml == 'NONE' and self.obs().models().is_empty():
  137 + msg = 'No model provided. Please specify an input XML model file'
  138 + raise RuntimeError(msg)
  139 + # Otherwise set the input model
  140 + else:
  141 + self.obs().models(modxml)
  142 +
  143 + # Unless all observations are On/Off query for mask definition
  144 + if n_onoff == n_cta:
  145 + pass
  146 + else:
  147 + self._mask = self['mask'].boolean()
  148 + if self._mask:
  149 + self['ra'].real()
  150 + self['dec'].real()
  151 + self['rad'].real()
114 152  
115 153 # Apply energy dispersion
116 154 self['edisp'].boolean()
... ... @@ -122,11 +160,56 @@ class csresspec(ctools.csobservation):
122 160 # Write input parameters into logger
123 161 self._log_parameters(gammalib.TERSE)
124 162  
125   - # Log census of input observations
126   - self._log_value(gammalib.TERSE, 'Unbinned CTA observations', n_unbinned)
127   - self._log_value(gammalib.TERSE, 'Binned CTA observations', n_binned)
128   - self._log_value(gammalib.TERSE, 'On/off CTA observations', n_onoff)
129   - self._log_value(gammalib.TERSE, 'Other observations', n_other)
  163 + # Return
  164 + return
  165 +
  166 + def _stack_observations(self):
  167 + # Use first observation to determine the type and apply the
  168 + # approptiate stacking method
  169 +
  170 + # If On/Off use the GCTAOnOffObservation constructor
  171 + if self.obs()[0].classname() == 'GCTAOnOffObservation':
  172 + msg = 'Stacking %s On/Off observations.'%(self.obs().size())
  173 + self._log_value(gammalib.TERSE, msg)
  174 + stacked_obs = gammalib.GCTAOnOffObservation(self.obs)
  175 + # reset observation container
  176 + self.obs = gammalib.GObservations()
  177 + # and append the new stacked observation
  178 + self.obs.append(stacked_obs)
  179 +
  180 + # If event list then bin observations
  181 + elif self.obs()[0].classname() == 'GCTAObservation':
  182 + msg = 'Stacking %s event lists.' % (self.obs().size())
  183 + self._log_value(gammalib.TERSE, msg)
  184 + stacked_obs = obsutils.get_stacked_obs(self, self.obs)
  185 + self.obs = stacked_obs
130 186  
131 187 # Return
132 188 return
  189 +
  190 + # Public methods
  191 + def run(self):
  192 + """
  193 + Run the script
  194 + """
  195 + # Switch screen logging on in debug mode
  196 + if self._logDebug():
  197 + self._log.cout(True)
  198 +
  199 + # Get parameters
  200 + self._get_parameters()
  201 +
  202 + # Write observation into logger
  203 + self._log_observations(gammalib.NORMAL, self.obs(), 'Observation')
  204 +
  205 + # Stack observations if requested
  206 + if self._stack:
  207 + self._stack_observations()
  208 +
  209 + # Loop over observations and calculate residuals
  210 + for obs in self.obs():
  211 +
  212 +
  213 +
  214 +
  215 +
... ...