Commit 3647c518d2885c9ab022d68340a659818ceab8ac

Authored by Jürgen Knödlseder
1 parent 48e0057e

Correct flux computation (was still based on the old interface and now leads to …

…excessive memory usage). This fixes bug #1425.
Showing 1 changed file with 15 additions and 13 deletions
scripts/cssens.py
@@ -154,7 +154,7 @@ class cssens(GApplication): @@ -154,7 +154,7 @@ class cssens(GApplication):
154 154
155 # Set some fixed parameters 155 # Set some fixed parameters
156 self.m_log = False # Logging in client tools 156 self.m_log = False # Logging in client tools
157 - self.m_debug = False # Debugging in client tools 157 + self.m_debug = False # self["debug"].boolean() # Debugging in client tools
158 158
159 # Return 159 # Return
160 return 160 return
@@ -424,26 +424,28 @@ class cssens(GApplication): @@ -424,26 +424,28 @@ class cssens(GApplication):
424 diffSens_value = [] 424 diffSens_value = []
425 iter = 0 425 iter = 0
426 test_flux = 0.1 # This is the initial test flux in Crab units 426 test_flux = 0.1 # This is the initial test flux in Crab units
427 - 427 +
428 # Loop until we break 428 # Loop until we break
429 while True: 429 while True:
430 - 430 +
431 # Update iteration counter 431 # Update iteration counter
432 iter += 1 432 iter += 1
433 - 433 +
434 # Write header 434 # Write header
435 if self.logExplicit(): 435 if self.logExplicit():
436 self.log.header2("Iteration "+str(iter)) 436 self.log.header2("Iteration "+str(iter))
437 437
438 # Set test flux 438 # Set test flux
439 - full_model["Test"]['Prefactor'].value(test_flux)  
440 - obs.models(full_model)  
441 - 439 + src_model = full_model.copy()
  440 + crab_flux = src_model["Test"]['Prefactor'].value()
  441 + src_model["Test"]['Prefactor'].value(crab_flux*test_flux)
  442 + obs.models(src_model)
  443 +
442 # Simulate events 444 # Simulate events
443 sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \ 445 sim = obsutils.sim(obs, nbins=self.m_nbins, seed=iter, \
444 binsz=self.m_binsz, npix=self.m_npix, \ 446 binsz=self.m_binsz, npix=self.m_npix, \
445 log=self.m_log, debug=self.m_debug) 447 log=self.m_log, debug=self.m_debug)
446 - 448 +
447 # Determine number of events in simulation 449 # Determine number of events in simulation
448 nevents = 0.0 450 nevents = 0.0
449 for run in sim: 451 for run in sim:
@@ -462,7 +464,7 @@ class cssens(GApplication): @@ -462,7 +464,7 @@ class cssens(GApplication):
462 result_bgm = like.obs().models().copy() 464 result_bgm = like.obs().models().copy()
463 LogL_bgm = like.opt().value() 465 LogL_bgm = like.opt().value()
464 npred_bgm = like.obs().npred() 466 npred_bgm = like.obs().npred()
465 - 467 +
466 # Assess quality based on a comparison between Npred and Nevents 468 # Assess quality based on a comparison between Npred and Nevents
467 quality_bgm = npred_bgm-nevents 469 quality_bgm = npred_bgm-nevents
468 470
@@ -484,7 +486,7 @@ class cssens(GApplication): @@ -484,7 +486,7 @@ class cssens(GApplication):
484 if self.logExplicit(): 486 if self.logExplicit():
485 self.log("Fit quality outside required range. Start over.\n") 487 self.log("Fit quality outside required range. Start over.\n")
486 continue 488 continue
487 - 489 +
488 # Write model fit results 490 # Write model fit results
489 if self.logExplicit(): 491 if self.logExplicit():
490 for model in result_bgm: 492 for model in result_bgm:
@@ -495,7 +497,7 @@ class cssens(GApplication): @@ -495,7 +497,7 @@ class cssens(GApplication):
495 self.log(str(par)+"\n") 497 self.log(str(par)+"\n")
496 498
497 # Fit background and test source 499 # Fit background and test source
498 - sim.models(full_model) 500 + sim.models(src_model)
499 like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug) 501 like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
500 result_all = like.obs().models().copy() 502 result_all = like.obs().models().copy()
501 LogL_all = like.opt().value() 503 LogL_all = like.opt().value()
@@ -541,8 +543,8 @@ class cssens(GApplication): @@ -541,8 +543,8 @@ class cssens(GApplication):
541 continue 543 continue
542 544
543 # Get fitted Crab, photon and energy fluxes 545 # Get fitted Crab, photon and energy fluxes
544 - cflux = result_all["Test"]['Prefactor'].value()  
545 - cflux_err = result_all["Test"]['Prefactor'].error() 546 + cflux = result_all["Test"]['Prefactor'].value() / crab_flux
  547 + cflux_err = result_all["Test"]['Prefactor'].error() / crab_flux
546 pflux = result_all["Test"].spectral().flux(emin, emax) 548 pflux = result_all["Test"].spectral().flux(emin, emax)
547 eflux = result_all["Test"].spectral().eflux(emin, emax) 549 eflux = result_all["Test"].spectral().eflux(emin, emax)
548 550