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

set auto_path [linsert $auto_path 0 $env(ELEGANT_LOCAL_LIB_PATH)]
set auto_path [linsert $auto_path 0 $env(ELEGANT_LOCAL_LIB_PATH)]

if 0 {
if [info exists env(LOCO_BINDIR)] {
    set locoBinDir $env(LOCO_BINDIR)
} else {
    puts stderr "Error: LOCO_BINDIR environment variable is not defined"
    exit
}
}
source fittingProcedures.tcl

set acceleratorCode elegant

set usage "calculateTwissCorrection -mode <correct|scan|plot|calcEmit|minEmit|runDqsTask> ..."
set args $argv

set argMode commandString
set scriptParametersFile ""
set verbose 0
set dispersionOnly 0
APSParseArguments {argMode scriptParametersFile verbose dispersionOnly}
if {[string compare $argMode file] == 0} {
    if {![string length $scriptParametersFile]} {
	puts stderr "Arguments: $argv"
	puts stderr $usage
	exit
    } else {
	if [catch {open $scriptParametersFile r} fid] {
	    return -code error "$fid"
	} else {
	    gets $fid args
	    close $fid
	}
    }
}

set mode ""
APSParseArguments {mode}
if ![string length $mode] {
    puts stderr "Arguments: $argv"
    puts stderr $usage
    exit
}

set outputSparsity 1
switch -exact -- $mode {
    correct {
	set requiredArgList [list calcMode lteMeasured beamlineMeasured lteDesired beamlineDesired \
				 outputFile elementListFile splitTasks workDir \
				 useTunes useDQS weightList calculateMatrix calculateInverse useRelative \
				 useInvSolution locoResultFile SVratio corFraction iterations \
				 abortFile useMeasDisp matrixFile]
	set addMeasurementNoise 0
	set doPerturbedMatrix 0
	set measDispFile ""
	APSParseArguments {addMeasurementNoise outputSparsity doPerturbedMatrix measDispFile}
	if $addMeasurementNoise {set requiredArgList [concat $requiredArgList betaNoise dispNoise rmNoise]}
    }
    scan {
	set requiredArgList [list calcMode lteMeasured beamlineMeasured lteDesired beamlineDesired workDir \
				 outputFile elementListFile weightList calculateMatrix useRelative \
				 SVstart SVend SVsteps matrixFile corFraction useDQS splitTasks iterations useTunes \
				 useMeasDisp measDispFile abortFile]
    }
    plot {
	set requiredArgList [list calcMode lteMeasured beamlineMeasured lteDesired beamlineDesired paramFileListAchieved \
				 weightList doneFile workDir twiMeasured twiDesired twiAchieved elementListFile \
				 useTunes useRelative]
    }
    minEmit {
	set requiredArgList [list calcMode workDir lteMeasured beamlineMeasured elementListFile emitArrayList optimTerm rootname]
    }
    calcEmit {
	set requiredArgList [list calcMode workDir lteMeasured beamlineMeasured lteDesired \
				 beamlineDesired paramFileListAchieved emitArrayList]
    }
    runDqsTask {
	set requiredArgList [list calcMode path eleTemplate varList useRelative acceleratorCode elementsUpdateFile]
    }
    default {
	puts stdout $usage
	exit
    }
}

foreach var $requiredArgList {set $var ""}
APSParseArguments $requiredArgList
set missingVarList ""
foreach var $requiredArgList {
    if ![string length [set $var]] { lappend missingVarList $var }
}
if [llength $missingVarList] {
    puts stderr "Arguments: $argv"
    puts stderr "The following variables are missing: $missingVarList"
    exit
}

#------ Not required variables:
set paramFileListMeasured ""
set paramFileListDesired ""
set tmpDir /tmp
set betaTargetFile ""
set beamEnergyMeV 7000
set elemByElemSR 0
APSParseArguments {paramFileListMeasured paramFileListDesired tmpDir betaTargetFile beamEnergyMeV elemByElemSR}

set outputStatusDevice stdout

#-----------------------------------------------------------------------------------------------------------------------

proc ReadListFromFile { args } {
    global quadCorrectors
    APSParseArguments {calcMode listFile}
    #------ See if correctors are in quads
    set quadCorrectors 0
    set paramList [join [exec sddsquery -para $listFile] ]
    if {[lsearch $paramList ElementType] != -1} {
	set elemType [lindex [join [exec sdds2stream $listFile -para=ElementType] ] 1]
	if {[string first QUAD $elemType] != -1} {set quadCorrectors 1}
    }
    #------ Read lists
    set mainList ""
    if {[string compare $calcMode beta] == 0} {
	set paramList [list Quad xBpm yBpm dBpm]
    } else {
	set paramList [list Quad hCorr vCorr xBpm yBpm dBpm]
    }
    foreach param $paramList {
	lappend mainList [join [exec sddsprocess $listFile -pipe=out -match=para,ListName=$param -nowarning \
				    | sdds2stream -pipe=in -col=List] ]
    }
    return $mainList
}

#-----------------------------------------------------------------------------------------------------------------------

proc TwissCorrection { args } {
    
    global workDir tmpDir outputStatusDevice outputStatusFile
    global measuredFiles achievedFiles desiredFiles
    global matrixFile matrixFileFiltered inverseMatrixFile measDifference
    global eleMatrixTemplate elementsUpdateFile measuredBpms measuredDispFile
    global xbpmList ybpmList dbpmList weightList doPerturbedMatrix
    global twissVectorDesired rmsIterationsFile errorMessage

    set adjustVectorValues 0
    set addMeasurementNoise 0
    set outputSparsity 1
    APSParseArguments { calcMode lteMeasured beamlineMeasured paramFileListMeasured lteDesired beamlineDesired paramFileListDesired \
			    outputFile quadList xbpmList ybpmList dbpmList hcorrList vcorrList splitTasks useTunes useDQS \
			    weightList calculateMatrix calculateInverse useRelative useInvSolution locoResultFile \
			    SVratio corFraction iterations outputStatusDevice outputStatusFile abortFile \
			    useMeasDisp measDispFile matrixFile adjustVectorValues addMeasurementNoise outputSparsity returnToBest}
    
    set deleteFiles ""
    OutputStatusMessage "Calculations started..."

    if [catch {InitialPreparations -quadList $quadList -lteMeasured $lteMeasured -beamlineMeasured $beamlineMeasured \
		   -paramFileListMeasured $paramFileListMeasured -lteDesired $lteDesired -useTunes $useTunes \
		   -beamlineDesired $beamlineDesired -paramFileListDesired $paramFileListDesired \
		   -hcorrList $hcorrList -vcorrList $vcorrList} result] {
	return -code error "InitialPreparations: $result"
    }

    #--- Calculate initial twiss vectors...
    OutputStatusMessage "Calculating initial vectors..."
    if [catch {BuildMeasurementVector -calcMode $calcMode -eleFile $desiredFiles(ele) -twissFile $desiredFiles(twiss) \
		   -hyFile $desiredFiles(hyrm) -vxFile $desiredFiles(vxrm)} measVectorDesired] {
	return -code error "BuildMeasurementVector (desired): $measVectorDesired"
    }
    if $adjustVectorValues {
	if [catch {AdjustVectorValues -vectorFile $measVectorDesired} result] {
	    return -code error "AdjustVectorValues: $result"
	}
    }
    if $useMeasDisp {
	if [catch {BuildMeasurementVector -calcMode $calcMode -eleFile $measuredFiles(ele) -dispFile $measDispFile \
		       -twissFile $measuredFiles(twiss) -hyFile $measuredFiles(hyrm) \
		       -vxFile $measuredFiles(vxrm) -addMeasurementNoise $addMeasurementNoise} measVectorMeasured] {
	    return -code error "BuildMeasurementVector: $measVectorMeasured"
	}
    } else {
	if [catch {BuildMeasurementVector -calcMode $calcMode -eleFile $measuredFiles(ele) -twissFile $measuredFiles(twiss) \
		       -hyFile $measuredFiles(hyrm) -vxFile $measuredFiles(vxrm) \
		       -addMeasurementNoise $addMeasurementNoise} measVectorMeasured] {
	    return -code error "BuildMeasurementVector (measured): $measVectorMeasured"
	}
    }
    set nameColumn Rootname
    set dataColumn Measurement
    if [catch {CalculateVectorDifference -vectorFile1 $measVectorDesired -vectorFile2 $measVectorMeasured \
		   -outputFile $measDifference -useRelative $useRelative -nameColumn $nameColumn -dataColumn $dataColumn} result ] {
	return -code error "CalculateVectorDifference: $result"
    }

    if [catch {rmsDisplayOutput -calcMode $calcMode -iteration -1 -diffFile $measDifference -twissFile $measuredFiles(twiss) \
		   -weightList $weightList -dataColumn MeasurementDiff} rmsList] {
	return -code error "rmsDisplayOutput: $rmsList"
    }
    set rmsTotal0 [lindex $rmsList 0]
    set tol [expr $rmsTotal0 / 10000]

    if !$useInvSolution {
	if $calculateMatrix {
	    OutputStatusMessage "Doing the $calcMode response matrix..."
	    if [catch {CalculateTwissResponseMatrix -calcMode $calcMode -matrixFile $matrixFile \
			   -quadList $quadList -useRelative $useRelative -splitTasks $splitTasks -useQsub $useDQS \
			   -doPerturbedMatrix $doPerturbedMatrix -abortFile $abortFile} result] {
		return -code error "CalculateTwissResponseMatrix: $result"
	    }
	}
	if $calculateInverse {
	    OutputStatusMessage "Doing the sddspseudoinverse..."
	    if [catch {CalculateInverse -matrixFile $matrixFile -matrixFileFiltered $matrixFileFiltered \
			   -inverseMatrixFile $inverseMatrixFile -SVratio $SVratio -quadList $quadList \
			   -diffFile $measDifference -useRelative $useRelative} result] {
		return -code error "CalculateInverse: $result"
	    }
	    OutputStatusMessage "sddspseudoinverse is done..."
	}
    }

    set tmpFile /tmp/[APSTmpString]-twissCorr
    set vectorFile $workDir/vector_out
    set rmsTotalList $rmsTotal0
    lappend deleteFiles $tmpFile $vectorFile $measVectorDesired $measVectorMeasured

#--- Main Iterations -------------------------------------------------------------------------------------
    
    if $useInvSolution {set iterations 1}
    set rmsTotalPrev $rmsTotal0
    set rmsTotalBest $rmsTotal0
    set iterBest -1
    set iterError 0
    if $returnToBest {file copy -force $workDir/tmp/elements_update.sdds $workDir/tmp/elements_update.sdds.best}
    for {set I 0} {$I < $iterations} {incr I} {
	
	if !$useInvSolution {
	    exec sddscombine $measDifference -pipe=out -merge \
		| sddsmatrixmult -pipe=in $inverseMatrixFile $vectorFile -commute
	}

	#--- Updating the elements_update file after the iteration...
	if !$useInvSolution {
	    exec sddsxref $elementsUpdateFile $vectorFile -pipe=out -take=MeasurementDiff \
		| sddsprocess -pipe "-redefine=col,ParameterValue,ParameterValue MeasurementDiff $corFraction * + " \
		| sddsconvert -pipe=in $tmpFile -delete=column,MeasurementDiff
	    file copy -force $tmpFile $elementsUpdateFile
	} else {
	    if [catch {exec sddscombine $locoResultFile -pipe=out -merge \
			   | sddsprocess -pipe -match=col,ElementParameter=K1 -nowarning \
			   | sddsprocess -pipe=in $elementsUpdateFile \
			   "-redef=col,ParameterValue,ParameterValue -1 *"} result] {
		OutputStatusMessage "Error updating: $result"; set iterError 1; break
	    }
	}

	#--- Create twissVector after the iteration...
	if [catch {BuildMeasurementVector -calcMode $calcMode -eleFile $achievedFiles(ele) -twissFile $achievedFiles(twiss) \
		       -hyFile $achievedFiles(hyrm) -vxFile $achievedFiles(vxrm)} measVectorAfterIter] {
	    OutputStatusMessage "BuildMeasurementVector: $measVectorAfterIter"
	    set iterError 1; break
	}
	if [catch {CalculateVectorDifference -vectorFile1 $measVectorDesired -vectorFile2 $measVectorAfterIter \
		       -outputFile $measDifference -useRelative $useRelative -nameColumn $nameColumn -dataColumn $dataColumn} result ] {
	    OutputStatusMessage "CalculateVectorDifference: $result"
	    set iterError 1; break
	}
	if {[expr $I - $I / $outputSparsity * $outputSparsity] == 0} {set verbose 1} else {set verbose 0}
	if [catch {rmsDisplayOutput -calcMode $calcMode -iteration $I -diffFile $measDifference -twissFile $achievedFiles(twiss) \
		       -weightList $weightList -dataColumn MeasurementDiff -verbose $verbose} rmsList] {
	    OutputStatusMessage "rmsDisplayOutput: $rmsList"
	    set iterError 1; break
	}
	file delete $measVectorAfterIter

	if [file exists $abortFile] {
	    if [catch { StopExecution } result] {
		return -code error "TwissCorrection: $result"
	    }
	    eval file delete $deleteFiles $measVectorAfterIter
	    return -code error "Interrupted by user."
	}
	file copy -force $elementsUpdateFile $outputFile
	if $returnToBest {
	    if {[lindex $rmsList 0] < $rmsTotalBest} {
		set rmsTotalBest [lindex $rmsList 0]
		set iterBest $I
		file copy -force $workDir/tmp/elements_update.sdds $workDir/tmp/elements_update.sdds.best
	    }
	}
	if {[expr abs($rmsTotalPrev - [lindex $rmsList 0])] < $tol} {
	    OutputStatusMessage "Iterations converged. Exiting."; break
	}
	set rmsTotalPrev [lindex $rmsList 0]
	lappend rmsTotalList $rmsTotalPrev
	exec sddsmakedataset -pipe=out -col=RMS,type=double -data=[join $rmsTotalList ,] \
	    | sddsprocess -pipe=in $rmsIterationsFile -def=col,Iteration,i_row
    }

    #------ Return to best if: iterBest != current iter or there was an error:
    if $returnToBest {
	if {$iterBest != $I || $iterError} {
	    OutputStatusMessage "Returning to best iteration: iteration $iterBest..."
	    file copy -force $workDir/tmp/elements_update.sdds.best $workDir/tmp/elements_update.sdds
	    if [catch {BuildMeasurementVector -calcMode $calcMode -eleFile $achievedFiles(ele) -twissFile $achievedFiles(twiss) \
			   -hyFile $achievedFiles(hyrm) -vxFile $achievedFiles(vxrm)} measVectorAfterIter] {
		eval file delete $deleteFiles
		return -code error "BuildMeasurementVector: $measVectorAfterIter"
	    }
	    if [catch {CalculateVectorDifference -vectorFile1 $measVectorDesired -vectorFile2 $measVectorAfterIter \
			   -outputFile $measDifference -useRelative $useRelative -nameColumn $nameColumn -dataColumn $dataColumn} result ] {
		eval file delete $deleteFiles $measVectorAfterIter
		return -code error "CalculateVectorDifference: $result"
	    }
	    if [catch {rmsDisplayOutput -calcMode $calcMode -iteration $I -diffFile $measDifference -twissFile $achievedFiles(twiss) \
			   -weightList $weightList -dataColumn MeasurementDiff -verbose $verbose} rmsList] {
		return -code error "rmsDisplayOutput: $rmsList"
	    }
	    file copy -force $elementsUpdateFile $outputFile
	}
    }
    eval file delete $deleteFiles $measVectorAfterIter
    return $rmsList
}

