Add support for energy dispersion.

Jürgen Knödlseder
1 parent 5f4af19a
Showing 2 changed files with 19 additions and 6 deletions
scripts/cspull.py
scripts/obsutils.py
@@ -94,11 +94,12 @@ class cspull(GApplication):
 			# Create default parfile
 			pars = GApplicationPars()
-			pars.append(GApplicationPar("srcmdl","f","a","$GAMMALIB/share/models/crab.xml","","","Source model"))
+			pars.append(GApplicationPar("srcmdl","f","a","$CTOOLS/share/models/crab.xml","","","Source model"))
 			pars.append(GApplicationPar("outfile","f","a","pull.dat","","","Output file name"))
 			pars.append(GApplicationPar("ntrials","i","a","10","","","Number of trials"))
 			pars.append(GApplicationPar("caldb","s","h","","","","Calibration database"))
 			pars.append(GApplicationPar("irf","s","a","cta_dummy_irf","","","Instrument response function"))
+			pars.append(GApplicationPar("edisp","b","h","no","","","Apply energy dispersion?"))
 			pars.append(GApplicationPar("ra","r","a","83.6331","0","360","RA of pointing (deg)"))
 			pars.append(GApplicationPar("dec","r","a","22.0145","-90","90","Dec of pointing (deg)"))
 			pars.append(GApplicationPar("emin","r","a","0.1","0.0","","Lower energy limit (TeV)"))
@@ -126,6 +127,7 @@ class cspull(GApplication):
 		self.m_ntrials  = self["ntrials"].integer()
 		self.m_caldb    = self["caldb"].string()
 		self.m_irf      = self["irf"].string()
+		self.m_edisp    = self["edisp"].boolean()
 		self.m_ra       = self["ra"].real()
 		self.m_dec      = self["dec"].real()
 		self.m_emin     = self["emin"].real()
@@ -268,6 +270,7 @@ class cspull(GApplication):
 		                   seed=seed, \
 		                   binsz=self.m_binsz, \
 		                   npix=self.m_npix, \
+                           edisp=self.m_edisp, \
 		                   log=self.m_log, debug=self.m_debug)
 		# Determine number of events in simulation
@@ -283,7 +286,8 @@ class cspull(GApplication):
 			self.log("\n")
 		# Fit model
-		like = obsutils.fit(obs, log=self.m_log, debug=self.m_debug)
+		like = obsutils.fit(obs, edisp=self.m_edisp, \
+                            log=self.m_log, debug=self.m_debug)
 		# Store results
 		logL   = like.opt().value()
@@ -25,7 +25,8 @@ import gammalib
 # ===================== #
 # Simulate observations #
 # ===================== #
-def sim(obs, log=False, debug=False, seed=0, nbins=0, binsz=0.05, npix=200):
+def sim(obs, log=False, debug=False, edisp=False, seed=0, nbins=0,
+        binsz=0.05, npix=200):
     """
     Simulate events for all observations in the container.
@@ -33,7 +34,8 @@ def sim(obs, log=False, debug=False, seed=0, nbins=0, binsz=0.05, npix=200):
      obs   - Observation container
     Keywords:
      log   - Create log file(s)
-     debug - Create screen dump
+     debug - Create console dump?
+     edisp - Apply energy dispersion?
      seed  - Seed value for simulations (default: 0)
      nbins - Number of energy bins (default: 0=unbinned)
      binsz - Pixel size for binned simulation (deg/pixel)
@@ -50,6 +52,9 @@ def sim(obs, log=False, debug=False, seed=0, nbins=0, binsz=0.05, npix=200):
     # Optionally switch-on debugging model
     if debug:
         sim["debug"].boolean(True)
+
+    # Optionally apply energy dispersion
+    sim["edisp"].boolean(edisp)
     # Run ctobssim application. This will loop over all observations in the
     # container and simulation the events for each observation. Note that
@@ -122,7 +127,7 @@ def sim(obs, log=False, debug=False, seed=0, nbins=0, binsz=0.05, npix=200):
 # ================ #
 # Fit observations #
 # ================ #
-def fit(obs, log=False, debug=False):
+def fit(obs, log=False, debug=False, edisp=False):
     """
     Perform maximum likelihood fitting of observations in the container.
@@ -131,6 +136,7 @@ def fit(obs, log=False, debug=False):
     Keywords:
      log   - Create log file(s)
      debug - Create screen dump
+     edisp - Apply energy dispersion?
     """
     # Allocate ctlike application
     like = ctools.ctlike(obs)
@@ -143,6 +149,9 @@ def fit(obs, log=False, debug=False):
     if debug:
         like["debug"].boolean(True)
+    # Optionally apply energy dispersion
+    like["edisp"].boolean(edisp)
+
     # Run ctlike application.
     like.run()
@@ -266,7 +275,7 @@ def modmap(obs, eref=0.1, proj=&quot;TAN&quot;, coord=&quot;GAL&quot;, xval=0.0, yval=0.0, \
 # ======================= #
 def set(pntdir, tstart=0.0, duration=1800.0, deadc=0.95, \
         emin=0.1, emax=100.0, rad=5.0, \
-        irf="cta_dummy_irf", caldb=""):
+        irf="cta_dummy_irf", caldb="$CTOOLS/share/caldb/data/cta/dummy"):
     """
     Returns a single CTA observation. By looping over this function we can
     add CTA observations to the observation container.