Andreas Specovius / ctools

Commit 8581bdb201326b5ce578202e8844edd80c3548c0

Authored by knodlseder
1 parent 0f5ab505

implement differential sensitivity estimate and improve convergence testing

Inline Side-by-side
Showing 1 changed file with 125 additions and 49 deletions
scripts/cssens.py
  View file @8581bdb
... ... @@ -57,7 +57,7 @@ class cssens(GApplication):
57 57 """
58 58 # Set name
59 59 self.name = "cssens"
60   - self.version = "0.2.0"
  60 + self.version = "0.3.0"
61 61  
62 62 # Make sure that parfile exists
63 63 file = self.parfile()
... ... @@ -117,7 +117,7 @@ class cssens(GApplication):
117 117 pars.append(GPar("npix","i","h","200","","","Number of pixels for binned"))
118 118 pars.append(GPar("binsz","r","h","0.05","","","Pixel size for binned (deg/pixel)"))
119 119 pars.append(GPar("sigma","r","h","5.0","","","Significance threshold"))
120   - pars.append(GPar("ts_use","b","h","no","","","Use TS instead of Sigma?"))
  120 + pars.append(GPar("ts_use","b","h","yes","","","Use TS?"))
121 121 pars.append(GPar("index","r","h","-2.48","","","Assumed spectral index"))
122 122 pars.append(GPar("max_iter","i","h","50","","","Maximum number of iterations"))
123 123 pars.append(GPar("num_avg","i","h","3","","","Number of iterations for sliding average"))
... ... @@ -183,7 +183,7 @@ class cssens(GApplication):
183 183 self.log("\n")
184 184  
185 185 # Initialise script
186   - colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux']
  186 + colnames = ['loge', 'emin', 'emax', 'crab', 'pflux', 'eflux', 'diffSens']
187 187 results = []
188 188  
189 189 # Initialise models. Note that we centre the point source at the Galactic
... ... @@ -389,12 +389,17 @@ class cssens(GApplication):
389 389 bkg_model - Background model
390 390 full_model - Source model
391 391 """
  392 + # Set TeV->erg conversion factor
  393 + tev2erg = 1.6021764
  394 +
392 395 # Determine energy boundaries from first observation in
393 396 # the container
394 397 for run in obs:
395   - emin = run.events().ebounds().emin()
396   - emax = run.events().ebounds().emax()
397   - loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
  398 + emin = run.events().ebounds().emin()
  399 + emax = run.events().ebounds().emax()
  400 + loge = math.log10(math.sqrt(emin.TeV()*emax.TeV()))
  401 + erg_mean = math.pow(10.0, loge) * tev2erg
  402 + erg_width = (emax.TeV()-emin.TeV()) * tev2erg
398 403 break
399 404  
400 405 # Write header
... ... @@ -403,11 +408,12 @@ class cssens(GApplication):
403 408 self.log.header2("Energies: "+str(emin)+" - "+str(emax))
404 409  
405 410 # Initialise loop
406   - flux_value = []
407   - pflux_value = []
408   - eflux_value = []
409   - iter = 0
410   - test_flux = 0.01 # This is the initial test flux in Crab units
  411 + flux_value = []
  412 + pflux_value = []
  413 + eflux_value = []
  414 + diffSens_value = []
  415 + iter = 0
  416 + test_flux = 0.01 # This is the initial test flux in Crab units
411 417  
412 418 # Loop until we break
413 419 while True:
... ... @@ -417,10 +423,10 @@ class cssens(GApplication):
417 423  
418 424 # Write header
419 425 if self.logExplicit():
420   - self.log.header2("Iteation "+str(iter))
  426 + self.log.header2("Iteration "+str(iter))
421 427  
422 428 # Set test flux
423   - full_model[1]['Prefactor'].value(test_flux)
  429 + full_model["Test"]['Prefactor'].value(test_flux)
424 430 obs.models(full_model)
425 431  
426 432 # Simulate events
... ... @@ -446,6 +452,9 @@ class cssens(GApplication):
446 452 result_bgm = like.obs().models().copy()
447 453 LogL_bgm = like.opt().value()
448 454 npred_bgm = like.obs().npred()
  455 +
  456 + # Assess quality based on a comparison between Npred and Nevents
  457 + quality_bgm = npred_bgm-nevents
