Luigi Tibaldo / ctools

Commit 4109dc8fbd6395a0936d0372d0d6878171082c06

Authored by Jürgen Knödlseder
1 parent b4137537

Remove obsolete code from BGDLIXE (#3738)

Inline Side-by-side
Showing 1 changed file with 47 additions and 107 deletions
modules/comscripts/comobsback.py
  View file @4109dc8
@@ -584,12 +584,6 @@ class comobsback(ctools.csobservation): @@ -584,12 +584,6 @@ class comobsback(ctools.csobservation):
584 scrat = drg.copy() 584 scrat = drg.copy()
585 drb = drg.copy() 585 drb = drg.copy()
586 586
587 - # Initialise Chi2 and nincl2 map for bookkeeping  
588 - chi2_map = drg.map().copy()  
589 - nincl2_map = drg.map().copy()  
590 - chi2_map.stack_maps()  
591 - nincl2_map.stack_maps()  
592 -  
593 # Do Phibar normalization of DRG and store result in SCRAT. 587 # Do Phibar normalization of DRG and store result in SCRAT.
594 for i3 in range(nphibar): 588 for i3 in range(nphibar):
595 sum_dre = 0.0 589 sum_dre = 0.0
@@ -614,105 +608,57 @@ class comobsback(ctools.csobservation): @@ -614,105 +608,57 @@ class comobsback(ctools.csobservation):
614 srcatfile = gammalib.GFilename(srcatname) 608 srcatfile = gammalib.GFilename(srcatname)
615 scrat.save(srcatfile, self['clobber'].boolean()) 609 scrat.save(srcatfile, self['clobber'].boolean())
616 610
617 - # Initialise array of  
618 - nincl2_hist = [0.0 for i in range(nincl2+1)]  
619 -  
620 # Make background model 611 # Make background model
621 for i1 in range(nchi): 612 for i1 in range(nchi):
622 for i2 in range(npsi): 613 for i2 in range(npsi):
623 614
624 - # Initialise nincl2 (add one since we reduce at beginning)  
625 - nincl2 = int(nincl/2) + 1  
626 -  
627 - # Loop while nincl2 is larger than 1  
628 - while nincl2 > 1:  
629 -  
630 - # Reduce nincl2  
631 - nincl2 -= 1  
632 -  
633 - # Initialise Chi-squared  
634 - chi2 = 0.0  
635 - ndof = 0.0  
636 -  
637 - # Loop over Phibar  
638 - for i3 in range(nphibar):  
639 -  
640 - # Determine index range for sum. There are two sums that  
641 - # go over [isel1, iex1[ and [ixe2, isel2[  
642 - iex1 = i3 - nexcl2  
643 - iex2 = i3 + nexcl2 + 1  
644 - if nexcl == 0:  
645 - iex2 -= 1  
646 - isel1 = max((i3 - nincl2), 0)  
647 - isel2 = min((i3 + nincl2 + 1), nphibar)  
648 - if isel1 == 0:  
649 - isel2 = min(nincl, nphibar)  
650 - if isel2 == nphibar:  
651 - isel1 = max(nphibar - nincl, 0)  
652 -  
653 - # Initialise sums  
654 - sum_dre = 0.0  
655 - sum_drm = 0.0  
656 - sum_scrat = 0.0  
657 -  
658 - # Compute average over small patch in Chi/Psi  
659 - for j1 in range(max(0, i1 - navgr2), min(nchi, i1 + navgr2 + 1)):  
660 - for j2 in range(max(0, i2 - navgr2), min(npsi, i2 + navgr2 + 1)):  
661 -  
662 - # Take sum over [isel1, iex1[  
663 - if iex1 >= 1:  
664 - for j3 in range(isel1, iex1):  
665 - jpix = j1 + j2*nchi + j3*npix  
666 - sum_dre += dre[jpix]  
667 - sum_drm += drm[jpix]  
668 - sum_scrat += scrat[jpix]  
669 -  
670 - # Take sum over [iex2, isel2[  
671 - if iex2 < nphibar:  
672 - for j3 in range(iex2, isel2):  
673 - jpix = j1 + j2*nchi + j3*npix  
674 - sum_dre += dre[jpix]  
675 - sum_drm += drm[jpix]  
676 - sum_scrat += scrat[jpix]  
677 -  
678 - # Re-normalise SCRAT to derive DRB  
679 - ipix = i1 + i2*nchi + i3*npix  
680 - if sum_scrat > 0.0:  
681 - drb[ipix] = scrat[ipix] * (sum_dre-sum_drm)/sum_scrat  
682 - else:  
683 - drb[ipix] = 0.0  
684 -  
685 - # Compute Chi2  
686 - if drb[ipix] > 0.0:  
687 - arg = dre[ipix] - drm[ipix] - drb[ipix]  
688 - chi2 += arg * arg / drb[ipix]  
689 - ndof += 1.0  
690 -  
691 - # Normalise Chi2  
692 - if ndof > 0.0:  
693 - chi2 /= ndof  
694 -  
695 - # Update maps  
696 - imap = i1 + i2*nchi  
697 - chi2_map[imap] = chi2  
698 - nincl2_map[imap] = nincl2  
699 -  
700 - # Break if Chi2 is below threshold  
701 - if chi2 < 1.1:  
702 - break  
703 -  
704 - # Log normalised Chi2  
705 - #msg = '(Chi,Psi)=(%3d,%3d)' % (i1, i2)  
706 - #value = '%.2f (%d)' % (chi2, nincl2)  
707 - #self._log_value(gammalib.NORMAL, msg, value)  
708 -  
709 - # Update nincl2 histogram  
710 - nincl2_hist[nincl2] += 1 615 + # Loop over Phibar
  616 + for i3 in range(nphibar):
  617 +
  618 + # Determine index range for sum. There are two sums that
  619 + # go over [isel1, iex1[ and [ixe2, isel2[
  620 + iex1 = i3 - nexcl2
  621 + iex2 = i3 + nexcl2 + 1
  622 + if nexcl == 0:
  623 + iex2 -= 1
  624 + isel1 = max((i3 - nincl2), 0)
  625 + isel2 = min((i3 + nincl2 + 1), nphibar)
  626 + if isel1 == 0:
  627 + isel2 = min(nincl, nphibar)
  628 + if isel2 == nphibar:
  629 + isel1 = max(nphibar - nincl, 0)
  630 +
  631 + # Initialise sums
  632 + sum_dre = 0.0
  633 + sum_drm = 0.0
  634 + sum_scrat = 0.0
711 635
712 - # Log nincl2 histogram  
713 - for i in range(nincl/2+1):  
714 - msg = '(Chi,Psi) with nincl2 = %d' % (i)  
715 - self._log_value(gammalib.NORMAL, msg, nincl2_hist[i]) 636 + # Compute average over small patch in Chi/Psi
  637 + for j1 in range(max(0, i1 - navgr2), min(nchi, i1 + navgr2 + 1)):
  638 + for j2 in range(max(0, i2 - navgr2), min(npsi, i2 + navgr2 + 1)):
  639 +
  640 + # Take sum over [isel1, iex1[
  641 + if iex1 >= 1:
  642 + for j3 in range(isel1, iex1):
  643 + jpix = j1 + j2*nchi + j3*npix
  644 + sum_dre += dre[jpix]
  645 + sum_drm += drm[jpix]
  646 + sum_scrat += scrat[jpix]
  647 +
  648 + # Take sum over [iex2, isel2[
  649 + if iex2 < nphibar:
  650 + for j3 in range(iex2, isel2):
  651 + jpix = j1 + j2*nchi + j3*npix
  652 + sum_dre += dre[jpix]
  653 + sum_drm += drm[jpix]
  654 + sum_scrat += scrat[jpix]
  655 +
  656 + # Re-normalise SCRAT to derive DRB
  657 + ipix = i1 + i2*nchi + i3*npix
  658 + if sum_scrat > 0.0:
  659 + drb[ipix] = scrat[ipix] * (sum_dre-sum_drm)/sum_scrat
  660 + else:
  661 + drb[ipix] = 0.0
716 662
717 # Get statistics 663 # Get statistics
718 sum_dre = 0.0 664 sum_dre = 0.0
@@ -727,13 +673,7 @@ class comobsback(ctools.csobservation): @@ -727,13 +673,7 @@ class comobsback(ctools.csobservation):
727 self._log_value(gammalib.NORMAL, 'Events in DRB', sum_drb) 673 self._log_value(gammalib.NORMAL, 'Events in DRB', sum_drb)
728 674
729 # Save DRB 675 # Save DRB
730 - fits = gammalib.GFits()  
731 - drb_hdu = drb.write(fits)  
732 - chi2_hdu = chi2_map.write(fits)  
733 - nincl2_hdu = nincl2_map.write(fits)  
734 - chi2_hdu.extname('Chi-squared')  
735 - nincl2_hdu.extname('nincl2')  
736 - fits.saveto(drbfile, self['clobber'].boolean()) 676 + drb.save(drbfile, self['clobber'].boolean())
737 677
738 # Log creation 678 # Log creation
739 self._log_value(gammalib.NORMAL, 'DRB file created', drbfile.url()) 679 self._log_value(gammalib.NORMAL, 'DRB file created', drbfile.url())