diff --git a/ChangeLog b/ChangeLog index 9772726..ff9f40f 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,9 +1,8 @@ -2021-11-29 +2021-12-01 * Version 2.0.0 released ======================== - Add phi_first and phi_last parameters to comlixfit and comlixmap (#3912) Add csadd2caldb script Add restart of bisection in case of bad convergence in ctulimit Add COMPTEL analysis scripts (#3738) diff --git a/NEWS b/NEWS index 14d15d8..06391c3 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 2.0.0 -29 November 2021 +1 December 2021 Introduction @@ -399,17 +399,11 @@ regions (#3248). comlixfit - COMPTEL maximum likelihood fitting using SRCLIX algorithm --------------------------------------------------------------------- -Added phi_first and phi_last parameters to allow selection of the Phibar layers -during the maximum likelihood fit (#3912). - Added script (#3738). comlixmap - Generation of COMPTEL Test Statistic map using SRCLIX algorithm --------------------------------------------------------------------------- -Added phi_first and phi_last parameters to allow selection of the Phibar layers -during the maximum likelihood fit (#3912). - Added script (#3738). diff --git a/README.md b/README.md index 9a0098d..d13bce0 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ctools information ================== -* Version: 2.0.0.dev (29 November 2021) +* Version: 2.0.0.dev (1 December 2021) * GammaLib dependency: 2.0.0.dev [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/) diff --git a/doc/source/users/reference_manual/comobsbin.rst b/doc/source/users/reference_manual/comobsbin.rst index 926b287..af7892a 100755 --- a/doc/source/users/reference_manual/comobsbin.rst +++ b/doc/source/users/reference_manual/comobsbin.rst @@ -66,6 +66,18 @@ General parameters ``(nphibar = 25) [integer]`` Number of bins in Phibar direction. +``(psdmin = 0) [integer]`` + Minimum PSD value. + +``(psdmax = 110) [integer]`` + Maximum PSD value. + +``(zetamin = 5.0) [real]`` + Minimum Earth horizon - Phibar (zeta) angle (deg). + +``(fpmtflag = 0) [integer]`` + Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT). + Standard parameters ------------------- diff --git a/doc/source/users/reference_manual/comobsmodel.rst b/doc/source/users/reference_manual/comobsmodel.rst index df3bfa9..f640cc4 100755 --- a/doc/source/users/reference_manual/comobsmodel.rst +++ b/doc/source/users/reference_manual/comobsmodel.rst @@ -41,9 +41,12 @@ General parameters ``iso [string]`` Isotropic component. -``modeltype [string]`` +``diffusetype [string]`` Diffuse model type. +``bkgtype [string]`` + Background model type. + ``(bremsmap = $COMDATA/../skymaps/galprop/map_comptel_bremsstrahlung.fits) [file]`` Bremsstrahlung map file name. diff --git a/modules/comscripts/comobsbin.par b/modules/comscripts/comobsbin.par index 3b159b7..527d5cc 100644 --- a/modules/comscripts/comobsbin.par +++ b/modules/comscripts/comobsbin.par @@ -39,7 +39,8 @@ npsi, i, h, 80,1,500, "Number of bins in Psi direction" nphibar, i, h, 25,1,100, "Number of bins in Phibar direction" psdmin, i, h, 0,0,110, "Minimum PSD value" psdmax, i, h, 110,0,110, "Maximum PSD value" -zetamin, r, h, 5.0,0.0,10.0, "Minimum Earth horizon - Phibar (zeta) angle" +zetamin, r, h, 5.0,0.0,10.0, "Minimum Earth horizon - Phibar (zeta) angle (deg)" +fpmtflag, i, h, 0,0,2, "Handling of D2 modules with failed PMT flag (0: exclude, 1: include, 2: exclude PMT)" # # Standard parameters diff --git a/modules/comscripts/comobsbin.py b/modules/comscripts/comobsbin.py index 7a654e1..beb5b9b 100755 --- a/modules/comscripts/comobsbin.py +++ b/modules/comscripts/comobsbin.py @@ -74,6 +74,7 @@ class comobsbin(ctools.csobservation): self['psdmin'].integer() self['psdmax'].integer() self['zetamin'].real() + self['fpmtflag'].integer() # Query ahead output model filename if self._read_ahead(): @@ -191,13 +192,14 @@ class comobsbin(ctools.csobservation): # Set phases self._select.phases(self._phases) - # Set suffix + # Set DRE, DRG and DRX suffix phases = self['phase'].string() - self._drx_suffix = '_phases'+phases.replace(';', '_') - self._drg_suffix = '_phases'+phases.replace(';', '_') self._dre_suffix = '_phases'+phases.replace(';', '_') + self._drg_suffix = '_phases'+phases.replace(';', '_') + self._drx_suffix = '_phases'+phases.replace(';', '_') # If PSD interval differs from standard interval then set the interval + # and append flag to suffix if self['psdmin'].integer() != 0 or self['psdmax'].integer() != 110: # Set PSD interval @@ -208,10 +210,22 @@ class comobsbin(ctools.csobservation): self._dre_suffix += '_psd%d-%d' % (self['psdmin'].integer(), self['psdmax'].integer()) - # If zeta angle differs from standard value then append zeta angle to suffix + # If zeta angle differs from standard value then append zeta angle to DRE + # and DRG suffix if self['zetamin'].real() != 5.0: - self._drg_suffix += '_zeta%.1f' % (self['zetamin'].real()) self._dre_suffix += '_zeta%.1f' % (self['zetamin'].real()) + self._drg_suffix += '_zeta%.1f' % (self['zetamin'].real()) + + # If handling of D2 modules with failed PMT flag differs from standard then + # set the flag and append flag to suffix + if self['fpmtflag'].integer() != 0: + + # Set handling of D2 modules with failed PMT flag + self._select.fpmtflag(self['fpmtflag'].integer()) + + # Set DRE and DRG suffix + self._dre_suffix += '_fpmt%1d' % (self['fpmtflag'].integer()) + self._drg_suffix += '_fpmt%1d' % (self['fpmtflag'].integer()) # Log input observations self._log_string(gammalib.NORMAL, str(self._select)) diff --git a/modules/comscripts/comobsmodel.par b/modules/comscripts/comobsmodel.par index 97fe421..b5a25f6 100644 --- a/modules/comscripts/comobsmodel.par +++ b/modules/comscripts/comobsmodel.par @@ -22,18 +22,19 @@ outmodel, f, a, models.xml,,, "Output model definition file" # # Script parameters #================== -ra, r, a, NONE,,, "Right Ascension of point source (deg)" -dec, r, a, NONE,,, "Declinaison of point source (deg)" -srcname, s, a, NONE,,, "Name of point source" -brems, s, a, NONE,NONE|MAP|CUBE,, "Bremsstrahlung component" -ic, s, a, NONE,NONE|MAP|CUBE,, "Inverse Compton component" -iso, s, a, NONE,NONE|CONST|CONSTFIX,, "Isotropic component" -modeltype, s, h, BINS,NODES|BINS,, "Diffuse model type" -ebinfile, f, h, NONE,,, "Name of the file containing the diffuse model energy bin definition" -bremsmap, f, h, $COMDATA/../skymaps/galprop/map_comptel_bremsstrahlung.fits,,, "Bremsstrahlungs map" -bremscube, f, h, $COMDATA/../skymaps/galprop/bremss_mapcube_54_77Xvarh7S.fits,,, "Bremsstrahlungs cube" -icmap, f, h, $COMDATA/../skymaps/galprop/map_comptel_ic.fits,,, "Inverse Compton map" -iccube, f, h, $COMDATA/../skymaps/galprop/ics_isotropic_mapcube_54_77Xvarh7S.fits,,, "Inverse Compton cube" +ra, r, a, NONE,,, "Right Ascension of point source (deg)" +dec, r, a, NONE,,, "Declinaison of point source (deg)" +srcname, s, a, NONE,,, "Name of point source" +brems, s, a, NONE,NONE|MAP|CUBE,, "Bremsstrahlung component" +ic, s, a, NONE,NONE|MAP|CUBE,, "Inverse Compton component" +iso, s, a, NONE,NONE|CONST|CONSTFIX,, "Isotropic component" +diffusetype, s, h, BINS,NODES|BINS,, "Diffuse model type" +bkgtype, s, h, BINS,NODES|BINS,, "Background model type" +ebinfile, f, h, NONE,,, "Name of the file containing the diffuse model energy bin definition" +bremsmap, f, h, $COMDATA/../skymaps/galprop/map_comptel_bremsstrahlung.fits,,, "Bremsstrahlungs map" +bremscube, f, h, $COMDATA/../skymaps/galprop/bremss_mapcube_54_77Xvarh7S.fits,,, "Bremsstrahlungs cube" +icmap, f, h, $COMDATA/../skymaps/galprop/map_comptel_ic.fits,,, "Inverse Compton map" +iccube, f, h, $COMDATA/../skymaps/galprop/ics_isotropic_mapcube_54_77Xvarh7S.fits,,, "Inverse Compton cube" # # Standard parameters diff --git a/modules/comscripts/comobsmodel.py b/modules/comscripts/comobsmodel.py index 973825b..86a56a6 100755 --- a/modules/comscripts/comobsmodel.py +++ b/modules/comscripts/comobsmodel.py @@ -64,7 +64,8 @@ class comobsmodel(ctools.csobservation): self['brems'].string() self['ic'].string() self['iso'].string() - self['modeltype'].string() + self['diffusetype'].string() + self['bkgtype'].string() if self['ebinfile'].is_valid(): self['ebinfile'].filename() self['bremsmap'].filename() @@ -111,9 +112,9 @@ class comobsmodel(ctools.csobservation): # Return energy boundaries return ebounds - def _generate_background_model(self, obs): + def _generate_background_nodes_model(self, obs): """ - Generate background model component for observation + Generate nodes background model component for observation Parameters ---------- @@ -164,6 +165,60 @@ class comobsmodel(ctools.csobservation): # Return model return model + def _generate_background_bins_model(self, obs): + """ + Generate bins background model component for observation + + Parameters + ---------- + obs : `~gammalib.GCOMObservation` + COMPTEL observation + + Returns + ------- + model : `~gammalib.GCOMModelDRBPhibarBins` + Background model + """ + # Set model attributes + name = 'Background_%s' % obs.id() + id = obs.id() + + # Extract attributes from observation + drb = obs.drb() + npix = drb.nchi() * drb.npsi() + nphibar = drb.nphibar() + dphibar = 2.0 # Fixed so far as we have no method to recover + + # Create XML model list object + xml = gammalib.GXml() + lib = xml.append('source_library title="source library"') + src = lib.append('source name="%s" type="DRBPhibarBins" instrument="COM" id="%s"' % (name, id)) + + # Loop over all phibar layers + for k in range(nphibar): + + # Compute DRB content for phibar layer + sum_drb = 0.0 + for i in range(npix): + index = i + k * npix + sum_drb += drb[index] + + # Free or fix parameter + if sum_drb > 0.0: + free = '1' + else: + free = '0' + + # Add bin + src.append('parameter name="Normalization" value="1.00" ' + 'scale="1" min="0" max="1000" free="%s"' % free) + + # Create background model from XML source element + model = gammalib.GCOMModelDRBPhibarBins(src) + + # Return model + return model + def _generate_source_model(self): """ Generate source model @@ -231,7 +286,7 @@ class comobsmodel(ctools.csobservation): ebds = ebounds # Re-sample model for 'NODES' model type - if self['modeltype'].string() == 'NODES': + if self['diffusetype'].string() == 'NODES': # Setup energies vector for bin mean energies energies = gammalib.GEnergies() @@ -485,7 +540,10 @@ class comobsmodel(ctools.csobservation): self._log_string(gammalib.NORMAL, str(obs)) # Generate model - model = self._generate_background_model(obs) + if self['bkgtype'].string() == 'NODES': + model = self._generate_background_nodes_model(obs) + else: + model = self._generate_background_bins_model(obs) # Log model self._log_string(gammalib.EXPLICIT, str(model))