From 4fb7bcac02c4f50e577b9fbbba2941423039d9f2 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 1 Feb 2022 11:21:44 +0100 Subject: [PATCH] Add "response" parameter to comobsbin --- doc/source/users/reference_manual/comobsbin.rst | 4 ++++ modules/comscripts/comobsbin.par | 1 + modules/comscripts/comobsbin.py | 111 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------- test/test_comobsbin.py | 5 +++-- 4 files changed, 82 insertions(+), 39 deletions(-) diff --git a/doc/source/users/reference_manual/comobsbin.rst b/doc/source/users/reference_manual/comobsbin.rst index 5694bce..a6f7f79 100755 --- a/doc/source/users/reference_manual/comobsbin.rst +++ b/doc/source/users/reference_manual/comobsbin.rst @@ -24,6 +24,10 @@ General parameters ``outfolder [string]`` Output folder for files. Put here the location of your COMPTEL data store. +``response [string]`` + Response type specifying whether ``MODEL`` IAQs or simulated IAQs provided + in the COMPTEL calibration database should be used. + ``ebinalg [string]`` Algorithm for defining energy bins. diff --git a/modules/comscripts/comobsbin.par b/modules/comscripts/comobsbin.par index 781288c..99e84f9 100644 --- a/modules/comscripts/comobsbin.par +++ b/modules/comscripts/comobsbin.par @@ -23,6 +23,7 @@ outfolder, s, a, $COMDATASTORE,,, "Output folder for files" # # Script parameters #================== +response, s, a, MODEL,MODEL|SIM2|SIM3,, "Response type" ebinalg, s, a, FILE,FILE|LIN|LOG,, "Algorithm for defining energy bins" emin, r, a, 0.75,0.05,50.0, "Minimum energy (MeV)" emax, r, a, 30.0,0.05,50.0, "Maximum energy (MeV)" diff --git a/modules/comscripts/comobsbin.py b/modules/comscripts/comobsbin.py index cc5378d..13eacba 100755 --- a/modules/comscripts/comobsbin.py +++ b/modules/comscripts/comobsbin.py @@ -60,6 +60,7 @@ class comobsbin(ctools.csobservation): self.obs().load(self['inobs'].filename()) # Query parameters + self['response'].string() self._get_ebounds() self._get_phases() if not self._phases.is_empty(): @@ -582,63 +583,99 @@ class comobsbin(ctools.csobservation): dpsi = self['dpsi'].real() dphibar = self['dphibar'].real() + # Get calibration database and response type + caldb = gammalib.GCaldb('cgro', 'comptel') + response = gammalib.toupper(self['response'].string()) + # Initialse list of IAQ names iaqnames = [] # Generate one IAQ for each energy boundary for i in range(ebounds.size()): - # Set IAQ filename - iaqname = '%s/iaq_%6.6d-%6.6dkeV.fits' % \ - (self['outfolder'].string(), - ebounds.emin(i).keV(), ebounds.emax(i).keV()) - iaqfile = gammalib.GFilename(iaqname) + # Handle simulated response + if response == 'SIM2' or response == 'SIM3': - # Set global IAQ filename - iaqname_global = '%s/iaq_%6.6d-%6.6dkeV.fits' % \ - (self._global_datastore, - ebounds.emin(i).keV(), ebounds.emax(i).keV()) - iaqfile_global = gammalib.GFilename(iaqname_global) + # Write header + self._log_header3(gammalib.NORMAL, + 'Use CALDB IAQ for %.3f - %.3f MeV' % + (ebounds.emin(i).MeV(), ebounds.emax(i).MeV())) - # Write header - self._log_header3(gammalib.NORMAL, 'Compute IAQ for %.3f - %.3f MeV' % \ - (ebounds.emin(i).MeV(), ebounds.emax(i).MeV())) + # Build response name + iaqname = '%s(%.2f-%.2f)MeV(%.0f)deg' % \ + (response, ebounds.emin(i).MeV(), ebounds.emax(i).MeV(), + dphibar) - # If IAQ file exists in global datastore then do nothing - if iaqfile_global.exists(): - iaqname = iaqname_global - self._log_value(gammalib.NORMAL, 'Global IAQ file exists', iaqfile_global.url()) + # Raise an exception if response does not exist + filename = caldb.filename('', '', 'IAQ', '', '', iaqname) + if filename.is_empty(): + msg = ('Could not find IAQ "'+iaqname+'" in calibration database. ' + 'Please select a different response type.') + raise RuntimeError(msg) - # If IAQ file exists then do nothing - elif iaqfile.exists(): - self._log_value(gammalib.NORMAL, 'Local IAQ file exists', iaqfile.url()) + # Log IAQ usage + self._log_value(gammalib.NORMAL, 'Use CALDB IAQ', iaqname) - # ... otherwise compute and save it + # Append IAQ name + iaqnames.append(iaqname) + + # ... otherwise handle modeled response else: - # Initialise IAQ - iaq = gammalib.GCOMIaq(55.0, dchi, 50.0, dphibar) + # Set IAQ filename + iaqname = '%s/iaq_%6.6d-%6.6dkeV.fits' % \ + (self['outfolder'].string(), + ebounds.emin(i).keV(), ebounds.emax(i).keV()) + iaqfile = gammalib.GFilename(iaqname) + + # Set global IAQ filename + iaqname_global = '%s/iaq_%6.6d-%6.6dkeV.fits' % \ + (self._global_datastore, + ebounds.emin(i).keV(), ebounds.emax(i).keV()) + iaqfile_global = gammalib.GFilename(iaqname_global) + + # Write header + self._log_header3(gammalib.NORMAL, + 'Compute IAQ for %.3f - %.3f MeV' % + (ebounds.emin(i).MeV(), ebounds.emax(i).MeV())) + + # If IAQ file exists in global datastore then do nothing + if iaqfile_global.exists(): + iaqname = iaqname_global + self._log_value(gammalib.NORMAL, + 'Global IAQ file exists', iaqfile_global.url()) + + # If IAQ file exists then do nothing + elif iaqfile.exists(): + self._log_value(gammalib.NORMAL, + 'Local IAQ file exists', iaqfile.url()) + + # ... otherwise compute and save it + else: - # Set IAQ energy range - ebds = gammalib.GEbounds(ebounds.emin(i), ebounds.emax(i)) + # Initialise IAQ + iaq = gammalib.GCOMIaq(55.0, dchi, 50.0, dphibar) - # Generate continuum IAQ - spectrum = gammalib.GModelSpectralPlaw(1.0,-2.0,gammalib.GEnergy(1.0,'MeV')) + # Set IAQ energy range + ebds = gammalib.GEbounds(ebounds.emin(i), ebounds.emax(i)) - # Compute IAQ - iaq.set(spectrum, ebds) + # Generate continuum IAQ + spectrum = gammalib.GModelSpectralPlaw(1.0,-2.0,gammalib.GEnergy(1.0,'MeV')) - # Save IAQ - iaq.save(iaqfile, True) + # Compute IAQ + iaq.set(spectrum, ebds) - # Log creation - self._log_value(gammalib.NORMAL, 'Local IAQ file created', iaqfile.url()) + # Save IAQ + iaq.save(iaqfile, True) + + # Log creation + self._log_value(gammalib.NORMAL, 'Local IAQ file created', iaqfile.url()) - # Log IAQ - self._log_string(gammalib.NORMAL, str(iaq)) + # Log IAQ + self._log_string(gammalib.NORMAL, str(iaq)) - # Append IAQ filename - iaqnames.append(iaqname) + # Append IAQ filename + iaqnames.append(iaqname) # Return return iaqnames diff --git a/test/test_comobsbin.py b/test/test_comobsbin.py index 1dd9d2f..4c8fae2 100755 --- a/test/test_comobsbin.py +++ b/test/test_comobsbin.py @@ -79,7 +79,7 @@ class Test(test): comobsbin = self._comscript('comobsbin') # Setup command - cmd = comobsbin+' inobs="'+self._obs+'" ebinalg="FILE" '+ \ + cmd = comobsbin+' inobs="'+self._obs+'" response="MODEL" ebinalg="FILE" '+ \ 'ebinfile="'+self._ebounds+'" phase="NONE" '+ \ 'coordsys="GAL" proj="TAN" '+ \ 'outfolder="comobsbin_cmd1" outobs="comobsbin_cmd1.xml" '+ \ @@ -93,7 +93,7 @@ class Test(test): self._check_result('comobsbin_cmd1.xml', 'comobsbin_cmd1') # Setup command - cmd = comobsbin+' inobs="file_that_does_not_exist.xml" ebinalg="FILE" '+ \ + cmd = comobsbin+' inobs="file_that_does_not_exist.xml" response="MODEL" ebinalg="FILE" '+ \ 'ebinfile="'+self._ebounds+'" phase="NONE" '+ \ 'coordsys="GAL" proj="TAN" '+ \ 'outfolder="comobsbin_cmd2" outobs="comobsbin_cmd2.xml" '+ \ @@ -117,6 +117,7 @@ class Test(test): # Set-up script bin = comscripts.comobsbin() bin['inobs'] = self._obs + bin['response'] = 'MODEL' bin['ebinalg'] = 'FILE' bin['ebinfile'] = self._ebounds bin['phase'] = 'NONE' -- libgit2 0.21.4