#!/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)]
set auto_path [linsert $auto_path 0 $env(LOCO_BINDIR)/tclLib]

APSStandardSetup

set usage "processLocalImpedanceBumps \[-mode <calcResponse|process>\] -workDir <dirname> -measRoot <rootname> -outputRoot <rootname> -bumpType <PositionBump|AngleBump> -bucketList <list> -idNumber <number> -bumpCurrentList <list> -ItoOrbitmm <number> -msSuff <ms|msAve> -rmFile <filename> \[-badBpmList0 <list>\] \[-correctOrbitDrift <0|1>\] \[-useOrbitDriftFilter <0|1>\]"

set args $argv
set mode process
APSParseArguments {mode}
switch -exact $mode {
    process {
	set requiredArgList [list workDir measRoot outputRoot bucketList idNumber bumpCurrentList ItoOrbitmm msSuff rmFile bumpType]
	set argList [concat $requiredArgList badBpmList0 correctOrbitDrift useOrbitDriftFilter]
    }
    calcResponse {
	set requiredArgList [list workDir lteFile idNumber rmFile]
	set argList [concat $requiredArgList paramFile]
    }
    default {return -code error "Wrong mode: $mode"}
}
foreach arg $argList {set $arg ""}
set correctOrbitDrift 0
set useOrbitDriftFilter 0
APSParseArguments $argList
set error 0
foreach arg $requiredArgList {if ![string length [set $arg]] {puts stdout "Argument $arg is missing."; set error 1}}
if $error {puts stderr $usage; exit}

