Commit ddbd70541db3cce08243b95a8ac946fe070723e5
Merge branch '1646-write-iact-docs' into integration
# Conflicts: # ChangeLog
Showing
14 changed files
with
1699 additions
and
60 deletions
ChangeLog
| @@ -4,6 +4,7 @@ | @@ -4,6 +4,7 @@ | ||
| 4 | ======================== | 4 | ======================== |
| 5 | 5 | ||
| 6 | Test and import IACT scripts only for Python 2.6+ (#1715) | 6 | Test and import IACT scripts only for Python 2.6+ (#1715) |
| 7 | + Add documentation for IACT analysis (#1646) | ||
| 7 | Adapt to revise GTime operators (#1714) | 8 | Adapt to revise GTime operators (#1714) |
| 8 | Add publish() methods to tools producing images (#1710) | 9 | Add publish() methods to tools producing images (#1710) |
| 9 | Read ahead output filename in cterror (#1700) | 10 | Read ahead output filename in cterror (#1700) |
cscripts/csiactdata.py
| @@ -40,9 +40,10 @@ class csiactdata(ctools.cscript): | @@ -40,9 +40,10 @@ class csiactdata(ctools.cscript): | ||
| 40 | Constructor. | 40 | Constructor. |
| 41 | """ | 41 | """ |
| 42 | # Set name | 42 | # Set name |
| 43 | - self.name = "csiactdata" | ||
| 44 | - self.version = "1.0.0" | ||
| 45 | - self.datapath = os.getenv("VHEFITS","") | 43 | + self.name = "csiactdata" |
| 44 | + self.version = "1.0.0" | ||
| 45 | + self.datapath = os.getenv("VHEFITS","") | ||
| 46 | + self.prodnames = [] | ||
| 46 | 47 | ||
| 47 | # Make sure that parfile exists | 48 | # Make sure that parfile exists |
| 48 | file = self.parfile() | 49 | file = self.parfile() |
| @@ -156,6 +157,9 @@ class csiactdata(ctools.cscript): | @@ -156,6 +157,9 @@ class csiactdata(ctools.cscript): | ||
| 156 | data = json.loads(json_data) | 157 | data = json.loads(json_data) |
| 157 | configs = data["datasets"] | 158 | configs = data["datasets"] |
| 158 | 159 | ||
| 160 | + # Initialise array for available names | ||
| 161 | + self.prodnames = [] | ||
| 162 | + | ||
| 159 | # Loop over configs and display unavailable storage first | 163 | # Loop over configs and display unavailable storage first |
| 160 | for config in configs: | 164 | for config in configs: |
| 161 | 165 | ||
| @@ -163,8 +167,11 @@ class csiactdata(ctools.cscript): | @@ -163,8 +167,11 @@ class csiactdata(ctools.cscript): | ||
| 163 | hdu = os.path.join(self.datapath, config["hduindx"]) | 167 | hdu = os.path.join(self.datapath, config["hduindx"]) |
| 164 | obs = os.path.join(self.datapath, config["obsindx"]) | 168 | obs = os.path.join(self.datapath, config["obsindx"]) |
| 165 | 169 | ||
| 170 | + filename_hdu = gammalib.GFilename(str(hdu)+"[HDU_INDEX]") | ||
| 171 | + filename_obs = gammalib.GFilename(str(obs)+"[OBS_INDEX]") | ||
| 172 | + | ||
| 166 | # Check if index files are available | 173 | # Check if index files are available |
| 167 | - if not (gammalib.is_fits(str(hdu)) and gammalib.is_fits(str(obs))): | 174 | + if not (filename_hdu.is_fits() and filename_obs.is_fits()): |
| 168 | self.log.parformat(str(config["name"])) | 175 | self.log.parformat(str(config["name"])) |
| 169 | self.log("Not available") | 176 | self.log("Not available") |
| 170 | self.log("\n") | 177 | self.log("\n") |
| @@ -180,13 +187,34 @@ class csiactdata(ctools.cscript): | @@ -180,13 +187,34 @@ class csiactdata(ctools.cscript): | ||
| 180 | hdu = os.path.join(self.datapath, config["hduindx"]) | 187 | hdu = os.path.join(self.datapath, config["hduindx"]) |
| 181 | obs = os.path.join(self.datapath, config["obsindx"]) | 188 | obs = os.path.join(self.datapath, config["obsindx"]) |
| 182 | 189 | ||
| 190 | + filename_hdu = gammalib.GFilename(str(hdu)+"[HDU_INDEX]") | ||
| 191 | + filename_obs = gammalib.GFilename(str(obs)+"[OBS_INDEX]") | ||
| 192 | + | ||
| 183 | # Check if index files are available | 193 | # Check if index files are available |
| 184 | - if gammalib.is_fits(str(hdu)) and gammalib.is_fits(str(obs)): | 194 | + if (filename_hdu.is_fits() and filename_obs.is_fits()): |
| 185 | 195 | ||
| 186 | # Log data information | 196 | # Log data information |
| 187 | self.log("\n") | 197 | self.log("\n") |
| 188 | self.log.header3(str(config["name"])) | 198 | self.log.header3(str(config["name"])) |
| 199 | + | ||
| 200 | + # Print important informations first | ||
| 201 | + self.log.parformat("Name") | ||
| 202 | + self.log(str(config["name"])) | ||
| 203 | + self.log("\n") | ||
| 204 | + self.log.parformat("Observation index") | ||
| 205 | + self.log(str(config["obsindx"])) | ||
| 206 | + self.log("\n") | ||
| 207 | + self.log.parformat("HDU index") | ||
| 208 | + self.log(str(config["hduindx"])) | ||
| 209 | + self.log("\n") | ||
| 210 | + | ||
| 211 | + # Append to available names | ||
| 212 | + self.prodnames.append(str(config["name"])) | ||
| 213 | + | ||
| 214 | + # Print additional information | ||
| 189 | for key in config: | 215 | for key in config: |
| 216 | + if key in ["name","hduindx","obsindx"]: | ||
| 217 | + continue | ||
| 190 | self.log.parformat(str(key)) | 218 | self.log.parformat(str(key)) |
| 191 | self.log(str(config[key])) | 219 | self.log(str(config[key])) |
| 192 | self.log("\n") | 220 | self.log("\n") |
| @@ -194,6 +222,14 @@ class csiactdata(ctools.cscript): | @@ -194,6 +222,14 @@ class csiactdata(ctools.cscript): | ||
| 194 | # Return | 222 | # Return |
| 195 | return | 223 | return |
| 196 | 224 | ||
| 225 | + def names(self): | ||
| 226 | + """ | ||
| 227 | + Return available FITS production names | ||
| 228 | + """ | ||
| 229 | + | ||
| 230 | + # Return | ||
| 231 | + return self.prodnames | ||
| 232 | + | ||
| 197 | # ======================== # | 233 | # ======================== # |
| 198 | # Main routine entry point # | 234 | # Main routine entry point # |
| 199 | # ======================== # | 235 | # ======================== # |
cscripts/csiactobs.py
| @@ -51,6 +51,7 @@ class csiactobs(ctools.cscript): | @@ -51,6 +51,7 @@ class csiactobs(ctools.cscript): | ||
| 51 | self.inmodels = None | 51 | self.inmodels = None |
| 52 | self.xml = gammalib.GXml() | 52 | self.xml = gammalib.GXml() |
| 53 | self.models = gammalib.GModels() | 53 | self.models = gammalib.GModels() |
| 54 | + self.m_runlist = [] | ||
| 54 | 55 | ||
| 55 | # Make sure that parfile exists | 56 | # Make sure that parfile exists |
| 56 | file = self.parfile() | 57 | file = self.parfile() |
| @@ -136,9 +137,27 @@ class csiactobs(ctools.cscript): | @@ -136,9 +137,27 @@ class csiactobs(ctools.cscript): | ||
| 136 | # Expand environment | 137 | # Expand environment |
| 137 | self.datapath = gammalib.expand_env(self.datapath) | 138 | self.datapath = gammalib.expand_env(self.datapath) |
| 138 | 139 | ||
| 139 | - # Read user parameters | 140 | + # Read FITS production |
| 140 | self.m_prodname = self["prodname"].string() | 141 | self.m_prodname = self["prodname"].string() |
| 141 | - self.m_runlistfile = self["infile"].filename() | 142 | + |
| 143 | + # Read runlist file if list not already filled | ||
| 144 | + if len(self.m_runlist) == 0: | ||
| 145 | + | ||
| 146 | + # Get file name | ||
| 147 | + self.m_runlistfile = self["infile"].filename() | ||
| 148 | + | ||
| 149 | + # Read runlist from file | ||
| 150 | + runfile = open(self.m_runlistfile.url()) | ||
| 151 | + for line in runfile.readlines(): | ||
| 152 | + if len(line) == 0: | ||
| 153 | + continue | ||
| 154 | + if line[0] == '#': | ||
| 155 | + continue | ||
| 156 | + if len(line.split()) > 0: | ||
| 157 | + self.m_runlist.append(line.split()[0]) | ||
| 158 | + runfile.close() | ||
| 159 | + | ||
| 160 | + # Read number of background parameters | ||
| 142 | self.m_bkgpars = self["bkgpars"].integer() | 161 | self.m_bkgpars = self["bkgpars"].integer() |
| 143 | 162 | ||
| 144 | # Output model file | 163 | # Output model file |
| @@ -220,6 +239,25 @@ class csiactobs(ctools.cscript): | @@ -220,6 +239,25 @@ class csiactobs(ctools.cscript): | ||
| 220 | # Return | 239 | # Return |
| 221 | return | 240 | return |
| 222 | 241 | ||
| 242 | + def runlist(self, runlist): | ||
| 243 | + """ | ||
| 244 | + Set observation list | ||
| 245 | + """ | ||
| 246 | + | ||
| 247 | + # Check if convertable to a python list | ||
| 248 | + try: | ||
| 249 | + | ||
| 250 | + # Store a copy of input argument | ||
| 251 | + self.m_runlist = list(runlist) | ||
| 252 | + | ||
| 253 | + except: | ||
| 254 | + | ||
| 255 | + # raise error if wrong type detected | ||
| 256 | + raise RuntimeError("Argument of csiactobs.runlist() must be of type 'list'") | ||
| 257 | + | ||
| 258 | + # Return | ||
| 259 | + return | ||
| 260 | + | ||
| 223 | def background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0): | 261 | def background_spectrum(self, run, prefactor, index, emin = 0.01, emax = 100.0): |
| 224 | 262 | ||
| 225 | # Handle constant spectral model | 263 | # Handle constant spectral model |
| @@ -390,30 +428,16 @@ class csiactobs(ctools.cscript): | @@ -390,30 +428,16 @@ class csiactobs(ctools.cscript): | ||
| 390 | self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\"")) | 428 | self.xml.append(gammalib.GXmlElement("observation_list title=\"observation list\"")) |
| 391 | lib = self.xml.element("observation_list", 0) | 429 | lib = self.xml.element("observation_list", 0) |
| 392 | 430 | ||
| 393 | - # Read runlist from file | ||
| 394 | - runlist = [] | ||
| 395 | - runfile = open(self.m_runlistfile.url()) | ||
| 396 | - for line in runfile.readlines(): | ||
| 397 | - if len(line) == 0: | ||
| 398 | - continue | ||
| 399 | - if line[0] == '#': | ||
| 400 | - continue | ||
| 401 | - if len(line.split()) > 0: | ||
| 402 | - runlist.append(line.split()[0]) | ||
| 403 | - else: | ||
| 404 | - runlist.append('') | ||
| 405 | - runfile.close() | ||
| 406 | - | ||
| 407 | # Log output | 431 | # Log output |
| 408 | if self.logTerse(): | 432 | if self.logTerse(): |
| 409 | self.log("\n") | 433 | self.log("\n") |
| 410 | - self.log.header1("Looping over "+str(len(runlist))+" runs") | 434 | + self.log.header1("Looping over "+str(len(self.m_runlist))+" runs") |
| 411 | 435 | ||
| 412 | # Initialise energy range values for logging | 436 | # Initialise energy range values for logging |
| 413 | self.m_ebounds.clear() | 437 | self.m_ebounds.clear() |
| 414 | 438 | ||
| 415 | # Loop over runs | 439 | # Loop over runs |
| 416 | - for obs_id in runlist: | 440 | + for obs_id in self.m_runlist: |
| 417 | 441 | ||
| 418 | # Create selection string | 442 | # Create selection string |
| 419 | obs_selection = "[OBS_ID=="+str(obs_id)+"]" | 443 | obs_selection = "[OBS_ID=="+str(obs_id)+"]" |
cscripts/csobsinfo.py
| @@ -40,6 +40,14 @@ class csobsinfo(ctools.cscript): | @@ -40,6 +40,14 @@ class csobsinfo(ctools.cscript): | ||
| 40 | self.name = "csobsinfo" | 40 | self.name = "csobsinfo" |
| 41 | self.version = "1.1.0" | 41 | self.version = "1.1.0" |
| 42 | 42 | ||
| 43 | + # Initialise some members | ||
| 44 | + if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): | ||
| 45 | + self.obs = argv[0] | ||
| 46 | + argv = argv[1:] | ||
| 47 | + else: | ||
| 48 | + self.obs = gammalib.GObservations() | ||
| 49 | + self.obs.clear() | ||
| 50 | + | ||
| 43 | # Make sure that parfile exists | 51 | # Make sure that parfile exists |
| 44 | file = self.parfile() | 52 | file = self.parfile() |
| 45 | 53 | ||
| @@ -99,9 +107,10 @@ class csobsinfo(ctools.cscript): | @@ -99,9 +107,10 @@ class csobsinfo(ctools.cscript): | ||
| 99 | """ | 107 | """ |
| 100 | Get parameters from parfile and setup the observation. | 108 | Get parameters from parfile and setup the observation. |
| 101 | """ | 109 | """ |
| 102 | - # Get parameters | ||
| 103 | - self.require_inobs("csobsinfo::get_parameters") | ||
| 104 | - self.obs = self.get_observations(False) | 110 | + # Get parameters |
| 111 | + if self.obs.size() == 0: | ||
| 112 | + self.require_inobs("csobsinfo::get_parameters") | ||
| 113 | + self.obs = self.get_observations(False) | ||
| 105 | 114 | ||
| 106 | # Initialise object position | 115 | # Initialise object position |
| 107 | self.obj_dir = gammalib.GSkyDir() | 116 | self.obj_dir = gammalib.GSkyDir() |
| @@ -137,16 +146,16 @@ class csobsinfo(ctools.cscript): | @@ -137,16 +146,16 @@ class csobsinfo(ctools.cscript): | ||
| 137 | 146 | ||
| 138 | # Initialise arrays to store certain values for reuse | 147 | # Initialise arrays to store certain values for reuse |
| 139 | # Todo, think about using a python dictionary | 148 | # Todo, think about using a python dictionary |
| 140 | - self.offsets = [] | ||
| 141 | - self.zeniths = [] | ||
| 142 | - self.azimuths = [] | 149 | + self.m_offsets = [] |
| 150 | + self.m_zeniths = [] | ||
| 151 | + self.m_azimuths = [] | ||
| 143 | self.pnt_ra = [] | 152 | self.pnt_ra = [] |
| 144 | self.pnt_dec = [] | 153 | self.pnt_dec = [] |
| 145 | obs_names = [] | 154 | obs_names = [] |
| 146 | 155 | ||
| 147 | # Initialise output to be filled | 156 | # Initialise output to be filled |
| 148 | - self.ebounds = gammalib.GEbounds() | ||
| 149 | - self.gti = gammalib.GGti() | 157 | + self.m_ebounds = gammalib.GEbounds() |
| 158 | + self.m_gti = gammalib.GGti() | ||
| 150 | ontime = 0.0 | 159 | ontime = 0.0 |
| 151 | livetime = 0.0 | 160 | livetime = 0.0 |
| 152 | n_events = 0 | 161 | n_events = 0 |
| @@ -195,10 +204,10 @@ class csobsinfo(ctools.cscript): | @@ -195,10 +204,10 @@ class csobsinfo(ctools.cscript): | ||
| 195 | 204 | ||
| 196 | # If avaliable append energy boundaries | 205 | # If avaliable append energy boundaries |
| 197 | if obs_bounds.size() > 0 : | 206 | if obs_bounds.size() > 0 : |
| 198 | - self.ebounds.merge(obs_bounds.emin(),obs_bounds.emax()) | 207 | + self.m_ebounds.append(obs_bounds.emin(),obs_bounds.emax()) |
| 199 | 208 | ||
| 200 | # Append time interval | 209 | # Append time interval |
| 201 | - self.gti.merge(obs_gti.tstart(), obs_gti.tstop()) | 210 | + self.m_gti.append(obs_gti.tstart(), obs_gti.tstop()) |
| 202 | 211 | ||
| 203 | # Increment global livetime and ontime | 212 | # Increment global livetime and ontime |
| 204 | ontime += obs.ontime() | 213 | ontime += obs.ontime() |
| @@ -269,18 +278,18 @@ class csobsinfo(ctools.cscript): | @@ -269,18 +278,18 @@ class csobsinfo(ctools.cscript): | ||
| 269 | # Log offset if possible | 278 | # Log offset if possible |
| 270 | if self.m_offset: | 279 | if self.m_offset: |
| 271 | offset = pnt_dir.dist_deg(self.obj_dir) | 280 | offset = pnt_dir.dist_deg(self.obj_dir) |
| 272 | - self.offsets.append(offset) | 281 | + self.m_offsets.append(offset) |
| 273 | self.log.parformat("Offset") | 282 | self.log.parformat("Offset") |
| 274 | self.log("%.2f" % (offset)) | 283 | self.log("%.2f" % (offset)) |
| 275 | self.log("\n") | 284 | self.log("\n") |
| 276 | else: | 285 | else: |
| 277 | - self.offsets.append(-1.0) | 286 | + self.m_offsets.append(-1.0) |
| 278 | 287 | ||
| 279 | # Retrieve zenith and azimuth and store for later use | 288 | # Retrieve zenith and azimuth and store for later use |
| 280 | zenith = obs.pointing().zenith() | 289 | zenith = obs.pointing().zenith() |
| 281 | azimuth = obs.pointing().azimuth() | 290 | azimuth = obs.pointing().azimuth() |
| 282 | - self.zeniths.append(zenith) | ||
| 283 | - self.azimuths.append(azimuth) | 291 | + self.m_zeniths.append(zenith) |
| 292 | + self.m_azimuths.append(azimuth) | ||
| 284 | 293 | ||
| 285 | # Log Zenith and Azimuth if required | 294 | # Log Zenith and Azimuth if required |
| 286 | if self.logExplicit(): | 295 | if self.logExplicit(): |
| @@ -315,15 +324,15 @@ class csobsinfo(ctools.cscript): | @@ -315,15 +324,15 @@ class csobsinfo(ctools.cscript): | ||
| 315 | # Log mean offset if possible | 324 | # Log mean offset if possible |
| 316 | if self.m_offset: | 325 | if self.m_offset: |
| 317 | self.log.parformat("Mean offset angle") | 326 | self.log.parformat("Mean offset angle") |
| 318 | - self.log("%.2f" % (sum(self.offsets)/len(self.offsets))) | 327 | + self.log("%.2f" % (sum(self.m_offsets)/len(self.m_offsets))) |
| 319 | self.log("\n") | 328 | self.log("\n") |
| 320 | 329 | ||
| 321 | # Log mean azimuth and zenith angle | 330 | # Log mean azimuth and zenith angle |
| 322 | self.log.parformat("Mean zenith angle") | 331 | self.log.parformat("Mean zenith angle") |
| 323 | - self.log("%.2f" % (sum(self.zeniths)/len(self.zeniths))) | 332 | + self.log("%.2f" % (sum(self.m_zeniths)/len(self.m_zeniths))) |
| 324 | self.log("\n") | 333 | self.log("\n") |
| 325 | self.log.parformat("Mean azimuth angle") | 334 | self.log.parformat("Mean azimuth angle") |
| 326 | - self.log("%.2f" % (sum(self.azimuths)/len(self.azimuths))) | 335 | + self.log("%.2f" % (sum(self.m_azimuths)/len(self.m_azimuths))) |
| 327 | self.log("\n") | 336 | self.log("\n") |
| 328 | 337 | ||
| 329 | # Log name of observations if requested | 338 | # Log name of observations if requested |
| @@ -338,25 +347,25 @@ class csobsinfo(ctools.cscript): | @@ -338,25 +347,25 @@ class csobsinfo(ctools.cscript): | ||
| 338 | # Log energy range | 347 | # Log energy range |
| 339 | self.log.header3("Energy range") | 348 | self.log.header3("Energy range") |
| 340 | self.log.parformat("Emin") | 349 | self.log.parformat("Emin") |
| 341 | - if self.ebounds.size() == 0: | 350 | + if self.m_ebounds.size() == 0: |
| 342 | self.log("undefined") | 351 | self.log("undefined") |
| 343 | else: | 352 | else: |
| 344 | - self.log(str(self.ebounds.emin())) | 353 | + self.log(str(self.m_ebounds.emin())) |
| 345 | self.log("\n") | 354 | self.log("\n") |
| 346 | self.log.parformat("Emax") | 355 | self.log.parformat("Emax") |
| 347 | - if self.ebounds.size() == 0: | 356 | + if self.m_ebounds.size() == 0: |
| 348 | self.log("undefined") | 357 | self.log("undefined") |
| 349 | else: | 358 | else: |
| 350 | - self.log(str(self.ebounds.emax())) | 359 | + self.log(str(self.m_ebounds.emax())) |
| 351 | self.log("\n\n") | 360 | self.log("\n\n") |
| 352 | 361 | ||
| 353 | # Log time range | 362 | # Log time range |
| 354 | self.log.header3("Time range") | 363 | self.log.header3("Time range") |
| 355 | self.log.parformat("Start [MJD]") | 364 | self.log.parformat("Start [MJD]") |
| 356 | - self.log(str(self.gti.tstart().mjd())) | 365 | + self.log(str(self.m_gti.tstart().mjd())) |
| 357 | self.log("\n") | 366 | self.log("\n") |
| 358 | self.log.parformat("Stop [MJD]") | 367 | self.log.parformat("Stop [MJD]") |
| 359 | - self.log(str(self.gti.tstop().mjd())) | 368 | + self.log(str(self.m_gti.tstop().mjd())) |
| 360 | self.log("\n") | 369 | self.log("\n") |
| 361 | 370 | ||
| 362 | # Log ontime and livetime in different units | 371 | # Log ontime and livetime in different units |
| @@ -406,31 +415,31 @@ class csobsinfo(ctools.cscript): | @@ -406,31 +415,31 @@ class csobsinfo(ctools.cscript): | ||
| 406 | """ | 415 | """ |
| 407 | Return zenith angles | 416 | Return zenith angles |
| 408 | """ | 417 | """ |
| 409 | - return self.zeniths | 418 | + return self.m_zeniths |
| 410 | 419 | ||
| 411 | def azimuths(self): | 420 | def azimuths(self): |
| 412 | """ | 421 | """ |
| 413 | Return azimuth angles | 422 | Return azimuth angles |
| 414 | """ | 423 | """ |
| 415 | - return self.azimuths | 424 | + return self.m_azimuths |
| 416 | 425 | ||
| 417 | def offsets(self): | 426 | def offsets(self): |
| 418 | """ | 427 | """ |
| 419 | Return offset angles | 428 | Return offset angles |
| 420 | """ | 429 | """ |
| 421 | - return self.offsets | 430 | + return self.m_offsets |
| 422 | 431 | ||
| 423 | def ebounds(self): | 432 | def ebounds(self): |
| 424 | """ | 433 | """ |
| 425 | Return energy boundaries | 434 | Return energy boundaries |
| 426 | """ | 435 | """ |
| 427 | - return self.ebounds | 436 | + return self.m_ebounds |
| 428 | 437 | ||
| 429 | def gti(self): | 438 | def gti(self): |
| 430 | """ | 439 | """ |
| 431 | Return good time intervals | 440 | Return good time intervals |
| 432 | """ | 441 | """ |
| 433 | - return self.gti | 442 | + return self.m_gti |
| 434 | 443 | ||
| 435 | def execute(self): | 444 | def execute(self): |
| 436 | """ | 445 | """ |
doc/source/reference_manual/csiactobs.rst
| @@ -67,7 +67,7 @@ General parameters | @@ -67,7 +67,7 @@ General parameters | ||
| 67 | ``(master_indx = master.json) [file]`` | 67 | ``(master_indx = master.json) [file]`` |
| 68 | Name of master index file. | 68 | Name of master index file. |
| 69 | 69 | ||
| 70 | -``(bkg_scale = yes) [boolean] | 70 | +``(bkg_scale = yes) [boolean]`` |
| 71 | Specifies whether the background scaling factor from the observation index | 71 | Specifies whether the background scaling factor from the observation index |
| 72 | file should be applied if available. | 72 | file should be applied if available. |
| 73 | 73 |
doc/source/reference_manual/csmodelmerge.rst
| @@ -12,10 +12,11 @@ Synopsis | @@ -12,10 +12,11 @@ Synopsis | ||
| 12 | This script merges several model XML file into a single model XML file. The | 12 | This script merges several model XML file into a single model XML file. The |
| 13 | files are passed via the ``inmodels`` parameter. There are 4 options to specify | 13 | files are passed via the ``inmodels`` parameter. There are 4 options to specify |
| 14 | a list of files that should be merged: | 14 | a list of files that should be merged: |
| 15 | -- a space-separated list of file names (e.g. model1.xml model2.xml) | ||
| 16 | -- a comma-separated list of file names (e.g. model1.xml,model2.xml) | ||
| 17 | -- a wildcard string (e.g. mymodel*.xml) | ||
| 18 | -- an ASCII file containing the file names, one file per line (e.g. @mymodels.txt) | 15 | + |
| 16 | +* a space-separated list of file names (e.g. model1.xml model2.xml) | ||
| 17 | +* a comma-separated list of file names (e.g. model1.xml,model2.xml) | ||
| 18 | +* a wildcard string (e.g. mymodel*.xml) | ||
| 19 | +* an ASCII file containing the file names, one file per line (e.g. @mymodels.txt) | ||
| 19 | 20 | ||
| 20 | Note that the "@" needs to be specified in case of an ASCII file name. It is | 21 | Note that the "@" needs to be specified in case of an ASCII file name. It is |
| 21 | also important to note that models with the same name cannot be appended. The | 22 | also important to note that models with the same name cannot be appended. The |
doc/source/reference_manual/cstsmapmerge.rst
| @@ -15,10 +15,11 @@ map again. Running :doc:`cttsmap` with the hidden parameters ``binmin`` and | @@ -15,10 +15,11 @@ map again. Running :doc:`cttsmap` with the hidden parameters ``binmin`` and | ||
| 15 | ``binmax`` specified will result in several sliced TS maps. This script merges | 15 | ``binmax`` specified will result in several sliced TS maps. This script merges |
| 16 | them into one FITS file. The input slices can be passed via the ``inmaps`` | 16 | them into one FITS file. The input slices can be passed via the ``inmaps`` |
| 17 | parameter. There are 4 options to specify a list of files that should be merged: | 17 | parameter. There are 4 options to specify a list of files that should be merged: |
| 18 | -- a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits) | ||
| 19 | -- a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits) | ||
| 20 | -- a wildcard string (e.g. tsmap*.xml) | ||
| 21 | -- an ASCII file containing the file names, one file per line (e.g. @mymaps.txt) | 18 | + |
| 19 | +* a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits) | ||
| 20 | +* a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits) | ||
| 21 | +* a wildcard string (e.g. tsmap*.xml) | ||
| 22 | +* an ASCII file containing the file names, one file per line (e.g. @mymaps.txt) | ||
| 22 | 23 | ||
| 23 | Note that the "@" needs to be specified in case of an ASCII file name. | 24 | Note that the "@" needs to be specified in case of an ASCII file name. |
| 24 | 25 |
doc/source/user_manual/iact_analysis/analysis.rst
0 → 100644
| 1 | +.. _sec_analysis: | ||
| 2 | + | ||
| 3 | +Analysing IACT data | ||
| 4 | +=================== | ||
| 5 | + | ||
| 6 | +In order to analyse data from current IACTs, the data have to be retrieved. Information how | ||
| 7 | +to copy IACT data to a local machine can be found here: :ref:`sec_copy` | ||
| 8 | + | ||
| 9 | +Check available FITS production | ||
| 10 | +------------------------------- | ||
| 11 | +Before starting an analysis, it is important to know what kind of data need to be analysed. | ||
| 12 | +In the current scheme of `storing IACT data <http://gamma-astro-data-formats.readthedocs.org/en/latest/data_storage/index.html>`_, | ||
| 13 | +every FITS data set is represented with a unique name, e.g. ``"hess-fits-pa-release-1.0-Prod26-MppStd"``. This string | ||
| 14 | +has to be passed to some set of ctools to find and use the data for analysis. | ||
| 15 | + | ||
| 16 | +To check the available names of FITS data sets, it is recommended to run :ref:`csiactdata`. | ||
| 17 | + | ||
| 18 | +.. code-block:: bash | ||
| 19 | + | ||
| 20 | + $ csiactdata | ||
| 21 | + Path were data is located [] /path/to/fits/data/ | ||
| 22 | + | ||
| 23 | +This will list names of available FITS data sets on the screen. In the | ||
| 24 | +following it is assumed the user has selected a data set called ``fits-data-name``. | ||
| 25 | +The query for the path where the data is located can be | ||
| 26 | +omitted by setting the environment variable ``VHEFITS`` to this locations. This might be convenient | ||
| 27 | +since some other tools also query for the same parameter and can use ``VHEFITS`` instead. | ||
| 28 | + | ||
| 29 | +Find observations | ||
| 30 | +----------------- | ||
| 31 | +To start an analysis, it is required to assemble a list of observations that should be used. For this purpose, | ||
| 32 | +:ref:`csfindobs` searches for observations according to user criteria. | ||
| 33 | + | ||
| 34 | +Search for sky pointings | ||
| 35 | +^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 36 | +The most common search for observations is by pointing position in the sky. | ||
| 37 | + | ||
| 38 | +.. code-block:: bash | ||
| 39 | + | ||
| 40 | + $ csfindobs | ||
| 41 | + Name of FITS production (Run csiactdata to view your options) [fits-data-name] | ||
| 42 | + Right ascension [83.6331] | ||
| 43 | + Declination [22.01] | ||
| 44 | + Search radius [2.5] | ||
| 45 | + Runlist outfile [runlist.lis] | ||
| 46 | + | ||
| 47 | + | ||
| 48 | +Specifying additional search criteria | ||
| 49 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 50 | +The search for observations can be further constrained arbitrarily. For this purpose, | ||
| 51 | +:ref:`csfindobs` contains a hidden parameter ``expression``. It is possible to use any | ||
| 52 | +FITS selection expression that is `supported by cfitsio <http://www.isdc.unige.ch/integral/download/osa/doc/10.1/osa_um_intro/node38.html>`_. | ||
| 53 | +Available properties (i.e. column names) for selection, can be found | ||
| 54 | +`here <http://gamma-astro-data-formats.readthedocs.org/en/latest/data_storage/obs_index/index.html>`_ | ||
| 55 | +or by simply browsing the observation index file with a FITS viewer. | ||
| 56 | + | ||
| 57 | +Example: | ||
| 58 | + | ||
| 59 | +.. code-block:: bash | ||
| 60 | + | ||
| 61 | + $ csfindobs expression="ZEN_PNT<30&&LIVETIME>1800" | ||
| 62 | + ... | ||
| 63 | + | ||
| 64 | +This command selects all observations that have a zenith angle less than | ||
| 65 | +30 degrees and a livetime larger than 1800 seconds (i.e. 30 minutes). | ||
| 66 | + | ||
| 67 | + | ||
| 68 | +Omitting the pointing selection | ||
| 69 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 70 | +You can also omit the pointing selection, for example in case where you | ||
| 71 | +don't know the sky coordinates of your object of interest. In that case, you | ||
| 72 | +may directly use the object name in the expression and pass ``NONE`` as | ||
| 73 | +coordinate to :ref:`csfindobs`: | ||
| 74 | + | ||
| 75 | +.. code-block:: bash | ||
| 76 | + | ||
| 77 | + $ csfindobs expression="OBJECT=='Crab Nebula'" | ||
| 78 | + Name of FITS production (Run csiactdata to view your options) [fits-data-name] | ||
| 79 | + Right ascension [NONE] | ||
| 80 | + Runlist outfile [runlist.lis] | ||
| 81 | + | ||
| 82 | +This will omit the selection of pointing positions. | ||
| 83 | + | ||
| 84 | +.. note:: | ||
| 85 | + | ||
| 86 | + By default, :ref:`csfindobs` only selects data of highest quality (i.e. QUALITY=0). | ||
| 87 | + You may overwrite this default by specifying the hidden parameter | ||
| 88 | + ``min_qual``. For example, ``min_qual=1`` selects all data with a | ||
| 89 | + looser quality criteria. | ||
| 90 | + | ||
| 91 | +Create an observation list | ||
| 92 | +-------------------------- | ||
| 93 | +The runlist ASCII file containing a list of selected observation IDs must now be converted to an | ||
| 94 | +observation definition XML file. This file contains information about the location of the files that are required | ||
| 95 | +for the analysis. The purpose of the tool :ref:`csiactobs` is to do this conversion. | ||
| 96 | + | ||
| 97 | +.. code-block:: bash | ||
| 98 | + | ||
| 99 | + $ csiactobs | ||
| 100 | + Data storage name (Run csiactdata to view your options) [] fits-data-name | ||
| 101 | + Runlist file [runlist.lis] | ||
| 102 | + Number of free parameters per background model [1] | ||
| 103 | + Output model XML file [bgmodels.xml] | ||
| 104 | + Observation XML outfile [obs.xml] | ||
| 105 | + | ||
| 106 | +This tools has now performed various tasks: | ||
| 107 | + | ||
| 108 | +* Use the IACT data storage to convert the observation list ``runlist.lis`` into an observation definition XML file ``obs.xml``. | ||
| 109 | +* Create a model XML file that contains the background models. In ctools, each observation has its own independent background model. The number of free parameters per model was set to 1, i.e. the normalisation of the background is left free for each observation. The models were saved in ``bgmodels.xml``. | ||
| 110 | +* The tool uses an internal hierarchy for assigning the background models. The hidden parameter ``bkg_mod_hiera`` (default=irf|aeff|gauss) steers the order how background models should be created. In case the IRF background model is not available, the script automatically falls back to the background model from effective area (GCTAModelAeffBackground). | ||
| 111 | +* There are some further hidden parameters to steer the start parameters for the Aeff and Gaussian background model. Have a look at :ref:`csiactobs` to view a full list of parameters. | ||
| 112 | +* In ``csiactobs.log`` (or on screen if ``debug=yes``), the script logs the complete energy range of the observation container. These values might be important for later usage (e.g. in binned analysis). | ||
| 113 | + | ||
| 114 | +In case a sky model is already prepared, it is possible to also provide the hidden parameter ``inmodel``. The output | ||
| 115 | +model XML file will then contain both, the background model and the input sky model: | ||
| 116 | + | ||
| 117 | +.. code-block:: bash | ||
| 118 | + | ||
| 119 | + $ csiactobs inmodel="mymodel.xml" | ||
| 120 | + | ||
| 121 | +Alternatively, models can be merged at any times using the simple tool :ref:`csmodelmerge`: | ||
| 122 | + | ||
| 123 | +.. code-block:: bash | ||
| 124 | + | ||
| 125 | + $ csmodelmerge | ||
| 126 | + Input model XML files [bgmodels.xml crab.xml] | ||
| 127 | + Output model file [crab_models.xml] | ||
| 128 | + | ||
| 129 | +Note that the number of files to merge is not limited to two. Detailled options how the input model XML file can | ||
| 130 | +be passed is given on the reference page of :ref:`csmodelmerge`. It is also important to know that each model in the container | ||
| 131 | +must have a unique name. This implies, for instance, merging the same XML model twice will result in an exception. | ||
| 132 | + | ||
| 133 | +A list of available sky models can be found `here <http://gammalib.sourceforge.net/user_manual/modules/model.html>`_. | ||
| 134 | +Of particular help to create sky models for your dataset is the section about :ref:`modelling CTA data <models>`. | ||
| 135 | + | ||
| 136 | + | ||
| 137 | +Example XML files | ||
| 138 | +----------------- | ||
| 139 | +To get familiar with the XML syntax and format, example files for an observation container and a model container | ||
| 140 | +are shown in the following. | ||
| 141 | + | ||
| 142 | +Observation XML file | ||
| 143 | +^^^^^^^^^^^^^^^^^^^^ | ||
| 144 | + | ||
| 145 | +.. code-block:: xml | ||
| 146 | + | ||
| 147 | + <?xml version="1.0" encoding="UTF-8" standalone="no"?> | ||
| 148 | + <observation_list title="observation list"> | ||
| 149 | + <observation name="Crab Nebula" id="11111" instrument="HESS"> | ||
| 150 | + <parameter name="EventList" file="/path/to/fits/file/events_11111.fits.gz" /> | ||
| 151 | + <parameter name="EffectiveArea" file="/path/to/fits/file/aeff_11111.fits.gz" /> | ||
| 152 | + <parameter name="PointSpreadFunction" file="/path/to/fits/file/psf_11111.fits.gz" /> | ||
| 153 | + <parameter name="EnergyDispersion" file="/path/to/fits/file/edisp_11111.fits.gz" /> | ||
| 154 | + <parameter name="Background" file="/path/to/fits/file/bgmodel_11111.fits.gz" /> | ||
| 155 | + </observation> | ||
| 156 | + <observation name="Crab Nebula" id="11112" instrument="HESS"> | ||
| 157 | + <parameter name="EventList" file="/path/to/fits/file/events_11112.fits.gz" /> | ||
| 158 | + <parameter name="EffectiveArea" file="/path/to/fits/file/aeff_11112.fits.gz" /> | ||
| 159 | + <parameter name="PointSpreadFunction" file="/path/to/fits/file/psf_11112.fits.gz" /> | ||
| 160 | + <parameter name="EnergyDispersion" file="/path/to/fits/file/edisp_11112.fits.gz" /> | ||
| 161 | + <parameter name="Background" file="/path/to/fits/file/bgmodel_11112.fits.gz" /> | ||
| 162 | + </observation> | ||
| 163 | + <observation name="Crab Nebula" id="11113" instrument="HESS"> | ||
| 164 | + <parameter name="EventList" file="/path/to/fits/file/events_11113.fits.gz" /> | ||
| 165 | + <parameter name="EffectiveArea" file="/path/to/fits/file/aeff_11113.fits.gz" /> | ||
| 166 | + <parameter name="PointSpreadFunction" file="/path/to/fits/file/psf_11113.fits.gz" /> | ||
| 167 | + <parameter name="EnergyDispersion" file="/path/to/fits/file/edisp_11113.fits.gz" /> | ||
| 168 | + <parameter name="Background" file="/path/to/fits/file/bgmodel_11113.fits.gz" /> | ||
| 169 | + </observation> | ||
| 170 | + <observation name="Crab Nebula" id="11114" instrument="HESS"> | ||
| 171 | + <parameter name="EventList" file="/path/to/fits/file/events_11114.fits.gz" /> | ||
| 172 | + <parameter name="EffectiveArea" file="/path/to/fits/file/aeff_11114.fits.gz" /> | ||
| 173 | + <parameter name="PointSpreadFunction" file="/path/to/fits/file/psf_11114.fits.gz" /> | ||
| 174 | + <parameter name="EnergyDispersion" file="/path/to/fits/file/edisp_11114.fits.gz" /> | ||
| 175 | + <parameter name="Background" file="/path/to/fits/file/bgmodel_11114.fits.gz" /> | ||
| 176 | + </observation> | ||
| 177 | + </observation_list> | ||
| 178 | + | ||
| 179 | +Model XML file | ||
| 180 | +^^^^^^^^^^^^^^ | ||
| 181 | + | ||
| 182 | +.. code-block:: xml | ||
| 183 | + | ||
| 184 | + <?xml version="1.0" encoding="UTF-8" standalone="no"?> | ||
| 185 | + <source_library title="source library"> | ||
| 186 | + <source name="bkg_11111" type="CTAIrfBackground" instrument="HESS" id="11111"> | ||
| 187 | + <spectrum type="ConstantValue"> | ||
| 188 | + <parameter name="Value" value="1" error="0" scale="1" min="0.01" max="100" free="1" /> | ||
| 189 | + </spectrum> | ||
| 190 | + </source> | ||
| 191 | + <source name="bkg_11112" type="CTAIrfBackground" instrument="HESS" id="11112"> | ||
| 192 | + <spectrum type="ConstantValue"> | ||
| 193 | + <parameter name="Value" value="1" error="0" scale="1" min="0.01" max="100" free="1" /> | ||
| 194 | + </spectrum> | ||
| 195 | + </source> | ||
| 196 | + <source name="bkg_11113" type="CTAIrfBackground" instrument="HESS" id="11113"> | ||
| 197 | + <spectrum type="ConstantValue"> | ||
| 198 | + <parameter name="Value" value="1" error="0" scale="1" min="0.01" max="100" free="1" /> | ||
| 199 | + </spectrum> | ||
| 200 | + </source> | ||
| 201 | + <source name="bkg_11114" type="CTAIrfBackground" instrument="HESS" id="11114"> | ||
| 202 | + <spectrum type="ConstantValue"> | ||
| 203 | + <parameter name="Value" value="1" error="0" scale="1" min="0.01" max="100" free="1" /> | ||
| 204 | + </spectrum> | ||
| 205 | + </source> | ||
| 206 | + <source name="Crab" type="PointSource"> | ||
| 207 | + <spectrum type="PowerLaw"> | ||
| 208 | + <parameter name="Prefactor" scale="1e-17" value="3.0" min="1e-07" max="1000.0" free="1"/> | ||
| 209 | + <parameter name="Index" scale="-1" value="2.5" min="0.0" max="+5.0" free="1"/> | ||
| 210 | + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/> | ||
| 211 | + </spectrum> | ||
| 212 | + <spatialModel type="SkyDirFunction"> | ||
| 213 | + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="1"/> | ||
| 214 | + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="1"/> | ||
| 215 | + </spatialModel> | ||
| 216 | + </source> | ||
| 217 | + </source_library> | ||
| 218 | + | ||
| 219 | +.. note:: | ||
| 220 | + | ||
| 221 | + It is important to ensure background models are properly linked to their respective observation. | ||
| 222 | + Therefore it is required to keep the attributes ``instrument`` and ``id`` the same for the observation | ||
| 223 | + and the corresponding background model. The tool :ref:`csiactobs` assures this automatically. | ||
| 224 | + | ||
| 225 | +Run ctselect | ||
| 226 | +------------ | ||
| 227 | +To prepare the data for analysis, cuts have to be applied to the event data. The selection is performed by :ref:`ctselect`. | ||
| 228 | +This tool writes out selected event lists into the local directory. If the observation XML file contains several runs, it is recommended | ||
| 229 | +to first create a separate folder and specify this folder in the hidden ``prefix`` argument. | ||
| 230 | + | ||
| 231 | +.. code-block:: bash | ||
| 232 | + | ||
| 233 | + $ mkdir selected | ||
| 234 | + | ||
| 235 | +.. code-block:: bash | ||
| 236 | + | ||
| 237 | + $ ctselect usethres=DEFAULT usepnt=yes prefix=selected/selected_ | ||
| 238 | + Input event list or observation definition file [events.fits] obs.xml | ||
| 239 | + Lower energy limit (TeV) [0.1] | ||
| 240 | + Upper energy limit (TeV) [100.0] | ||
| 241 | + Radius of ROI (degrees) (0-180) [3.0] 2.5 | ||
| 242 | + Start time (CTA MET in seconds) [0.0] | ||
| 243 | + End time (CTA MET in seconds) [0.0] | ||
| 244 | + Output event list or observation definition file [selected_events.fits] selected_obs.xml | ||
| 245 | + | ||
| 246 | +For IACT analysis, it is recommended to use the hidden parameter ``usethres="DEFAULT"``. This instructs :ref:`ctselect` | ||
| 247 | +to extract the safe energy range from the instrument response functions and apply them to the data. This safe energy range | ||
| 248 | +is thus superior to the energy limit passed via the user parameters. In addition, to analyse the complete field of view, | ||
| 249 | +the parameter ``usepnt=yes`` uses, for each observation, the pointing position as centre for the selection radius. | ||
| 250 | +The radius parameter is dependent on the intrument, for an instrument with a 5 degree field of view, a radius of 2.5 degrees | ||
| 251 | +seems reasonable. The time selection is not applied in the above example; | ||
| 252 | +specifying 0 as start and end time skips the time selection. | ||
| 253 | +For time-resolved analysis, it is important to know the MET time that is required to extract. | ||
| 254 | +The result of the selection step is written into the observation XML file ``selected_obs.xml``, which now contains references | ||
| 255 | +to the new selected event FITS files. | ||
| 256 | + | ||
| 257 | +Unbinned analysis | ||
| 258 | +----------------- | ||
| 259 | +Once the data is selected, the easiest way to analyse is an unbinned analysis. Note that the input model XML | ||
| 260 | +file must now contain the background models and source components to describe the field of view. | ||
| 261 | + | ||
| 262 | +.. code-block:: bash | ||
| 263 | + | ||
| 264 | + $ ctlike | ||
| 265 | + Event list, counts cube or observation definition file [events.fits] selected_obs.xml | ||
| 266 | + Source model [$CTOOLS/share/models/crab.xml] crab_models.xml | ||
| 267 | + Source model output file [crab_results.xml] | ||
| 268 | + | ||
| 269 | +The result of the fit was stored in ``crab_results.xml``. Note that fitted parameters, ``Prefactors`` in particular, | ||
| 270 | +typically use MeV as energy unit. To monitor the progress of the fit on the screen, one can simply run with the option ``debug=yes``. | ||
| 271 | +Alternatively, the logfile ``ctlike.log`` can be inspected after the fit. | ||
| 272 | + | ||
| 273 | +On default, energy dispersion is not considered in the fit. To switch on the usage of the energy migration matrix, | ||
| 274 | +the hidden parameter ``edisp=yes`` can be provided. Note that this will cause a significant reduction of the computing | ||
| 275 | +speed. | ||
| 276 | + | ||
| 277 | +Stacked (binned) analysis | ||
| 278 | +------------------------- | ||
| 279 | +The stacked analysis mode is using a binned analysis where all observations are included and stacked into one event cube. | ||
| 280 | +This analysis mode is much faster than unbinned analysis when having a large dataset (e.g. > 100 hours). | ||
| 281 | +For this type of binned analysis, some intermediate data products have to be produced. The products are a binned data cube, | ||
| 282 | +an exposure cube, a PSF cube, and a background cube. | ||
| 283 | + | ||
| 284 | +Bin events | ||
| 285 | +^^^^^^^^^^ | ||
| 286 | + | ||
| 287 | +.. code-block:: bash | ||
| 288 | + | ||
| 289 | + $ ctbin | ||
| 290 | + Event list or observation definition file [events.fits] selected_obs.xml | ||
| 291 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 292 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 293 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 294 | + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] | ||
| 295 | + Image scale (in degrees/pixel) [0.02] | ||
| 296 | + Size of the X axis in pixels [200] | ||
| 297 | + Size of the Y axis in pixels [200] | ||
| 298 | + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] | ||
| 299 | + Start value for first energy bin in TeV [0.1] 0.5 | ||
| 300 | + Stop value for last energy bin in TeV [100.0] 50 | ||
| 301 | + Number of energy bins [20] | ||
| 302 | + Output counts cube [cntcube.fits] | ||
| 303 | + | ||
| 304 | +Note that bins only get filled if the bin of the cube is fully contained in the energy range and RoI of a considered observation. | ||
| 305 | +It is therefore useful to provide the energy range given by :ref:`csiactobs` above. This ensures a maximum agreement between | ||
| 306 | +observations and binning and reduces the loss of data. | ||
| 307 | + | ||
| 308 | +Create exposure cube | ||
| 309 | +^^^^^^^^^^^^^^^^^^^^ | ||
| 310 | +After binning the events into a three-dimensional cube, an exposure cube has to be computed. | ||
| 311 | +The exposure is defined as the effective area times the dead-time corrected observation time. | ||
| 312 | +Each observation from the input container gets stacked in the resulting cube. The exposure is stored | ||
| 313 | +in units of :math:`cm^2 s`. The exposure cube does not have to contain the same binning as the event cube | ||
| 314 | +but for simplicity, the event cube can be passed to adopt the binning parameters. Note, however, that the | ||
| 315 | +exposure cube is defined in true sky coordinates and energy while the counts cube is defined in reconstructed | ||
| 316 | +sky coordinates and energy. Consequently, the sky area and energy range covered by the exposure cube should be | ||
| 317 | +slightly larger than that of the counts cube to accommodate for spill over of events due to the point spread function | ||
| 318 | +and energy dispersion. By computing the exposure cube on the same grid as the counts cube, the spill over of events | ||
| 319 | +from sources at the edge of cube will not be handled correctly. In the example below, however, no source at the edge of the | ||
| 320 | +field of view is present. Therefore, for simplicity, the count cube is used as input cube to extract the binning. | ||
| 321 | + | ||
| 322 | +This task of computing the exposure cube is is performed by :ref:`ctexpcube`. | ||
| 323 | + | ||
| 324 | +.. code-block:: bash | ||
| 325 | + | ||
| 326 | + $ ctexpcube | ||
| 327 | + Event list or observation definition file [NONE] selected_obs.xml | ||
| 328 | + Input counts cube file to extract exposure cube definition [NONE] cntcube.fits | ||
| 329 | + Output exposure cube file [expcube.fits] | ||
| 330 | + | ||
| 331 | +Alternatively, the exposure cube can be created with different binning than the event cube: | ||
| 332 | + | ||
| 333 | +.. code-block:: bash | ||
| 334 | + | ||
| 335 | + $ ctexpcube | ||
| 336 | + Input event list or observation definition XML file [NONE] selected_obs.xml | ||
| 337 | + Input counts cube file to extract exposure cube definition [NONE] | ||
| 338 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 339 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 340 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 341 | + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] | ||
| 342 | + Image scale (in degrees/pixel) [0.02] 0.04 | ||
| 343 | + Size of the X axis in pixels [200] 100 | ||
| 344 | + Size of the Y axis in pixels [200] 100 | ||
| 345 | + Lower energy limit (TeV) [0.5] | ||
| 346 | + Upper energy limit (TeV) [50.0] | ||
| 347 | + Number of energy bins [20] 30 | ||
| 348 | + Output exposure cube file [expcube.fits] | ||
| 349 | + | ||
| 350 | +Create PSF cube | ||
| 351 | +^^^^^^^^^^^^^^^ | ||
| 352 | +As a next step for the binned analysis, a cube containing the point spread function (PSF) must be computed. | ||
| 353 | +Since the PSF cannot be stored by one single parameter, the PSF cube computed by :ref:`ctpsfcube` has a fourth dimension. | ||
| 354 | +In each bin of the cube, the PSF is stored as a function of offset from source. The granularity of the | ||
| 355 | +PSF histogram is determined by the hidden parameter ``anumbins`` (default: 200). Therefore, when passing | ||
| 356 | +the event cube to adopt the sky binning for the PSF cube, the resulting FITS file can become quite large due to | ||
| 357 | +the fourth dimension. Usually in IACT analysis, the PSF doesn't change too dramatically across the field of view. | ||
| 358 | +Therefore the user can think about reducing the spatial binning of the PSF cube: | ||
| 359 | + | ||
| 360 | +.. code-block:: bash | ||
| 361 | + | ||
| 362 | + $ ctpsfcube | ||
| 363 | + Input event list or observation definition XML file [NONE] selected_obs.xml | ||
| 364 | + Input counts cube file to extract PSF cube definition [NONE] | ||
| 365 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 366 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 367 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 368 | + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] | ||
| 369 | + Image scale (in degrees/pixel) [1.0] 0.2 | ||
| 370 | + Size of the X axis in pixels [10] 20 | ||
| 371 | + Size of the Y axis in pixels [10] 20 | ||
| 372 | + Lower energy limit (TeV) [0.1] 0.5 | ||
| 373 | + Upper energy limit (TeV) [100.0] 50 | ||
| 374 | + Number of energy bins [20] | ||
| 375 | + Output PSF cube file [psfcube.fits] | ||
| 376 | + | ||
| 377 | +Depending on the required PSF precision, one could reduce the number of offset bins via the hidden parameter | ||
| 378 | +``anumbins``: | ||
| 379 | + | ||
| 380 | +.. code-block:: bash | ||
| 381 | + | ||
| 382 | + $ ctpsfcube anumbins=100 | ||
| 383 | + | ||
| 384 | +Create background cube | ||
| 385 | +^^^^^^^^^^^^^^^^^^^^^^ | ||
| 386 | +Last but not least, for a binned IACT analysis a cube containing the background rate in sky coordinates and | ||
| 387 | +reconstructed energy has to be computed. This task is performed by :ref:`ctbkgcube`. The binning here can also differ | ||
| 388 | +from the event cube. For simplicity, however, the example below uses the event cube to adopt the binning. | ||
| 389 | + | ||
| 390 | +.. code-block:: bash | ||
| 391 | + | ||
| 392 | + $ ctbkgcube debug=yes | ||
| 393 | + Input event list or observation definition XML file [NONE] selected_obs.xml | ||
| 394 | + Input counts cube file to extract background cube definition [NONE] cntcube.fits | ||
| 395 | + Input model XML file [NONE] crab_models.xml | ||
| 396 | + Output background cube file [bkgcube.fits] | ||
| 397 | + Output model XML file [NONE] binned_models.xml | ||
| 398 | + | ||
| 399 | +Note that this tool also requires the parameters of an input and output model. In the model XML file that came out of | ||
| 400 | +:ref:`csiactobs`, one background model per observation is included. This models get merged and averaged in the background | ||
| 401 | +sky cube. In the output model the background models per obseervation will be removed. Instead, a global background model | ||
| 402 | +for the newly created background cube is included. Sky models present in the input model XML file will also be included in the | ||
| 403 | +new XML file, which subsequently can be used for binned :ref:`ctlike`. | ||
| 404 | + | ||
| 405 | +Example for stacked model XML file | ||
| 406 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 407 | +The output model of :ref:`ctbkgcube` looks the following: | ||
| 408 | + | ||
| 409 | +.. code-block:: xml | ||
| 410 | + | ||
| 411 | + <?xml version="1.0" encoding="UTF-8" standalone="no"?> | ||
| 412 | + <source_library title="source library"> | ||
| 413 | + <source name="BackgroundModel" type="CTACubeBackground" instrument="CTA,HESS,MAGIC,VERITAS"> | ||
| 414 | + <spectrum type="PowerLaw"> | ||
| 415 | + <parameter name="Prefactor" value="1" error="0" scale="1" min="0.01" max="100" free="1" /> | ||
| 416 | + <parameter name="Index" value="0" error="0" scale="1" min="-5" max="5" free="1" /> | ||
| 417 | + <parameter name="Scale" value="1" scale="1e+06" free="0" /> | ||
| 418 | + </spectrum> | ||
| 419 | + </source> | ||
| 420 | + <source name="Crab" type="PointSource"> | ||
| 421 | + <spectrum type="PowerLaw"> | ||
| 422 | + <parameter name="Prefactor" scale="1e-17" value="3.0" min="1e-07" max="1000.0" free="1"/> | ||
| 423 | + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/> | ||
| 424 | + <parameter name="Scale" scale="1e6" value="1." min="0.01" max="1000.0" free="0"/> | ||
| 425 | + </spectrum> | ||
| 426 | + <spatialModel type="SkyDirFunction"> | ||
| 427 | + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/> | ||
| 428 | + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/> | ||
| 429 | + </spatialModel> | ||
| 430 | + </source> | ||
| 431 | + </source_library> | ||
| 432 | + | ||
| 433 | +The background model with ``type=CTACubeBackground`` is used to scale the background cube stored in the FITS file | ||
| 434 | +created by :ref:`ctbkgcube`. | ||
| 435 | + | ||
| 436 | +Run ctlike | ||
| 437 | +^^^^^^^^^^ | ||
| 438 | +Having all the intermediate data products ready, a binned analysis can be conducted using :ref:`ctlike`. | ||
| 439 | + | ||
| 440 | +.. code-block:: bash | ||
| 441 | + | ||
| 442 | + $ ctlike | ||
| 443 | + Input event list, counts cube or observation definition XML file [events.fits] cntcube.fits | ||
| 444 | + Input exposure cube file (only needed for stacked analysis) [NONE] expcube.fits | ||
| 445 | + Input PSF cube file (only needed for stacked analysis) [NONE] psfcube.fits | ||
| 446 | + Input background cube file (only needed for stacked analysis) [NONE] bkgcube.fits | ||
| 447 | + Input model XML file [binned_models.xml] | ||
| 448 | + Output model XML file [binned_results.xml] | ||
| 449 | + | ||
| 450 | +Note that when passing an event cube to :ref:`ctlike`, the tool behaves differently than in unbinned mode. | ||
| 451 | +It queries directly for the additional ingredients for the binned analysis. It is important to pass the background model | ||
| 452 | +generated by :ref:`ctbkgcube` here to ensure the proper modelling of the background in the fit. | ||
| 453 | +The usage of energy dispersion cannot be included in binned analyses yet. Therefore, the hidden parameter ``edisp=yes`` | ||
| 454 | +has no effect in this case. | ||
| 455 | + | ||
| 456 | + | ||
| 457 | + | ||
| 458 | + |
doc/source/user_manual/iact_analysis/copy.rst
0 → 100644
| 1 | +.. _sec_copy: | ||
| 2 | + | ||
| 3 | +Copying IACT data | ||
| 4 | +================== | ||
| 5 | + | ||
| 6 | +Before you start | ||
| 7 | +---------------- | ||
| 8 | +Currently the way how IACT data is stored is quite prototypical. The ctools can work | ||
| 9 | +with IACT data that is stored following the `specifications here <http://gamma-astro-data-formats.readthedocs.org/en/latest/data_storage/index.html>`_. | ||
| 10 | +The script :ref:`csiactcopy` can be used to transfer these data from one system to another. | ||
| 11 | +IACT data is the property of collaborations (such as HESS, MAGIC or VERITAS) and not publicly available. If you are member of one of these | ||
| 12 | +collaborations, you should have access to a central place where the FITS data is stored and maintained. | ||
| 13 | + | ||
| 14 | +Mount remote file system | ||
| 15 | +------------------------ | ||
| 16 | +Note that :ref:`csiactcopy` does not include functionality to establish a conection to a remote file system. | ||
| 17 | +For the moment, it is capable of copying data from one location to another. Therefore it is necessary to | ||
| 18 | +mount the remote file system where the data is located on the local machine. This can be achieved using e.g. | ||
| 19 | +`sshfs <https://github.com/libfuse/sshfs>`_: | ||
| 20 | + | ||
| 21 | +.. code-block:: bash | ||
| 22 | + | ||
| 23 | + $ sshfs user@remote.server.com:/path/to/remote/fits/data/ /path/to/mountpoint/ | ||
| 24 | + | ||
| 25 | +The success of this operation can be verified by checking if there is a file called ``master.json`` | ||
| 26 | +in the mounted directory. | ||
| 27 | + | ||
| 28 | +Copy data | ||
| 29 | +--------- | ||
| 30 | +After mounting the file system to a folder of your choice (``/path/to/mountpoint/`` in this case), the copying of the data can be initiated. | ||
| 31 | +To monitor the progress on the screen, the hidden parameter ``debug=yes`` can be used. To increase the chattiness of the output, | ||
| 32 | +the hidden parameter ``chatter`` can e.g. be set to 3. | ||
| 33 | + | ||
| 34 | +.. code-block:: bash | ||
| 35 | + | ||
| 36 | + $ csiactcopy debug=yes chatter=3 | ||
| 37 | + Location of remote master file [/path/to/mountpoint/master.json] | ||
| 38 | + Name of FITS production to download [iact-fits] | ||
| 39 | + Destination path of FITS data [/path/to/local/fits/data/] | ||
| 40 | + | ||
| 41 | +The name of the FITS production must be available in the remote | ||
| 42 | +``master.json`` file. In case the names of the available productions are unkown | ||
| 43 | +there are two options to determine an appropriate name: | ||
| 44 | + | ||
| 45 | +* Run :ref:`csiactdata` and set as ``datapath`` the mountpoint | ||
| 46 | +* Try any name, :ref:`csiactcopy` will throw a ``RuntimeError`` with the available production names if the input name is not found | ||
| 47 | + | ||
| 48 | +Copy only a subset of the data | ||
| 49 | +------------------------------ | ||
| 50 | +If the user is interested only in a small subset of the remote FITS production, | ||
| 51 | +it is also possible to download only parts of the data. The user simply has to provide an | ||
| 52 | +ASCII file containing a list of observation IDs to be downloaded (one ID per line). The hidden | ||
| 53 | +parameter ``runlist`` then needs to be specified. | ||
| 54 | + | ||
| 55 | +.. code-block:: bash | ||
| 56 | + | ||
| 57 | + $ csiactcopy runlist=myrunlist.txt debug=yes | ||
| 58 | + | ||
| 59 | +Troubleshooting | ||
| 60 | +--------------- | ||
| 61 | +In case a download was interupted (e.g. the remote file system was disconnected), the script can simply be rerun. | ||
| 62 | +The index files which list the files available on the system are refreshed after the copying process. | ||
| 63 | +The hidden boolean parameter ``clobber`` (default=yes) specifies if the user wishes that local content should be overwritten with | ||
| 64 | +remote content. In turn, specifying ``clobber=no`` can speed up the process since remote content does not need to be copied | ||
| 65 | +if it is already available on the user machine. | ||
| 66 | + | ||
| 67 | +To verify the copying process was successful, it is recommended to run :ref:`csiactdata` after the process in order to list available productions | ||
| 68 | +on the local machine. The output should now include the copied FITS production name. |
doc/source/user_manual/iact_analysis/high_level.rst
0 → 100644
| 1 | +.. _sec_high_level: | ||
| 2 | + | ||
| 3 | +High level analysis tools for IACT data | ||
| 4 | +======================================= | ||
| 5 | + | ||
| 6 | +This section describes and explains how to create high level output. The following tools are explained: | ||
| 7 | + | ||
| 8 | +* :ref:`csobsinfo` (inspect and observation XML files) | ||
| 9 | +* :ref:`csmodelinfo` (inspect model XML files) | ||
| 10 | +* :ref:`ctulimit` (compute upper limits) | ||
| 11 | +* :ref:`cterror` (compute asymmetric errors) | ||
| 12 | +* :ref:`csspec` (compute spectral points) | ||
| 13 | +* :ref:`cslightcrv` (compute a light curve) | ||
| 14 | +* :ref:`csresmap` (compute a residual map) | ||
| 15 | +* :ref:`cttsmap` (compute a test statistics map) | ||
| 16 | + | ||
| 17 | +Inspecting observation definition files | ||
| 18 | +--------------------------------------- | ||
| 19 | +After having created an observation definition file, the user might want to know a summary of its content. | ||
| 20 | +Important quantities usually are the number of observations, the total livetime, or the mean zenith angle. | ||
| 21 | +For this purpose, :ref:`csobsinfo` inspects an observation XML container for its content. | ||
| 22 | + | ||
| 23 | +.. code-block:: bash | ||
| 24 | + | ||
| 25 | + $ csobsinfo | ||
| 26 | + Event list, counts cube, or observation definition file [obs.xml] selected_obs.xml | ||
| 27 | + | ||
| 28 | +This dumps general information of the observation container in a logfile ``csobsinfo.log``. | ||
| 29 | +For a direct output on the screen, pass the hidden parameter ``debug=yes`` on the command line: | ||
| 30 | + | ||
| 31 | +.. code-block:: bash | ||
| 32 | + | ||
| 33 | + $ csobsinfo debug=yes | ||
| 34 | + Event list, counts cube, or observation definition file [obs.xml] selected_obs.xml | ||
| 35 | + ... | ||
| 36 | + 2016-02-19T15:43:08: +=========+ | ||
| 37 | + 2016-02-19T15:43:08: | Summary | | ||
| 38 | + 2016-02-19T15:43:08: +=========+ | ||
| 39 | + 2016-02-19T15:43:08: Unbinned observations .....: 4 | ||
| 40 | + 2016-02-19T15:43:08: Number of events ..........: 6085 | ||
| 41 | + 2016-02-19T15:43:08: Binned observations .......: 0 | ||
| 42 | + 2016-02-19T15:43:08: Number of bins ............: 0 | ||
| 43 | + 2016-02-19T15:43:08: | ||
| 44 | + 2016-02-19T15:43:08: === Pointings === | ||
| 45 | + 2016-02-19T15:43:08: Mean zenith angle .........: 46.89 | ||
| 46 | + 2016-02-19T15:43:08: Mean azimuth angle ........: 13.63 | ||
| 47 | + 2016-02-19T15:43:08: | ||
| 48 | + 2016-02-19T15:43:08: === Energy range === | ||
| 49 | + 2016-02-19T15:43:08: Emin ......................: 461.48 GeV | ||
| 50 | + 2016-02-19T15:43:08: Emax ......................: 37.8832 TeV | ||
| 51 | + 2016-02-19T15:43:08: | ||
| 52 | + 2016-02-19T15:43:08: === Time range === | ||
| 53 | + 2016-02-19T15:43:08: Start [MJD] ...............: 53343.9212153 | ||
| 54 | + 2016-02-19T15:43:08: Stop [MJD] ................: 53347.9321065 | ||
| 55 | + 2016-02-19T15:43:08: Total ontime ..............: 6877.00 s = 114.62 min = 1.91 h | ||
| 56 | + 2016-02-19T15:43:08: Total livetime ............: 6427.52 s = 107.13 min = 1.79 h | ||
| 57 | + ... | ||
| 58 | + | ||
| 59 | +Further options are described on the reference page :ref:`csobsinfo`. | ||
| 60 | + | ||
| 61 | +Inspecting model XML files | ||
| 62 | +-------------------------- | ||
| 63 | +Similar to the observation containers, model XML file can be inspected, too. | ||
| 64 | +The tool :ref:`csmodelinfo` gives a summary of the model container. In particular, | ||
| 65 | +the number of free parameters, number of sky and background models and also number of | ||
| 66 | +parameters that are at their limits might be of interest. | ||
| 67 | + | ||
| 68 | +.. code-block:: bash | ||
| 69 | + | ||
| 70 | + $ csmodelinfo debug=yes | ||
| 71 | + Input model XML file [$CTOOLS/share/models/crab.xml] | ||
| 72 | + 2016-02-19T15:50:40: +=========+ | ||
| 73 | + 2016-02-19T15:50:40: | Summary | | ||
| 74 | + 2016-02-19T15:50:40: +=========+ | ||
| 75 | + 2016-02-19T15:50:40: === Instrument specific models === | ||
| 76 | + 2016-02-19T15:50:40: All .......................: 1 | ||
| 77 | + 2016-02-19T15:50:40: CTA .......................: 1 | ||
| 78 | + 2016-02-19T15:50:40: === Model types === | ||
| 79 | + 2016-02-19T15:50:40: PointSource ...............: 1 | ||
| 80 | + 2016-02-19T15:50:40: CTAIrfBackground ..........: 1 | ||
| 81 | + 2016-02-19T15:50:40: | ||
| 82 | + 2016-02-19T15:50:40: +=======================+ | ||
| 83 | + 2016-02-19T15:50:40: | Parameter information | | ||
| 84 | + 2016-02-19T15:50:40: +=======================+ | ||
| 85 | + 2016-02-19T15:50:40: All parameters ............: 10 | ||
| 86 | + 2016-02-19T15:50:40: Fixed parameters ..........: 6 | ||
| 87 | + 2016-02-19T15:50:40: Free parameters (total) ...: 4 | ||
| 88 | + 2016-02-19T15:50:40: Free background parameters : 2 | ||
| 89 | + 2016-02-19T15:50:40: Free source parameters ....: 2 | ||
| 90 | + 2016-02-19T15:50:40: Free spectral parameters ..: 2 | ||
| 91 | + 2016-02-19T15:50:40: Parameters at limit .......: 0 | ||
| 92 | + | ||
| 93 | + | ||
| 94 | +Compute upper limit | ||
| 95 | +------------------- | ||
| 96 | +Very often in gamma-ray astronomy sources are at the verge of detection or even not detectable. | ||
| 97 | +In such cases, it is useful to derive an upper limit using :ref:`ctulimit`. | ||
| 98 | + | ||
| 99 | +.. code-block:: bash | ||
| 100 | + | ||
| 101 | + $ ctulimit | ||
| 102 | + Input event list, counts cube or observation definition XML file [selected_obs.xml] | ||
| 103 | + Input model XML file [crab_models.xml] | ||
| 104 | + Source of interest [Crab] | ||
| 105 | + | ||
| 106 | +The upper limit will be stored in the log file. To get the limit printed on screen, use the hidden parameter | ||
| 107 | +``debug=yes``. | ||
| 108 | + | ||
| 109 | +Compute asymmetric errors | ||
| 110 | +------------------------- | ||
| 111 | +When an analysis approaches its final state, it makes sense to have asymmetric errors on the parameters of the | ||
| 112 | +source of interest. For this purpose, the tool :ref:`cterror` can be used: | ||
| 113 | + | ||
| 114 | +.. code-block:: bash | ||
| 115 | + | ||
| 116 | + $ cterror debug=yes | ||
| 117 | + Input event list, counts cube or observation definition XML file [selected_obs.xml] | ||
| 118 | + Input model XML file [crab_models.xml] | ||
| 119 | + Source of interest [Crab] | ||
| 120 | + Output model XML file [cterror_results.xml] | ||
| 121 | + | ||
| 122 | +The output model does not contain asymmetric errors yet. The positive and negative uncertainties can be read from the | ||
| 123 | +logfile (or from screen) if ``debug=yes`` was specified. | ||
| 124 | + | ||
| 125 | +Compute spectral points | ||
| 126 | +----------------------- | ||
| 127 | +A very common task in astronomy is to compute spectral data points. To determine a spectral data point, a small energy | ||
| 128 | +range is considered and the model prefactor and its uncertainty is evaluated. The tool :ref:`csspec` works on both, binned | ||
| 129 | +and unbinned data. | ||
| 130 | + | ||
| 131 | +Unbinned | ||
| 132 | +^^^^^^^^ | ||
| 133 | + | ||
| 134 | +.. code-block:: bash | ||
| 135 | + | ||
| 136 | + $ csspec debug=yes | ||
| 137 | + Input event list, counts cube, or observation definition XML file [events.fits] selected_obs.xml | ||
| 138 | + Input model XML file [$CTOOLS/share/models/crab.xml] crab_models.xml | ||
| 139 | + Source name [Crab] | ||
| 140 | + Algorithm for defining energy bins (FILE|LIN|LOG) [LOG] | ||
| 141 | + Lower energy limit for spectral points (TeV) [0.1] 0.5 | ||
| 142 | + Upper energy limit for spectral points (TeV) [100.0] 50.0 | ||
| 143 | + Number of spectral points (1-10000) [20] 10 | ||
| 144 | + Output spectrum file [spectrum.fits] | ||
| 145 | + | ||
| 146 | +Binned | ||
| 147 | +^^^^^^ | ||
| 148 | + | ||
| 149 | +.. code-block:: bash | ||
| 150 | + | ||
| 151 | + $ csspec debug=yes | ||
| 152 | + Input event list, counts cube, or observation definition XML file [cntcube.fits] | ||
| 153 | + Input exposure cube file (only needed for stacked analysis) [expcube.fits] | ||
| 154 | + Input PSF cube file (only needed for stacked analysis) [psfcube.fits] | ||
| 155 | + Input background cube file (only needed for stacked analysis) [bkgcube.fits] | ||
| 156 | + Input model XML file [binned_models.xml] | ||
| 157 | + Source name [Crab] | ||
| 158 | + Number of spectral points (1-10000) [10] | ||
| 159 | + Lower energy limit for spectral points (TeV) [0.5] | ||
| 160 | + Upper energy limit for spectral points (TeV) [50.0] | ||
| 161 | + Output spectrum file [spectrum.fits] | ||
| 162 | + | ||
| 163 | + | ||
| 164 | +Plot spectral points | ||
| 165 | +^^^^^^^^^^^^^^^^^^^^ | ||
| 166 | +Instead of supporting a plotting library, simple example scripts are available to visualise data products. | ||
| 167 | +Have a look at ``$CTOOLS/examples/show_spectrum.py`` how a spectrum can be plotted using ``matplotlib``. | ||
| 168 | +To have first glance at the above computed spectrum one can use this script in the following way: | ||
| 169 | + | ||
| 170 | +.. code-block:: bash | ||
| 171 | + | ||
| 172 | + $ python $CTOOLS/examples/show_spectrum.py spectrum.fits | ||
| 173 | + | ||
| 174 | + | ||
| 175 | +Compute light curves | ||
| 176 | +-------------------- | ||
| 177 | +Euqally common in astronomy are light curves from a time-series analysis. The tool :ref:`cslightcrv` takes as input only | ||
| 178 | +an unbinned observation container (or a single event list). Counts cube are not possible as inout since the time informstion | ||
| 179 | +is lost during the binning procedure. Nevertheless, in case of large time bins, it is possible to require a binned analysis | ||
| 180 | +in each time bin. The tool will accordingly slice the event list and create the data products for the stacked analysis in each | ||
| 181 | +time span. | ||
| 182 | + | ||
| 183 | +.. code-block:: bash | ||
| 184 | + | ||
| 185 | + $ cslightcrv | ||
| 186 | + Input event list, counts cube, or observation definition XML file [selected_obs.xml] | ||
| 187 | + Input model XML file [crab_models.xml] | ||
| 188 | + Source name [Crab] | ||
| 189 | + Algorithm for defining time bins (FILE|LIN|GTI) [GTI] | ||
| 190 | + Number of energy bins per light curve bin (0=unbinned) [0] | ||
| 191 | + Lower energy limit of events (TeV) [0.1] 0.5 | ||
| 192 | + Upper energy limit of events (TeV) [100.0] 50 | ||
| 193 | + Output light curve file [lightcurve.fits] | ||
| 194 | + | ||
| 195 | +Compute a residual map | ||
| 196 | +---------------------- | ||
| 197 | +A frequent means to visually inspect the fitted model with respect to the input data is to create residual maps. | ||
| 198 | +Using :ref:`csresmap` such a map can easily be computed. The tool internally bins the data according to user parameters. | ||
| 199 | +(if data is not already provided in a binned state). Taking into account the corresponding IRFs, | ||
| 200 | +:reF:`csresmap` computes a model map (running :ref:`ctmodel`). Subsequently, there are several choices how data | ||
| 201 | +and model should be compared. There are three options: | ||
| 202 | + | ||
| 203 | +* ``SUB``: the subtraction of the model from the counts. The resulting map will display differences in absolute counts. | ||
| 204 | +* ``SUBDIV``: the subtraction and division by the model. The resulting map will display relative differences of the data with respect to the model. | ||
| 205 | +* ``SUBDIVSQRT``: the subtraction and division by the square root of the model. The resulting map will display an approximation of a residual significance map. In case of sufficient count statistic per bin, the bin value represents the significance. | ||
| 206 | + | ||
| 207 | +Example unbinned observation container | ||
| 208 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 209 | + | ||
| 210 | +.. code-block:: bash | ||
| 211 | + | ||
| 212 | + $ csresmap | ||
| 213 | + Input event list, counts cube, or observation definition XML file [selected_obs.xml] | ||
| 214 | + Input model XML file [crab_results.xml] | ||
| 215 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 216 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 217 | + Coordinate System (CEL|GAL) [CEL] | ||
| 218 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 219 | + Size of the X axis in pixels [100] | ||
| 220 | + Size of the Y axis in pixels [100] | ||
| 221 | + Pixel size (deg/pixel) [0.02] | ||
| 222 | + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV] | ||
| 223 | + Output residual map file [resmap.fits] | ||
| 224 | + | ||
| 225 | +Example binned/stacked observation container | ||
| 226 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 227 | +The tool behaves differently if the input observation parameter is a count cube. It will query for the additional binned reponse components. | ||
| 228 | + | ||
| 229 | +.. code-block:: bash | ||
| 230 | + | ||
| 231 | + $ csresmap | ||
| 232 | + Input event list, counts cube, or observation definition XML file [cntcube.fits] | ||
| 233 | + Input model cube file (generated with ctmodel) [NONE] | ||
| 234 | + Input exposure cube file (only needed for stacked analysis) [expcube.fits] | ||
| 235 | + Input PSF cube file (only needed for stacked analysis) [psfcube.fits] | ||
| 236 | + Input background cube file (only needed for stacked analysis) [bkgcube.fits] | ||
| 237 | + Input model XML file [binned_results.xml] | ||
| 238 | + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUBDIV] | ||
| 239 | + Output residual map file [resmap.fits] | ||
| 240 | + | ||
| 241 | +Note that the tool also queries for an input model cube file (which we set to ``NONE`` here). This is very convenient in case the model cube has already been precomputed | ||
| 242 | +using :ref:`ctmodel`. | ||
| 243 | + | ||
| 244 | +Example with input model cube | ||
| 245 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 246 | +In case a model cube can be provided, :ref:`csresmap` simply collapses both cube into one skymap and applies the specified algorithm. | ||
| 247 | + | ||
| 248 | +.. code-block:: bash | ||
| 249 | + | ||
| 250 | + $ csresmap | ||
| 251 | + Input event list, counts cube, or observation definition XML file [cntcube.fits] | ||
| 252 | + Input model cube file (generated with ctmodel) [modcube.fits] | ||
| 253 | + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB] | ||
| 254 | + Output residual map file [resmap.fits] | ||
| 255 | + | ||
| 256 | + | ||
| 257 | +Compute a test statistics (TS) map | ||
| 258 | +---------------------------------- | ||
| 259 | +The test statistics quantity is a more precise measurement of statistic than the ``SUBDIVSQRT`` option in :ref:`csresmap`. | ||
| 260 | +To get familiar with the concept of this procedure, :ref:`read more about TS maps <sec_tsmap>`. | ||
| 261 | + | ||
| 262 | +.. code-block:: bash | ||
| 263 | + | ||
| 264 | + $ cttsmap | ||
| 265 | + Input event list, counts cube or observation definition XML file [selected_obs.xml] | ||
| 266 | + Input model XML file [crab_results.xml] | ||
| 267 | + Test source name [Crab] | ||
| 268 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 269 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 270 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 271 | + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] | ||
| 272 | + Image scale (in degrees/pixel) [0.05] | ||
| 273 | + Size of the X axis in pixels [20] | ||
| 274 | + Size of the Y axis in pixels [20] | ||
| 275 | + Output Test Statistic map file [tsmap.fits] | ||
| 276 | + | ||
| 277 | +The output file will contain a sky map holding the TS map. In complementary extensions, for each free parameter of the source model, a map displaying the parameter | ||
| 278 | +value is added to the output file. For instance a map of the prefactor can be visualised with `DS9 <http://ds9.si.edu/>`_ as follows: | ||
| 279 | + | ||
| 280 | +.. code-block:: bash | ||
| 281 | + | ||
| 282 | + $ ds9 tsmap.fits[Prefactor] | ||
| 283 | + | ||
| 284 | +Since :ref:`cttsmap` has to run a parameter optimisation in each skymap bin, it is very time comsuming to | ||
| 285 | +compute a fine-granulated TS sky map. Read here how to split up the :ref:`computation into several jobs <sec_tips>`. | ||
| 286 | +Similar as many other tools :ref:`cttsmap` can also work on binned observation input: | ||
| 287 | + | ||
| 288 | +.. code-block:: bash | ||
| 289 | + | ||
| 290 | + $ cttsmap debug=yes | ||
| 291 | + Input event list, counts cube or observation definition XML file [cntcube.fits] | ||
| 292 | + Input exposure cube file (only needed for stacked analysis) [expcube.fits] | ||
| 293 | + Input PSF cube file (only needed for stacked analysis) [psfcube.fits] | ||
| 294 | + Input background cube file (only needed for stacked analysis) [bkgcube.fits] | ||
| 295 | + Input model XML file [binned_results.xml] | ||
| 296 | + Test source name [Crab] | ||
| 297 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 298 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 299 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 300 | + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] | ||
| 301 | + Image scale (in degrees/pixel) [0.05] | ||
| 302 | + Size of the X axis in pixels [20] | ||
| 303 | + Size of the Y axis in pixels [20] | ||
| 304 | + Output Test Statistic map file [tsmap.fits] | ||
| 305 | + | ||
| 306 | + | ||
| 0 | \ No newline at end of file | 307 | \ No newline at end of file |
doc/source/user_manual/iact_analysis/index.rst
0 → 100644
doc/source/user_manual/iact_analysis/simulation.rst
0 → 100644
| 1 | +.. _sec_simulation: | ||
| 2 | + | ||
| 3 | +Simulating IACT data | ||
| 4 | +==================== | ||
| 5 | +This section describes how to run simulations with current IACT data. In particular, a guide how to simulate data for a given set of IRFs, or a given set of | ||
| 6 | +observations is provided. In addition, tools to compute sensitivity curves or pull distributions are presented. | ||
| 7 | +All tools are based on the event list simulator :ref:`ctobssim`. | ||
| 8 | + | ||
| 9 | +Simulations for a specific observation container | ||
| 10 | +------------------------------------------------ | ||
| 11 | +On default, :ref:`ctobssim` builds a new observation container by querying relevant user parameters. In case | ||
| 12 | +an observation XML container is already in place, the hidden parameter ``inobs`` can be specified (:ref:`see here how to create an | ||
| 13 | +observation container <sec_analysis>`). It might be useful to run :ref:`ctselect` beforehand to ensure that the simulation knows about | ||
| 14 | +energy thresholds, field of view, etc. If an unselected observation container is used, the simulation tool will query for global energy | ||
| 15 | +boundaries and a global RoI radius. The following example illustrates the usage of an already selected observation container. | ||
| 16 | + | ||
| 17 | +.. code-block:: bash | ||
| 18 | + | ||
| 19 | + $ ctobssim inobs=selected_obs.xml | ||
| 20 | + Input model XML file [crab_models.xml] | ||
| 21 | + Output event data file or observation definition XML file [sim_obs.xml] | ||
| 22 | + | ||
| 23 | +Note that :ref:`ctobssim` writes a new event list file per observation into the current directory. If the observation container | ||
| 24 | +is large, it might be convenient to create a new folder and specify the ``prefix`` keyword (similar to :ref:`ctselect`): | ||
| 25 | + | ||
| 26 | +.. code-block:: bash | ||
| 27 | + | ||
| 28 | + $ mkdir simulated | ||
| 29 | + $ ctobssim inobs=selected_obs.xml prefix=simulated/sim_events_ | ||
| 30 | + ... | ||
| 31 | + | ||
| 32 | +The newly created observation container can be used in the same way as a normal observation container, i.e. it can be analysed with | ||
| 33 | +:ref:`ctlike` etc. | ||
| 34 | + | ||
| 35 | +.. note:: | ||
| 36 | + | ||
| 37 | + On default, :ref:`ctobssim` uses the same seed in every execution, i.e. the tool will always return the same simulated event list, regardless | ||
| 38 | + when, where and how it was executed. To steer this the user can provide the hidden integer parameter ``seed``: | ||
| 39 | + | ||
| 40 | + .. code-block:: bash | ||
| 41 | + | ||
| 42 | + $ ctobssim inobs=selected_obs.xml seed=2 | ||
| 43 | + | ||
| 44 | + This might be important if simulations are run in parallel. | ||
| 45 | + | ||
| 46 | +Simulating observation time for a given set of IRFs | ||
| 47 | +--------------------------------------------------- | ||
| 48 | +There are other use cases to run simulations, for instance when there is no observation container available yet. This might be the case e.g | ||
| 49 | +for an observation proposal: one would want to use a certain set of IRFs to simulate e.g. 50 hours of data to check whether a source | ||
| 50 | +of interest can be detected. For this purpose it is necessary to first find a set of suitable IRFs to be used for | ||
| 51 | +simulations. This can be easily done using :ref:`csfindobs` and using the ``expression`` parameter. The following example extracts | ||
| 52 | +observations with at a certain zenith angle range (20-25 degrees) and azimuth angle range (90-100 degrees): | ||
| 53 | + | ||
| 54 | +.. code-block:: bash | ||
| 55 | + | ||
| 56 | + $ csfindobs expression="AZ_PNT>90.0&&AZ_PNT<100.0&&ZEN_PNT>20.0&&ZEN_PNT<25.0" | ||
| 57 | + Name of FITS production (Run csiactdata to view your options) [fits-prod-name] | ||
| 58 | + Right ascension [83.6331] NONE | ||
| 59 | + Runlist outfile [zen20_az90.lis] | ||
| 60 | + | ||
| 61 | +Of course, further details could be specified, e.g. one might want to check for a certain range of muon efficieny, or any other parameter which impacts IRFs. | ||
| 62 | +The output runlist subsequently has to be converted to an observation XML file using :ref:`csiactobs` (:ref:`see here how to execute this tool <sec_analysis>`). | ||
| 63 | +In the following it is assumed that the observation XML file of the selected representative observations is produced and called ``myobs.xml``. | ||
| 64 | + | ||
| 65 | +Create caldb entry | ||
| 66 | +^^^^^^^^^^^^^^^^^^ | ||
| 67 | +Now, one of the observations in this container have to be added to the local calibration database. This database is used | ||
| 68 | +by ctools to extract IRFs for simulations. | ||
| 69 | + | ||
| 70 | +.. code-block:: bash | ||
| 71 | + | ||
| 72 | + $ csobs2caldb | ||
| 73 | + Input observation definition file [myobs.xml] | ||
| 74 | + Response output name (e.g. Zenith50) [zen20_az90] | ||
| 75 | + | ||
| 76 | +Note that the response name has been specified to "zen20_az90". Any name can be chosen which represents the observation which was inserted into the database. | ||
| 77 | +On default, :ref:`csobs2caldb` uses the first observation in the provided XML file. This can be modifyied by specifying the hidden integer parameter index, e.g.: | ||
| 78 | + | ||
| 79 | +.. code-block:: bash | ||
| 80 | + | ||
| 81 | + $ csobs2caldb index=1 | ||
| 82 | + | ||
| 83 | +This would use the second entry instead of the first from the observation container to extract the IRFs. | ||
| 84 | + | ||
| 85 | +.. note:: | ||
| 86 | + | ||
| 87 | + The ``index`` parameter starts counting from zero. | ||
| 88 | + | ||
| 89 | +Check caldb entries | ||
| 90 | +^^^^^^^^^^^^^^^^^^^ | ||
| 91 | +In order to check if the calibration database has successfully been added, the tool :ref:`cscaldb` can be used to inspect the calibration database: | ||
| 92 | + | ||
| 93 | +.. code-block:: bash | ||
| 94 | + | ||
| 95 | + $ cscaldb debug=yes | ||
| 96 | + ... | ||
| 97 | + 2016-02-25T15:57:11: +==============+ | ||
| 98 | + 2016-02-25T15:57:11: | Mission: cta | | ||
| 99 | + 2016-02-25T15:57:11: +==============+ | ||
| 100 | + 2016-02-25T15:57:11: === Response functions in database "hess" === | ||
| 101 | + 2016-02-25T15:57:11: zen20_az90 | ||
| 102 | + 2016-02-25T15:57:11: | ||
| 103 | + 2016-02-25T15:57:11: === Response functions in database "prod2" === | ||
| 104 | + 2016-02-25T15:57:11: North_0.5h | ||
| 105 | + 2016-02-25T15:57:11: North_50h | ||
| 106 | + 2016-02-25T15:57:11: North_5h | ||
| 107 | + 2016-02-25T15:57:11: South_0.5h | ||
| 108 | + 2016-02-25T15:57:11: South_50h | ||
| 109 | + 2016-02-25T15:57:11: South_5h | ||
| 110 | + ... | ||
| 111 | + | ||
| 112 | +Note that the observation of the actual instrument ("hess" in this example) was added to the CTA mission. This is mainly for structural and simplicity reasons | ||
| 113 | +and might be improved eventually. Nevertheless, this caldb entry can now be used to simulate any observation duration. | ||
| 114 | + | ||
| 115 | +Simulate | ||
| 116 | +^^^^^^^^ | ||
| 117 | + | ||
| 118 | +.. code-block:: bash | ||
| 119 | + | ||
| 120 | + $ ctobssim | ||
| 121 | + RA of pointing (degrees) (0-360) [83.63] | ||
| 122 | + Dec of pointing (degrees) (-90-90) [22.01] | ||
| 123 | + Radius of FOV (degrees) (0-180) [2.5] | ||
| 124 | + Start time (MET in s) [0.0] | ||
| 125 | + End time (MET in s) [18000.0] | ||
| 126 | + Lower energy limit (TeV) [0.5] | ||
| 127 | + Upper energy limit (TeV) [50] | ||
| 128 | + Calibration database [hess] | ||
| 129 | + Instrument response function [zen20_az90] | ||
| 130 | + Input model XML file [input_model.xml] | ||
| 131 | + Output event data file or observation definition XML file [simulated_events.fits] | ||
| 132 | + | ||
| 133 | +Here, a total duration of 5 hours at the crab position was simulated using the newly created caldb entry. Note that the choice of the energy range is up to the user. | ||
| 134 | +One should bear in mind that simulations only make sense in energy ranges where the IRFs are well-defined. The simulation | ||
| 135 | +assumes one observation of defined duration here. Accordingly, only one FITS event list is written as output. The input XML model | ||
| 136 | +must now contain only one background model and (of course) all sky models that are needed. Here is an example of an input model | ||
| 137 | +e.g. for the Crab Nebula: | ||
| 138 | + | ||
| 139 | +.. code-block:: xml | ||
| 140 | + | ||
| 141 | + <?xml version="1.0" standalone="no"?> | ||
| 142 | + <source_library title="source library"> | ||
| 143 | + <source name="Crab" type="PointSource"> | ||
| 144 | + <spectrum type="PowerLaw"> | ||
| 145 | + <parameter name="Prefactor" scale="1e-17" value="3.5" min="1e-07" max="1000.0" free="1"/> | ||
| 146 | + <parameter name="Index" scale="-1" value="2.5" min="0.0" max="+5.0" free="1"/> | ||
| 147 | + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/> | ||
| 148 | + </spectrum> | ||
| 149 | + <spatialModel type="SkyDirFunction"> | ||
| 150 | + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/> | ||
| 151 | + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/> | ||
| 152 | + </spatialModel> | ||
| 153 | + </source> | ||
| 154 | + <source name="CTABackgroundModel" type="CTAIrfBackground" instrument="CTA"> | ||
| 155 | + <spectrum type="PowerLaw"> | ||
| 156 | + <parameter name="Prefactor" scale="1.0" value="1.0" min="1e-3" max="1e+3" free="1"/> | ||
| 157 | + <parameter name="Index" scale="1.0" value="0.0" min="-5.0" max="+5.0" free="1"/> | ||
| 158 | + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/> | ||
| 159 | + </spectrum> | ||
| 160 | + </source> | ||
| 161 | + </source_library> | ||
| 162 | + | ||
| 163 | +Note that since the caldb entry (:ref:`csobs2caldb`) has been produced under the mission "CTA", the event list is now considered a CTA observation (regardless of the input instrument). | ||
| 164 | +Therefore, it is important to have "CTA" specified in the instrument attribute of the background model XML component. Of course the addition of any sky model can be | ||
| 165 | +done without a problem. | ||
| 166 | + | ||
| 167 | +The output of this kind of simulation is a single FITS event list file that can be used for further analysis. | ||
| 168 | +For instance one could run :ref:`ctlike` in the following way: | ||
| 169 | + | ||
| 170 | +.. code-block:: bash | ||
| 171 | + | ||
| 172 | + $ ctlike | ||
| 173 | + Input event list, counts cube or observation definition XML file [simulated_events.fits] | ||
| 174 | + Calibration database [hess] | ||
| 175 | + Instrument response function [zen20_az90] | ||
| 176 | + Input model XML file [input_model.xml] | ||
| 177 | + Output model XML file [simulated_fit_results.xml] | ||
| 178 | + | ||
| 179 | +Note that :ref:`ctlike` automatically queries for the IRF database in case no full observation container is provided. Any other tool taking | ||
| 180 | +unbinned observations as input will behave the same way. Now all means to run simulations with current IACTs are at hand. | ||
| 181 | +The following sections will describe tools for running higher level simulation analysis. | ||
| 182 | + | ||
| 183 | +Create pull distributions | ||
| 184 | +------------------------- | ||
| 185 | +Pull distributions inspect the capability of the science tools to find source parameters in the data that were injected into the simulations. | ||
| 186 | +The tool :ref:`cspull` will run several times the sequence of simulating and fitting the data. After each sequence, the pull value of the | ||
| 187 | +fitted parameters is written into an ASCII file. The pull value ``p`` of parameter ``X`` is defined as follows: | ||
| 188 | + | ||
| 189 | +.. math:: | ||
| 190 | + p = \frac{X_{\rm fitted} - X_{\rm true}}{\Delta X_{\rm fitted}} | ||
| 191 | + | ||
| 192 | +Where :math:`\Delta X_{\rm fitted}` is the derived error of the fitted parameter :math:`X_{\rm fitted}`. After several trials, the pull value should follow | ||
| 193 | +a normal distribution. Knowing this, the tool may help to distinguish if a source is detectable and can be normally reconstructed. | ||
| 194 | +The tool can either work on an input observation container or on the caldb entry that was created above. The number of trials is queried by the tool, too. | ||
| 195 | + | ||
| 196 | +Example: observation container | ||
| 197 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 198 | + | ||
| 199 | +.. code-block:: bash | ||
| 200 | + | ||
| 201 | + $ cspull inobs=selected_obs.xml | ||
| 202 | + Number of energy bins (0=unbinned) [0] | ||
| 203 | + Input model XML file [crab_results.xml] | ||
| 204 | + Output pull distribution file [pull.dat] | ||
| 205 | + Number of trials [100] | ||
| 206 | + | ||
| 207 | +Note that the tool asks for the energy binning. It will use an unbinned analysis if zero bins are provided. Otherwise | ||
| 208 | +it will query for binning parameters and will convert the simulated observations into a binned/stacked observation by running | ||
| 209 | +:ref:`ctbin`, :ref:`ctexpcube`, :ref:`ctpsfcube` and :ref:`ctbkgcube` in order before fitting the simulated data with :ref:`ctlike`. | ||
| 210 | + | ||
| 211 | +Example: caldb entry | ||
| 212 | +^^^^^^^^^^^^^^^^^^^^ | ||
| 213 | + | ||
| 214 | +.. code-block:: bash | ||
| 215 | + | ||
| 216 | + $ cspull | ||
| 217 | + RA of pointing (deg) (0-360) [83.6331] | ||
| 218 | + Dec of pointing (deg) (-90-90) [22.0145] | ||
| 219 | + Duration (in s) [18000.0] | ||
| 220 | + Lower energy limit (TeV) [0.5] | ||
| 221 | + Upper energy limit (TeV) [50] | ||
| 222 | + Calibration database [hess] | ||
| 223 | + Instrument response function [zen20_az90] | ||
| 224 | + Number of energy bins (0=unbinned) [0] | ||
| 225 | + Input model XML file [input_models.xml] | ||
| 226 | + Output pull distribution file [pull.dat] | ||
| 227 | + Number of trials [100] | ||
| 228 | + | ||
| 229 | +Again, the caldb entry ``zen20_az90`` which was created above is used to simulate 5 hours of data. The data is simulated and analysed 100 | ||
| 230 | +times in a row and the output is written into the file ``pull.dat``. | ||
| 231 | + | ||
| 232 | +Visualise pull distributions | ||
| 233 | +---------------------------- | ||
| 234 | +The output ASCII file of :ref:`cspull` can be inspected by the user using own scripts. As a starting point, the ``example`` folder | ||
| 235 | +contains a script that can be used to visualise the pull distribution histogram: | ||
| 236 | + | ||
| 237 | +.. code-block:: bash | ||
| 238 | + | ||
| 239 | + $ python $CTOOLS/examples/show_pull_histogram.py pull.dat Pull_Crab_Prefactor 50 | ||
| 240 | + | ||
| 241 | +Run the script without arguments to see the usage. In case a wrong parameter name was provided, the tool will print available parameters on the screen. | ||
| 242 | +The last integer parameter defines the binning of the histogram. In addition, the tool shows the standard normal distribution which should | ||
| 243 | +match the pull distribution for a high number of trials if everything went well. | ||
| 244 | + | ||
| 245 | +Create sensitivity curve | ||
| 246 | +------------------------ | ||
| 247 | +A very important means to illustrate if a source can be detected with a certain significance is a sensitivity curve. | ||
| 248 | +The tool :ref:`cssens` can be used to compute either an integral, or a differential sensitivity curve for a given observation container. | ||
| 249 | +Analogous to :ref:`cspull`, the tool works either on an unbinned observation container or using a caldb entry. | ||
| 250 | +The script :ref:`cssens` computes the differential or integrated CTA sensitivity using maximum likelihood fitting of a test source. | ||
| 251 | +The differential sensitivity is determined for a number of energy bins, the integral sensitivity is determined for a number | ||
| 252 | +of energy thresholds. The test source is fitted to simulated data using ctlike to determine it’s detection significance as a | ||
| 253 | +function of source flux. The source flux is then varied until the source significance achieves a given level, specified by the | ||
| 254 | +(hidden) significance parameter ``sigma``. To damp variations between individual Monte Carlo simulations, a sliding average is | ||
| 255 | +applied in the significance computation (controlled by hidden ``num_avg parameter``). Note that this procedure might take a while to compute. | ||
| 256 | +To switch from diffferential to integral sensitivity, the hidden parameter ``type="Integral"`` has to be specified. It should be noted that | ||
| 257 | +the results are written to a file ``sensitivity.dat``. To change the output file name, the hidden parameter ``outfile=mysensitivity.dat`` | ||
| 258 | +can be provided. | ||
| 259 | + | ||
| 260 | + | ||
| 261 | +Example: observation container | ||
| 262 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 263 | + | ||
| 264 | +.. code-block:: bash | ||
| 265 | + | ||
| 266 | + $ cssens inobs=selected_obs.xml | ||
| 267 | + Input model XML file [crab_results.xml] | ||
| 268 | + Source name [Crab] | ||
| 269 | + Lower energy limit (TeV) [0.5] | ||
| 270 | + Upper energy limit (TeV) [50] | ||
| 271 | + Number of energy bins for differential sensitivity computation [21] | ||
| 272 | + | ||
| 273 | +Example: caldb entry | ||
| 274 | +^^^^^^^^^^^^^^^^^^^^ | ||
| 275 | + | ||
| 276 | +.. code-block:: bash | ||
| 277 | + | ||
| 278 | + $ cssens | ||
| 279 | + Calibration database [hess] | ||
| 280 | + Instrument response function [zen20_az90] | ||
| 281 | + Effective exposure time (s) [18000.0] | ||
| 282 | + Radius of ROI (deg) [2.5] | ||
| 283 | + Input model XML file [input_models.xml] | ||
| 284 | + Source name [Crab] | ||
| 285 | + Lower energy limit (TeV) [0.5] | ||
| 286 | + Upper energy limit (TeV) [50] | ||
| 287 | + Number of energy bins for differential sensitivity computation [21] | ||
| 288 | + | ||
| 289 | +Note that the input model in the latter case only requires to contain one background model (see section above). Using this setup, | ||
| 290 | +the user can specify the observation time for the sensitivity curve. In this example, again, 5 hours are used. | ||
| 291 | + | ||
| 292 | +Visualise sensitivity curves | ||
| 293 | +---------------------------- | ||
| 294 | +In order to show a sensitivity curve, there is also a script in the example folder: | ||
| 295 | + | ||
| 296 | +.. code-block:: bash | ||
| 297 | + | ||
| 298 | + $ python $CTOOLS/examples/show_sensitivity.py sensitivity.dat | ||
| 299 | + | ||
| 300 | +This will display a sensitivity curve for the applied settings. |
doc/source/user_manual/iact_analysis/tips.rst
0 → 100644
| 1 | +.. _sec_tips: | ||
| 2 | + | ||
| 3 | +Tips, Tricks and FAQs | ||
| 4 | +===================== | ||
| 5 | +This section shows some of the wrinkles that could fascilitate the users' life when running IACT analyses. | ||
| 6 | + | ||
| 7 | +Visualise observations | ||
| 8 | +---------------------- | ||
| 9 | +A very common task is to inspect the observation container and visualise several quantities. | ||
| 10 | + | ||
| 11 | +Create a skymap of counts | ||
| 12 | +^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 13 | +The tool :ref:`ctskymap` allows to create a single map of counts contained inside an observation container. | ||
| 14 | +In contrast to :ref:`ctbin`, ctskymap takes into account all events that fall into the skymap energy range and spatial range. | ||
| 15 | +Since :ref:`ctbin` is intended for binned/stacked analysis, it respects energy thresholds, RoI selections etc. | ||
| 16 | +Accordingly, :ref:`ctskymap` should be used for a simple visualisation of the counts. | ||
| 17 | + | ||
| 18 | + | ||
| 19 | +.. code-block:: bash | ||
| 20 | + | ||
| 21 | + $ ctskymap debug=yes | ||
| 22 | + Input event list or observation definition XML file [events.fits] obs.xml | ||
| 23 | + Lower energy limit (TeV) [0.1] | ||
| 24 | + Upper energy limit (TeV) [100.0] | ||
| 25 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 26 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 27 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 28 | + Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL] | ||
| 29 | + Image scale (in degrees/pixel) [0.02] | ||
| 30 | + Size of the X axis in pixels [200] | ||
| 31 | + Size of the Y axis in pixels [200] | ||
| 32 | + Output skymap file [skymap.fits] | ||
| 33 | + | ||
| 34 | +Create region file of pointings | ||
| 35 | +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 36 | +In particular, one might want to know how the pointings of the observation container distribute across the sky. | ||
| 37 | +For this purpose, :ref:`csobsinfo` provides the hidden parameter ``ds9file`` which signals to produce a region file that | ||
| 38 | +contains these information: | ||
| 39 | + | ||
| 40 | +.. code-block:: bash | ||
| 41 | + | ||
| 42 | + $ csobsinfo ds9file=pointings.reg debug=yes | ||
| 43 | + Event list, counts cube, or observation definition file [obs.xml] | ||
| 44 | + | ||
| 45 | +Together with the previously produced skymap, one can plot the pointings using e.g. ds9: | ||
| 46 | + | ||
| 47 | +.. code-block:: bash | ||
| 48 | + | ||
| 49 | + $ ds9 skymap.fits -regions pointings.reg | ||
| 50 | + | ||
| 51 | +Plot zenith angle distribution of observations | ||
| 52 | +---------------------------------------------- | ||
| 53 | +The tool :ref:`csobsinfo` dumps more information into the log file. Some of them might be of interest to catch in a script. | ||
| 54 | +Therefore, one could simply use this tool via python and plot some distributions: | ||
| 55 | + | ||
| 56 | +.. code-block:: python | ||
| 57 | + | ||
| 58 | + # Run script | ||
| 59 | + import cscripts | ||
| 60 | + info = cscripts.csobsinfo() | ||
| 61 | + info["inobs"] = "selected_obs.xml" | ||
| 62 | + info.run() | ||
| 63 | + | ||
| 64 | + # Plot data | ||
| 65 | + import matplotlib.pyplot as plt | ||
| 66 | + plt.hist(info.zeniths(), bins=10, range=[0,90]) | ||
| 67 | + plt.show() | ||
| 68 | + | ||
| 69 | +This will plot a histogram showing the zenith angle distribution of the observation container. In addition there are more methods | ||
| 70 | +to access information in python: | ||
| 71 | + | ||
| 72 | +.. code-block:: python | ||
| 73 | + | ||
| 74 | + # csobsinfo can also return: | ||
| 75 | + info.azimuths() # list of azimuth values | ||
| 76 | + info.offsets() # list of offset values (only computed if hidden parameter offset=yes was specified) | ||
| 77 | + info.ebounds() # gammalib.GEbounds object of energy ranges | ||
| 78 | + info.gti() # gammalib.GGti object containing good time intervals | ||
| 79 | + | ||
| 80 | + | ||
| 81 | +Visualise models | ||
| 82 | +---------------- | ||
| 83 | +Similar to :ref:`csobsinfo`, there is :ref:`csmodelinfo` to find out more about the content of a model XML file. In order to show the | ||
| 84 | +position and sizes of the model on top of a skymap, this tool has the hidden parameter ``ds9file``, too. | ||
| 85 | + | ||
| 86 | +.. code-block:: bash | ||
| 87 | + | ||
| 88 | + csmodelinfo ds9file=models.reg debug=yes | ||
| 89 | + Input model XML file [$CTOOLS/share/models/crab.xml] | ||
| 90 | + | ||
| 91 | +There are several options regarding the color, text and other attributes of the region file. To see a full list, | ||
| 92 | +:ref:`visit the reference page <csmodelinfo>`. | ||
| 93 | +Analogous to the pointings of the observation container, the models can be visualised using e.g. ds9: | ||
| 94 | + | ||
| 95 | +.. code-block:: bash | ||
| 96 | + | ||
| 97 | + $ ds9 skymap.fits -regions models.reg | ||
| 98 | + | ||
| 99 | +Manipulating models in python | ||
| 100 | +----------------------------- | ||
| 101 | +This step will give some guidance how to work with model XML files in python and how to manipulate their content | ||
| 102 | + | ||
| 103 | +.. code-block:: python | ||
| 104 | + | ||
| 105 | + import gammalib | ||
| 106 | + | ||
| 107 | + # Open model file | ||
| 108 | + models = gammalib.GModels("$CTOOLS/share/models/crab.xml") | ||
| 109 | + | ||
| 110 | + # Access a model component | ||
| 111 | + src = models["Crab"] | ||
| 112 | + print(src) | ||
| 113 | + | ||
| 114 | + # Retrieve the spectral or spatial component | ||
| 115 | + print(src.spectral()) | ||
| 116 | + print(src.spatial()) | ||
| 117 | + | ||
| 118 | + # Print the spectral parameter "Prefactor" | ||
| 119 | + prefactor = src.spectral()["Prefactor"] | ||
| 120 | + print("Prefactor: "+str(prefactor.value())+" +- "+str(prefactor.error())) | ||
| 121 | + | ||
| 122 | + # Set the prefactor value | ||
| 123 | + prefactor.value(3.5e-17) | ||
| 124 | + | ||
| 125 | + # Loop over models and fix all parameters | ||
| 126 | + for model in models: | ||
| 127 | + | ||
| 128 | + # Loop over model parameters | ||
| 129 | + for par in model: | ||
| 130 | + | ||
| 131 | + # Fix parameter | ||
| 132 | + par.fix() | ||
| 133 | + | ||
| 134 | + # Release one specific parameter | ||
| 135 | + models["Crab"]["Prefactor"].free() | ||
| 136 | + | ||
| 137 | + # Set the parameter range | ||
| 138 | + models["Crab"]["Prefactor"].min(1e-18) | ||
| 139 | + models["Crab"]["Prefactor"].max(1e-16) | ||
| 140 | + | ||
| 141 | + # ... or in one step: | ||
| 142 | + models["Crab"]["Prefactor"].range(1e-18, 1e-16) | ||
| 143 | + | ||
| 144 | + # Remove model from container | ||
| 145 | + models.remove("Crab") | ||
| 146 | + | ||
| 147 | + # Save model to another XML file | ||
| 148 | + models.save("mymodels.xml") | ||
| 149 | + | ||
| 150 | + | ||
| 151 | +Retrieve likelihood values from :ref:`ctlike` | ||
| 152 | +--------------------------------------------- | ||
| 153 | +For some purposes, it might be useful to retrieve the fitted likelihood value and other results of the fit with :ref:`ctlike`. | ||
| 154 | + | ||
| 155 | +.. code-block:: python | ||
| 156 | + | ||
| 157 | + import gammalib | ||
| 158 | + import ctools | ||
| 159 | + | ||
| 160 | + # Create ctlike | ||
| 161 | + like = ctools.ctlike() | ||
| 162 | + like["inobs"] = "selected_obs.xml" | ||
| 163 | + like["inmodel"] = "mymodels.xml" | ||
| 164 | + like.run() | ||
| 165 | + | ||
| 166 | + # Get optimizer from ctlike | ||
| 167 | + opt = like.opt() | ||
| 168 | + print(opt) | ||
| 169 | + | ||
| 170 | + # Get likelihood value from optimizer | ||
| 171 | + print(opt.value()) | ||
| 172 | + | ||
| 173 | + # Get fit status: | ||
| 174 | + print(opt.status()) | ||
| 175 | + # 0: converged | ||
| 176 | + # 1: stalled | ||
| 177 | + # 2: singular curvature matrix encountered | ||
| 178 | + # 3: curvature matrix not positive definite | ||
| 179 | + # 4: errors are inaccurate | ||
| 180 | + | ||
| 181 | +Speed up analysis | ||
| 182 | +----------------- | ||
| 183 | +In some cases it may occur that the analysis takes very long. There are several reasons that can slow down the analysis quite severly. | ||
| 184 | +Not all issues can be tackled. Nevertheless, here is a list of actions that have proven to speed up the fit. | ||
| 185 | + | ||
| 186 | +* Consider switching to binned analysis if the observation time exceeds 50-100 hours. | ||
| 187 | +* Diffuse and extended sources in the RoI are computational expensive. However, in unbinned analysis they impact the speed much more. | ||
| 188 | +* Usage of energy dispersion in unbinned analysis could also cause a long fit. Try switching it off if speed is important. | ||
| 189 | +* Switch to a machine that supports OpenMP. This allows to parallelise the fit onto several cores. | ||
| 190 | + | ||
| 191 | +Compute excess maps | ||
| 192 | +------------------- | ||
| 193 | +A very important means to visualise images of sources are excess maps. The tool :ref:`csresmap` can also be used to create such a map. | ||
| 194 | +The user, however, has to manually remove the source(s) of interest from the input model XML file. | ||
| 195 | + | ||
| 196 | +.. code-block:: python | ||
| 197 | + | ||
| 198 | + import gammalib | ||
| 199 | + models = gammalib.GModels("crab_models.xml") | ||
| 200 | + models.remove("Crab") | ||
| 201 | + models.save("models_without_crab.xml") | ||
| 202 | + | ||
| 203 | +Of course, this can also be done by editing the XML file with an editor. | ||
| 204 | + | ||
| 205 | +Subsequently, :ref:`csresmap` can be executed using the ``algorithm=SUB`` parameter: | ||
| 206 | + | ||
| 207 | +.. code-block:: bash | ||
| 208 | + | ||
| 209 | + $ csresmap | ||
| 210 | + Input event list, counts cube, or observation definition XML file [selected_obs.xml] | ||
| 211 | + Input model XML file [models_without_crab.xml] | ||
| 212 | + First coordinate of image center in degrees (RA or galactic l) (0-360) [83.63] | ||
| 213 | + Second coordinate of image center in degrees (DEC or galactic b) (-90-90) [22.01] | ||
| 214 | + Coordinate System (CEL|GAL) [CEL] | ||
| 215 | + Projection method (AIT|AZP|CAR|MER|MOL|STG|TAN) [CAR] | ||
| 216 | + Size of the X axis in pixels [100] | ||
| 217 | + Size of the Y axis in pixels [100] | ||
| 218 | + Pixel size (deg/pixel) [0.02] | ||
| 219 | + Residual map computation algorithm (SUB|SUBDIV|SUBDIVSQRT) [SUB] | ||
| 220 | + Output residual map file [excessmap.fits] | ||
| 221 | + | ||
| 222 | + | ||
| 223 | +Compute flux maps | ||
| 224 | +----------------- | ||
| 225 | +The creation of a flux map is much more expensive in terms of computing than an excess map. In the section about | ||
| 226 | +:ref:`high level analysis tools <sec_high_level>`, the usage of :ref:`cttsmap` was demonstrated. This tool can also | ||
| 227 | +be used to generate flux maps. For this, the user has to incoorporate a test source in the XML model file. This source | ||
| 228 | +should consist of a spectral component using a PowerLaw2 model: | ||
| 229 | + | ||
| 230 | +.. code-block:: xml | ||
| 231 | + | ||
| 232 | + <source name="TestSource" type="PointSource"> | ||
| 233 | + <spectrum type="PowerLaw2"> | ||
| 234 | + <parameter scale="1e-07" name="Integral" min="1e-07" max="1000.0" value="1.0" free="1"/> | ||
| 235 | + <parameter scale="1.0" name="Index" min="-5.0" max="+5.0" value="-2.0" free="1"/> | ||
| 236 | + <parameter scale="1.0" name="LowerLimit" min="10.0" max="1000000.0" value="100.0" free="0"/> | ||
| 237 | + <parameter scale="1.0" name="UpperLimit" min="10.0" max="1000000.0" value="500000.0" free="0"/> | ||
| 238 | + </spectrum> | ||
| 239 | + <spatialModel type="SkyDirFunction"> | ||
| 240 | + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" /> | ||
| 241 | + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" /> | ||
| 242 | + </spatialModel> | ||
| 243 | + </source> | ||
| 244 | + | ||
| 245 | +The further content of the rest of the XML model file depends on the user requirements: | ||
| 246 | + | ||
| 247 | +* For an absolute flux map, all other sky models should be removed | ||
| 248 | +* For a residual flux map, all other sky models should be kept | ||
| 249 | + | ||
| 250 | +The tool :ref:`cttsmap` will create one skymap per free spectral parameter in the model. In the result FITS file, | ||
| 251 | +there will be an extension called "Integral" that contains the flux map. For the above XML example, the source name | ||
| 252 | +"TestSource" should be specified to :ref:`cttsmap`. | ||
| 253 | + | ||
| 254 | +Speed up TS map computation | ||
| 255 | +--------------------------- | ||
| 256 | +Since the computation of a TS map can be extremely time consuming, the option of splitting the computation into several jobs | ||
| 257 | +is supported. This might be of particular interest if the user has access to a batch farm. | ||
| 258 | + | ||
| 259 | +Split TS map computation | ||
| 260 | +^^^^^^^^^^^^^^^^^^^^^^^^ | ||
| 261 | +For the purpose of job splitting, the hidden parameter ``binmin`` and ``binmax`` were included in the tool. These are integer parameter | ||
| 262 | +that specify which bins should be computed. For instance, if the map should consist of 30x30(=900) pixels, the user could e.g. execute | ||
| 263 | + | ||
| 264 | +.. code-block:: bash | ||
| 265 | + | ||
| 266 | + $ cttsmap binmin=0 binmax=299 outmap=tsmap_0_299.fits | ||
| 267 | + ... | ||
| 268 | + $ cttsmap binmin=300 binmax=599 outmap=tsmap_300_599.fits | ||
| 269 | + ... | ||
| 270 | + $ cttsmap binmin=600 binmax=899 outmap=tsmap_600_899.fits | ||
| 271 | + ... | ||
| 272 | + | ||
| 273 | +Each command could run on a different machine. | ||
| 274 | + | ||
| 275 | +.. note:: | ||
| 276 | + | ||
| 277 | + The output file name of individual jobs should be different. Otherwise files could overwrite each other. The naming of the individual slices | ||
| 278 | + is up to the user. | ||
| 279 | + | ||
| 280 | +Merge sliced TS maps | ||
| 281 | +^^^^^^^^^^^^^^^^^^^^ | ||
| 282 | +The tool :ref:`cstsmapmerge` is intended to take care of merging the files that were produced while splitting the TS map computation. | ||
| 283 | +There are several options to pass the files to be merged as arguments: | ||
| 284 | + | ||
| 285 | +* a space-separated list of file names (e.g. tsmap1.fits tsmap2.fits) | ||
| 286 | +* a comma-separated list of file names (e.g. tsmap1.fits,tsmap2.fits) | ||
| 287 | +* a wildcard string (e.g. tsmap*.xml) | ||
| 288 | +* an ASCII file containing the file names, one file per line (e.g. @mymaps.txt) | ||
| 289 | + | ||
| 290 | +In this example the ASCII file method is presented: | ||
| 291 | + | ||
| 292 | +.. code-block:: bash | ||
| 293 | + | ||
| 294 | + # Put slice files into an ascii file | ||
| 295 | + $ ls tsmap_*.fits > tsmapfiles.txt | ||
| 296 | + | ||
| 297 | + # Run cstsmapmerge | ||
| 298 | + $ cstsmapmerge | ||
| 299 | + Input TS map FITS files [@tsmapfiles.txt] | ||
| 300 | + Output TS map FITS file [mytsmap.fits] | ||
| 301 | + | ||
| 302 | + | ||
| 303 | +Creating a python analysis pipeline | ||
| 304 | +----------------------------------- | ||
| 305 | +It is easily possible to build an own analysis workflow with a simple python script. The following source code example shows a python script | ||
| 306 | +running from gathering observations until fitting spectral points without storing intermediate data products on disk. It assumes that the | ||
| 307 | +environment variable ``$VHEFITS`` is set to the path where IACT FITS data is located. | ||
| 308 | + | ||
| 309 | +.. code-block:: python | ||
| 310 | + | ||
| 311 | + import gammalib | ||
| 312 | + import ctools | ||
| 313 | + import cscripts | ||
| 314 | + | ||
| 315 | + # Set debug flag | ||
| 316 | + debug = True | ||
| 317 | + | ||
| 318 | + # Set flag to use energy dispersion | ||
| 319 | + edisp = False | ||
| 320 | + | ||
| 321 | + # Set inmodel file name | ||
| 322 | + inmodel = "$GAMMALIB/test/data/model_point_plaw.xml" | ||
| 323 | + | ||
| 324 | + # Expand environment variable | ||
| 325 | + inmodel = gammalib.expand_env(inmodel) | ||
| 326 | + | ||
| 327 | + # Extract coordinates and model properties | ||
| 328 | + models = gammalib.GModels(inmodel) | ||
| 329 | + srcname = models[0].name() | ||
| 330 | + ra = models[0]["RA"].value() | ||
| 331 | + dec = models[0]["DEC"].value() | ||
| 332 | + | ||
| 333 | + # Find FITS production name | ||
| 334 | + iactdata = cscripts.csiactdata() | ||
| 335 | + iactdata["debug"] = debug | ||
| 336 | + iactdata.run() | ||
| 337 | + | ||
| 338 | + # Use first available production | ||
| 339 | + prodname = iactdata.names()[0] | ||
| 340 | + | ||
| 341 | + # Run csfindobs | ||
| 342 | + findobs = cscripts.csfindobs() | ||
| 343 | + findobs["ra"] = ra | ||
| 344 | + findobs["dec"] = dec | ||
| 345 | + findobs["rad"] = 2.5 | ||
| 346 | + findobs["prodname"] = prodname | ||
| 347 | + findobs["debug"] = debug | ||
| 348 | + findobs["outfile"] = "NONE" | ||
| 349 | + findobs.run() | ||
| 350 | + | ||
| 351 | + # Retrieve obervation IDs (runlist) | ||
| 352 | + obs_ids = findobs.obs_ids() | ||
| 353 | + | ||
| 354 | + # Build observation container | ||
| 355 | + iactobs = cscripts.csiactobs() | ||
| 356 | + iactobs["prodname"] = prodname | ||
| 357 | + iactobs["inmodel"] = inmodel | ||
| 358 | + iactobs["bkgpars"] = 1 | ||
| 359 | + iactobs["outobs"] = "NONE" | ||
| 360 | + iactobs["outmodel"] = "NONE" | ||
| 361 | + iactobs["debug"] = debug | ||
| 362 | + iactobs.runlist(obs_ids) | ||
| 363 | + iactobs.run() | ||
| 364 | + | ||
| 365 | + # Retrieve observation container and energy boundaries | ||
| 366 | + obs = iactobs.obs() | ||
| 367 | + ebounds = iactobs.ebounds() | ||
| 368 | + | ||
| 369 | + # Run ctselect | ||
| 370 | + select = ctools.ctselect(obs) | ||
| 371 | + select["usepnt"] = True | ||
| 372 | + select["rad"] = 2.5 | ||
| 373 | + select["usethres"] = "DEFAULT" | ||
| 374 | + select["tmin"] = 0.0 | ||
| 375 | + select["tmax"] = 0.0 | ||
| 376 | + select["emin"] = 0.1 | ||
| 377 | + select["emax"] = 100.0 | ||
| 378 | + select["debug"] = debug | ||
| 379 | + select.run() | ||
| 380 | + | ||
| 381 | + # Pass selected observations to ctlike | ||
| 382 | + like = ctools.ctlike(select.obs()) | ||
| 383 | + like["debug"] = debug | ||
| 384 | + like["edisp"] = edisp | ||
| 385 | + like.run() | ||
| 386 | + | ||
| 387 | + # Compute a spectrum and save | ||
| 388 | + spec = cscripts.csspec(like.obs()) | ||
| 389 | + spec["srcname"] = srcname | ||
| 390 | + spec["emin"] = ebounds.emin().TeV() | ||
| 391 | + spec["emax"] = ebounds.emax().TeV() | ||
| 392 | + spec["enumbins"] = 10 | ||
| 393 | + spec["edisp"] = edisp | ||
| 394 | + spec["ebinalg"] = "LOG" | ||
| 395 | + spec["debug"] = debug | ||
| 396 | + spec["outfile"] = "spectrum.fits" | ||
| 397 | + spec.execute() | ||
| 398 | + | ||
| 399 | + # Remove model for excess map computation | ||
| 400 | + like.obs().models().remove(srcname) | ||
| 401 | + | ||
| 402 | + # Compute an excess map | ||
| 403 | + resmap = cscripts.csresmap(like.obs()) | ||
| 404 | + resmap["xref"] = ra | ||
| 405 | + resmap["yref"] = dec | ||
| 406 | + resmap["proj"] = "CAR" | ||
| 407 | + resmap["coordsys"] = "CEL" | ||
| 408 | + resmap["emin"] = ebounds.emin().TeV() | ||
| 409 | + resmap["emax"] = ebounds.emax().TeV() | ||
| 410 | + resmap["nxpix"] = 100 | ||
| 411 | + resmap["nypix"] = 100 | ||
| 412 | + resmap["binsz"] = 0.02 | ||
| 413 | + resmap["outmap"] = "excessmap.fits" | ||
| 414 | + resmap.execute() | ||
| 415 | + | ||
| 416 | +The results of this analysis workflow can be inspected using an example python script and e.g. ds9: | ||
| 417 | + | ||
| 418 | +.. code-block:: bash | ||
| 419 | + | ||
| 420 | + $ python $CTOOLS/examples/show_spectrum.py spectrum.fits | ||
| 421 | + $ ds9 excessmap.fits | ||
| 422 | + | ||
| 0 | \ No newline at end of file | 423 | \ No newline at end of file |