#-------------------------------------------------------------------------------------------------------------------------
proc AdjustVectorValues {args} {
    global tmpDir betaTargetFile
    APSParseArguments {vectorFile}
    set tmpRoot $tmpDir/[APSTmpString]
    if [catch {exec sdds2stream $betaTargetFile -para=Plane} result] {
	return -code error "Error reading $betaTargetFile"
    } else {
	if {[llength $result] != 4} {
	    return -code error "Error: $betaTargetFile must be 4 pages long."
	}
    }
    exec sddsxref $vectorFile $betaTargetFile -pipe=out -take=TargetBeta -match=Rootname -nowarning -fillin \
	| sddsprocess -pipe "-redef=col,Measurement,TargetBeta 0 == ? Measurement : TargetBeta \$" \
	| sddsconvert -pipe=in $tmpRoot -del=col,TargetBeta
    file copy -force $tmpRoot $vectorFile
    file delete $tmpRoot
}
#-------------------------------------------------------------------------------------------------------------------------

proc InvertLOCOSolution { args } {
    global outputStatusDevice
    APSParseArguments { outputStatusDevice }
    global locoResultFile

    set locoResultsFile ""
    set dbwin [APSUniqueName .]
    APSDialogBox $dbwin -name "LOCO results file" -width 60 -contextHelp " " -cancelCommand {set locoResultFile ""}
    APSLabeledEntry .le1 -parent $dbwin.userFrame -width 60 \
	-label "LOCO result file:" \
	-textVariable locoResultFile \
	-contextHelp "Name of the file with quadrupoles calculated by LOCO."
    update
    tkwait window $dbwin
    if ![file exists $locoResultFile] {
	return -code error "Error: File $locoResultFile does not exist."
    }
    OutputStatusMessage "LOCO result file has been chosen."
    return $locoResultFile
}

#-------------------------------------------------------------------------------------------------------------------------
proc MakeElementNameFromTagName {args} {
    set tagNameColumn TagName
    set output ""
    APSParseArguments {input output tagNameColumn}
    if {[string length $output] == 0} {set localOutput $input.tmp} else {set localOutput $output}
    #------ elegant can create a mixture of names when some names have #N in it and others don't.
    #------ To handle it, add #1 at the end, then delete second # symbol and after it.
    if [catch {exec sddsprocess $input -pipe=out "-edit=col,ElementName,$tagNameColumn,S?@#@ K" \
		   "-edit=col,TagNameCopy,$tagNameColumn,e i@#1@ a 2S?@#@ 2d" \
		   "-scan=col,ElementOccurence,TagNameCopy,%ld,type=long,edit=Z#" \
		   | sddsconvert -pipe=in $localOutput -del=col,TagNameCopy} result] {
	return -code error "Error making ElementName from TagName ($input -- $tagNameColumn): $result"
    }
    if {[string length $output] == 0} {file copy -force $localOutput $input; file delete $localOutput}
}
#-------------------------------------------------------------------------------------------------------------------------