#---------------------------------------------------------------------------------------------------------
proc CalculateOrbitDiff {args} {
    global rmFile msSuff
    APSParseArguments {rootname bucketList nBumps correctOrbitDrift}
    for {set i 0} {$i < $nBumps} {incr i} {
	foreach bucket $bucketList {
	    if $correctOrbitDrift {
		set filename0 $rootname-bump-0-bucket-$bucket.sdds.driftCorrected
		set filename1 $rootname-bump-$i-bucket-$bucket.sdds.driftCorrected
	    } else {
		set filename0 $rootname-bump-0-bucket-$bucket.sdds
		set filename1 $rootname-bump-$i-bucket-$bucket.sdds
	    }
	    exec sddstranspose $filename1 -pipe=out -oldcolumn=BPMName -root=Data \
		| sddsprocess -pipe -match=col,BPMName=S*:$msSuff:y "-reedit=col,BPMName,%@:$msSuff:y@@" \
		| sddsrowstats -pipe -mean=Orbit,Data* \
		| sddsconvert -pipe=in $filename1.orbit -retain=col,BPMName,Orbit
	    if {$i != 0} {
		set filename $rootname-bump-$i-bucket-$bucket.diff
		exec sddsxref $filename0.orbit $filename1.orbit -pipe=out -take=Orbit -rename=col,Orbit=Orbit1 \
		    | sddsxref -pipe $rmFile -take=s -match=BPMName -nowarning \
		    | sddsprocess -pipe=in $filename "-redef=col,OrbitDiff,Orbit1 Orbit -" "-match=col,BPMName=*P3,BPMName=*P4,|"
	    }
	}
    }
}
#---------------------------------------------------------------------------------------------------------
proc RemoveZeroSumBpms {args} {
    global bumpCurrentList bucketList msSuff
    APSParseArguments {workDir measRoot outputRoot badBpmList0}

    set bumpCounter 0
    set fileList ""
    set outputFileList ""
    foreach bumpCurrent $bumpCurrentList {
	foreach bucket $bucketList {
	    lappend fileList $measRoot-bump-$bumpCounter-bucket-$bucket.sdds
#	    exec sddsprocess $measRoot-bump-$bumpCounter-bucket-$bucket.sdds -nowarning -redef=para,BumpCurrent,$bumpCurrent
#	    file delete $measRoot-bump-$bumpCounter-bucket-$bucket.sdds~
	    lappend outputFileList $outputRoot-bump-$bumpCounter-bucket-$bucket.sdds
	}
	incr bumpCounter
    }

    set badBpmList ""
    foreach filename $fileList {
	exec sddstranspose $filename -pipe=out -oldcolumn=BPMpv -root=Data \
	    | sddsprocess -pipe -match=col,BPMpv=S*:sy \
	    | sddsrowstats -pipe -mean=Sum,Data* \
	    | sddsprocess -pipe -filter=col,Sum,0,1e10 \
	    | sddsoutlier -pipe=in $filename.tmp -col=Sum -stdev=4
	#--- We record sum signal only with msAve suffix:
	set badBpms [join [exec sddstranspose $filename -pipe=out -oldcolumn=BPMpv -root=Data \
			       | sddsprocess -pipe -match=col,BPMpv=S*:sy \
			       | sddsselect -pipe $filename.tmp -match=BPMpv -inv \
			       | sddsprocess -pipe "-reedit=col,BPMpv,%@:msAve:sy@@" \
			       | sdds2stream -pipe=in -col=BPMpv] ]
	file delete $filename.tmp
	foreach badBpm $badBpms {
	    if {[lsearch -exact $badBpmList $badBpm] == -1} {lappend badBpmList $badBpm}
	}
    }
    puts stdout "List of bad BPMs (based on the sum signal): $badBpmList"
    set badBpmListPattern ""
    foreach badBpm $badBpmList {lappend badBpmListPattern ${badBpm}*}
    set badBpmListPattern0 ""
    foreach badBpm $badBpmList0 {lappend badBpmListPattern0 ${badBpm}*}
    foreach filename1 $fileList filename2 $outputFileList {
	exec sddsconvert $filename1 $filename2 -del=col,[join [concat $badBpmListPattern $badBpmListPattern0] ,]
    }
}
#---------------------------------------------------------------------------------------------------------
proc RemoveOrbitDrift {args} {
    global bucketList ItoOrbitmm bumpCurrentList msSuff
    APSParseArguments {rootname useOrbitDriftFilter correctOrbitDrift}
    set nBumps [llength $bumpCurrentList]
    set tmpRoot /tmp/[APSTmpString]
    #--- Subtract orbit drift:
###--- To remove ringing in msAve that we had on 2016/04/19
    set clipPoints 1
###---
    foreach bucket $bucketList {
	eval exec sddscombine [glob $rootname-bump-*-bucket-$bucket.sdds] -pipe=out \
	    | sddsprocess -pipe -def=col,CorrectorCurrent,BumpCurrent \"-redef=col,Ybump,BumpCurrent $ItoOrbitmm *,units=mm\" \
	    | sddsprocess -pipe=in $rootname-bucket-$bucket.raw -clip=$clipPoints,0
    }
    if {$correctOrbitDrift == 2} {
	#--- For bunch 0 based drift:
	if $useOrbitDriftFilter {
	    #--- Finding bad BPMs based on bucket 0 std:
	    exec sddscombine $rootname-bucket-0.raw -pipe=out -merge \
		| sddsconvert -pipe -retain=col,S*P\[34\]*$msSuff:y \
		| sddstranspose -pipe -oldcolumn=PVName -root=Orbit \
		| sddsprocess -pipe "-reedit=col,PVName,2S@:@ 20d" \
		| sddsrowstats -pipe -stand=StdDev,Orbit* \
		| tee $tmpRoot.0 \
		| sddsoutlier -pipe=in $tmpRoot.1 -col=StdDev -stdev=3 -pass=3
	    if [catch {join [exec sddsselect $tmpRoot.0 $tmpRoot.1 -pipe=out -match=PVName -inv \
				 | sdds2stream -pipe=in -col=PVName] } badBpmList] {
		return -code error "Error reading bad bpms: $badBpmList"
	    }
	    file delete $tmpRoot.0 $tmpRoot.1
	    if [llength $badBpmList] {
		puts stdout "List of bad BPMs (based on bucket 0 drift): $badBpmList"
		set bpmListPattern ""
		foreach bpmName $badBpmList {lappend bpmListPattern ${bpmName}*}
		foreach bucket $bucketList {
		    exec sddsconvert $rootname-bucket-$bucket.raw $rootname-bucket-$bucket.raw.filtered -del=col,[join $bpmListPattern ,]
		    file copy -force $rootname-bucket-$bucket.raw.filtered $rootname-bucket-$bucket.raw
		    file delete $rootname-bucket-$bucket.raw.filtered
		}
	    }
	}

	#--- Commented command below generates drift in absolute value, so when subtracting this drift, the resulting values are close 
	#--- to zero since average position is also subtracted. It makes comparison plot before and after drift subtraction look bad.
	#--- So instead of this simple command, the two command below create drift relative to the first point of the measurement,
	#--- then the plotting is good.
	#exec sddsprocess $rootname-bucket-0.raw $rootname-bucket-0.raw.aver -process=S*$msSuff:y,average,%sAver
	exec sddsconvert $rootname-bucket-0.raw $rootname-bucket-0.raw.page0 -keep=1
	exec sddscombine $rootname-bucket-0.raw.page0 $rootname-bucket-0.raw -pipe=out \
	    | sddsprocess -pipe -process=S*$msSuff:y,average,%sAver \
	    | sddscollapse -pipe \
	    | sddsbreak -pipe -rowlimit=1 \
	    | sddschanges -pipe -change=*Aver \
	    | sddsconvert -pipe "-edit=col,*Aver,%@ChangeIn@@" \
	    | sddscombine -pipe -merge \
	    | sddsexpand -pipe=in $rootname-bucket-0.raw.aver

	foreach bucket $bucketList {
	    #--- Non-zero buckets give opposite orbit reading. Something about commutation????
	    if {$bucket == 0} {set sign -1} else {set sign 1}
	    exec sddsxref $rootname-bucket-$bucket.raw $rootname-bucket-0.raw.aver -pipe=out -leave=* -transfer=para,S*Aver \
		| sddsprocess -pipe "-redef=col,%s,%s %sAver $sign * +,select=S*$msSuff:y" \
		| sddsconvert -pipe=in $rootname-bucket-$bucket.driftCorrected -del=para,S*Aver
	}
	eval file delete $rootname-bucket-0.raw.aver $rootname-bucket-0.raw.page0
    } else {
	#--- For zero-bump based drift:
	foreach bucket $bucketList {
	    exec sddscombine $rootname-bucket-$bucket.raw -pipe=out -merge \
		| sddsprocess -pipe -process=Time,aver,Time0 "-redef=col,Time,Time Time0 -" \
		| tee $rootname-drift-bucket-$bucket.sdds \
		| sddsprocess -pipe -filter=col,CorrectorCurrent,0,0 \
		| sddsmpfit -pipe=in $rootname-drift-bucket-$bucket.mpfit -indep=Time -dep=S*:$msSuff:y
	    exec sddsxref $rootname-drift-bucket-$bucket.sdds $rootname-drift-bucket-$bucket.mpfit -pipe=out \
		-leave=* -transfer=para,S*Slope,S*Intercept \
		| sddsprocess -pipe "-redef=col,%s,%s %sSlope Time * -,select=S*:$msSuff:y" \
		| sddssort -pipe -col=Time \
		| sddsbreak -pipe -change=CorrectorCurrent \
		| sddsconvert -pipe -del=para,S*:y?*,BumpCurrent \
		| sddsprocess -pipe=in $rootname-bucket-$bucket.driftCorrected -process=CorrectorCurrent,first,BumpCurrent
	    file delete $rootname-drift-bucket-$bucket.sdds
	}

	#--- Find and remove BPMs based of smoothness of the drift:
	if $useOrbitDriftFilter {
	    exec sddscollapse $rootname-drift-bucket-0.mpfit -pipe=out \
		| sddstranspose -pipe -oldcolumn=ColumnName -root=Data \
		| sddsprocess -pipe -match=col,ColumnName=S*:* \
		"-edit=col,Suff,ColumnName,3Z:" "-edit=col,BPMName,ColumnName,S@:$msSuff@ 50d" \
		| sddssort -pipe -col=Suff \
		| sddsbreak -pipe -change=Suff \
		| sddsconvert -pipe -del=para,* \
		| sddsprocess -pipe -process=Suff,first,SuffParam \
		| sddstranspose -pipe -newcol=BPMName \
		| sddsprocess -pipe -print=col,Suff,%s,SuffParam \
		| sddscombine -pipe -merge \
		| sddsconvert -pipe -del=para,* -del=col,OldColumnNames \
		| sddstranspose -pipe=in $tmpRoot.1 -newcol=Suff -oldcolumn=BPMName
	    #	exec sddsoutlier $tmpRoot.1 $tmpRoot.2 -col=yReducedChiSquared -std=4 -pass=20
	    # make it easy to pass for now
	    exec sddsoutlier $tmpRoot.1 $tmpRoot.2 -col=yReducedChiSquared -std=4 -pass=2
	    if [catch {join [exec sddsselect $tmpRoot.1 $tmpRoot.2 -pipe=out -match=BPMName -inv \
				 | sdds2stream -pipe=in -col=BPMName] } bpmList] {
		return -code error "Error selecting bpms: $bpmList"
	    }
	    file delete $tmpRoot.1 $tmpRoot.2
	    puts stdout "List of bad BPMs (based on orbit drift fit): $bpmList"
	    set bpmListPattern ""
	    foreach bpmName $bpmList {lappend bpmListPattern ${bpmName}*}
	    for {set i 0} {$i < $nBumps} {incr i} {
		foreach bucket $bucketList {
		    exec sddsconvert $rootname-bump-$i-bucket-$bucket.sdds $tmpRoot.1 -del=col,[join $bpmListPattern ,]
		    exec sddsconvert $rootname-bump-$i-bucket-$bucket.sdds.driftCorrected $tmpRoot.2 -del=col,[join $bpmListPattern ,]
		    file copy -force $tmpRoot.1 $rootname-bump-$i-bucket-$bucket.sdds
		    file copy -force $tmpRoot.2 $rootname-bump-$i-bucket-$bucket.sdds.driftCorrected
		    file delete $tmpRoot.1 $tmpRoot.2
		}
	    }
	    file delete $rootname-drift-bucket-$bucket.mpfit
	}
    }
    for {set i 0} {$i < $nBumps} {incr i} {
	foreach bucket $bucketList {
	    exec sddsconvert $rootname-bucket-$bucket.driftCorrected \
		$rootname-bump-$i-bucket-$bucket.sdds.driftCorrected -keep=[expr $i + 1]
	}
    }
    foreach bucket $bucketList {file delete $rootname-drift-bucket-$bucket.mpfit}
}
#-------------------------------------------------------------------------------------------------------------------
proc CalculateBunchLength {args} {
    set current ""
    APSParseArguments {current}
    if {$current == ""} {
	return -code error "CalculateBunchLength: current not specified"
    }
    set bunchLength [exec rpnl "$current $current ln 0.0346 * 0.1484 + pow 25.1e-12 * c_mks *"]
    return $bunchLength
}
#---------------------------------------------------------------------------------------------------------------
proc CalculateCorrector {args} {
    APSParseArguments {diffVector rmFile correctorList rootname}
    set tmpRoot /tmp/[APSTmpString]
    exec sddsmakedataset $tmpRoot.corlist -col=Corrector,type=string -data=[join $correctorList ,]

    #--- Getting corrector strength from inverted response matrix
    exec sddsselect $rmFile $diffVector -pipe=out -match=BPMName \
	| sddsconvert -pipe -retain=col,BPMName,[join $correctorList ,] \
	| sddsmatrixop -pipe -inv -push=$diffVector -multiply \
	| sddsconvert -pipe -rename=col,doubleColumn0=CorrectorStrength \
	| sddsxref -pipe=in $tmpRoot.corlist $rootname.cor -take=Corrector

    #--- Getting fitted orbit for plotting:
    exec sddsconvert $rmFile -pipe=out -retain=col,BPMName,[join $correctorList ,] \
	| sddsmatrixop -pipe -push=$rootname.cor -mult \
	| sddsxref -pipe $rmFile -take=s,BPMName \
	| sddsconvert -pipe=in $tmpRoot.fit -rename=col,doubleColumn0=OrbitFit

    #--- Calculating residual orbit error:
    exec sddsxref $diffVector $rmFile -pipe=out -take=s -match=BPMName -nowarning \
	| sddsxref -pipe $tmpRoot.fit -take=OrbitFit -match=BPMName -nowarning \
	| sddsprocess -pipe "-redef=col,OrbitDiff,Orbit OrbitFit -" -process=OrbitDiff,rms,OrbitDiffRms \
	-process=Orbit,rms,OrbitRms \
	| sddsprocess -pipe=in $rootname.matrFit

    file delete $tmpRoot.corlist $tmpRoot.fit
}
#-------------------------------------------------------------------------------------------------------------------
proc CalculateResponse {args} {
    APSParseArguments {lteFile paramFile idNumber outputFile}
    set beamEnergy 7000
    set localRoot localImpResponse
    set eleFile $localRoot.ele
    set twiFile $localRoot.twi
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -run_setup [list "lattice $lteFile use_beamline RING p_central_mev $beamEnergy rootname $localRoot concat_order 0"] \
		   -load_parameters [list "filename_list $paramFile"] \
		   -twiss_output [list "filename $twiFile"] \
		   -run_control [list "n_steps 1"] \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE1: $result"
    }
    set command "elegant $eleFile"
    puts stdout "Running elegant to calculate orbit response: $command"
    if [catch {eval exec $command > $eleFile.log} result] {
	return -code error "Error running elegant (see $eleFile.log): $result"
    }
    lappend deleteFiles $eleFile $eleFile.log $twiFile
    set idMidName S${idNumber}END
    set corUpName S${idNumber}B:V1
    set corDownName S[expr ${idNumber} + 1]A:V1
    foreach elementName [list $corUpName $idMidName $corDownName] {
	set paramList [join [exec sddsprocess $twiFile -pipe=out -match=col,ElementName=$elementName \
				 | sdds2stream -pipe=in -col=s,betay,psiy] ]
	foreach [list s0 betay0 psiy0] $paramList {}
	exec sddsprocess $twiFile $twiFile.$elementName -match=col,ElementType=MONI \
	    "-redef=col,Mult,$betay0 betay * sqrt nuy pi * sin / 2 /" \
	    "-redef=col,Cos,s $s0 < ? psiy $psiy0 - pi nuy * + cos : psiy $psiy0 - pi nuy * - cos \$" \
	    "-redef=col,Response,Mult Cos *"
	lappend deleteFiles $twiFile.$elementName
    }
    exec sddsxref $twiFile.$corUpName $twiFile.$corDownName -pipe=out -take=Response -rename=col,Response=Response1 \
	| sddsprocess -pipe "-redef=col,AngleBump,Response Response1 -" \
	| sddsconvert -pipe=in $twiFile.angle -rename=col,Response=$corUpName -rename=col,Response1=$corDownName 
    lappend deleteFiles $twiFile.angle
    exec sddsxref $twiFile.$idMidName $twiFile.angle -pipe=out -take=AngleBump,$corUpName,$corDownName \
	| sddsprocess -pipe -redef=col,S${idNumber}B:V0,Response \
	| sddsconvert -pipe -rename=col,Response=PositionBump -rename=col,ElementName=BPMName \
	| sddsconvert -pipe=in $outputFile -retain=col,BPMName,s,PositionBump,AngleBump,$corUpName,$corDownName,S${idNumber}B:V0 -del=para,*
    eval file delete $deleteFiles
}
#-------------------------------------------------------------------------------------------------------------------
proc CalculateSlopesVsBump {args} {
    global msSuff OAGGlobal
    APSParseArguments {rootname bucketList}
    set twiFile $OAGGlobal(SRLatticesDirectory)/default/aps.twi
    foreach bucket $bucketList {
	set rawFile $rootname-bucket-$bucket.raw
	set rawFileDC $rootname-bucket-$bucket.driftCorrected
	foreach driftFile [list $rawFile $rawFileDC] {
	    exec sddscombine $driftFile -pipe=out -merge \
		| sddsconvert -pipe -retain=col,S*$msSuff:y,Ybump \
		| sddstranspose -pipe -oldcolumn=BPMName \
		| sddsbreak -pipe -rowlimit=1 \
		| sddsprocess -pipe -process=BPMName,first,BPMNameParam \
		| sddstranspose -pipe -root=Orbit \
		| tee $driftFile.byPage.tmp \
		| sddsconvert -pipe -del=col,OldColumnNames -del=para,*Aver \
		| sddsprocess -pipe=in $driftFile.byPage "-match=para,BPMNameParam=Ybump,!" -match=para,BPMNameParam=S*:P\[34\]*
	    exec sddsprocess $driftFile.byPage.tmp $driftFile.byPage.ybump -match=para,BPMNameParam=Ybump
	    exec sddsxref $driftFile.byPage $driftFile.byPage.ybump -pipe=out -take=Orbit -rename=col,Orbit=Ybump -reuse=page \
		| sddspfit -pipe -col=Ybump,Orbit -terms=2 -generate -copyParameters \
		| sddsprocess -pipe "-print=para,PlotTitle,%s: %8.2e +- %8.2e,BPMNameParam,Slope,SlopeSigma" \
		| tee $driftFile.slopeVsBump.fit \
		| sddscollapse -pipe \
		| sddsprocess -pipe "-edit=col,BPMName,BPMNameParam,%@:$msSuff:y@@" \
		| sddsxref -pipe $twiFile -take=s -match=BPMName=ElementName -nowarning \
		| sddsprocess -pipe "-redef=col,RmsRel,RmsResidual Slope / abs" "-redef=col,SlopeRel,SlopeSigma Slope / abs" \
		| sddssort -pipe=in $driftFile.slopeVsBump.summ -col=s
	    file delete $driftFile.byPage.tmp $driftFile.byPage $driftFile.byPage.ybump
	}
    }
}
#-------------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------------
#-------------------------------------------------------------------------------------------------------------------

