#!/bin/sh
# \
exec oagtclsh "$0" "$@"


set auto_path [linsert $auto_path 0 /usr/local/oag/apps/lib/$env(HOST_ARCH)]
set auto_path [linsert $auto_path 0 /usr/local/oag/lib_patch/$env(HOST_ARCH)]
APSStandardSetup

set CVSRevisionAuthor "\$Revision: 1.12 $ \$Author: borland $"

set usage "usage: processCTRIFScan -fileRoot <string> -totalCharge <pC> -lowFreqFitLimit <1/microns> -highFreqFitLimit <1/microns> -plotOnly <1 | 0> -baseSubtract <1 | 0>"
set fileRoot ""
set totalCharge 3000.0
set lowFreqFitLimit 0.002
set highFreqFitLimit 0.004
set totalMacroPulseLength 12.0
set baseSubtract 1
set plotOnly 1
set args $argv

if {[APSStrictParseArguments {fileRoot totalCharge lowFreqFitLimit highFreqFitLimit totalMacroPulseLength plotOnly baseSubtract}] || ![string length $fileRoot] || ![string length ${totalCharge}] || ![string length ${totalMacroPulseLength}] || ![string length ${plotOnly}] || ![string length ${baseSubtract}]} {
    puts stderr $usage
    exit 1
}

if [catch {exec sddsconvert ${fileRoot}.sdds -pipe=out -topage=1 \
             | sdds2stream -pipe -parameter=CTRSystem} CTRSystem] {
    set CTRSystem LI
}

# Initial processing, obtain the FFT, take sqrt(FFT) and deriv(FFT).

    if $baseSubtract {
        if [catch {exec sddsbaseline ${fileRoot}.sdds ${fileRoot}.base -col=${CTRSystem}:CTR:scanDataWF \
                     -select=endpoints=5 -method=fit} result] {
            return -code error "error (1): $result"}
        set baseFile ${fileRoot}.base
    } else {
        set baseFile ${fileRoot}.sdds
    }

if [catch {exec sddsprocess -pipe=out $baseFile -filter=col,Index,0,160 \
             "-redefine=col,${CTRSystem}:CTR:scanPositionsWF,${CTRSystem}:CTR:scanPositionsWF 2.0 *" \
             | sddsfft -pipe -col=${CTRSystem}:CTR:scanPositionsWF,${CTRSystem}:CTR:scanDataWF \
             | sddsprocess -pipe=in ${fileRoot}.fft \
             "-define=col,rho,FFT${CTRSystem}:CTR:scanDataWF sqrt,symbol=\$gr\$r\$n"
    exec sddsderiv ${fileRoot}.fft ${fileRoot}.deriv \
             -differentiate=rho -versus=f \
             -SavitzkyGolay=1,1,1
    exec sddsxref ${fileRoot}.fft ${fileRoot}.deriv -take=rhoDeriv -nowarning
    exec rm ${fileRoot}.deriv ${fileRoot}.fft~} result] {
    return -code error "error (1): $result"
}


# Initial processing, obtain the FFT, Remove all data for wavelengths above 1 mm (< 0.001 micron^-1), 
# Repair spectrum by quadratic extrapolation to DC, take sqrt(FFT) and deriv(FFT).

if [catch {exec sddsprocess ${fileRoot}.fft ${fileRoot}_HighFreq.fft \
             "-filter=col,f,-0.001,${lowFreqFitLimit},!"
    exec sddsprocess ${fileRoot}.fft ${fileRoot}_LowFreq.fft \
             "-filter=col,f,-0.001,${lowFreqFitLimit}"
    exec sddspfit ${fileRoot}.fft ${fileRoot}_FitFreqs.fft \
             -col=f,rho -orders=0,2 \
             -range=${lowFreqFitLimit},${highFreqFitLimit}
    exec sddsxref -pipe=out ${fileRoot}_LowFreq.fft ${fileRoot}_FitFreqs.fft \
             -transfer=param,Intercept,Curvature "-leave=*" \
             | sddsprocess -pipe=in ${fileRoot}_LowFreqEvalFit.fft \
             "-redefine=col,rho,Curvature f sqr * Intercept +"
    exec sddscombine -pipe=out ${fileRoot}_HighFreq.fft ${fileRoot}_LowFreqEvalFit.fft \
             -merge \
             | sddssort -pipe=in ${fileRoot}_DCExt.fft \
             -col=f,increasing
    exec sddsderiv ${fileRoot}_DCExt.fft ${fileRoot}_DCExt.deriv \
             -differentiate=rho -versus=f \
             -SavitzkyGolay=1,1,1
    exec sddsxref ${fileRoot}_DCExt.fft ${fileRoot}_DCExt.deriv -take=rhoDeriv -nowarning
    exec rm ${fileRoot}_DCExt.deriv ${fileRoot}_DCExt.fft~} result] {
    return -code error "error (2): $result"
}