proc InitialPreparations { args } {

    global workDir tmpDir beamEnergyMeV elemByElemSR acceleratorCode
    global matrixFile matrixFileFiltered inverseMatrixFile measDifference
    global eleMatrixTemplate elementsUpdateFile measuredBpms
    global measuredFiles achievedFiles desiredFiles
    global xbpmList ybpmList dbpmList weightList calcMode correctorFile quadCorrectors

    APSParseArguments {quadList lteMeasured beamlineMeasured paramFileListMeasured lteDesired beamlineDesired \
			    paramFileListDesired useTunes hcorrList vcorrList}

    set tmpRoot $tmpDir/[APSTmpString]-iniPrep
    catch {exec mkdir $workDir/tmp}
    # ----- Required internal files -----------------------
    set matrixFileFiltered $workDir/tmp/twissCorrectionMatrixFiltered.sdds
    set inverseMatrixFile  $workDir/tmp/twissCorrectionMatrixInverse.sdds
    set measDifference     $workDir/tmp/measDifference.sdds
    set eleMatrixTemplate  $workDir/tmp/twissCorrectionMatrixTemplate.ele
    set elementsUpdateFile $workDir/tmp/elements_update.sdds
    set measuredBpms       $workDir/tmp/measuredBpms.sdds
    
    set  desiredFiles(lte)       $lteDesired
    set  desiredFiles(beamline)  $beamlineDesired
    set  desiredFiles(paramList) $paramFileListDesired
    set  desiredFiles(ele)       $workDir/tmp/desired.ele
    set  desiredFiles(outParam)  $workDir/tmp/desired.param
    set  desiredFiles(twiss)     $workDir/tmp/desired.twi
    set  desiredFiles(hyrm)      $workDir/tmp/desired.hyrm
    set  desiredFiles(vxrm)      $workDir/tmp/desired.vxrm

    set measuredFiles(lte)       $lteMeasured
    set measuredFiles(beamline)  $beamlineMeasured
    set measuredFiles(paramList) $paramFileListMeasured
    set measuredFiles(ele)       $workDir/tmp/measured.ele
    set measuredFiles(outParam)  $workDir/tmp/measured.param
    set measuredFiles(twiss)     $workDir/tmp/measured.twi
    set measuredFiles(hyrm)      $workDir/tmp/measured.hyrm
    set measuredFiles(vxrm)      $workDir/tmp/measured.vxrm

    set achievedFiles(lte)       $lteMeasured
    set achievedFiles(beamline)  $beamlineMeasured
    set achievedFiles(paramList) [concat $paramFileListMeasured $elementsUpdateFile]
    set achievedFiles(ele)       $workDir/tmp/achieved.ele
    set achievedFiles(outParam)  $workDir/tmp/achieved.param
    set achievedFiles(twiss)     $workDir/tmp/achieved.twi
    set achievedFiles(hyrm)      $workDir/tmp/achieved.hyrm
    set achievedFiles(vxrm)      $workDir/tmp/achieved.vxrm

    #--- Make twiss file to get element types:
    set eleFile $workDir/tmp/twiss.ele
    set twissFile $workDir/tmp/twiss.twi
    set elegantOutFile $workDir/tmp/twiss.log
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -run_setup [list "lattice $lteDesired use_beamline $beamlineDesired default_order 2 p_central_mev $beamEnergyMeV"] \
		   -load_parameters [list "filename_list \"$paramFileListDesired\" force_occurence_data 1"] \
		   -twiss_output [list "filename $twissFile"] \
		   -run_control [list "n_steps 1"] \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleFile): $result"
    }
    if [catch {exec $acceleratorCode $eleFile > $elegantOutFile} result] {
	return -code error "elegant (ele file: $eleFile; output file: $elegantOutFile): $result"
    }
    if [catch {AddTagNameColumn -filename $twissFile} result] {return -code error "AddTagNameColumn: $result"}
    file delete $eleFile $elegantOutFile

    #--- Create elements_update.sdds file for iterations...
    if [llength $quadList] {
	if [catch {exec sddsmakedataset -pipe=out -col=TagName,type=string -data=[join $quadList ,] \
		       | sddsprocess -pipe=in $elementsUpdateFile -print=col,ParameterMode,differential -def=col,ParameterValue,0} result] {
	    return -code error "Error making elementsUpdateFile: $result"
	}
    } else {
	return -code error "Error: List quadList has zero length."
    }
    if [catch {MakeElementNameFromTagName -input $elementsUpdateFile} result] {
	return -code error "MakeElementNameFromTagName: $result"
    }

    #--- Add ElementType column to elements_update.sdds and make ElementParameter column based on ElementType:
    exec sddsmakedataset $tmpRoot.types -col=ElementType,type=string -data=QUAD,KQUAD,SBEN,CSBEND,CCBEND,MULT -col=ElementParameter,type=string -data=K1,K1,K1,K1,K1,KNL
    exec sddsxref $elementsUpdateFile $twissFile -pipe=out -take=ElementType -match=TagName -nowarning \
	| sddsxref -pipe=in $tmpRoot.types $elementsUpdateFile.tmp -take=ElementParameter -match=ElementType -reuse -nowarning
    file copy -force $elementsUpdateFile.tmp $elementsUpdateFile
    file delete $elementsUpdateFile.tmp $tmpRoot.types
    
    #--- Create files with lists of bpms for later filtering...
    if $useTunes { set tuneList [list nux nuy] } else { set tuneList "" }
    set bpmList     [list $xbpmList $ybpmList $dbpmList $tuneList]
    set fileExtList [list xbpm ybpm dbpm nu]
    set planeList   [list X Y D Q]
    if {[string compare $calcMode beta] == 0} {
	set weightList1 $weightList
    } else {
	set weightList1 [list [lindex $weightList 0] [lindex $weightList 0] [lindex $weightList 1]]
    }
    foreach bpmList $bpmList fileExt $fileExtList plane $planeList weight $weightList1 {
	if [llength $bpmList] {
	    set dataOption "-data=[join $bpmList ,]"
	} else {
	    set dataOption "-data"
	}
	if [catch {exec sddsmakedataset -pipe=out -col=TagName,type=string $dataOption \
		       | sddsprocess -nowarning -pipe=in $measuredBpms.$fileExt \
		       -print=para,Plane,$plane -def=para,Weight,$weight} result] {
	    return -code error "Error making measuredBpm file for plane $plane: $result"
	}
	if [catch {MakeElementNameFromTagName -input $measuredBpms.$fileExt} result] {
	    return -code error "MakeElementNameFromTagName: $result"
	}
    }
    if {[string compare $calcMode coupling] == 0} {
	set correctorFile $workDir/correctors.param
	if $quadCorrectors {
	    set hSteeringString HSTEERING
	    set vSteeringString VSTEERING
	} else {
	    set hSteeringString STEERING
	    set vSteeringString STEERING
	}
	exec sddsmakedataset -pipe=out -col=TagName,type=string -data=[join $hcorrList ,] \
	    | sddsprocess -pipe=in $tmpRoot.hcorr -print=col,ElementParameter,$hSteeringString -def=col,ParameterValue,1
	exec sddsmakedataset -pipe=out -col=TagName,type=string -data=[join $vcorrList ,] \
	    | sddsprocess -pipe=in $tmpRoot.vcorr -print=col,ElementParameter,$vSteeringString -def=col,ParameterValue,1
	exec sddscombine $tmpRoot.hcorr $tmpRoot.vcorr $correctorFile -merge -overWrite
	if [catch {MakeElementNameFromTagName -input $correctorFile} result] {
	    return -code error "MakeElementNameFromTagName: $result"
	}
	file delete $tmpRoot.hcorr $tmpRoot.vcorr
    }

    #--- Create elegant files...
    #--- RFC CHANGE_T is required if the file contains RF and closed orbit with fixed_length required (it has allow_missing_elements=1)
    foreach arrayName [list measuredFiles achievedFiles desiredFiles] {
	array set filesArray [array get $arrayName]
	if $elemByElemSR {
	    set alterElemString [list "name RFC item CHANGE_T value 0" "name * type CSBEN* item SYNCH_RAD value 1"]
	} else {
	    set alterElemString [list "name RFC item CHANGE_T value 0"]
	}
	switch $calcMode {
	    beta {
		if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $filesArray(ele) \
			       -run_setup [list "lattice $filesArray(lte) use_beamline $filesArray(beamline) parameters $filesArray(outParam) default_order 2 p_central_mev $beamEnergyMeV"] \
			       -load_parameters [list "filename_list \"$filesArray(paramList)\" force_occurence_data 1"] \
			       -alter_elements $alterElemString \
			       -closed_orbit [list "fixed_length 1 output $filesArray(twiss).orb iteration_fraction 0.1 closed_orbit_accuracy 1e-8 closed_orbit_iterations 400"] \
			       -twiss_output [list "filename $filesArray(twiss)"] \
			       -run_control [list "n_steps 1"] \
			       -bunched_beam [list "n_particles_per_bunch 1"] \
			   } result] {
		    return -code error "Fit_GenerateElegantFileFromLTE (ele file $filesArray(ele)): $result"
		}
	    }
	    coupling {
		lappend alterElemString "name * type *KICK item STEERING value 0" "name * type *QUAD item HSTEERING value 0" "name * type *QUAD item VSTEERING value 0"
		if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $filesArray(ele) \
			       -run_setup [list "lattice $filesArray(lte) use_beamline $filesArray(beamline) parameters $filesArray(outParam) default_order 2 p_central_mev $beamEnergyMeV"] \
			       -load_parameters [list "filename_list \"$filesArray(paramList)\" force_occurence_data 1"] \
			       -alter_elements $alterElemString \
			       -load_parameters2 [list "filename $correctorFile force_occurence_data 1"] \
			       -closed_orbit [list "fixed_length 1 output $filesArray(twiss).orb closed_orbit_accuracy 1e-8 closed_orbit_iterations 400 iteration_fraction 0.9"] \
			       -twiss_output [list "filename $filesArray(twiss)"] \
			       -run_control [list "n_steps 1"] \
			       -correct [list "mode orbit method global n_xy_cycles 1 fixed_length 1 closed_orbit_accuracy 1e-7 closed_orbit_iterations 400"] \
			       -correction_matrix_output [list "response2 $filesArray(hyrm) response3 $filesArray(vxrm) coupled 1 fixed_length 1"] \
			       -bunched_beam [list "n_particles_per_bunch 1"] \
			   } result] {
		    return -code error "Fit_GenerateElegantFileFromLTE (ele file $filesArray(ele)): $result"
		}
	    }
	}
    }
}

#-------------------------------------------------------------------------------------------------------------------------
proc AddMeasurementNoise {args} {
    global betaNoise dispNoise rmNoise
    APSParseArguments {calcMode vectorFile}
    #--- Beta noise is in relative units
    #--- Dispersion noise is in meters
    #--- RM noise is in meters, but RM is in m/rad. Assume 5e-5 rad orbit measurement kick
    set rmKick 5e-5
    set glimit 2
    switch $calcMode {
	beta {
	    #--- First 2 pages are beta functions, 3rd page is dispersion
	    exec sddsprocess $vectorFile -pipe=out -def=param,PageIndex,i_page \
		"-redef=col,Measurement,PageIndex 3 < ? $glimit grndl $betaNoise * 1 + Measurement * : Measurement \$" \
		"-redef=col,Measurement,PageIndex 3 == ? $glimit grndl $dispNoise * Measurement + : Measurement \$" \
		| sddsconvert -pipe=in $vectorFile.tmp -del=para,PageIndex
	}
	coupling {
	    #--- First 2 pages are response, 3rd page is dispersion
	    exec sddsprocess $vectorFile -pipe=out -def=param,PageIndex,i_page \
		"-redef=col,Measurement,PageIndex 3 < ? $glimit grndl $rmNoise * $rmKick / Measurement + : Measurement \$" \
		"-redef=col,Measurement,PageIndex 3 == ? $glimit grndl $dispNoise * Measurement + : Measurement \$" \
		| sddsconvert -pipe=in $vectorFile.tmp -del=para,PageIndex
	}
    }
    file copy -force $vectorFile.tmp $vectorFile
    file delete $vectorFile.tmp
}
#-------------------------------------------------------------------------------------------------------------------------

proc BuildMeasurementVector {args} {
    set dispFile ""
    set addMeasurementNoise 0
    APSParseArguments {calcMode eleFile twissFile dispFile hyFile vxFile addMeasurementNoise}
    switch $calcMode {
	beta {
	    if [catch {BuildTwissVector -eleFile $eleFile -twissFile $twissFile -dispFile $dispFile} outputFile] {
		return -code error "BuildTwissVector: $outputFile"
	    }
	}
	coupling {
	    if [catch {BuildRMVector -eleFile $eleFile -hyFile $hyFile -vxFile $vxFile -twissFile $twissFile} outputFile] {
		return -code error "BuildRMVector: $outputFile"
	    }
	}
    }
    if $addMeasurementNoise {
	if [catch {AddMeasurementNoise -calcMode $calcMode -vectorFile $outputFile} result] {
	    return -code error "AddMeasurementNoise: $result"
	}
    }
    return $outputFile
}

#-------------------------------------------------------------------------------------------------------------------------

proc BuildRMVector {args} {
    # h and v stand for corrector plane, and x and y stand for bpm plane.
    APSParseArguments {eleFile hyFile vxFile twissFile}
    global tmpDir measuredBpms acceleratorCode weightList dispersionOnly

    set deleteFiles ""
    set tmpRoot $tmpDir/[APSTmpString]-rmVector
    set outputFile $tmpRoot.out
    set elegantOutFile $tmpRoot.ele.log
    file delete $elegantOutFile
    lappend deleteFiles $elegantOutFile
    if [catch {exec $acceleratorCode $eleFile > $elegantOutFile} result] {
	return -code error "elegant (ele file: $eleFile; output file: $elegantOutFile): $result"
    }
    if [catch {AddTagNameColumn -filename $twissFile} result] {return -code error "AddTagNameColumn: $result"}

    #--- Making one-column file from the RM:
    set fileList ""
    foreach filename [list $hyFile $vxFile] ext [list hyrm vxrm] bpmExt [list ybpm xbpm] {
	if [catch {AddHashtags2rmFile -rmFile $filename} result] {return -code error "AddHashtags2rmFile: $result"}
	if [catch {exec sddsselect $filename $measuredBpms.$bpmExt -pipe=out -match=TagName \
		       | sddsconvert -pipe -del=col,s \
		       | sddsmatrix2column -pipe -data=Measurement -rowname=TagName -major=column \
		       | sddsxref -pipe $measuredBpms.$bpmExt -leave=* -transfer=para,Plane \
		       | sddsprocess -pipe=in $tmpRoot.$ext -redef=para,Weight,[lindex $weightList 0]} result] {
	    return -code error "Error making one-column file (rm file: $filename; ext: $ext): $result"
	}
	if $dispersionOnly {exec sddsprocess $tmpRoot.$ext -nowarning -match=col,Rootname=DUMMY; file delete $tmpRoot.${ext}~}
	lappend fileList $tmpRoot.$ext
	lappend deleteFiles $tmpRoot.$ext
    }
    #--- Add dispersion from twiss file 
    if [catch {exec sddsconvert $twissFile -pipe=out -retain=col,TagName,etay \
		   | sddsselect -pipe $measuredBpms.dbpm -match=TagName \
		   | sddsprocess -pipe -redef=para,Weight,[lindex $weightList 1] \
		   | sddsconvert -pipe=in $tmpRoot.etay -rename=col,etay=Measurement -rename=col,TagName=Rootname} result] {
	return -code error "Error adding dispersion (twissFile: $twissFile): $result"
    }
    lappend fileList $tmpRoot.etay
    lappend deleteFiles $tmpRoot.etay
    eval exec sddscombine $fileList $outputFile
    eval file delete $deleteFiles
    return $outputFile
}