if {[string first calcResponse $mode] != -1} {
    if [catch {CalculateResponse -lteFile $lteFile -paramFile $paramFile -idNumber $idNumber -outputFile $rmFile} result] {
	puts stdout "CalculateResponse: $result"
	exit 1
    }
    puts stdout "Done. Results are in file $rmFile."
    exit 
} else {
    if ![file exists $rmFile] {puts stdout "Error: No response file $rmFile. Check the file name or calculate the response."; exit}
}

set E 7e9
set R [expr 1104 / 2 / 3.14159]
set beta_simQuad 3 
set beta_imped 4.3 
set idNumber1 [expr $idNumber + 1]
set excludeBPMList [list S${idNumber}B:P3 S${idNumber1}A:P3]
set localRoot $workDir/$outputRoot
set nBumps [llength $bumpCurrentList]

#--- For position bump, correctorList should be either S${idNumber}B:V0 or PositionBump
#--- For angle bump, should be AngleBump
#--- Could be a combination like "S4B:V1 S4B:V0 S5A:V1", but I didn't test it
set correctorList [list $bumpType]
#set correctorList [list S${idNumber}B:V0]

set currentFile $workDir/$outputRoot.current.summary
if ![file exists $currentFile] {
    return -code error "File $currentFile does not exist. Process bucket currents first."
}
#--- Read bucket currents:
foreach bucketNumber $bucketList {
    set bucketCurrentList($bucketNumber) [join [exec sddsprocess $currentFile -pipe=out \
						    -filter=col,BucketNumber,$bucketNumber,$bucketNumber \
						    | sdds2stream -pipe=in -col=Current] ]
}

