Maria Haupt / ctools

Commit 4a212ecef8f6651aba7ae5d800786a8d9f92ba03

Authored by Rolf Buehler
1 parent 30f87636

Added handling of stacked-binned ctlike - Issue 1503

Inline Side-by-side
Showing 1 changed file with 78 additions and 12 deletions
scripts/cspull.py
  View file @4a212ec
@@ -44,11 +44,15 @@ class cspull(ctools.cscript): @@ -44,11 +44,15 @@ class cspull(ctools.cscript):
44 self.version = "1.0.0" 44 self.version = "1.0.0"
45 45
46 # Initialise some members 46 # Initialise some members
47 - self.obs = None  
48 - self.model = None  
49 - self.m_inmodel = None  
50 - self.m_edisp = False  
51 - self.m_profile = False 47 + self.obs = None
  48 + self.model = None
  49 + self.m_inmodel = None
  50 + self.m_edisp = False
  51 + self.m_profile = False
  52 + self.m_exposure = None
  53 + self.m_psfcube = None
  54 + self.m_bckcube = None
  55 + self.m_stackmodels = None
52 56
53 # Make sure that parfile exists 57 # Make sure that parfile exists
54 file = self.parfile() 58 file = self.parfile()
@@ -194,8 +198,63 @@ class cspull(ctools.cscript): @@ -194,8 +198,63 @@ class cspull(ctools.cscript):
194 198
195 # Return 199 # Return
196 return 200 return
197 -  
198 - 201 +
  202 + def set_stackedirf(self):
  203 + """
  204 + For stacked analysis prepare stacked irfs
  205 + """
  206 +
  207 + #Get stacked exposure
  208 + expcube = ctools.ctexpcube(self.obs)
  209 + expcube["incube"].filename("NONE")
  210 + expcube["usepnt"].boolean(True)
  211 + expcube["ebinalg"].string("LOG")
  212 + expcube["binsz"].real(self.m_binsz)
  213 + expcube["nxpix"].integer(self.m_npix)
  214 + expcube["nypix"].integer(self.m_npix)
  215 + expcube["enumbins"].integer(self.m_enumbins)
  216 + expcube["emin"].real(self["emin"].real())
  217 + expcube["emax"].real(self["emax"].real())
  218 + expcube["coordsys"].string("GAL")
  219 + expcube["proj"].string("TAN")
  220 + expcube.run()
  221 +
  222 + #Get stacked psf
  223 + psfcube = ctools.ctpsfcube(self.obs)
  224 + psfcube["incube"].filename("NONE")
  225 + psfcube["usepnt"].boolean(True)
  226 + psfcube["ebinalg"].string("LOG")
  227 + psfcube["binsz"].real(self.m_binsz)
  228 + psfcube["nxpix"].integer(self.m_npix)
  229 + psfcube["nypix"].integer(self.m_npix)
  230 + psfcube["enumbins"].integer(self.m_enumbins)
  231 + psfcube["emin"].real(self["emin"].real())
  232 + psfcube["emax"].real(self["emax"].real())
  233 + psfcube["coordsys"].string("GAL")
  234 + psfcube["proj"].string("TAN")
  235 + psfcube.run()
  236 +
  237 + #Get stacked background
  238 + bkgcube = ctools.ctbkgcube(self.obs)
  239 + bkgcube["incube"].filename("NONE")
  240 + bkgcube["usepnt"].boolean(True)
  241 + bkgcube["ebinalg"].string("LOG")
  242 + bkgcube["binsz"].real(self.m_binsz)
  243 + bkgcube["nxpix"].integer(self.m_npix)
  244 + bkgcube["nypix"].integer(self.m_npix)
  245 + bkgcube["enumbins"].integer(self.m_enumbins)
  246 + bkgcube["emin"].real(self["emin"].real())
  247 + bkgcube["emax"].real(self["emax"].real())
  248 + bkgcube["coordsys"].string("GAL")
  249 + bkgcube["proj"].string("TAN")
  250 + bkgcube.run()
  251 +
  252 + #Store results
  253 + self.m_exposure = expcube.expcube().copy()
  254 + self.m_psfcube = psfcube.psfcube().copy()
  255 + self.m_bckcube = bkgcube.bkgcube().copy()
  256 + self.m_stackmodels = bkgcube.models().copy()
  257 +
199 def run(self): 258 def run(self):
200 """ 259 """
201 Run the script. 260 Run the script.
@@ -207,6 +266,10 @@ class cspull(ctools.cscript): @@ -207,6 +266,10 @@ class cspull(ctools.cscript):
207 # Get parameters 266 # Get parameters
208 self.get_parameters() 267 self.get_parameters()
209 268
  269 + #If several observations and binned: prepare stacked irfs
  270 + if self.obs.size()>0 and self.m_enumbins>0:
  271 + self.set_stackedirf()
  272 +
210 # Write input parameters into logger 273 # Write input parameters into logger
211 if self.logTerse(): 274 if self.logTerse():
212 self.log_parameters() 275 self.log_parameters()
@@ -276,6 +339,7 @@ class cspull(ctools.cscript): @@ -276,6 +339,7 @@ class cspull(ctools.cscript):
276 Parameters: 339 Parameters:
277 seed - Random number generator seed 340 seed - Random number generator seed
278 """ 341 """
  342 +
279 # Write header 343 # Write header
280 if self.logNormal(): 344 if self.logNormal():
281 self.log.header2("Trial "+str(seed+1)) 345 self.log.header2("Trial "+str(seed+1))
@@ -290,7 +354,12 @@ class cspull(ctools.cscript): @@ -290,7 +354,12 @@ class cspull(ctools.cscript):
290 log=self.m_log, \ 354 log=self.m_log, \
291 debug=self.m_debug, \ 355 debug=self.m_debug, \
292 chatter=self.m_chatter) 356 chatter=self.m_chatter)
293 - 357 +
  358 + #If stacked, add stacked responses and model
  359 + if self.obs.size()>0 and self.m_enumbins>0:
  360 + obs[0].response(self.m_exposure,self.m_psfcube,self.m_bckcube)
  361 + obs.models(self.m_stackmodels)
  362 +
294 # Determine number of events in simulation 363 # Determine number of events in simulation
295 nevents = 0.0 364 nevents = 0.0
296 for run in obs: 365 for run in obs:
@@ -305,16 +374,13 @@ class cspull(ctools.cscript): @@ -305,16 +374,13 @@ class cspull(ctools.cscript):
305 374
306 # Fit model 375 # Fit model
307 if self.m_profile: 376 if self.m_profile:
308 - models = obs.models()  
309 - new_models = gammalib.GModels() 377 + models = obs.models()
310 for i in range(models.size()): 378 for i in range(models.size()):
311 model_name = models[i].name() 379 model_name = models[i].name()
312 like = obsutils.cterror(obs, model_name, \ 380 like = obsutils.cterror(obs, model_name, \
313 log=self.m_log, \ 381 log=self.m_log, \
314 debug=self.m_debug, 382 debug=self.m_debug,
315 chatter=self.m_chatter) 383 chatter=self.m_chatter)
316 - new_models.append(like.obs().models()[i])  
317 - like.obs().models(new_models)  
318 else: 384 else:
319 like = obsutils.fit(obs, edisp=self.m_edisp, \ 385 like = obsutils.fit(obs, edisp=self.m_edisp, \
320 log=self.m_log, \ 386 log=self.m_log, \