#-------------------------------------------------------------------------------------------------------------------------
proc AddTagNameColumn {args} {
    APSParseArguments {filename}
    set columnList [exec sddsquery $filename -col]
    if {[lsearch -exact $columnList BPMName] != -1 } {
	#--- BPMName already has #occurence, we only add #1 if there is no one
	exec sddsprocess $filename $filename.tagname "-edit=col,TagName,BPMName,e i@#1@ a 2S?@#@ 5d"
    } else {
	exec sddsprocess $filename $filename.tagname "-print=col,TagName,%s#%1d,ElementName,ElementOccurence"
    }
    file copy -force $filename.tagname $filename
    file delete $filename.tagname
}
#-------------------------------------------------------------------------------------------------------------------------
proc BuildTwissVector { args } {

    global tmpDir measuredBpms acceleratorCode
    set eleFile ""
    set dispFile ""
    APSStrictParseArguments {eleFile twissFile dispFile}
    
    set tmpRoot $tmpDir/[APSTmpString].twissVector
    set outputFile $tmpRoot.out
    set elegantOutFile $tmpRoot.ele.log

    file delete $elegantOutFile
    if [catch {exec $acceleratorCode $eleFile > $elegantOutFile} result] {
	return -code error "elegant (ele file: $eleFile; output file: $elegantOutFile): $result"
    } else {
	file delete $elegantOutFile
    }
    if [catch {AddTagNameColumn -filename $twissFile} result] {
	return -code error "AddTagNameColumn: $result"
    }

    #--- Substituting etax column in twissFile with the measured dispersion.
    if [string length $dispFile] {
	exec sddsxref -nowarning $twissFile -pipe=out $dispFile -take=xDispersion -match=TagName=Rootname \
	    | sddsselect -pipe $measuredBpms.dbpm -match=TagName \
	    | sddsprocess -pipe "-redef=col,etax,xDispersion 1000.0 /,units=m" \
	    | sddsconvert -pipe=in $tmpRoot.twiss -del=col,xDispersion
	file copy -force $tmpRoot.twiss $twissFile
	file delete $tmpRoot.twiss
    }

    set fileList ""
    foreach fileExt [list xbpm ybpm dbpm] column [list betax betay etax] {
	exec sddsxref $measuredBpms.$fileExt $twissFile -pipe=out -take=$column -match=TagName -nowarning \
	    | sddsconvert -pipe=in $tmpRoot.$fileExt -rename=col,$column=Measurement -rename=col,TagName=Rootname 
	lappend fileList $tmpRoot.$fileExt
    }

    exec sddscollapse $twissFile -pipe=out \
	| sddsconvert -pipe -retain=col,nu? \
	| sddscollect -pipe=in $tmpRoot.tmp "-collect=match=nu*,col=Measurement,edit="
    set tuneList [exec sdds2stream $tmpRoot.tmp -col=Measurement]
    if {[string first "NaN" [lindex $tuneList 0]] > -1 || [string first "NaN" [lindex $tuneList ]] > -1} {
	file delete $tmpRoot.xbpm $tmpRoot.ybpm $tmpRoot.dbpm $tmpRoot.tmp
	return -code error "Error: Tunes are unstable!"
    }

    exec sddsxref $measuredBpms.nu $tmpRoot.tmp -pipe=out -take=Measurement -match=TagName=Rootname -nowarning \
	| sddsconvert -pipe=in $tmpRoot.nu -rename=col,TagName=Rootname
    lappend fileList $tmpRoot.nu
    
    eval exec sddscombine $fileList $outputFile -overWrite
    eval file delete $fileList $tmpRoot.tmp
    return $outputFile
}

#-------------------------------------------------------------------------------------------------------------------------

proc CalculateVectorDifference {args} {
    APSStrictParseArguments {vectorFile1 vectorFile2 outputFile useRelative nameColumn dataColumn}
    #--- vectorFile2 is current, vectorFile1 is desired
    if $useRelative {
	#--- Dispersion is not relative because it can be zero:
	if [catch {exec sddsxref -nowarning $vectorFile2 $vectorFile1 -pipe=out -take=$dataColumn -rename=col,$dataColumn=${nameColumn}Value \
		       -match=$nameColumn \
		       | sddsconvert -pipe -del=col,Element* \
		       | sddsprocess -pipe "-redef=col,${nameColumn}DivColumn,Plane \"D\" streq ? 1 : ${nameColumn}Value \$" \
		       "-redef=col,${nameColumn}DivColumn,Plane \"Q\" streq ? ${nameColumn}DivColumn : ${nameColumn}DivColumn \$" \
		       | sddsprocess -nowarning -pipe \
		       "-def=col,%sDiff,%s ${nameColumn}Value - ${nameColumn}DivColumn / Weight *,select=*,exclude=${nameColumn}*" \
		       | sddsconvert -pipe=in $outputFile -retain=col,*Diff,$nameColumn} result] {
	    return -code error "Error calculating difference (vectorFile1 $vectorFile1, vectorFile2 $vectorFile2, dataColumn $dataColumn, \                                                   nameColumn $nameColumn): $result"
	}
    } else {
	if [catch {exec sddsxref -nowarning $vectorFile2 $vectorFile1 -pipe=out -take=$dataColumn -rename=col,$dataColumn=${nameColumn}Value \
		       -match=$nameColumn \
		       | sddsconvert -pipe -del=col,Element* \
		       | sddsprocess -nowarning -pipe "-def=col,%sDiff,%s ${nameColumn}Value - Weight *,select=*,exclude=${nameColumn}*" \
		       | sddsconvert -pipe=in $outputFile -retain=col,*Diff,$nameColumn} result] {
	    return -code error "Error calculating difference: $result"
	}
    }
    #--- Calculate RMS if only one Diff column (which means we are not processing response matrices here):
    if [catch {exec sddsconvert $outputFile -pipe=out -retain=col,*Diff | sddsquery -pipe -col -readAll | wc} result] {
	return -code error "Error reading number of columns ($outputFile): $result"
    }
    set diffColumnNumber [lindex $result 0]
    if {$diffColumnNumber == 1} {
	exec sddsprocess $outputFile $outputFile.1 -process=${dataColumn}Diff,rms,Rms -nowarning
	file copy -force $outputFile.1 $outputFile
	file delete $outputFile.1
    }
}

#-------------------------------------------------------------------------------------------------------------------------

proc rmsDisplayOutput {args} {

    set verbose 1
    APSStrictParseArguments {calcMode iteration diffFile twissFile weightList dataColumn verbose}
    
    if [catch {exec sddsprocess $diffFile -pipe=out "-redef=col,$dataColumn,$dataColumn Weight /" -nowarning \
		   | sddscombine -pipe -merge \
		   | sddsprocess -pipe -process=$dataColumn,standardDeviation,rmsTotal \
		   | sdds2stream -pipe=in -para=rmsTotal} rmsTotal] {
	return -code error "Error calculating total rms (diffFile $diffFile): $rmsTotal"
    }
    set returnList $rmsTotal

    set rowList [join [exec sdds2stream $diffFile -rows=bare] ]
    set rmsList [join [exec sdds2stream $diffFile -parameter=Rms] ]
    if {[string compare $calcMode beta] == 0} {
	set calculatedTuneX [exec sdds2stream -para=nux $twissFile]
	set calculatedTuneY [exec sdds2stream -para=nuy $twissFile]
	if {[string first "NaN" $calculatedTuneX] > -1 || [string first "NaN" $calculatedTuneY] > -1 \
		||[string first "nan" $calculatedTuneX] > -1 || [string first "nan" $calculatedTuneY] > -1} {
	    return -code error "Error: Tunes are unstable!"
	}
	set TUNEX [format "%10.4f" $calculatedTuneX]
	set TUNEY [format "%10.4f" $calculatedTuneY]
    }
    switch $calcMode {
	beta {
	    foreach rms [list rmsX rmsY rmsD rmsN] index [list 0 1 2 3] {
		if [lindex $rowList $index] {
		    set scaledRms [expr [lindex $rmsList $index] / [lindex $weightList $index]]
		    set $rms [format "%7.4f" $scaledRms]
		    lappend returnList $scaledRms
		} else { 
		    set $rms NaN 
		}
	    }
	    set outputLine "Iteration $iteration. \nQx = $TUNEX,   Qy = $TUNEY \n"
	    append outputLine "RMS total: [format "%7.4f" $rmsTotal]; --> "
	    append outputLine "x: $rmsX; y: $rmsY; disp: $rmsD; nu: $rmsN"
	}
	coupling {
	    foreach rms [list rmsX rmsY rmsD] index [list 0 1 2] \
		weight [list [lindex $weightList 0] [lindex $weightList 0] [lindex $weightList 1]] {
		    if [lindex $rowList $index] {
			set scaledRms [expr [lindex $rmsList $index] / $weight]
			set $rms [format "%7.4f" $scaledRms]
			lappend returnList $scaledRms
		    } else { 
			set $rms NaN 
		    }
		}
	    set outputLine "Iteration $iteration.\n"
	    append outputLine "RMS total: [format "%7.4f" $rmsTotal]; --> "
	    append outputLine "x: $rmsX; y: $rmsY; disp: $rmsD"
	    
	}
    }
    if $verbose {OutputStatusMessage $outputLine}
    return $returnList
}


#-------------------------------------------------------------------------------------------------------------------------