exec sddsmakedataset $localRoot.excludeBPMs -col=BPMName,type=string -data=[join $excludeBPMList ,]

#--- Removing BPMs with zero sum signal:
if [catch {RemoveZeroSumBpms -measRoot $workDir/$measRoot -outputRoot $localRoot \
	       -badBpmList0 $badBpmList0} result] {
    return -code error "RemoveZeroSumBpms: $result"
}

#--- Remove orbit drifs:
if {$correctOrbitDrift != 0} {
    if [catch {RemoveOrbitDrift -rootname $localRoot -useOrbitDriftFilter $useOrbitDriftFilter -correctOrbitDrift $correctOrbitDrift} result] {
	return -code error "RemoveOrbitDrift: $result"
    }
} else {
    foreach bucket $bucketList {
	eval exec sddscombine [glob $localRoot-bump-*-bucket-$bucket.sdds] -pipe=out \
	    | sddsprocess -pipe=in $localRoot-bucket-$bucket.raw -def=col,CorrectorCurrent,BumpCurrent \
	    \"-redef=col,Ybump,BumpCurrent $ItoOrbitmm *,units=mm\"
    }
}

#--- Calculate orbit differences:
if [catch {CalculateOrbitDiff -rootname $localRoot -bucketList $bucketList \
	       -nBumps $nBumps -correctOrbitDrift $correctOrbitDrift} result] {
    return -code error "CalculateOrbitDiff: $result"
}