449 458  
450 459 # Write background fit results
451 460 if self.logExplicit():
... ... @@ -456,13 +465,25 @@ class cssens(GApplication):
456 465 self.log.parformat("Number of predicted events")
457 466 self.log(npred_bgm)
458 467 self.log("\n")
  468 + self.log.parformat("Fit quality")
  469 + self.log(quality_bgm)
  470 + self.log("\n")
  471 +
  472 + # Start over if the fit quality was bad
  473 + if abs(quality_bgm) > 3.0:
  474 + if self.logExplicit():
  475 + self.log("Fit quality outside required range. Start over.\n")
  476 + continue
  477 +
  478 + # Write model fit results
  479 + if self.logExplicit():
459 480 for model in result_bgm:
460 481 self.log.parformat("Model")
461 482 self.log(model.name())
462 483 self.log("\n")
463 484 for par in model:
464 485 self.log(str(par)+"\n")
465   -
  486 +
466 487 # Fit background and test source
467 488 sim.models(full_model)
468 489 like = obsutils.fit(sim, log=self.m_log, debug=self.m_debug)
... ... @@ -471,11 +492,8 @@ class cssens(GApplication):
471 492 npred_all = like.obs().npred()
472 493 ts = 2.0*(LogL_bgm-LogL_all)
473 494  
474   - # Get fitted Crab, photon and energy fluxes
475   - cflux = result_all["Test"]['Prefactor'].value()
476   - cflux_err = result_all["Test"]['Prefactor'].error()
477   - pflux = cast_GModelSky(result_all["Test"]).spectral().flux(emin, emax)
478   - eflux = cast_GModelSky(result_all["Test"]).spectral().eflux(emin, emax)
  495 + # Assess quality based on a comparison between Npred and Nevents
  496 + quality_all = npred_all-nevents
479 497  
480 498 # Write background and test source fit results
481 499 if self.logExplicit():
... ... @@ -489,6 +507,59 @@ class cssens(GApplication):
489 507 self.log.parformat("Number of predicted events")
490 508 self.log(npred_all)
491 509 self.log("\n")
  510 + self.log.parformat("Fit quality")
  511 + self.log(quality_all)
  512 + self.log("\n")
  513 + #
  514 + for model in result_all:
  515 + self.log.parformat("Model")
  516 + self.log(model.name())
  517 + self.log("\n")
  518 + for par in model:
  519 + self.log(str(par)+"\n")
  520 +
  521 + # Start over if the fit quality was bad
  522 + if abs(quality_all) > 3.0:
  523 + if self.logExplicit():
  524 + self.log("Fit quality outside required range. Start over.\n")
  525 + continue
  526 +
  527 + # Start over if TS was non-positive
  528 + if ts <= 0.0:
  529 + if self.logExplicit():
  530 + self.log("Non positive TS. Start over.\n")
  531 + continue
  532 +
  533 + # Get fitted Crab, photon and energy fluxes
  534 + cflux = result_all["Test"]['Prefactor'].value()
  535 + cflux_err = result_all["Test"]['Prefactor'].error()
  536 + pflux = cast_GModelSky(result_all["Test"]).spectral().flux(emin, emax)
  537 + eflux = cast_GModelSky(result_all["Test"]).spectral().eflux(emin, emax)
  538 +
  539 + # Compute differential sensitivity
  540 + diffSens = pflux / erg_width * erg_mean*erg_mean
  541 +
  542 + # Compute flux correction factor based on average TS
  543 + correct = 1.0
  544 + if self.m_use_ts:
  545 + if ts > 0:
  546 + correct = math.sqrt(self.m_ts_thres/ts)
  547 + else:
  548 + if cflux > 0 and cflux_err > 0:
  549 + correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
  550 +
  551 + # Compute extrapolated fluxes
  552 + flux = correct * cflux
  553 + pflux = correct * pflux
  554 + eflux = correct * eflux
  555 + diffSens = correct * diffSens
  556 + flux_value.append(flux)
  557 + pflux_value.append(pflux)
  558 + eflux_value.append(eflux)
  559 + diffSens_value.append(diffSens)
  560 +
  561 + # Write background and test source fit results
  562 + if self.logExplicit():
492 563 self.log.parformat("Photon flux")
493 564 self.log(pflux)
494 565 self.log(" ph/cm2/s\n")
... ... @@ -498,6 +569,9 @@ class cssens(GApplication):
498 569 self.log.parformat("Crab flux")
499 570 self.log(cflux*1000.0)
500 571 self.log(" mCrab\n")
  572 + self.log.parformat("Differential sensitivity")
  573 + self.log(diffSens)
  574 + self.log(" erg/cm2/s\n")