proc CalculateInverse { args } {

    global tmpDir measuredBpms
    APSStrictParseArguments {matrixFile matrixFileFiltered inverseMatrixFile SVratio quadList diffFile useRelative}

    #--- Prepare matrix for inversion (initial matrixFile is not (F-F0)/delta but just F).
    if ![llength $quadList] {return -code error "Error: quadList is zero length."}
    set takeColumnsOption [join $quadList ","]
    if [catch {exec sddsxref $diffFile $matrixFile -pipe=out -take=$takeColumnsOption -match=Rootname -nowarning \
		   | sddsconvert -pipe=in $matrixFileFiltered -del=col,MeasurementDiff} result] {
	return -code error "Error filtering out non-used bpms (1) (diffFile: $diffFile; matrixFile $matrixFile): $result"
    }
    if [catch {exec sddsxref $diffFile $matrixFile -pipe=out -take=OriginalValue -match=Rootname -nowarning \
		   | sddsconvert -pipe=in $matrixFileFiltered.orig -del=col,MeasurementDiff} result] {
	return -code error "Error filtering out non-used bpms (2) (diffFile: $diffFile; matrixFile $matrixFile): $result"
    }
    lappend deleteFiles $matrixFileFiltered.orig
    #--- Calculate F-F0 (also pages are multiplied by Weight parameter)
    if [catch {CalculateVectorDifference -vectorFile1 $matrixFileFiltered.orig -vectorFile2 $matrixFileFiltered \
		   -outputFile $matrixFileFiltered.diff -useRelative $useRelative -nameColumn Rootname -dataColumn OriginalValue} result] {
	return -code error "CalculateVectorDifference: $result"
    }
    lappend deleteFiles $matrixFileFiltered.diff

    #--- To avoid problem with the column naming in pseudoinverse, reprint the Rootname with unique names for each page
    if {$SVratio < 1} {set svratioOption "-minimumSingularValueRatio=$SVratio"} else {set svratioOption "-largestSingularValues=$SVratio"}
    set delta [lindex [exec sdds2stream $matrixFile -para=QuadDelta] 0]
    if [catch {exec sddsprocess $matrixFileFiltered.diff -pipe=out -nowarning "-reprint=col,Rootname,%s%s,Rootname,Plane" \
		   "-redef=col,%s,%s $delta / -1 *,select=*Diff" \
		   | sddscombine -pipe -merge \
		   | sddspseudoinverse -pipe=in $inverseMatrixFile $svratioOption \
		   -sFile=$matrixFileFiltered.SV -economy} result] {
	return -code error "Error inverting matrix: $result"
    }
    eval file delete $deleteFiles
}

#-------------------------------------------------------------------------------------------------------------------------

proc CalculateFilesForPlotting {args} {
    global measuredBpms beamEnergyMeV
    APSParseArguments {calcMode lteMeasured beamlineMeasured paramFileListMeasured lteDesired beamlineDesired paramFileListDesired \
			   paramFileListAchieved weightList doneFile workDir twiMeasured twiDesired twiAchieved \
			   xbpmList ybpmList dbpmList useTunes useRelative}
    
    set deleteFiles ""
    OutputStatusMessage "Calculations started..."

    #--- Create files with lists of bpms for later filtering...
    set measuredBpms $workDir/measuredBpms.sdds
    if $useTunes { set tuneList [list nux nuy] } else { set tuneList "" }
    set bpmList     [list $xbpmList $ybpmList $dbpmList $tuneList]
    set fileExtList [list xbpm ybpm dbpm nu]
    set planeList   [list X Y D Q]
    foreach bpmList $bpmList fileExt $fileExtList plane $planeList weight $weightList {
	if [llength $bpmList] {set dataOption "-data=[join $bpmList ,]"} else {set dataOption "-data"}
	if [catch {exec sddsmakedataset -pipe=out -col=ElementName,type=string $dataOption \
		       | sddsprocess -nowarning -pipe=in $measuredBpms.$fileExt \
		       -print=para,Plane,$plane -def=para,Weight,$weight} result] {
	    return -code error "Error making measuredBpm file for plane $plane: $result"
	}
	lappend deleteFiles $measuredBpms.$fileExt
    }

    set lteFileList [list $lteMeasured $lteDesired $lteMeasured]
    set beamlineList [list $beamlineMeasured $beamlineDesired $beamlineMeasured]
    set paramFileListList [list $paramFileListMeasured $paramFileListDesired $paramFileListAchieved]
    set twiFileList [list $twiMeasured $twiDesired $twiAchieved]
    set rootList [list measured desired achieved]
    foreach lteFile $lteFileList beamline $beamlineList paramFileList $paramFileListList \
	twissFile $twiFileList root $rootList {
	    #--- Create elegant files...
	    set eleFile $workDir/$root.ele
	    set vectorName ${root}Vector
	    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
			   -run_setup [list "lattice $lteFile use_beamline $beamline default_order 2 p_central_mev $beamEnergyMeV"] \
			   -load_parameters [list "filename_list \"$paramFileList\" force_occurence_data 1"] \
			   -twiss_output [list "filename $twissFile"] \
			   -run_control [list "n_steps 1"] \
			   -bunched_beam [list "n_particles_per_bunch 1"] \
		       } result] {
		return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleFile): $result"
	    }
	    #------ Calculating beta functions:
	    if [catch {BuildMeasurementVector -calcMode $calcMode -eleFile $eleFile -twissFile $twissFile} $vectorName] {
		return -code error "BuildTwissVector: [set $vectorName]"
	    }
	}

    #------ Calculating difference:
    set measDifference $workDir/diff.twi
    lappend deleteFiles $measDifference
    if [catch {CalculateVectorDifference -vectorFile1 $desiredVector -vectorFile2 $achievedVector \
		   -outputFile $measDifference -useRelative $useRelative -nameColumn Rootname \
		   -dataColumn Measurement} result ] {
	return -code error "CalculateVectorDifference: $result"
    }
    if [catch {rmsDisplayOutput -calcMode $calcMode -iteration -1 -diffFile $measDifference -twissFile $workDir/achieved.twi \
		   -weightList $weightList -dataColumn MeasurementDiff} result] {
	return -code error "rmsDisplayOutput: $result"
    }
    eval file delete $deleteFiles
    exec touch $doneFile
}

#-----------------------------------------------------------------------------------------------------------------------

proc CalculateEmittances {args} {
    global acceleratorCode desiredFiles beamEnergyMeV
    APSParseArguments {calcMode workDir lteMeasured beamlineMeasured paramFileListMeasured lteDesired \
			   beamlineDesired paramFileListDesired paramFileListAchieved emitArrayList}
    set deleteFiles ""
    OutputStatusMessage "Calculations started (make sure the .lte file contains kick elements and a single rf cavity named RF)..."

    array set emitArray $emitArrayList
    #------ Run elegant to calculate beam losses and rf frequency:
    if [catch {GetRFParameters -lteFile $lteDesired -beamline $beamlineDesired -paramFileList $paramFileListDesired \
		   -rfVoltage $emitArray(rfVoltage) -rfH $emitArray(rfH) -workDir $workDir} result] {
	return -code error "GetRFParameters: $result"
    }
    set rfVolt $emitArray(rfVoltage)
    set rfPhi [lindex $result 0]
    set rfFreq [lindex $result 1]

    #------ Make elegant template for moments calculations:
    OutputStatusMessage "Calculating beam moments..."
    set lteFileList [list $lteMeasured $lteDesired $lteMeasured]
    set beamlineList [list $beamlineMeasured $beamlineDesired $beamlineMeasured]
    set paramFileListList [list $paramFileListMeasured $paramFileListDesired $paramFileListAchieved]
    set rootList [list measured desired achieved]
    set eleFile $workDir/moments.ele
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -divide_elements [list "name S35B:M divisions 11"] \
		   -run_setup [list "lattice <LTEFILE> use_beamline <BEAMLINE> default_order 2 parameters [file rootname $eleFile].param semaphore_file <DONEFILE> p_central_mev $beamEnergyMeV"] \
		   -load_parameters [list "filename_list \"<PARAMFILELIST>\" allow_missing_files 1 force_occurence_data 1"] \
		   -alter_elements [list "name RF item VOLT value $rfVolt" \
					"name RF item PHASE value $rfPhi" \
					"name RF item FREQ value $rfFreq" \
					"name * type KQUAD item N_KICKS value 50" \
					"name * type KSEXT item N_KICKS value 3" \
					"name * type CSBEND item N_KICKS value 100" \
					"name * type CSBEND item EDGE_ORDER value 1" \
					"name * type CSBEND item INTEGRATION_ORDER value 4" \
					"name * type KQUAD item ISR value 1" \
					"name * type KSEXT item ISR value 1" \
					"name * type CSBEND item ISR value 1" \
					"name * type KQUAD item SYNCH_RAD value 1" \
					"name * type KSEXT item SYNCH_RAD value 1" \
					"name * type CSBEND item SYNCH_RAD value 1"] \
		   -run_control [list "n_steps 1"] \
		   -closed_orbit [list "closed_orbit_accuracy 1e-7 closed_orbit_iterations 500 iteration_fraction 0.25"] \
		   -twiss_output [list "output_at_each_step 0 radiation_integrals 1" "filename <TWISSFILE> radiation_integrals 1"] \
		   -moments_output [list "filename <MOMENTSFILE> radiation 1 n_slices 100 matched 1 verbosity 1"] \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
		   -track [list "center_on_orbit 1"] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleFile): $result"
    }
    #------ Run elegant calculations:
    foreach lteFile $lteFileList beamline $beamlineList paramFileList $paramFileListList root $rootList {
	set twissFile $workDir/moments.$root.twi
	set momentsFile $workDir/moments.$root.mom
	set doneFile $workDir/moments.$root.done
	catch {file delete $eleFile.log}
	if {[string length $paramFileList] == 0} {set paramFileList dummy}
	if [catch {exec $acceleratorCode $eleFile \
		       -macro=LTEFILE=$lteFile,BEAMLINE=$beamline,PARAMFILELIST="$paramFileList",TWISSFILE=$twissFile,MOMENTSFILE=$momentsFile,DONEFILE=$doneFile \
		       > $eleFile.$root.log} result] {
	    return -code error "Error running elegant (see $eleFile.log\ncommand: \"$acceleratorCode $eleFile -macro=LTEFILE=$lteFile,BEAMLINE=$beamline,PARAMFILELIST=\"$paramFileList\",TWISSFILE=$twissFile,MOMENTSFILE=$momentsFile\"): $result"
	}
	set e$root [exec sdds2stream $momentsFile -para=e1,e2]
	set s35ph$root [exec sddsprocess $momentsFile -pipe=out -match=col,ElementName=S35B:M -clip=1,0,inv \
	    | sdds2stream -pipe=in -col=s3]
    }
    #------ Check is there is rf voltage in the lattice:
    if [catch {exec sddsprocess [file rootname $eleFile].param -pipe=out -match=col,ElementType=RFCA* \
		   -match=col,ElementParameter=VOLT -process=ParameterValue,sum,RFVOLTAGE \
		   | sdds2stream -pipe=in -para=RFVOLTAGE} rfVolt] {
	return -code error "Error calculating total rf voltage: $result"
    } else {
	if {$rfVolt < 6e6 || $rfVolt > 12e6} {
	    OutputStatusMessage "----> Warning: total RF voltage is unusual: $rfVolt"
	}
    }
    OutputStatusMessage "Results of moments calculations:"
    OutputStatusMessage "Ideal emittances:     ex = [format %10.2e [lindex $edesired 0]]; ey = [format %10.2e [lindex $edesired 1]]; S35B:M.s3 = [format %10.2e $s35phdesired]."
    OutputStatusMessage "Measured  emittances: ex = [format %10.2e [lindex $emeasured 0]]; ey = [format %10.2e [lindex $emeasured 1]]; S35B:M.s3 = [format %10.2e $s35phmeasured]."
    OutputStatusMessage "Corrected emittances: ex = [format %10.2e [lindex $eachieved 0]]; ey = [format %10.2e [lindex $eachieved 1]]; S35B:M.s3 = [format %10.2e $s35phachieved]."
}