#--- Calculate slopes as a function of s:
if [catch {CalculateSlopesVsBump -rootname $localRoot -bucketList $bucketList} result] {
    return -code error "CalculateSlopesVsBump: $result"
}

#--- Calculation of Z for individual bumps:
set zList ""
set fitErrorList ""
for {set bump 1} {$bump < $nBumps} {incr bump} {
    set bucketCurrent0 [lindex $bucketCurrentList(0) $bump]
    set bucket_sigmas0 [CalculateBunchLength -current $bucketCurrent0]
    set bumpAmplitude [expr [lindex $bumpCurrentList $bump] * $ItoOrbitmm]
    if {$bumpAmplitude ==0} {
	puts stdout "Skipping Bump $bump -- Orbit [format %.2f $bumpAmplitude]mm -- Current_0 [format %.2f $bucketCurrent0]mA -- Bunch length_0 [format %.3e $bucket_sigmas0] m"
	continue
    } else {
	puts stdout "Bump $bump -- Orbit [format %.2f $bumpAmplitude]mm -- Current_0 [format %.2f $bucketCurrent0]mA -- Bunch length_0 [format %.3e $bucket_sigmas0] m"
    }
    for {set i 1} {$i < [llength $bucketList]} {incr i} {
	set bucket [lindex $bucketList $i]
	set bucketCurrent [lindex $bucketCurrentList($bucket) $bump]
	set bucket_sigmas [CalculateBunchLength -current $bucketCurrent]
	set filename $localRoot-bump-$bump-bucket-$bucket.diff
	
	exec sddsselect $filename -pipe=out $localRoot.excludeBPMs -match=BPMName -inv \
	    | sddsprocess -pipe "-reedit=col,BPMName,%@:$msSuff:y@@" \
	    | sddsconvert -pipe=in $outputRoot.diffVector -retain=col,BPMName,OrbitDiff \
	    -rename=col,OrbitDiff=Orbit
	
	set diffVector $outputRoot.diffVector
	#--- Run fit three times to allow for outlier removal
	for {set k 0} {$k < 3} {incr k} {
	    if [catch {CalculateCorrector -diffVector $diffVector -rmFile $rmFile -correctorList $correctorList \
			   -rootname $localRoot} result] {
		return -code error "CalculateCorrector: $result"
	    }
	    set corrStrength [exec sdds2stream $localRoot.cor -col=CorrectorStrength]
	    set fitError [exec sdds2stream $localRoot.matrFit -para=OrbitDiffRms]
	    set orbitRms [exec sdds2stream $localRoot.matrFit -para=OrbitRms]
	    exec sddsoutlier $localRoot.matrFit $localRoot.outl -col=OrbitDiff -stdev=3 -passes=3
	    set badBpmList [join [exec sddsselect $localRoot.matrFit $localRoot.outl -pipe=out -match=BPMName -inv \
				      | sdds2stream -pipe=in -col=BPMName] ]
	    exec sddsselect $diffVector $localRoot.outl $localRoot.diffVector.$k -match=BPMName
	    set diffVector $localRoot.diffVector.$k
	    file delete $localRoot.outl
	    lappend deleteFiles $localRoot.diffVector.$k
	}
	eval file delete $deleteFiles
	exec sddsprocess $localRoot.matrFit $localRoot-bump-$bump-bucket-$bucket.matrFit \
	    -def=para,BumpNumber,$bump,type=long -def=para,BucketNumber,$bucket,type=long \
	    -def=para,CorrKick,$corrStrength,units=mrad -def=para,BumpAmplitude,$bumpAmplitude,units=mm \
	    -def=para,BunchCurrentDelta,[expr $bucketCurrent0 - $bucketCurrent],units=mA \
	    -def=para,BunchCurrentOverSigmaSDelta,[expr $bucketCurrent / $bucket_sigmas - $bucketCurrent0 / $bucket_sigmas0],units=mA/m \
	    -def=para,BunchCurrent,$bucketCurrent,units=mA -def=para,Sigma_s,$bucket_sigmas,units=m \
	    "-print=para,PlotTopline,Bucket [format %1d $bucket] -- Bump [format %1d $bump]" \
	    "-print=para,PlotTitle,Fit error: [format %10.2e $fitError]"
	lappend fitFileList $localRoot-bump-$bump-bucket-$bucket.matrFit
	
	if $i {
	    set outList [exec sddsprocess $localRoot.cor -pipe=out \
			     "-redef=col,Z,CorrectorStrength $E * 2 / $R / $beta_simQuad * $beta_imped / $bumpAmplitude / $bucketCurrent $bucket_sigmas / $bucketCurrent0 $bucket_sigmas0 / - / abs" \
			     | sdds2stream -pipe=in -col=CorrectorStrength,Z]
	    puts stdout "Bucket [format %4d $bucket]; I= [format %5.2f $bucketCurrent]mA; Kick [format %9.2e [lindex $outList 0]](mrad); Z= [format %5.1f [lindex $outList 1]](kOhm/m); Orbit rms: [format %8.2e $orbitRms]mm; Fit error: [format %8.2e $fitError](mm); Ratio: [format %5.3f [expr $fitError / $orbitRms]]"
	    lappend zList [lindex $outList 1]
	    lappend fitErrorList $fitError
	}
	file delete $localRoot.cor $localRoot.corlist $localRoot.diffVector $localRoot.matrFit 
    }
}
set outList [join [exec sddsmakedataset -pipe=out -col=Z,type=double -data=[join $zList ,] \
		       -col=FitError,type=double -data=[join $fitErrorList ,] \
		       | sddsprocess -pipe "-redef=col,Numerator,Z FitError sqr /" "-redef=col,Denominator,1 FitError sqr /" \
		       -process=Numerator,sum,NumSum -process=Denominator,sum,DenomSum \
		       "-redef=para,AverW,NumSum DenomSum /" \
		       -process=Z,aver,Aver -process=Z,stand,Std \
		       | sdds2stream -pipe=in -para=Aver,Std,AverW] ]
