sics/peaksearch.tcl

#---------------------------------------------------------------------------
# peaksearch a peak search utility for TRICS using the PSD detectors.
#
# Mark Koennecke, November 2001
#---------------------------------------------------------------------------

proc initPeakSearch {} {
    VarMake ps.phiStart Float User
    ps.phiStart 0
    VarMake ps.phiEnd Float User
    ps.phiEnd 180
    VarMake ps.phiStep Float User
    ps.phiStep 3.

    VarMake ps.chiStart Float User
    ps.chiStart 0
    VarMake ps.chiEnd Float User
    ps.chiEnd 180
    VarMake ps.chiStep Float User
    ps.chiStep 12.

    VarMake ps.omStart Float User
    ps.omStart 0
    VarMake ps.omEnd Float User
    ps.omEnd 30
    VarMake ps.omStep Float User
    ps.omStep 3.

    VarMake ps.sttStart Float User
    ps.sttStart 5
    VarMake ps.sttEnd Float User
    ps.sttEnd 70
    VarMake ps.sttStep Float User
    ps.sttStep 3.

    VarMake ps.threshold Int User
    ps.threshold 30
    VarMake ps.steepness Int User
    ps.steepness 3
    VarMake ps.window Int User
    ps.window 7
    VarMake ps.cogwindow Int User 
    ps.cogwindow 60
    VarMake ps.cogcontour Float User
    ps.cogcontour .2

    VarMake ps.countmode Text User
    ps.countmode monitor

    VarMake ps.preset Float User
    ps.preset 1000

    Publish ps.phirange User
    Publish ps.chirange User
    Publish ps.omrange User
    Publish ps.sttrange User
    Publish ps.countpar User
    Publish ps.maxpar   User
#------- these are for debugging only!
    Publish checkomega User
    Publish optimizepeak User
}
#--------------------------------------------------------------------------
proc ps.phirange args {
    if { [llength $args] >= 3 } {
	ps.phiStart [lindex $args 0]
	ps.phiEnd  [lindex $args 1]
        ps.phiStep [lindex $args 2]
    }
    clientput "Peak Search Phi Range:"
    return [format "   Start = %6.2f, End = %6.2f, Step = %6.2f" \
           [SplitReply [ps.phiStart]] [SplitReply [ps.phiEnd]] \
           [SplitReply [ps.phiStep]]] 	
}
#--------------------------------------------------------------------------
proc ps.chirange args {
    if { [llength $args] >= 3 } {
	ps.chiStart [lindex $args 0]
	ps.chiEnd  [lindex $args 1]
        ps.chiStep [lindex $args 2]
    }
    clientput "Peak Search Chi Range:"
    return [format "   Start = %6.2f, End = %6.2f, Step = %6.2f" \
           [SplitReply [ps.chiStart]] [SplitReply [ps.chiEnd]] \
           [SplitReply [ps.chiStep]]] 	
}
#--------------------------------------------------------------------------
proc ps.omrange args {
    if { [llength $args] >= 3 } {
	ps.omStart [lindex $args 0]
	ps.omEnd  [lindex $args 1]
        ps.omStep [lindex $args 2]
    }
    clientput "Peak Search Omega Range:"
    return [format "   Start = %6.2f, End = %6.2f, Step = %6.2f" \
           [SplitReply [ps.omStart]] [SplitReply [ps.omEnd]] \
           [SplitReply [ps.omStep]]] 	
}
#--------------------------------------------------------------------------
proc ps.sttrange args {
    if { [llength $args] >= 3 } {
	ps.sttStart [lindex $args 0]
	ps.sttEnd  [lindex $args 1]
        ps.sttStep [lindex $args 2]
    }
    clientput "Peak Search Two Theta Range:"
    return [format "   Start = %6.2f, End = %6.2f, Step = %6.2f" \
           [SplitReply [ps.sttStart]] [SplitReply [ps.sttEnd]] \
           [SplitReply [ps.sttStep]]] 	
}
#-------------------------------------------------------------------------
proc ps.countpar args {
    if { [llength $args] >= 2 } {
        if { [catch {counter setmode [lindex $args 0]} msg] != 0} {
	    error "ERROR: Invalid countmode specified"
        } 
	ps.countmode [lindex $args 0]
	ps.preset  [lindex $args 1]
    }
    clientput "Peak Search Count Parameters:"
    return [format "   Mode = %s, Preset = %12.2f" \
           [SplitReply [ps.countmode]] [SplitReply [ps.preset]]] 
}
#-------------------------------------------------------------------------
proc ps.maxpar args {
    if { [llength $args] >= 5 } {
	ps.window [lindex $args 0]
	ps.threshold  [lindex $args 1]
	ps.steepness [lindex $args 2]
	ps.cogwindow  [lindex $args 3]
	ps.cogcontour  [lindex $args 4]
    }
    clientput "Peak Search Maximum Detection Parameters:"
    set t1 [format "   Window = %d, Threshold = %d, Steepness = %d" \
           [SplitReply [ps.window]] [SplitReply [ps.threshold]] \
           [SplitReply [ps.steepness] ]]
    set t2 [format "   COGWindow = %d, COGcontour = %f " \
           [SplitReply [ps.cogwindow]] [SplitReply [ps.cogcontour]]] 
    return [format "%s\n%s" $t1 $t2]
}
#------------------------------------------------------------------------
proc checknewomega {omega maxIntensity} {
  if {[catch {drive $omega} msg] != 0} {
     error $msg
  }
  if {[catch {hmc start [SplitReply [ps.preset]] [SplitReply \
            [ps.countmode]]} msg] != 0} {
     error $msg
  }
  if {[catch {Success}  msg] != 0} {
     error $msg
  }
  if { [catch {lomax cog $hm $x $y} result] != 0} {
	     error "Failed to calculate COG"
  }
  if {[lindex $result 2] > $maxIntensity } {
      return $result
  } else {
      return 0
  }
}
#------------------------------------------------------------------------
proc optimizepeak {hm x y} {
    if { [catch {lomax cog $hm $x $y} result] != 0} {
	error "Failed to calculate COG"
    }
    set xMax [lindex $result 0]
    set ymax [lindex $result 1]
    set maxIntensity [lindex $result 2]
    set maxOmega [SplitReply [om]]
    set startOmega $maxOmega
#--------- move to positive omega until maximum found
    while {1} {
	set newOm [expr [SplitReply [om]] + [SplitReply [ps.omStep]]]
        if {[catch {checknewomega $newOm $maxIntensity} result] != 0} {
	    error $result
        }
        if {$result != 0} {
	    set xMax [lindex $result 0]
            set yMax [lindex $result 1]
            set maxIntensity [lindex $result 2]
            set maxOmega [SplitReply [om]]
        } else {
	    break
        }
    }
#--------- if maxOmega is still startOmega then we were on the right
#          side of the peak. In this case try to find the maximum in
#          negative direction
    if {$maxOmega == $startOmega} {
       while {1} {
	   set newOm [expr [SplitReply [om]] - [SplitReply [ps.omStep]]]
           if {[catch {checknewomega $newOm $maxIntensity} result] != 0} {
	      error $result
           }
           if {$result != 0} {
	       set xMax [lindex $result 0]
	       set yMax [lindex $result 1]
	       set maxIntensity [lindex $result 2]
	       set maxOmega [SplitReply [om]]
	   } else {
	       break
	   }
       }
    } 
#----------- print the results we have found
}