#-----------------------------------------------------------------------------------------------------------------------

proc GetRFParameters {args} {
    global acceleratorCode beamEnergyMeV
    set outputParamFile ""
    APSParseArguments {workDir lteFile beamline paramFileList rfVoltage rfH outputParamFile}
    OutputStatusMessage "Calculating rf parameters..."
    set eleFile $workDir/rfparam.ele
    set twissFile $workDir/rfparam.twi
    if {[string length $paramFileList] == 0} {set paramFileList dummy}
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -run_setup [list "lattice $lteFile use_beamline $beamline default_order 2 \
                                        parameters [file rootname $eleFile].param centroid [file rootname $eleFile].cen p_central_mev $beamEnergyMeV"] \
		   -load_parameters [list "filename_list \"$paramFileList\" allow_missing_files 1 force_occurence_data 1"] \
		   -alter_elements [list "name * type KQUAD item N_KICKS value 50" \
					"name * type KSEXT item N_KICKS value 3" \
					"name * type CSBEND item N_KICKS value 100" \
					"name * type CSBEND item EDGE_ORDER value 1" \
					"name * type CSBEND item INTEGRATION_ORDER value 4"] \
		   -run_control [list "n_steps 1"] \
		   -closed_orbit [list "closed_orbit_accuracy 1e-7 closed_orbit_iterations 500 iteration_fraction 0.25"] \
		   -twiss_output [list "output_at_each_step 0 radiation_integrals 1" "filename $twissFile radiation_integrals 1"] \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
		   -track [list "center_on_orbit 1"] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleFile): $result"
    }
    catch {file delete $eleFile.log}
    if [catch {exec $acceleratorCode $eleFile > $eleFile.log} result] {
	return -code error "Error running elegant (see $eleFile.log\ncommand: \"$acceleratorCode $eleFile\"): $result"
    }
    set U0 [expr [exec sdds2stream $twissFile -par=U0] * 1e6]
    set phi [exec rpnl "180 $U0 $rfVoltage / dasin -"]
    set C [exec sddsprocess [file rootname $eleFile].cen -pipe=out -process=Cs,last,C \
	       | sdds2stream -pipe=in -para=C]
    set freq [exec rpnl "c_mks $C / $rfH *"]
    if [string length $outputParamFile] {
	exec sddsmakedataset $outputParamFile \
	    -col=ElementName,type=string -data=RF,RF,RF \
	    -col=ElementParameter,type=string -data=VOLT,PHASE,FREQ \
	    -col=ParameterValue,type=double -data=$rfVoltage,$phi,$freq
    }
    return [list $phi $freq]
}

#-----------------------------------------------------------------------------------------------------------------------

proc MinimizeEmittance {args} {
    global acceleratorCode beamEnergyMeV
    APSParseArguments {workDir quadList lteMeasured beamlineMeasured paramFileListMeasured emitArrayList optimTerm rootname}
    set deleteFiles ""
    OutputStatusMessage "Calculations started (make sure the .lte file contains kick elements and a single rf cavity named RF)..."
    array set emitArray $emitArrayList
    set rfParamFile $workDir/$rootname.rfparam
    #------ Run elegant to calculate beam losses and rf frequency:
    if [catch {GetRFParameters -lteFile $lteMeasured -beamline $beamlineMeasured -paramFileList $paramFileListMeasured \
		   -rfVoltage $emitArray(rfVoltage) -rfH $emitArray(rfH) -workDir $workDir \
		   -outputParamFile $rfParamFile} result] {
	return -code error "GetRFParameters: $result"
    }
    set paramFileList [concat $paramFileListMeasured $rfParamFile]

    #------ Insert element M1 into .lte file after every *END element:
    set lteFile $workDir/$rootname.lte
    if [catch {open $lteMeasured r} fid1] {return -code error "Error opening file: $fid1"}
    if [catch {open $lteFile w} fid2] {return -code error "Error opening file: $fid2"}
    for {set i 1} {$i < 41} {incr i} {lappend nameList S${i}END}
    puts $fid2 "M1: MARK, FIT_POINT=1"
    while {[gets $fid1 line] >= 0} {
	foreach name $nameList {
	    if {[string first $name $line] != -1 && [string first MARK $line] == -1} {
		set index [expr [string first $name $line] + [string length $name]]
		set string1 [string range $line 0 [expr $index - 1]]
		set string2 [string range $line $index 999]
		set line ""
		append line $string1 ,M1 $string2
		break
	    }
	}
	puts $fid2 $line
    }
    close $fid1
    close $fid2

    #------ Make elegant file for emittance minimization:
    set eleFile $workDir/$rootname.ele
    set stepSize 0.001
    set lowerLimit -0.12
    set upperLimit 0.12
    foreach skewQuad $quadList {
	lappend optimizationVariableList "name $skewQuad item K1 step_size $stepSize lower_limit $lowerLimit upper_limit $upperLimit"
    }
    #------ Here optimization term is made of form "\"term\"", that's why I need so many \\\.
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -run_setup [list "lattice $lteFile use_beamline $beamlineMeasured default_order 2 parameters [file rootname $eleFile].param final [file rootname $eleFile].fin p_central_mev $beamEnergyMeV"] \
		   -load_parameters [list "filename_list \"$paramFileList\" allow_missing_files 1 force_occurence_data 1"] \
		   -alter_elements [list "name * type KQUAD item N_KICKS value 10" \
					"name * type KSEXT item N_KICKS value 3" \
					"name * type CSBEND item N_KICKS value 10" \
					"name * type CSBEND item EDGE_ORDER value 1" \
					"name * type CSBEND item INTEGRATION_ORDER value 4" \
					"name * type KQUAD item ISR value 1" \
					"name * type KSEXT item ISR value 1" \
					"name * type CSBEND item ISR value 1" \
					"name * type KQUAD item SYNCH_RAD value 1" \
					"name * type KSEXT item SYNCH_RAD value 1" \
					"name * type CSBEND item SYNCH_RAD value 1"] \
		   -run_control [list "n_steps 1"] \
		   -moments_output [list "filename [file rootname $eleFile].mom radiation 1 n_slices 10 matched 1 verbosity 0"] \
		   -optimization_setup [list "mode minimize method simplex target 1e-14 tolerance 1e-14 n_passes 5 n_restarts 1 n_evaluations 200 log_file /dev/tty verbose 1 output_sparsing_factor 10"] \
		   -optimization_term [list "term \"\\\"$optimTerm\\\"\" field_string @ field_initial_value 1 field_final_value 39 field_interval 1"] \
		   -optimization_variable $optimizationVariableList \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
		   -optimize [list ""] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleFile): $result"
    }
    
    #------ Running optimization:
    exec $acceleratorCode $eleFile >&@ stdout
}

#-----------------------------------------------------------------------------------------------------------------------

proc CalculateTwissResponseMatrix { args } {

    global workDir desiredFiles measuredFiles eleMatrixTemplate correctorFile acceleratorCode elementsUpdateFile
    global qsubRespProcCommand queueSystemName beamEnergyMeV
    APSStrictParseArguments {calcMode matrixFile quadList useRelative abortFile splitTasks useQsub doPerturbedMatrix}

    #------ For submitting several jobs without queue:
    if $useQsub {set qsubCommand "exec \$scriptFile &"} else {set qsubCommand "exec \$scriptFile"}
    set queueSystemName PS
    set qsubRespProcCommand "lindex \$qsubResponse 0"
    set remotePath ""

    set waitTime 3600
    set waitInterval 2
    set submissionPause 0
    set continuePrevious 0 
    set usePopupWindow 0

    switch $calcMode {
	beta {
	    if $doPerturbedMatrix {set paramFiles $measuredFiles(paramList)} else {set paramFiles $desiredFiles(paramList)}
	    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleMatrixTemplate \
			   -run_setup [list "lattice $desiredFiles(lte) use_beamline $desiredFiles(beamline) default_order 2 p_central_mev $beamEnergyMeV"] \
			   -load_parameters [list "filename_list \"$paramFiles\" force_occurence_data 1" "filename <PARAMFILE> change_defined_values 0"] \
			   -closed_orbit [list "fixed_length 1 iteration_fraction 0.1 closed_orbit_accuracy 1e-6 closed_orbit_iterations 400"] \
			   -twiss_output [list "filename <TWIFILE>"] \
			   -run_control [list "n_steps <NSTEPS>"] \
			   -bunched_beam [list "n_particles_per_bunch 1"] \
		       } result] {
		return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleMatrixTemplate): $result"
	    }
	} 
	coupling {
	    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleMatrixTemplate \
			   -run_setup [list "lattice $desiredFiles(lte) use_beamline $desiredFiles(beamline) p_central_mev $beamEnergyMeV"] \
			   -load_parameters [list "filename_list \"$desiredFiles(paramList)\" force_occurence_data 1" "filename <PARAMFILE> change_defined_values 0"] \
			   -alter_elements [list "name * type *KICK item STEERING value 0" "name * type *QUAD item HSTEERING value 0" "name * type *QUAD item VSTEERING value 0"] \
			   -load_parameters2 [list "filename $correctorFile force_occurence_data 1"] \
			   -closed_orbit [list "fixed_length 0 closed_orbit_accuracy 1e-7"] \
			   -twiss_output [list "filename <TWIFILE>"] \
			   -run_control [list "n_steps <NSTEPS>"] \
			   -correct [list "mode orbit method global n_xy_cycles 0 fixed_length 1 closed_orbit_accuracy 1e-7 closed_orbit_iterations 400"] \
			   -correction_matrix_output [list "response2 <HYRMFILE> response3 <VXRMFILE> coupled 1 fixed_length 1"] \
			   -bunched_beam [list "n_particles_per_bunch 1"] \
		       } result] {
		return -code error "Fit_GenerateElegantFileFromLTE (ele file $eleMatrixTemplate): $result"
	    }
	}
    }

    set varList [concat OriginalValue $quadList]
    #--- Make elementsUpdateFile that has "OriginalValue" in it:
    exec sddsprocess $elementsUpdateFile $elementsUpdateFile.1 -clip=1,0,inv -reprint=col,TagName,OriginalValue
    exec sddscombine $elementsUpdateFile.1 $elementsUpdateFile $elementsUpdateFile.origValue -merge -over
    file delete $elementsUpdateFile.1
    set scriptName calculateTwissCorrection
    set scriptParameters " -mode runDqsTask -calcMode $calcMode -eleTemplate $eleMatrixTemplate -useRelative $useRelative -acceleratorCode $acceleratorCode -elementsUpdateFile $elementsUpdateFile.origValue"

    set rootTaskName twiDeriv
    set tmpDirName $workDir/dqs/twiss
    if ![file exists $workDir/dqs] {exec mkdir $workDir/dqs}
    if ![file exists $workDir/dqs/twiss] {exec mkdir $workDir/dqs/twiss}
    if [llength $quadList] {
	OutputStatusMessage "Doing quadrupoles..."
	set verbose 0
	if [catch {Fit_CalculateResponseMatrixDerivative -varList $varList \
		       -splitTasks $splitTasks -scriptName $scriptName \
		       -scriptParameters $scriptParameters -matrixFile $matrixFile -tmpDirName $tmpDirName \
		       -useQsub $useQsub -qsubCommand $qsubCommand -verbose $verbose \
		       -rootTaskName $rootTaskName -continuePrevious $continuePrevious -usePopupWindow $usePopupWindow \
		       -waitTime $waitTime -waitInterval $waitInterval -submissionPause $submissionPause \
		       -abortFile $abortFile} result] {
	    return -code error "Fit_CalculateResponseMatrixDerivative: $result"
	}
    } else {
	return -code error "Error: quadList has zero length."
    }
    file delete $eleMatrixTemplate

}