puts stdout "Z= [lindex $outList 0] +- [lindex $outList 1]"
puts stdout "Zweighted = [lindex $outList 2]"

eval exec sddscombine $fitFileList -pipe=out \
    | sddssort -pipe=in $localRoot.orbit.fit -para=BucketNumber

exec sddscollapse $localRoot.orbit.fit -pipe=out \
    | sddsprocess -pipe=in $localRoot.summary "-redef=col,KickPerMA,BunchCurrentDelta 0 == ? 0 : CorrKick BunchCurrentDelta / \$"

#--- 2D fit over all points:
#--- Fit assuming equation
#                             ( I_bucket_i       I_bucket_0 )            2R
#   theta_ij = C * y_bump_j * ( ----------  --   ---------- ), where C = -- * Z
#                             ( sigma_s_i        sigma_s_0  )            E
#
#   Using least square method:
#
#        SUM (theta_ij * y_bump_j * IoverSigmaDelta)
#   C = --------------------------------------------
#        SUM (y_bump_j^2 * IoverSigmaDelta^2)
#
if [catch {exec sddsprocess $localRoot.summary -pipe=out \
	       "-redef=col,thetaxy,BumpAmplitude BunchCurrentOverSigmaSDelta * CorrKick *" \
	       "-redef=col,xy2,BumpAmplitude sqr BunchCurrentOverSigmaSDelta sqr *" \
	       -process=thetaxy,sum,thetaxySum -process=xy2,sum,xy2Sum \
	       "-redef=para,C,thetaxySum xy2Sum /" \
	       | sdds2stream -pipe=in -para=C} coeffC] {
    return -code error "Error calculating 2D fit: $coeffC"
}
set Z2D [expr $coeffC / 2 / $R * $E]
puts stdout "From 2D fit: Z = [format %5.1f $Z2D](kOhm/m)"