501 575 for model in result_all:
502 576 self.log.parformat("Model")
503 577 self.log(model.name())
... ... @@ -508,52 +582,50 @@ class cssens(GApplication):
508 582 self.log.parformat("Iteration "+str(iter))
509 583 self.log("TS=")
510 584 self.log(ts)
  585 + self.log(" ")
  586 + self.log("corr=")
  587 + self.log(correct)
511 588 self.log(" ")
512 589 self.log(pflux)
513 590 self.log(" ph/cm2/s = ")
514 591 self.log(eflux)
515 592 self.log(" erg/cm2/s = ")
516 593 self.log(cflux*1000.0)
517   - self.log(" mCrab\n")
518   -
519   - # Compute extrapolated fitted prefactor, photon and energy fluxes
520   - correct = 0.0
521   - if self.m_use_ts:
522   - if ts > 0:
523   - correct = math.sqrt(self.m_ts_thres/ts)
524   - else:
525   - if cflux > 0 and cflux_err > 0:
526   - correct = math.sqrt(self.m_ts_thres)/(cflux/cflux_err)
527   - flux = correct * cflux
528   - pflux = correct * pflux
529   - eflux = correct * eflux
530   - flux_value.append(flux)
531   - pflux_value.append(pflux)
532   - eflux_value.append(eflux)
  594 + self.log(" mCrab = ")
  595 + self.log(diffSens)
  596 + self.log(" erg/cm2/s\n")
533 597  
534 598 # Compute sliding average of extrapolated fitted prefactor,
535 599 # photon and energy flux. This damps out fluctuations and
536 600 # improves convergence
537   - flux = 0.0
538   - pflux = 0.0
539   - eflux = 0.0
540   - num = 0.0
  601 + flux = 0.0
  602 + pflux = 0.0
  603 + eflux = 0.0
  604 + diffSens = 0.0
  605 + num = 0.0
541 606 for k in range(self.m_num_avg):
542   - inx = iter - k - 1
  607 + inx = len(flux_value) - k - 1
543 608 if inx >= 0:
544   - flux += flux_value[inx]
545   - pflux += pflux_value[inx]
546   - eflux += eflux_value[inx]
547   - num += 1.0
548   - flux /= num
549   - pflux /= num
550   - eflux /= num
  609 + flux += flux_value[inx]
  610 + pflux += pflux_value[inx]
  611 + eflux += eflux_value[inx]
  612 + diffSens += diffSens_value[inx]
  613 + num += 1.0
  614 + flux /= num
  615 + pflux /= num
  616 + eflux /= num
  617 + diffSens /= num
551 618  
552 619 # Compare average flux to last average
553 620 if iter > self.m_num_avg:
554 621 if test_flux > 0:
555 622 ratio = flux/test_flux
556   - if ratio >= 0.99 and ratio <= 1.01:
  623 +
  624 + # We have 2 convergence criteria:
  625 + # 1. The average flux does not change
  626 + # 2. The flux correction factor is small
  627 + if ratio >= 0.99 and ratio <= 1.01 and \
  628 + correct >= 0.9 and correct <= 1.1:
557 629 if self.logTerse():
558 630 self.log(" Converged ("+str(ratio)+")\n")
559 631 break
... ... @@ -586,6 +658,9 @@ class cssens(GApplication):
586 658 self.log.parformat("Crab flux")
587 659 self.log(cflux*1000.0)
588 660 self.log(" mCrab\n")
  661 + self.log.parformat("Differential sensitivity")
  662 + self.log(diffSens)
  663 + self.log(" erg/cm2/s\n")
589 664 self.log.parformat("Number of simulated events")
590 665 self.log(nevents)
591 666 self.log("\n")
... ... @@ -618,7 +693,8 @@ class cssens(GApplication):
618 693  
619 694 # Store result
620 695 result = {'loge': loge, 'emin': emin.TeV(), 'emax': emax.TeV(), \
621   - 'crab': flux, 'pflux': pflux, 'eflux': eflux}
  696 + 'crab': flux, 'pflux': pflux, 'eflux': eflux, \
  697 + 'diffSens': diffSens}
622 698  
623 699 # Return result
624 700 return result
... ...