#-------------------------------------------------------------------------------------------------------------------------

proc OutputStatusMessage {text} {
    global outputStatusDevice outputStatusFile errorMessage
    set errorMessage $text
    switch -exact $outputStatusDevice {
        statusScreen { APSSetVarAndUpdate status $text }
        stdout { puts stdout $text }
        file {
            if {![info exists outputStatusFile] || ![string length $outputStatusFile]} {
                return -code error "OutputStatusMessage: outputStatusFile variable is not defined."
            }
            if [file exists $outputStatusFile] {
                set fid [open $outputStatusFile a]
                puts $fid $text
                close $fid
            } else {
                set fid [open $outputStatusFile w]
                puts $fid $text
                close $fid
            }
        }
        default { return -code error "OutputStatusMessage: wrong outputStatusDevice: $outputStatusDevice" }
    }
}

#-------------------------------------------------------------------------------------------------------------------------

proc StopExecution { } {
}

#-----------------------------------------------------------------------------------------------------------------

proc ScanSVnumber {args} {

    global outputStatusDevice outputStatusFile

    APSParseArguments {calcMode lteMeasured beamlineMeasured paramFileListMeasured lteDesired beamlineDesired paramFileListDesired \
			   outputFile quadList xbpmList ybpmList dbpmList weightList calculateMatrix useRelative \
			   SVstart SVend SVsteps matrixFile corFraction useDQS splitTasks iterations useTunes \
			   useMeasDisp measDispFile outputStatusDevice abortFile hcorrList vcorrList}

    set outputFileList ""
    set svList ""
    set k1List ""
    set rmsList ""
    set rmsXList ""
    set rmsYList ""
    set rmsDList ""
    set rmsNList ""
    
    set SVstep [expr ($SVend - $SVstart) / ($SVsteps - 1)]
    for {set I 0} {$I < $SVsteps} {incr I} {
	OutputStatusMessage "Step: $I."

	if [file exists $abortFile] {
	    return -code error "Scan interrupted by user."
	}

	set SVnumber [expr int($SVstart + $SVstep * $I)]
	if [catch {TwissCorrection \
		       -calcMode $calcMode \
		       -lteMeasured $lteMeasured \
		       -beamlineMeasured $beamlineMeasured \
		       -paramFileListMeasured $paramFileListMeasured \
		       -lteDesired $lteDesired \
		       -beamlineDesired $beamlineDesired \
		       -paramFileListDesired $paramFileListDesired \
		       -outputFile $outputFile \
		       -quadList $quadList \
		       -xbpmList $xbpmList \
		       -ybpmList $ybpmList \
		       -dbpmList $dbpmList \
		       -hcorrList $hcorrList \
		       -vcorrList $vcorrList \
		       -weightList $weightList \
		       -calculateMatrix 0 \
		       -calculateInverse 1 \
		       -useRelative $useRelative \
		       -SVratio $SVnumber \
		       -matrixFile $matrixFile \
		       -corFraction $corFraction \
		       -useDQS $useDQS \
		       -splitTasks $splitTasks \
		       -iterations $iterations \
		       -useTunes $useTunes \
		       -useInvSolution 0 \
		       -locoResultFile dummy \
		       -useMeasDisp $useMeasDisp \
		       -measDispFile $measDispFile \
		       -outputStatusDevice $outputStatusDevice \
		       -abortFile abortFile} resultList] {
	    OutputStatusMessage "TwissCorrection: $resultList"
	} else {
	    exec sddsprocess $outputFile $outputFile.$I -nowarning \
		-def=para,SV,$SVnumber \
		-process=ParameterValue,rms,K1rms \
		-def=para,TotalRms,[lindex $resultList 0]
	    lappend outputFileList $outputFile.$I
	    lappend svList $SVnumber
	    lappend k1List [exec sdds2stream $outputFile.$I -param=K1rms]
	    lappend rmsList [lindex  $resultList 0]
	    lappend rmsXList [lindex $resultList 1]
	    lappend rmsYList [lindex $resultList 2]
	    lappend rmsDList [lindex $resultList 3]
	    lappend rmsNList [lindex $resultList 4]
	}
    }
    exec sddsmakedataset twissScan.sdds \
	-col=SV,type=double -data=[join $svList ,] \
	-col=rmsK1,type=double -data=[join $k1List ,] \
	-col=rmsTotal,type=double -data=[join $rmsList ,] \
	-col=rmsX,type=double -data=[join $rmsXList ,] \
	-col=rmsY,type=double -data=[join $rmsYList ,] \
	-col=rmsD,type=double -data=[join $rmsDList ,] \
	-col=rmsN,type=double -data=[join $rmsNList ,]
    eval exec sddscombine $outputFileList -pipe=out \
	| sddsprocess -pipe=in $outputFile.scan "\"-print=para,BottomLine,Krms=%8.2e; Meas rms=%8.2e,K1rms,TotalRms\"" \
	"\"-print=para,TopLine,Number of SV = %3.0f,SV\""
    eval file delete $outputFileList

}

#-------------------------------------------------------------------------------------------------------------------------

proc RunDqsTask {args} {

    set verbose 0
    APSParseArguments {calcMode path eleTemplate varList useRelative acceleratorCode verbose elementsUpdateFile}

    cd $path
    set elegantOutFile $path/elegant.out
    set twiFile        $path/aps.twi
    set hyFile         $path/aps.hyrm
    set vxFile         $path/aps.vxrm
    set matrixFile     $path/responseMatrix.sdds
    set paramFile      $path/variables.param
    set delta 0.0001
    set Nruns [llength $varList]

    #------ Run elegant to calculate the multipage response file:
    exec sddsmakedataset -pipe=out -col=TagName,type=string -data=[join $varList ,] \
	| sddsprocess -pipe -print=col,ParameterMode,DIFFERENTIAL -def=col,ParameterValue,$delta \
	| sddsxref -pipe $elementsUpdateFile -take=ElementParameter -match=TagName -nowarning \
	| sddsbreak -pipe=in $paramFile -rowlimit=1
    if [catch {MakeElementNameFromTagName -input $paramFile} result] {
	return -code error "MakeElementNameFromTagName ($paramFile): $result"
    }
 
    catch {file delete $elegantOutFile}
    set macroCommand PARAMFILE=$paramFile,TWIFILE=$twiFile,HYRMFILE=$hyFile,VXRMFILE=$vxFile,NSTEPS=$Nruns
    if $verbose {
	puts stdout "Running elegant..."
	puts stdout "$acceleratorCode $eleTemplate -macro=$macroCommand > $elegantOutFile"
    }
    if [catch {exec $acceleratorCode $eleTemplate -macro=$macroCommand > $elegantOutFile} result] {
	puts stderr "elegant: $result"
	return -code error
    }
    if [catch {AddTagNameColumn -filename $twiFile} result] {
	return -code error "AddTagNameColumn ($twiFile): $result"
    }

    #----- Postprocessing the response file...

    set tmpRoot /tmp/[APSTmpString]-twissMatrix
    if $verbose {puts stdout "Postprocessing..."}

    #--- Make file from columns of twiss file (betax,betay,etax for beta correction or etay for coupling correction)
    if {[string compare $calcMode beta] == 0} {set colList [list betax betay etax]} else {set colList etay}
    if [catch {exec sddsprocess $twiFile -pipe=out -match=col,ElementType=MONI \
		   | sddsconvert -pipe -retain=col,TagName,[join $colList ,] -del=para,* \
		   | sddssortcolumn -pipe -sortList=TagName,[join $colList ,] \
		   | sddsmatrix2column -pipe -rowNameColumn=TagName -dataColumnName=Data -rootnameColumnName=Rootname -majorOrder=column \
		   | sddstranspose -pipe \
		   | sddsxref -pipe $paramFile -take=TagName \
		   | sddscombine -pipe -merge \
		   | sddstranspose -pipe -newcolumn=TagName -oldcolumn=Rootname \
		   | sddsprocess -pipe -del=para,* -def=para,QuadDelta,$delta \
		   "-edit=col,TwissFunctionName,Rootname,s@eta@ 1f K" "-reedit=col,Rootname,%@betax@@ %@betay@@ %@etax@@ %@etay@@" \
		   | sddsbreak -pipe -change=TwissFunctionName \
		   | sddsprocess -pipe=in $tmpRoot.twiss -process=TwissFunctionName,first,ResponseQuantity} result] {
	return -code error "Error processing columns in $twiFile: $result"
    }

    if {[string compare $calcMode beta] == 0} {
	#------ Beta mode calculations:
	#--- Make file from parameters of twiss file (nux,nuy)
	if [catch {exec sddsconvert $twiFile -pipe=out -retain=para,nux,nuy \
		       | sddscollapse -pipe \
		       | sddsbreak -pipe -rowlimit=1 \
		       | sddsxref -pipe $paramFile -take=TagName \
		       | sddscombine -pipe -merge \
		       | sddsconvert -pipe -del=col,PageNumber \
		       | sddstranspose -pipe -oldcolumnnames=Rootname \
		       | sddsconvert -pipe -del=para,* \
		       | sddsprocess -pipe=in $tmpRoot.nu -print=col,TwissFunctionName,nu -print=para,ResponseQuantity,nu} result] {
	    return -code error "Error processing parameters in $twiFile: $result"
	}
	set fileList [list $tmpRoot.twiss $tmpRoot.nu]
    } else {
	#------ Coupling mode calculations:
	foreach ext [list hyrm vxrm] rmFile [list $hyFile $vxFile] {
	    if [catch {AddHashtags2rmFile -rmFile $rmFile} result] {return -code error "AddHashtags2rmFile: $result"}
	    if [catch {exec sddsconvert $rmFile -pipe=out -del=col,s \
			   | sddsmatrix2column -pipe -rowNameColumn=TagName -dataColumnName=Data -rootnameColumnName=Rootname -majorOrder=column \
			   | sddstranspose -pipe \
			   | sddsxref -pipe $paramFile -take=TagName \
			   | sddscombine -pipe -merge \
			   | sddstranspose -pipe -newcolumn=TagName -oldcolumn=Rootname \
			   | sddsprocess -pipe=in $tmpRoot.$ext -del=para,* -def=para,QuadDelta,$delta -print=para,ResponseQuantity,$ext} result] {
		return -code error "Error processing file $rmFile: $result"
	    }
	}
	set fileList [list $tmpRoot.hyrm $tmpRoot.vxrm $tmpRoot.twiss]
    }
    eval exec sddscombine $fileList -pipe=out \
	| sddsconvert -pipe=in $matrixFile -del=col,TwissFunctionName
    eval file delete $fileList
    exec echo END > $path/task.completed
    exit
}