set aList [join [exec sddsprocess $localRoot.summary -pipe=out \
		     "-redef=col,thetaxy,BumpAmplitude BunchCurrentOverSigmaSDelta * CorrKick *" \
		     "-redef=col,thetax3y,BumpAmplitude 3 pow BunchCurrentOverSigmaSDelta * CorrKick *" \
		     "-redef=col,x2y2,BumpAmplitude 2 pow BunchCurrentOverSigmaSDelta 2 pow *" \
		     "-redef=col,x4y2,BumpAmplitude 4 pow BunchCurrentOverSigmaSDelta 2 pow *" \
		     "-redef=col,x6y2,BumpAmplitude 6 pow BunchCurrentOverSigmaSDelta 2 pow *" \
		     -process=thetaxy,sum,thetaxySum -process=thetax3y,sum,thetax3ySum \
		     -process=x2y2,sum,x2y2Sum -process=x4y2,sum,x4y2Sum -process=x6y2,sum,x6y2Sum \
		     | sdds2stream -pipe=in -para=thetaxySum,thetax3ySum,x2y2Sum,x4y2Sum,x6y2Sum] ]
exec sddsmakedataset aaa.matr -col=C1,type=double -data=[lindex $aList 2],[lindex $aList 3] \
-col=C2,type=double -data=[lindex $aList 3],[lindex $aList 4]
exec sddsmakedataset aaa.vector -col=C1,type=double -data=[lindex $aList 0],[lindex $aList 1]
set cList [join [exec sddsmatrixop aaa.matr -pipe=out -inv -push=aaa.vector -mult \
		     | sdds2stream -pipe=in -col=doubleColumn0] ]
puts stdout "cList: $cList"
set Z2D [expr [lindex $cList 0] / 2 / $R * $E]
puts stdout "From 2D fit: Z = [format %5.1f $Z2D](kOhm/m)"

eval file delete $fitFileList $localRoot.excludeBPMs
set deleteFiles ""
for {set i 0} {$i < $nBumps} {incr i} {
    foreach bucket $bucketList {
	if $correctOrbitDrift {
	    lappend deleteFiles $localRoot-bump-$i-bucket-$bucket.sdds $localRoot-bump-$i-bucket-$bucket.sdds.driftCorrected \
		$localRoot-bump-$i-bucket-$bucket.sdds.driftCorrected.orbit
	} else {
	    lappend deleteFiles $localRoot-bump-$i-bucket-$bucket.sdds $localRoot-bump-$i-bucket-$bucket.sdds.orbit
	}
    }
}
catch [eval file delete $deleteFiles]