proc plotBunchSpectrumAndFormFactor {rootName} {

    global lowFreqFitLimit highFreqFitLimit CTRSystem baseSubtract

    if $baseSubtract {
        set baseFile ${rootName}.base
    } else {
        set baseFile ${rootName}.sdds
    }
    if [string compare [file dirname $rootName] .]==0 {
        set title [pwd]/$rootName
    } else {
        set title $rootName
    }
    exec sddsplot -title=$title -yscales=names \
             -col=Frequency,rhoAtFrequency ${rootName}_Phase.sdds -graph=line,type=0 "-legend=specified=DC Uncorrected" \
             -col=Frequency,rhoAtFrequency ${rootName}_Phase.sdds -graph=symbol,subtype=0,scale=2 \
             -col=Frequency,rhoAtFrequency ${rootName}_DCExt_Phase.sdds -graph=line,type=1 "-legend=specified=DC Corrected" \
             -col=Frequency,rhoAtFrequency ${rootName}_DCExt_Phase.sdds -graph=symbol,subtype=1,scale=2 \
             "-string=\$s5\$e,x=$lowFreqFitLimit,y=0,scale=1,justify=rc,angle=90,linetype=4" \
             "-string=\$s5\$e,x=$highFreqFitLimit,y=0,scale=1,justify=rc,angle=90,linetype=6" \
             "-topline=Amplitude Spectrum Corrected and Uncorrected at DC" -end \
             -col=Frequency,Psim ${rootName}_Phase.sdds -graph=line,type=0 "-legend=specified=DC Uncorrected" \
             -col=Frequency,Psim ${rootName}_Phase.sdds -graph=symbol,subtype=0,scale=2 \
             -col=Frequency,Psim ${rootName}_DCExt_Phase.sdds -graph=line,type=1 "-legend=specified=DC Corrected" \
             -col=Frequency,Psim ${rootName}_DCExt_Phase.sdds -graph=symbol,subtype=1,scale=2 \
             "-string=\$s5\$e,x=$lowFreqFitLimit,y=0,scale=1,justify=rc,angle=90,linetype=4" \
             "-string=\$s5\$e,x=$highFreqFitLimit,y=0,scale=1,justify=rc,angle=90,linetype=6" \
             "-topline=Minimal Phase Spectrum Corrected and Uncorrected at DC" -end \
             -col=z,Sz ${rootName}_DistFinal.sdds -graph=line,type=0 "-legend=specified=DC Uncorrected" \
             -col=z,Sz ${rootName}_DistFinal.sdds -graph=symbol,scale=2,subtype=0 \
             -col=z,Sz ${rootName}_DCExt_DistFinal.sdds -graph=line,type=1 "-legend=specified=DC Corrected" \
             -col=z,Sz ${rootName}_DCExt_DistFinal.sdds -graph=symbol,scale=2,subtype=1 \
             "-topline=Bunch Distribution Derived From Autocorrelation Using the Minimal Phase Technique" -end \
             -col=t,St ${rootName}_DistFinal.sdds -graph=line,type=0 "-legend=specified=DC Uncorrected" \
             -col=t,St ${rootName}_DistFinal.sdds -graph=symbol,scale=2,subtype=0 \
             -col=t,St ${rootName}_DCExt_DistFinal.sdds -graph=line,type=1 "-legend=specified=DC Corrected" \
             -col=t,St ${rootName}_DCExt_DistFinal.sdds -graph=symbol,scale=2,subtype=1 \
             "-topline=Bunch Distribution Derived From Autocorrelation Using the Minimal Phase Technique" -end \
             -col=t,It -graph=symbol,scale=2,subtype=1 ${rootName}_DCExt_DistFinal.sdds \
             -col=t,It -graph=line ${rootName}_DCExt_DistFinal.sdds \
             -string=@QoString,p=0.7,q=0.9 \
             "-topline=Microbunch Current Distribution From DC Corrected Raw Distribution" -end \
             -col=t,It -graph=symbol,scale=2,subtype=1 ${rootName}_DistFinal.sdds \
             -col=t,It -graph=line ${rootName}_DistFinal.sdds \
             -string=@QoString,p=0.7,q=0.9 \
             "-topline=Microbunch Current Distribution From Un-DC Corrected Raw Distribution" -end \
             -col=${CTRSystem}:CTR:scanPositionsWF,${CTRSystem}:CTR:scanDataWF $baseFile -factor=xMultiplier=2.0 \
             -graph=line,type=0 "-legend=specified=Measured" \
             -col=${CTRSystem}:CTR:scanPositionsWF,${CTRSystem}:CTR:scanDataWF $baseFile -factor=xMultiplier=2.0 \
             -graph=symbol,scale=2,subtype=0 \
             "-yLabel=Autocorrelation (Arb Units)" \
             -col=z,SzConv ${rootName}_DistFinalConv.sdds \
             -graph=line,type=1 "-legend=specified=DC Uncorrected" \
             -col=z,SzConv ${rootName}_DistFinalConv.sdds \
             -graph=symbol,subtype=1 "-yLabel=Autocorrelation (Arb Units)" \
             -col=z,SzConv ${rootName}_DCExt_DistFinalConv.sdds \
             -graph=line,type=2 "-legend=specified=DC Corrected" \
             -col=z,SzConv ${rootName}_DCExt_DistFinalConv.sdds \
             -graph=symbol,subtype=2 \
             "-topline=Measured Autocorrelation Compared to Derived Autocorrelation" \
             "-xLabel=Optical Path Difference (microns)" \&
}