#-----------------------------------------------------------------------------------------------------------------
proc AddHashtags2rmFile {args} {
    global tmpDir
    APSParseArguments {rmFile}
    set tmpRoot $tmpDir/[APSTmpString]-hashtag
    if [catch {AddTagNameColumn -filename $rmFile} result] {
	return -code error "AddTagNameColumn: $result"
    }
    #--- Adding #1 at the end of column names (elegant does not always put #Occurence at the end)
    exec sddsconvert $rmFile -pipe=out "-edit=col,*,e i@#1@ a 2S?@#@ 2d" \
	| sddsconvert -pipe=in $tmpRoot -rename=col,BPMName#1=BPMName -rename=col,TagName#1=TagName -rename=col,s#1=s
    file copy -force $tmpRoot $rmFile
    file delete $tmpRoot
}
#-----------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------
#-----------------------------------------------------------------------------------------------------------------

if {[string compare $mode calcEmit] != 0 && [string compare $mode runDqsTask] != 0} {
    if [catch {ReadListFromFile -calcMode $calcMode -listFile $elementListFile} elementList] {
	puts stderr "ReadListFromFile: $elementList"
	exit
    }
    if {[string compare $calcMode beta] == 0} {
	set quadList [lindex $elementList 0]
	set xbpmList [lindex $elementList 1]
	set ybpmList [lindex $elementList 2]
	set dbpmList [lindex $elementList 3]
	set hcorrList ""
	set vcorrList ""
    } else {
	set quadList [lindex $elementList 0]
	set hcorrList [lindex $elementList 1]
	set vcorrList [lindex $elementList 2]
	set xbpmList [lindex $elementList 3]
	set ybpmList [lindex $elementList 4]
	set dbpmList [lindex $elementList 5]
    }
}
if [string length $betaTargetFile] {
    set adjustVectorValues 1
} else {
    set adjustVectorValues 0
}
switch -exact $mode {
    correct {
	set rmsIterationsFile $workDir/tmp/iterations.rms
	catch {file delete $rmsIterationsFile}
	if [catch {TwissCorrection \
		       -calcMode $calcMode \
		       -lteMeasured $lteMeasured \
		       -beamlineMeasured $beamlineMeasured \
		       -paramFileListMeasured $paramFileListMeasured \
		       -lteDesired $lteDesired \
		       -beamlineDesired $beamlineDesired \
		       -paramFileListDesired $paramFileListDesired \
		       -outputFile $outputFile \
		       -quadList $quadList \
		       -xbpmList $xbpmList \
		       -ybpmList $ybpmList \
		       -dbpmList $dbpmList \
		       -hcorrList $hcorrList \
		       -vcorrList $vcorrList \
		       -weightList $weightList \
		       -calculateMatrix $calculateMatrix \
		       -calculateInverse $calculateInverse \
		       -useRelative $useRelative \
		       -SVratio $SVratio \
		       -matrixFile $matrixFile \
		       -corFraction $corFraction \
		       -useDQS $useDQS \
		       -splitTasks $splitTasks \
		       -iterations $iterations \
		       -useTunes $useTunes \
		       -useInvSolution $useInvSolution \
		       -locoResultFile $locoResultFile \
		       -useMeasDisp $useMeasDisp \
		       -measDispFile $measDispFile \
		       -outputStatusDevice $outputStatusDevice \
		       -adjustVectorValues $adjustVectorValues \
		       -addMeasurementNoise $addMeasurementNoise \
		       -outputSparsity $outputSparsity \
		       -returnToBest 1 \
		       -abortFile $abortFile} result] {
	    OutputStatusMessage "TwissCorrection: $result \nInterrupt iterations."
	    if 0 {
	    OutputStatusMessage "TwissCorrection exited abnormally, will re-run to the best iteration..."
	    if [file exists $rmsIterationsFile] {
		set minIteration [exec sddsprocess $rmsIterationsFile -pipe=out -process=RMS,min,IterMin,position,functionof=Iteration \
				      | sdds2stream -pipe=in -para=IterMin]
		if {$minIteration != 0} {
		    set iterations $minIteration
		    set calculateMatrix 0
		    set calculateInverse 0
		    catch {file delete $rmsIterationsFile}
		    if [catch {TwissCorrection \
				   -calcMode $calcMode \
				   -lteMeasured $lteMeasured \
				   -beamlineMeasured $beamlineMeasured \
				   -paramFileListMeasured $paramFileListMeasured \
				   -lteDesired $lteDesired \
				   -beamlineDesired $beamlineDesired \
				   -paramFileListDesired $paramFileListDesired \
				   -outputFile $outputFile \
				   -quadList $quadList \
				   -xbpmList $xbpmList \
				   -ybpmList $ybpmList \
				   -dbpmList $dbpmList \
				   -hcorrList $hcorrList \
				   -vcorrList $vcorrList \
				   -weightList $weightList \
				   -calculateMatrix $calculateMatrix \
				   -calculateInverse $calculateInverse \
				   -useRelative $useRelative \
				   -SVratio $SVratio \
				   -matrixFile $matrixFile \
				   -corFraction $corFraction \
				   -useDQS $useDQS \
				   -splitTasks $splitTasks \
				   -iterations $iterations \
				   -useTunes $useTunes \
				   -useInvSolution $useInvSolution \
				   -locoResultFile $locoResultFile \
				   -useMeasDisp $useMeasDisp \
				   -measDispFile $measDispFile \
				   -outputStatusDevice $outputStatusDevice \
				   -adjustVectorValues $adjustVectorValues \
				   -addMeasurementNoise $addMeasurementNoise \
				   -outputSparsity $outputSparsity \
				   -returnToBest 0 \
				   -abortFile $abortFile} result] {
			OutputStatusMessage "TwissCorrection: $result \nInterrupt iterations."
			exit 1
		    }
		    set finalOutput "Final results: "
		    foreach rms $result {append finalOutput [format %8.4f $rms]; append finalOutput "  "}
		    OutputStatusMessage $finalOutput
		    OutputStatusMessage "Iterations are done."
		} else {
		    #--- No good non-zero iterations
		    puts stdout "No good iterations were found. Exiting..."
		    exit 1
		}
	    } else {
		#--- No rmsIterationsFile found:
		puts stdout "Error running script."
		exit 1
	    }
	    }
	    exit 1
	} else {
	    set finalOutput "Final results: "
	    foreach rms $result {append finalOutput [format %8.4f $rms]; append finalOutput "  "}
	    OutputStatusMessage $finalOutput
	    OutputStatusMessage "Iterations are done."
	}
    }
    scan {
	if [catch {ScanSVnumber \
		       -calcMode $calcMode \
		       -lteMeasured $lteMeasured \
		       -beamlineMeasured $beamlineMeasured \
		       -paramFileListMeasured $paramFileListMeasured \
		       -lteDesired $lteDesired \
		       -beamlineDesired $beamlineDesired \
		       -paramFileListDesired $paramFileListDesired \
		       -outputFile $outputFile \
		       -quadList $quadList \
		       -xbpmList $xbpmList \
		       -ybpmList $ybpmList \
		       -dbpmList $dbpmList \
		       -hcorrList $hcorrList \
		       -vcorrList $vcorrList \
		       -weightList $weightList \
		       -calculateMatrix $calculateMatrix \
		       -useRelative $useRelative \
		       -SVstart $SVstart \
		       -SVend $SVend \
		       -SVsteps $SVsteps \
		       -matrixFile $matrixFile \
		       -corFraction $corFraction \
		       -useDQS $useDQS \
		       -splitTasks $splitTasks \
		       -iterations $iterations \
		       -useTunes $useTunes \
		       -useMeasDisp $useMeasDisp \
		       -measDispFile $measDispFile \
		       -outputStatusDevice $outputStatusDevice \
		       -abortFile $abortFile} result] {
	    OutputStatusMessage "ScanSVnumber: $result \nInterrupt iterations."
	    exit 1
	} else {
	    OutputStatusMessage "Scan is done."
	}
    }
    plot {
	if [catch {CalculateFilesForPlotting \
		       -calcMode $calcMode \
		       -lteMeasured $lteMeasured \
		       -beamlineMeasured $beamlineMeasured \
		       -paramFileListMeasured $paramFileListMeasured \
		       -lteDesired $lteDesired \
		       -beamlineDesired $beamlineDesired \
		       -paramFileListDesired $paramFileListDesired \
		       -paramFileListAchieved $paramFileListAchieved \
		       -xbpmList $xbpmList \
		       -ybpmList $ybpmList \
		       -dbpmList $dbpmList \
		       -weightList $weightList \
		       -doneFile $doneFile \
		       -workDir $workDir \
		       -useTunes $useTunes \
		       -useRelative $useRelative \
		       -twiMeasured $twiMeasured \
		       -twiDesired $twiDesired \
		       -twiAchieved $twiAchieved} result] {
	    OutputStatusMessage "CalculateFilesForPlotting: $result"
	    exit 1
	}
    }
    minEmit {
	if [catch {MinimizeEmittance \
		       -workDir $workDir \
		       -quadList $quadList \
		       -lteMeasured $lteMeasured \
		       -beamlineMeasured $beamlineMeasured \
		       -paramFileListMeasured $paramFileListMeasured \
		       -optimTerm $optimTerm \
		       -rootname $rootname \
		       -emitArrayList $emitArrayList} result] {
	     OutputStatusMessage "MinimizeEmittance: $result"
	    exit 1
	}
    }
    calcEmit {
	if [catch {CalculateEmittances \
		       -calcMode $calcMode \
		       -workDir $workDir \
		       -lteMeasured $lteMeasured \
		       -beamlineMeasured $beamlineMeasured \
		       -paramFileListMeasured $paramFileListMeasured \
		       -lteDesired $lteDesired \
		       -beamlineDesired $beamlineDesired \
		       -paramFileListDesired $paramFileListDesired \
		       -paramFileListAchieved $paramFileListAchieved \
		       -emitArrayList $emitArrayList} result] {
	    OutputStatusMessage "CalculateEmittances: $result"
	    exit 1
	}
    }
    runDqsTask {
	if [catch {RunDqsTask -calcMode $calcMode -path $path -eleTemplate $eleTemplate -varList $varList -elementsUpdateFile $elementsUpdateFile \
		       -useRelative $useRelative -acceleratorCode $acceleratorCode -verbose $verbose} result] {
	    OutputStatusMessage "RunDqsTaks: $result"
	    exit 1
	}
    }
}
