IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41720


Ignore:
Timestamp:
Jul 14, 2021, 3:30:26 PM (5 years ago)
Author:
eugene
Message:

improved plots, option to restrict to camera-confirmed sources

Location:
trunk/tools/eam/teleff
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/tools/eam/teleff/cammask.pro

    r41719 r41720  
    9292 sprintf output "%s/%s.%s.cam.stf" $OUTDIR $expname $OUTVERSION
    9393 write -fits DATA $output $OUT_FIELDS
     94
     95 # save this name for warp and diff analysis, if desired
     96 $CAMSOURCE = $output
    9497end
    9598

  • trunk/tools/eam/teleff/camstats.sh

    r41716 r41720  
    11
    2 macro camstats.glob
     2macro load.camstats
    33 if ($0 != 2)
    4    echo "USAGE: camstats (glob)"
     4   echo "USAGE: load.camstats (glob)"
    55   break
    66 end
    77 
    8  delete -q Rref Dref Xref Yref Q1ref Xpsf Ypsf Q1psf
     8 # accumulate these fields for comparison plots
     9 $FIELDS = RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF Rref Dref Mref
     10
     11 # delete accumulation vectors
     12 delete -q $FIELDS
    913
    1014 list name -glob $1
    1115 for i 0 $name:n
    1216
     17   echo "reading from $name:$i"
    1318   data $name:$i
    1419   read -fits DATA -sizes
    … …  
    2025   calc.camstats
    2126   if ($i == 0)
    22      set magFound_sum = magFound
     27     set numTotal_sum = numTotal
    2328     set numFound_sum = numFound
    24      set numTotal_sum = numTotal
    2529     set numPoss_sum  = numPoss
    26      set numOK_sum    = numOK
     30     set numGood_sum  = numGood
    2731   else
    28      set numOK_sum    = numOK    + numOK_sum
     32     set numTotal_sum = numTotal + numTotal_sum
     33     set numFound_sum = numFound + numFound_sum
    2934     set numPoss_sum  = numPoss  + numPoss_sum
    30      set numFound_sum = numFound + numFound_sum
    31      set numTotal_sum = numTotal + numTotal_sum
     35     set numGood_sum  = numGood  + numGood_sum
    3236   end
    33    concat RrefSS Rref
    34    concat DrefSS Dref
    35    concat XrefSS Xref
    36    concat YrefSS Yref
    37    concat Q1refSS Q1ref
    38    concat X_PSF_found Xpsf
    39    concat Y_PSF_found Ypsf
    40    concat Q1_PSF_found Q1psf
     37   foreach field $FIELDS
     38     concat $field\_out $field
     39   end
     40
     41   # cleanup by deleting the loaded vectors
     42   delete -q $allfields
    4143 end
    4244   
    43  set magFound = magFound_sum
     45 set numTotal = numTotal_sum
    4446 set numFound = numFound_sum
    45  set numTotal = numTotal_sum
    46  set numOK    = numOK_sum
    47 
    48  set fracFound = numFound_sum / numTotal_sum
    49  set fracPoss  = numPoss_sum  / numTotal_sum
    50  set fracOK    = numOK_sum    / numTotal_sum
    51  set dfracFound = sqrt(numFound_sum) / numTotal_sum
    52  set dfracOK    = sqrt(numOK_sum)    / numTotal_sum
    53 
    54  break -auto off
    55  threshold -q -r magFound fracFound 0.25 -range 25 {magFound[]-1}
    56  break -auto on
    57  $MagLim = $threshold
    58 
    59  plot.camstats
    60 end
    61 
    62 macro camstats.single
    63  if ($0 != 2)
    64    echo "USAGE: camstats (file)"
    65    break
    66  end
    67  
    68  data $1
    69  read -fits DATA -list-fields -q
    70  list tfields -join allfields
    71  read -fits DATA $allfields
    72 
    73  # assumes data has been loaded
    74  # results in magFound, fracFound, numFound
    75  calc.camstats
    76 
    77  break -auto off
    78  threshold -q -r magFound fracFound 0.25  -range 25 {magFound[]-1}
    79  break -auto on
    80  $MagLim = $threshold
    81 
    82  plot.camstats
    83 end
    84 
    85 macro plot.camstats
     47 set numPoss  = numPoss_sum
     48 set numGood  = numGood_sum
     49 delete numTotal_sum numFound_sum numPoss_sum numGood_sum
     50
     51 set fracFound = (numTotal > 0) ? numFound / numTotal : zero(numTotal)
     52 set fracPoss  = (numTotal > 0) ? numPoss  / numTotal : zero(numTotal)
     53 set fracGood  = (numTotal > 0) ? numGood  / numTotal : zero(numTotal)
     54end
     55
     56if (not($?HISTBIN)) set HISTBIN = 0.01
     57macro plot.camstats.resid
     58
     59  clear -s
     60  resize 1800 1800
     61  label -fn helvetica 18
     62  $BST1 = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 0.5 +labelpady 0.0
     63  $BST2 = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 0.5 -labelpady 0.0 +ypad 4.5 +labelpady 4.0 -labels 1001
     64 
     65  section aM 0.00 0.66 0.65 0.33
     66  section aR 0.00 0.33 0.65 0.33
     67  section aD 0.00 0.00 0.65 0.33
     68
     69  section bM 0.65 0.66 0.35 0.33
     70  section bR 0.65 0.33 0.35 0.33
     71  section bD 0.65 0.00 0.35 0.33
     72
     73  set dR = 3600*(RA_PSF - Rref)*dcos(Dref)
     74  set dD = 3600*(DEC_PSF - Dref)
     75  set dM = CAL_PSF_MAG - Mref
     76
     77  set keepGood = (PSF_QF > 0.85)
     78  set keepBright = keepGood && (Mref > 15) && (Mref < 17)
     79
     80  subset dRs = dR where keepBright
     81  subset dDs = dD where keepBright
     82  subset dMs = dM where keepBright
     83
     84  # allow mag and pos histograms to have different bin sizes?
     85  histogram dRs NdRs -0.5 0.5 $HISTBIN -range dx
     86  histogram dDs NdDs -0.5 0.5 $HISTBIN -range dx
     87  histogram dMs NdMs -0.5 0.5 $HISTBIN -range dx
     88
     89  NdRs[0] = 0; NdRs[-1] = 0
     90  NdDs[0] = 0; NdDs[-1] = 0
     91  NdMs[0] = 0; NdMs[-1] = 0
     92
     93  $C3 = 0.0
     94
     95  peak -q dx NdRs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 3*$HISTBIN; vgauss dx NdRs con NdRf
     96  peak -q dx NdDs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 3*$HISTBIN; vgauss dx NdDs con NdDf
     97  peak -q dx NdMs; $C0 = $peakpos; $C2 = $peaknum; $C1 = 3*$HISTBIN; vgauss dx NdMs con NdMf
     98
     99  section aM; lim 13 23 -0.5 0.5; box $BST1
     100  plot CAL_PSF_MAG dM -pt box -sz 0.5 -op 0.2 where keepGood
     101  label -x mag -y "&sd&h mag"
     102 
     103  section aR; lim 13 23 -0.5 0.5; box $BST1
     104  plot CAL_PSF_MAG dR -pt box -sz 0.5 -op 0.2 where keepGood
     105  label -x mag -y "&sd&h R.A."
     106 
     107  section aD; lim 13 23 -0.5 0.5; box $BST1
     108  plot CAL_PSF_MAG dD -pt box -sz 0.5 -op 0.2 where keepGood
     109  label -x mag -y "&sd&h Dec"
     110 
     111  section bM; lim dx NdMs; box $BST2
     112  plot dx NdMs -x hist -lw 3
     113  plot dx NdMf -x line -c red -lw 2
     114  label -x "&sd&h mag" +y "N (15 < mag < 17)"
     115 
     116  section bR; lim dx NdRs; box $BST2
     117  plot dx NdRs -x hist -lw 3
     118  plot dx NdRf -x line -c red -lw 2
     119  label -x "&sd&h R.A" +y "N (15 < mag < 17)"
     120 
     121  section bD; lim dx NdDs; box $BST2
     122  plot dx NdDs -x hist -lw 3
     123  plot dx NdDf -x line -c red -lw 2
     124  label -x "&sd&h Dec" +y "N (15 < mag < 17)"
     125end
     126
     127# assumes we have loaded values from the stf tables
     128# magBin, numTotal, numFound, numPoss, numGood
     129macro plot.camstats.fracs
     130
    86131 clear -s
    87132 resize 1800 1200
    88133 label -fn helvetica 18
     134 $BSTY = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 4.5 +labelpady 4.0
     135
     136 ## calculate the rough detection threshold (50%)
     137 # find the bin with the most detections
     138 peak -q magBin numFound
     139
     140 # search backwards along the fracFound curve until it crosses 50%
     141 break -auto off
     142 threshold -q -r magBin fracFound 0.5 -range $peaknum {magBin[]-1}
     143 break -auto on
     144 $MagLim = $threshold
     145
     146 # fraction error bars
     147 set dfracFound = sqrt(numFound) / numTotal
     148 set dfracGood  = sqrt(numGood)  / numTotal
     149
     150 # plot the detection fractions, using left axis
    89151 section a 0 0 1 1
    90  $BSTY = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 4.5 +labelpady 4.0
    91  lim magFound -0.02 1.02; box $BSTY -ticks 1110;
    92  plot -x hist magFound fracFound -lw 3 -c red70
    93  plot -x hist magFound fracPoss  -lw 3 -c blue70
    94  plot -x hist magFound fracOK    -lw 3 -c red
     152 lim magBin -0.02 1.02; box $BSTY -ticks 1110;
     153
     154 plot -x hist magBin fracFound -lw 3 -c red70
     155 plot -x hist magBin fracPoss  -lw 3 -c blue70
     156 plot -x hist magBin fracGood  -lw 3 -c red
     157
     158 subset tmp = fracGood  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightGood  = $MEDIAN
     159 subset tmp = fracPoss  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightPoss  = $MEDIAN
     160 subset tmp = fracFound where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightFound = $MEDIAN
     161
     162 line -c red  -lt dot -lw 2 100 $FracBrightGood to 0 $FracBrightGood
     163 line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
     164
     165 fprintf "%.1f-pct of all bright sources were detected unmasked" {100*$FracBrightGood}
     166 fprintf "%.1f-pct of unmasked bright sources were found" {100*$FracBrightFound / $FracBrightPoss}
     167 fprintf "%.1f-pct of unmasked bright sources were detected unmasked" {100*$FracBrightGood / $FracBrightPoss}
     168 fprintf "mag limit %.1f" $MagLim
     169
    95170 label -x mag -y "detected fraction"
    96  subset tmp = fracOK    where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightOK    = $MEDIAN
    97  subset tmp = fracPoss  where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightPoss  = $MEDIAN
    98  subset tmp = fracFound where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightFound = $MEDIAN
    99  line -c red  -lt dot -lw 2 100 $FracBrightOK to 0 $FracBrightOK
    100  line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
    101  fprintf "%.1f-pct of all bright sources were detected unmasked" {100*$FracBrightOK}
    102  fprintf "%.1f-pct of unmasked bright sources were found" {100*$FracBrightFound / $FracBrightPoss}
    103  fprintf "%.1f-pct of unmasked bright sources were detected unmasked" {100*$FracBrightOK / $FracBrightPoss}
    104  fprintf "mag limit %.1f" $MagLim
    105 
     171
     172 # plot the total number of detected sources, using the right axis
    106173 section b 0 0 1 1
    107  lim magFound numFound; box $BSTY -ticks 1011 -labels 0001; plot -x hist magFound numFound -c grey70 -lw 2
     174 lim magBin numTotal; box $BSTY -ticks 1011 -labels 0001
     175
     176 plot -x hist magBin numFound -c grey50 -lw 2
     177 plot -x hist magBin numTotal -c grey80 -lw 2
     178
    108179 line -c red -lt dot -lw 2 $MagLim -100 to $MagLim 10000
    109180 label +y "detected number"
    110181
    111182 vstat -q numFound; $Nfound = $TOTAL
    112  vstat -q numOK;    $Nok    = $TOTAL
    113  
    114  fprintf "%d total sources detected, %d unmasked" $Nfound $Nok
     183 vstat -q numGood;  $Ngood  = $TOTAL
     184 
     185 fprintf "%d total sources detected, %d unmasked" $Nfound $Ngood
    115186
    116187 $Xlegend = 0.02; $dXlegend = 0.03
    … …  
    121192end
    122193
     194# generate histograms as a function of magnitude for the following:
     195# number of sources (numTotal)
     196# number of sources which were not masked, possibly detected (numPoss)
     197# number detected (numFound)
     198# number of detections with good quality (numGood)
     199
     200# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
    123201macro calc.camstats
    124202 $Mmin = 10
    125203 $Mmax = 25
    126  $Mdel = 0.5
    127  
    128  delete -q fracFound magFound numFound numTotal numPoss fracPoss numOK fracOK
     204 $Mdel = 0.25
     205 
     206 delete -q  magBin numTotal numPoss numFound numGood
    129207 for mag $Mmin $Mmax $Mdel
    130    concat $mag magFound
    131 
    132    # what fraction of the sources were actually detected?
    133    set keep = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel)
    134    # X_PSF_found is NAN if the source is not found (detected in smf), finite otherwise
    135    subset XinBin = X_PSF_found if keep
    136    vstat -q XinBin
    137    concat {$NPTS/XinBin[]} fracFound
    138    concat XinBin[] numTotal
     208   concat $mag magBin
     209
     210   # sources in mag bin
     211   set keepSrc = (Mref_out > $mag - 0.5*$Mdel) && (Mref_out <= $mag + 0.5*$Mdel)
     212   subset tmp = X_PSF_out if keepSrc
     213   concat tmp[] numTotal
     214
     215   # NOTE: X_PSF_out is NAN if the source is not found (detected in smf), finite otherwise
     216   # vstat ignores NAN values, so NPTS is count of non-NAN values
     217   vstat -q tmp
    139218   concat $NPTS numFound
    140219
     220   # what fraction of the sources could have been detected (mask fraction > 0.85 at position)?
     221   set keepMsk = keepSrc && (Q1ref_out > 0.85)
     222   subset tmp = X_PSF_out if keepMsk
     223   concat tmp[] numPoss
     224
    141225   # what fraction of the sources were actually detected with good quality?
    142    set keep = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel) && (Q1_PSF_found > 0.85) && not(isnan(Q1_PSF_found))
    143    # X_PSF_found is NAN if the source is not found (detected in smf), finite otherwise
    144    subset XinBinOK = X_PSF_found if keep
    145    concat {XinBinOK[]/XinBin[]} fracOK
    146    concat XinBinOK[] numOK
    147 
    148    # what fraction of the sources could have been detected (mask fraction > 0.85 at position)?
    149    set keep = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel) && (Q1refSS > 0.85)
    150    subset XinBinGood = X_PSF_found if keep
    151    concat {XinBinGood[]/XinBin[]} fracPoss
    152    concat XinBinGood[] numPoss
    153  end
     226   set keepDet = keepSrc && (PSF_QF_out > 0.85) && not(isnan(PSF_QF_out))
     227   subset tmp = X_PSF_out if keepDet
     228   concat tmp[] numGood
     229 end
     230
     231 delete -q tmp keepSrc keepMsk keepDet
    154232end
    155233
    … …  
    234312end
    235313
    236 macro load.camstats.glob
     314macro load.camstats.glob.old
    237315 if ($0 != 2)
    238316   echo "USAGE: load.camstats.glob (glob)"
    … …  
    265343end
    266344
     345# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
     346macro test.camstats
     347 $Mmin = 10
     348 $Mmax = 25
     349 $Mdel = 0.25
     350 
     351 delete -q  magBin numTotal numFound
     352 for mag $Mmin $Mmax $Mdel
     353   concat $mag magBin
     354
     355   # sources in mag bin
     356   set keepSrc = (MrefDiff_onchp > $mag - 0.5*$Mdel) && (MrefDiff_onchp <= $mag + 0.5*$Mdel)
     357   subset tmp = RpsfDiff_onchp if keepSrc
     358   concat tmp[] numTotal
     359
     360   # NOTE: X_PSF_out is NAN if the source is not found (detected in smf), finite otherwise
     361   # vstat ignores NAN values, so NPTS is count of non-NAN values
     362   vstat -q tmp
     363   concat $NPTS numFound
     364 end
     365
     366 delete -q tmp keepSrc
     367end

  • trunk/tools/eam/teleff/diffmask.pro

    r41719 r41720  
    1313if (not($?MATCH_RADIUS)) set MATCH_RADIUS = 2.0
    1414
     15if (not($?USE_CAMSOURCE)) set USE_CAMSOURCE = 0
     16
    1517$SMF_FIELDS_DIF = X_PSF Y_PSF RA_PSF DEC_PSF CAL_PSF_MAG PSF_QF PSF_QF_PERFECT
    1618$REF_FIELDS_DIF = Rref Dref Mref
    … …  
    5052  # load the reference source table
    5153  # how to choose the name of the filter?
    52   data $refsrcs; read -fits DATA $REF_FIELDS_DIF
     54  data $refsrcs;
     55  if ($USE_CAMSOURCE)
     56   read -fits DATA Rref_out Dref_out Mref_out
     57   set Rref = Rref_out
     58   set Dref = Dref_out
     59   set Mref = Mref_out
     60  else
     61    read -fits DATA $REF_FIELDS_DIF
     62  end
     63
    5364  set Nref = ramp(Rref)
    5465  if ($VERBOSE) echo "loaded Rref[] sources from reference $refsrcs"

  • trunk/tools/eam/teleff/teleff.pro

    r41719 r41720  
    7070
    7171  # use a larger match radius for the asteroids
    72   $MATCH_RADIUS = 3.0
     72  # MATCH_RADIUS is in pixels : this should probably be even larger
     73  $MATCH_RADIUS = 8.0
    7374  $OUTVERSION = MPC
    7475  $COPYMASKS = 0
    7576
    76   make.truth.sources gpc1 $expname $mpcfile MPC
     77  make.truth.sources gpc1 $expname $mpcfile MPC; # sets the value of REFSOURCE
    7778  get.camera.eff     gpc1 $expname EU.nightlyscience $REFSOURCE
    7879  get.warp.eff       gpc1 $expname $camID            $REFSOURCE
    7980  get.diff.eff       gpc1 $expname $warpID           $REFSOURCE
     81  # load data from cam stage output:
     82  $USE_CAMSOURCE = 1
     83  get.warp.eff       gpc1 $expname $camID            $CAMSOURCE
     84  get.diff.eff       gpc1 $expname $warpID           $CAMSOURCE
    8085
    8186  return

  • trunk/tools/eam/teleff/warpmask.pro

    r41719 r41720  
    55$OUTDIR     = .
    66$OUTVERSION = XXX
     7
     8if (not($?USE_CAMSOURCE)) set USE_CAMSOURCE = 0
    79
    810$REF_FIELDS_WRP = Rref Dref Mref
    … …  
    4042  # load the reference source table
    4143  # how to choose the name of the filter?
    42   data $refsrcs
    4344  # NOTE: Mref is chosen to match filter for this exposure
    44   data $refsrcs; read -fits DATA $REF_FIELDS_WRP
     45  data $refsrcs                 ;
     46  if ($USE_CAMSOURCE)
     47   read -fits DATA Rref_out Dref_out Mref_out
     48   set Rref = Rref_out
     49   set Dref = Dref_out
     50   set Mref = Mref_out
     51  else
     52    read -fits DATA $REF_FIELDS_WRP
     53  end
     54
    4555  set Nref = ramp(Rref)
    4656  if ($VERBOSE) echo "loaded Rref[] sources from reference $refsrcs"

  • trunk/tools/eam/teleff/warpstats.sh

    r41712 r41720  
    11
    2 macro warpstats.glob
     2macro load.warpstats
    33 if ($0 != 2)
    4    echo "USAGE: warpstats (glob)"
     4   echo "USAGE: load.warpstats (glob)"
    55   break
    66 end
    77 
    8  delete -q Rref Dref Xref Yref Q1ref
     8 # accumulate these fields for comparison plots
     9 $FIELDS = Rref Dref Mref
     10
     11 # delete accumulation vectors
     12 delete -q $FIELDS
    913
    1014 list name -glob $1
    1115 for i 0 $name:n
    1216
     17   echo "reading from $name:$i"
    1318   data $name:$i
    1419   read -fits DATA -sizes
    … …  
    2025   calc.warpstats
    2126   if ($i == 0)
    22      set magFound_sum = magFound
    2327     set numTotal_sum = numTotal
    2428     set numPoss_sum  = numPoss
    2529   else
     30     set numTotal_sum = numTotal + numTotal_sum
    2631     set numPoss_sum  = numPoss  + numPoss_sum
    27      set numTotal_sum = numTotal + numTotal_sum
    2832   end
    29    concat RrefSS Rref
    30    concat DrefSS Dref
    31    concat XrefSS Xref
    32    concat YrefSS Yref
    33    concat Q1refSS Q1ref
     33   foreach field $FIELDS
     34     concat $field\_out $field
     35   end
     36
     37   # cleanup by deleting the loaded vectors
     38   delete -q $allfields
    3439 end
    3540   
    36  set magFound = magFound_sum
    3741 set numTotal = numTotal_sum
    38  set fracPoss  = numPoss_sum  / numTotal_sum
     42 set numPoss  = numPoss_sum
     43 delete numTotal_sum numPoss_sum
    3944
    40  plot.warpstats
    41 
    42  # plot.warpsky
     45 set fracPoss  = (numTotal > 0) ? numPoss  / numTotal : zero(numTotal)
    4346end
    4447
    45 macro warpstats.single
    46  if ($0 != 2)
    47    echo "USAGE: warpstats (file)"
    48    break
    49  end
    50  
    51  data $1
    52  read -fits DATA -list-fields -q
    53  list tfields -join allfields
    54  read -fits DATA $allfields
     48# assumes we have loaded values from the stf tables
     49# magBin, numTotal, numFound, numPoss, numGood
     50macro plot.warpstats.fracs
    5551
    56  # assumes data has been loaded
    57  # results in magFound, fracFound, numFound
    58  calc.warpstats
    59 
    60  plot.warpstats
    61 end
    62 
    63 macro plot.warpstats
    6452 clear -s
    6553 resize 1800 1200
    6654 label -fn helvetica 18
     55 $BSTY = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 4.5 +labelpady 4.0
     56
     57 # plot the detection fractions, using left axis
    6758 section a 0 0 1 1
    68  $BSTY = -tickpad 0.2 -xpad 3.5 -labelpadx 3.0 -ypad 4.5 -labelpady 4.0 +ypad 5.5 +labelpady 5.0
    69  lim magFound -0.02 1.02; box $BSTY -ticks 1110;
    70  plot -x hist magFound fracPoss  -lw 3 -c blue70
    71  label -x mag -y "detectable fraction"
    72  subset tmp = fracPoss where (magFound > 15) && (magFound < 20); vstat -q tmp; $FracBrightPoss = $MEDIAN
     59 lim magBin -0.02 1.02; box $BSTY -ticks 1110;
     60
     61 plot -x hist magBin fracPoss  -lw 3 -c blue70
     62
     63 subset tmp = fracPoss  where (magBin > 15) && (magBin < 20); vstat -q tmp; $FracBrightPoss = $MEDIAN
    7364 line -c blue -lt dot -lw 2 100 $FracBrightPoss to 0 $FracBrightPoss
    74  fprintf "%.1f-pct of sources are unmasked" {100*$FracBrightPoss}
    7565
     66 fprintf "%.1f-pct of bright sources were detectable (unmasked)" {100*$FracBrightPoss}
     67
     68 label -x mag -y "detected fraction"
     69
     70 # plot the total number of detected sources, using the right axis
    7671 section b 0 0 1 1
    77  lim magFound numTotal; box $BSTY -ticks 1011 -labels 0001; plot -x hist magFound numTotal -c grey70 -lw 2
    78  label +y "reference catalog numbers"
     72 lim magBin numTotal; box $BSTY -ticks 1011 -labels 0001
     73
     74 plot -x hist magBin numTotal -c grey80 -lw 2
     75 label +y "detected number"
     76
     77 $Xlegend = 0.02; $dXlegend = 0.03
     78 $Ylegend = 0.95; $dYlegend = 0.02
     79 line -frac $Xlegend $Ylegend to {$Xlegend + $dXlegend} $Ylegend -c blue  -lw 3; textline -frac -justify 5 {$Xlegend + $dXlegend + 0.01} $Ylegend "unmasked source positions (psfqf > 0.85)"; $Ylegend -= $dYlegend
    7980end
    8081
    81 macro plot.warpsky
    82 
    83   dev -n sky
    84   vstat -q Rref; $Ro = $MEDIAN
    85   vstat -q Dref; $Do = $MEDIAN; $PlotRadius = 0.6*($MAX - $MIN)
    86   region -n sky $Ro $Do $PlotRadius
    87   cplot -pt box -sz 0.3 -op 0.1 -c black Rref Dref
    88 
    89   # create t 0 360 0.01
    90   # set dR = 3*dcos(t)/2.0
    91   # set dD = 3*dsin(t)/2.0
    92   # set Dc = $Do + dD
    93   # set Rc = $Ro + dR/dcos(Dc)
    94   # cplot -x line Rc Dc -c blue -lw 2
    95 end
    96 
    97 
     82# relies on the following vectors: Mref_out, Q1ref_out, PSF_QF_out, X_PSF_out
    9883macro calc.warpstats
    9984 $Mmin = 10
    … …  
    10186 $Mdel = 0.25
    10287 
    103  delete -q magFound numTotal numPoss fracPoss
     88 delete -q  magBin numTotal numPoss
    10489 for mag $Mmin $Mmax $Mdel
    105    concat $mag magFound
     90   concat $mag magBin
    10691
    107    # what fraction of the sources were actually detected?
    108    set keepmag = (MrefSS > $mag - 0.5*$Mdel) && (MrefSS <= $mag + 0.5*$Mdel)
    109    subset XinBin = XrefSS if keepmag
    110    concat XinBin[] numTotal
     92   # sources in mag bin
     93   set keepSrc = (Mref_out > $mag - 0.5*$Mdel) && (Mref_out <= $mag + 0.5*$Mdel)
     94   subset tmp = Mref_out if keepSrc
     95   concat tmp[] numTotal
    11196
    11297   # what fraction of the sources could have been detected (mask fraction > 0.85 at position)?
    113    set keepgood = keepmag && (Q1refSS > 0.85)
    114    subset XinBinGood = XrefSS if keepgood
    115    concat {XinBinGood[]/XinBin[]} fracPoss
    116    concat XinBinGood[] numPoss
     98   set keepMsk = keepSrc && (Q1ref_out > 0.85)
     99   subset tmp = Mref_out if keepMsk
     100   concat tmp[] numPoss
    117101 end
     102
     103 delete -q tmp keepSrc keepMsk keepDet
    118104end
    119 
Note: See TracChangeset for help on using the changeset viewer.