proc computeBunchFormFactor {rootName totalCharge totalMacroPulseLength} { 

    if {![file exists CTRIFScan0_DCExt_DistFinal_Int.sdds]} {

    set frequencyList [APSGetSDDSColumn -fileName ${rootName}.fft -column f]
    set frequencySpacing [APSGetSDDSParameter -fileName ${rootName}.fft -parameter fftFrequencySpacing]
    set index 0
    set fileList ""

# Compute the minimal phase vs frequency.

    foreach frequency $frequencyList {
        set upperLimit [expr $frequency + ($frequencySpacing / 10.0)]
        set lowerLimit [expr $frequency - ($frequencySpacing / 10.0)]
        exec sddsprocess ${rootName}.fft ${rootName}_${index}.fft \
          "-define=param,Frequency,$frequency,units=1/microns" \
          -process=rho,maximum,rhoAtFrequency,functionOf=f,lowerLimit=${lowerLimit},upperLimit=${upperLimit},symbol=\$gr\$r\$n\(\$gw\$r\$n\) \
          -process=rhoDeriv,maximum,rhoDerivAtFrequency,functionOf=f,lowerLimit=${lowerLimit},upperLimit=${upperLimit} \
          "-define=col,integrand,f Frequency \== \? rhoDerivAtFrequency rhoAtFrequency / pi / -1.0 * : rho rhoAtFrequency / ln f sqr Frequency sqr - / -2.0 * Frequency * pi / \$ "
        lappend fileList ${rootName}_${index}.fft
        incr index
    }
    eval exec sddscombine ${fileList} ${rootName}_Proc.fft -overWrite
    eval exec rm ${fileList}
    
    puts ""
    puts "Computing the minimal phase for each frequency."
    puts ""

    exec sddsinteg -pipe=out ${rootName}_Proc.fft \
      -integrate=integrand -versus=f \
      | sddsprocess -pipe \
      "-process=integrandInteg,last,Psim,symbol=\$gY\$r\$bm\$n\(\$gw\$r\$n\)" \
      -convertunits=parameter,rhoAtFrequency,microns,,1 \
      -convertunits=parameter,Psim,radians,,1 \
      | sddscollapse -pipe \
      | sddsprocess -pipe=in ${rootName}_Phase.sdds \
      "-define=col,rhoAtFrequencySqr,rhoAtFrequency sqr"

    exec sddsfft ${rootName}_Phase.sdds ${rootName}_Phase.fft \
      -col=Frequency,rhoAtFrequencySqr
    
# Compute the bunch distribution function using the minimal phase and rho vs frequency.

    set fileList ""
    for {set z -100} {$z <= 1000} {incr z 5} {
        exec sddsprocess ${rootName}_Phase.sdds ${rootName}_$z.sdds \
          "-define=param,z,$z,units=microns,symbol=z" \
          "-define=col,integrand,Psim 2.0 pi * Frequency * $z * - cos rhoAtFrequency * pi /"
        lappend fileList ${rootName}_${z}.sdds
    }
    eval exec sddscombine ${fileList} ${rootName}_Dist.sdds -overWrite
    eval exec rm ${fileList}

    puts ""
    puts "Computing the bunch distribution function."
    puts ""

    exec sddsinteg -pipe=out ${rootName}_Dist.sdds \
      -integrate=integrand -versus=Frequency \
      | sddsprocess -pipe \
      "-process=integrandInteg,last,Sz,symbol=S\(z\)" \
      "-convertunits=parameter,Sz,,1/microns,1" \
      -define=param,St,Sz,symbol=S\(t\) \
      "-define=param,t,z 300 /,units=ps,symbol=t" \
      | sddscollapse -pipe=in ${rootName}_DistFinal.sdds

    exec sddsconvolve ${rootName}_DistFinal.sdds ${rootName}_DistFinal.sdds ${rootName}_DistFinalConv.sdds \
      -signalColumns=z,Sz -responseColumns=z,Sz -outputColumns=z,SzConv

}
    puts ""
    puts "Computing current distribution function."
    puts ""

    if {[catch {exec sddsinteg -pipe=out ${rootName}_DistFinal.sdds \
                  -integrate=St -versus=t \
                  | sddsprocess -pipe=in ${rootName}_DistFinal_Int.sdds -process=StInteg,maximum,Norm} result]} {
        puts $result
    }

    if {[catch {exec sddsxref ${rootName}_DistFinal.sdds ${rootName}_DistFinal_Int.sdds \
                  "-transfer=param,Norm" "-leave=*" -nowarning} result]} {
        puts $result
    }

    if [exec rpnl "$totalMacroPulseLength 0.0 == ? 1 : 0 \$ "] {
        puts "Total macropulse length cannot be zero.  Setting to default 12 ns."
        set totalMacroPulseLength 12.0
    } elseif [exec rpnl "$totalMacroPulseLength 0.0 < ? 1 : 0 \$ "] {
        puts "Total macropulse length is positive.  Setting to its absolute value."
        set totalMacroPulseLength [expr abs(${totalMacroPulseLength})]
    }
              
    set totalMicroPulses [expr int($totalMacroPulseLength * 2.856)]
    
    if {[catch {exec sddsprocess ${rootName}_DistFinal.sdds -nowarning \
                  "-redefine=col,It,${totalCharge} Norm / ${totalMicroPulses} / St *,units=Amperes,symbol=I(t)" \
                  "-redefine=param,Qo,${totalCharge} ${totalMicroPulses} /,units=pC" \
                  "-reprint=param,QoString,Q\$bo\$n = %0.1 f %s,Qo,Qo.units"} result]} {
        puts $result
    }

    puts ""
    puts "Done."
    puts ""
}

if {!$plotOnly} {
    computeBunchFormFactor ${fileRoot}_DCExt $totalCharge $totalMacroPulseLength
    computeBunchFormFactor $fileRoot $totalCharge $totalMacroPulseLength
}

plotBunchSpectrumAndFormFactor $fileRoot

exit
