#!/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)]

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

if [info exists env(LOCO_TMPDIR)] {
    set tmpDir $env(LOCO_TMPDIR)
} else {
    set tmpDir $env(HOME)/tmp
}
if ![file exists $tmpDir] {exec mkdir -p $tmpDir}
source $locoBinDir/locoCommonProcedures

set args $argv
set rememberArgs $args

set runMode ""
APSParseArguments {runMode}
if {![string length $runMode]} {
    puts stderr "Error: variable runMode is empty."
    puts stderr "Calling arguments (LOCO_BINDIR: $env(LOCO_BINDIR)): $rememberArgs"
    exit
}

if ![info exists env(ELEGANT_BINDIR)] {
    puts stderr "   Error: ELEGANT_BINDIR environment variable is not defined"; exit
} else {
    set accelCodeBinDir $env(ELEGANT_BINDIR)
}

#
#--- Argument checking first (too many arguments, need to check them all - good for debugging)...
#

switch -exact $runMode {
    getSDDSLattice {
	set fullVarList [list lteFile beamlineName latticeParamFileList workDir sddsLatticeFile bRho] 
	set requiredVarList [list lteFile beamlineName workDir sddsLatticeFile bRho]
	set defaultVarList ""
	set defaultValueList ""
    }
    getBetaFunctions {
	set fullVarList [list lteFile beamlineName workDir latticeParamFileList twiFile closedOrbit \
			     twissRefFile twissRefElement malignElement fixedOrbitLength \
			     verbose bRho]
	set requiredVarList [list lteFile beamlineName workDir twiFile bRho]
	set defaultVarList [list closedOrbit fixedOrbitLength fourD verbose elemByElemSR]
	set defaultValueList [list 1 0 1 0 0]
    }
    getResponseMatrix {
	set varList [list lteFile beamlineName workDir latticeParamFileList twiFile \
			 hrmFile vrmFile hCorrList vCorrList kickX kickY useKickFiles kickFileX kickFileY \
			 fitDispersion dispFitWeight dispColumn \
			 malignElement twissRefFile twissRefElement \
			 verbose sddsLatticeFile bRho \
			 bpmGainXFile bpmGainYFile bpmTiltFile bpmCoefFile bpmNoise bpmNoiseSeed dispNoise usePelegant]
	set calcList [list closedOrbit fixedOrbitLength directCalc coupledMatrix useDoubleOrbits elemByElemSR useElegant4bpmErrors \
			  directDispersionCalc dispDirectParamList additionalErrorSources additionalErrorSourcesSeed doneFile averageOrbits]
	set qsubList [list useQsub numberOfQsubTasks waitInterval waitTime qsubCommand qsubRespProcCommand \
			  submissionPause queueSystemName]
	set fullVarList [concat $varList $calcList $qsubList]
	set requiredVarList [list lteFile beamlineName workDir twiFile \
					hrmFile vrmFile hCorrList vCorrList directCalc kickX kickY coupledMatrix \
					fitDispersion dispFitWeight dispColumn useQsub \
					useDoubleOrbits bRho]
	set defaultVarList [list closedOrbit fixedOrbitLength directCalc fourD verbose usQsub numberOfQsubTasks submissionPause bpmNoise useKickFiles \
				elemByElemSR useElegant4bpmErrors dispNoise directDispersionCalc additionalErrorSourcesSeed averageOrbits usePelegant]
	set defaultValueList [list 1 0 0 1 0 0 1 0 0 0 0 0 0 0 1234567 0 0]
    }
    runCalculateOrbits {
	set fullVarList [list locoBinDir accelCodeBinDir latticeFile orbitFile workDir nonlin correctorFile corrList kick \
			     sddsLatticeFile bRho tmpDir fixedOrbitLength \
			     closedOrbit refCorr useDoubleOrbits useKickFiles kickFile verbose \
			     useElegant4bpmErrors bpmNoise bpmNoiseSeed additionalErrorSources additionalErrorSourcesSeed averageOrbits]
	set requiredVarList [list locoBinDir accelCodeBinDir latticeFile orbitFile workDir nonlin correctorFile corrList kick \
					sddsLatticeFile bRho tmpDir fixedOrbitLength \
					closedOrbit refCorr useDoubleOrbits]
	set defaultVarList [list useKickFiles elemByElemSR verbose useElegant4bpmErrors bpmNoise additionalErrorSourcesSeed averageOrbits]
	set defaultValueList [list 0 0 0 0 0 1234567 0]
    }
}
set error 0
foreach var $fullVarList {set $var ""}
if [llength $defaultVarList] {foreach varName $defaultVarList varValue $defaultValueList {set $varName $varValue}}
APSParseArguments $fullVarList
foreach var $requiredVarList {
    if ![string length [set $var]] {set error 1; puts stderr "Error: variable $var is empty."}
}
if $error {exit}

if {[string compare $runMode "getResponseMatrix"] == 0 && $directDispersionCalc} {
    if [string length $dispDirectParamList] {
	array set testArray $dispDirectParamList
	set argList [array names testArray]
	set error 0
	foreach arg [list dp0 dp1 Steps malignElementName] {
	    if {[lsearch $argList $arg] == -1} {puts stdout "$arg missing from \"$dispDirectParamList\""; set error 1}
	}
	if $error {exit}
    } else {puts stdout "Argument dispDirectParamList is missing. Should be \"dp0 <number> dp1 <number> Steps <number> malignElementName <name>\""; exit}

}

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

#--- Nondirect uses elegant response_matrix_output command. It can use beta functions (directCalc=0) or from orbits (directCalc=2)

proc NondirectElegantResponse {args} {
    global bRho
    set elemByElemSR 0
    APSParseArguments {accelCodeBinDir workDir lteFile beamlineName latticeParamFileList hCorrList vCorrList hrmFile vrmFile \
			   twiFile fixedOrbitLength coupledMatrix useResponseFromOrbits includeTwissCalc kickX kickY \
			   useQsub qsubCommand numberOfQsubTasks waitInterval waitTime elemByElemSR verbose}

    set deleteFiles ""
    set commandList ""
    set jobNameList ""
    set doneFileList ""
    set eleLogFileList ""

    set alterElemList [list "name * type RFCA item CHANGE_T value 0 allow_missing_elements 1"]
    if $elemByElemSR {
	set alterElemList [concat [list "name * type *BEN* item SYNCH_RAD value 1" "name * type *BEN* item ISR value 1" \
				       "name * type *BEN* item ISR1PART value 0"] $alterElemList]
    }

    #------ Prepare job for twiss calculation:
    if $includeTwissCalc {
	set tmpRoot $workDir/[APSTmpString]-calculateElegantBetas
	set eleFile $tmpRoot.ele
	set localTwiFile $tmpRoot.twi
	set localCloFile $tmpRoot.clo
	set doneFile $tmpRoot.done
	set elegantOutFile $tmpRoot.log
	lappend deleteFiles $eleFile $doneFile $elegantOutFile $localTwiFile $localCloFile
	set energy [format %10.4e [expr $bRho * 299792458e-6]]
	if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		       -run_setup [list "lattice $lteFile use_beamline $beamlineName default_order 2 semaphore_file $doneFile p_central_mev $energy"] \
		       -load_parameters [list "filename_list \"$latticeParamFileList\""] \
		       -alter_elements $alterElemList \
		       -run_control [list "n_steps 1"] \
		       -closed_orbit [list "fixed_length $fixedOrbitLength output $localCloFile"] \
		       -twiss_output [list "filename $localTwiFile matched 1 output_at_each_step 1"] \
		       -bunched_beam [list "n_particles_per_bunch 1"] \
		   } result] {
	    return -code error "Fit_GenerateElegantFileFromLTE1: $result"
	}
	lappend commandList "$accelCodeBinDir/elegant $eleFile > $elegantOutFile"
	lappend eleLogFileList $elegantOutFile
	lappend jobNameList beta
	lappend doneFileList $doneFile
    }

    #------ Prepare jobs for RM calculations:
    if [catch {Fit_SplitList -List [concat $hCorrList $vCorrList] -splitTasks $numberOfQsubTasks} listAfterSplitting] {
	return -code error "Fit_SplitList: $listAfterSplitting"
    }
    set rootList ""
    set tmpRoot0 $workDir/[APSTmpString]-singleElegResp
    set counter 0
    foreach taskList $listAfterSplitting {
	set tmpRoot ${tmpRoot0}$counter
	if [catch {SingleElegantResponse -accelCodeBinDir $accelCodeBinDir -elemByElemSR $elemByElemSR \
		       -lteFile $lteFile -beamlineName $beamlineName -kickX $kickX -kickY $kickY \
		       -latticeParamFileList $latticeParamFileList -corrList $taskList -twiFile $twiFile \
		       -coupledMatrix $coupledMatrix -fixedOrbitLength $fixedOrbitLength \
		       -useResponseFromOrbits $useResponseFromOrbits -tmpRoot $tmpRoot} result] {
	    return -code error "SingleElegantResponse: $result"
	}
	lappend commandList [lindex $result 0]
	lappend eleLogFileList [lindex $result 1]
	lappend jobNameList orbit$counter
	lappend doneFileList $tmpRoot.done
	lappend rootList $tmpRoot
	set deleteFiles [concat $deleteFiles [lindex $result 2]]
	incr counter
    }

    #------ If RM measurement orbit kick was too large, elegant will report diverging orbits. 
    #------ If that happens, reduce the kick and repeat upto 4 3 times using this loop:
    set kickReductionFactor 4
    set reduceKickValueIterMax 3
    set kickXactual $kickX
    set kickYactual $kickY
    for {set reduceKickValueIter 0} {$reduceKickValueIter < $reduceKickValueIterMax} {incr reduceKickValueIter} {
	#------ Submit and wait
	if $useQsub {
	    if [catch {Fit_SubmitJobs -commandList $commandList -jobNameList $jobNameList -doneFileList $doneFileList \
			   -verbose $verbose} returnList] {
		return -code error "Fit_SubmitJobs: $returnList"
	    }
	    set pidList [lindex $returnList 0]
	    set logFileList [lindex $returnList 1]
	    set scriptFileList [lindex $returnList 2]
	    set dotOFileList [lindex $returnList 3]
	    set deleteFiles [concat $deleteFiles $scriptFileList $dotOFileList]
	    if [catch {Fit_WaitForTasks -pidList $pidList -jobNameList $jobNameList -commandList $commandList \
			   -logFileList $logFileList -doneFileList $doneFileList -waitTime $waitTime \
			   -waitInterval $waitInterval -usePopupWindow 0 \
			   -useQsub $useQsub -abortFile [APSTmpString] \
			   -verbose $verbose} result] {
		return -code error "Fit_WaitForTasks: $result"
	    }
	} else {
	    foreach command $commandList {
		if [catch {eval exec $command} result] {
		    return -code error "$command: $result"
		}
	    }
	}
	#------ Analyze if elegant printed error in orbit convergence (elegant does not abort in case of orbit divergence):
	set repeatCalcs 0
	foreach logFile $eleLogFileList {
	    #--- grep gives error if no match
	    if ![catch {exec grep "error: closed orbit did not converge" $logFile} result] {set repeatCalcs 1}
	}
	if {$repeatCalcs == 0} {
	    break
	} else {
	    #--- Edit elegant files to reduce orbit kick for RM calculation:
	    if {$reduceKickValueIter == [expr $reduceKickValueIterMax-1]} {
		return -code error "======>>> Error:: elegant reported not converging orbits even after reducing orbit kicks $reduceKickValueIterMax times, 
                                   see $logFile."
	    }
	    if $verbose {
		puts stdout "======>>> WARNING:: "
		puts stdout "======>>> WARNING:: elegant reported not converging orbits, see $logFile."
		puts stdout "======>>> WARNING:: Reducing measurement kick by a factor of $kickReductionFactor and repeating calculations..."
	    }
	    set kickXactual [expr $kickXactual / $kickReductionFactor]
	    set kickYactual [expr $kickYactual / $kickReductionFactor]
	    foreach command $commandList {
		if [catch {ReduceOrbitKick -eleFile [lindex $command 1] -divideFactor $kickReductionFactor} result] {
		    return -code error "ReduceOrbitKick: $result"
		}
	    }
	}
    }

    if $includeTwissCalc {
	if [catch {AddTagNameColumn -filename $localTwiFile} result] {return -code error "AddTagNameColumn: $result"}
	file copy -force $localTwiFile $twiFile
	set orbitOutput [exec editstring -stream $twiFile -edit=%@.twi@.clo@]
	file copy -force $localCloFile $orbitOutput
    }

    #------ Analyze if elegant printed error in orbit convergence (elegant does not abort in case of orbit divergence):
    foreach logFile $eleLogFileList {
	if ![catch {exec grep "error: closed orbit did not converge" $logFile} result] {
	    #--- grep gives error if no match
	    return -code error "elegant reported not converging orbits, see $logFile."
	}
    }

    #------ Sometimes elegant mixes "name#N" and "name" columns. Fix it here:
    if [catch {FixColumnNames -rootList $rootList} result] {
	return -code error "FixColumnNames: $result"
    }

    #------ Assemble rm files:
    set tmpRoot $workDir/[APSTmpString]-nondirectResponse
    if {[llength $rootList] == 1} {
	file copy -force $rootList.hrm $tmpRoot.hrm
	file copy -force $rootList.vrm $tmpRoot.vrm
	file copy -force $rootList.hvrm $tmpRoot.hvrm
	file copy -force $rootList.vhrm $tmpRoot.vhrm
    } else {
	foreach root $rootList {
	    lappend hrmFileList $root.hrm
	    lappend vrmFileList $root.vrm
	    lappend hvrmFileList $root.hvrm
	    lappend vhrmFileList $root.vhrm
	}
	eval exec sddsxref $hrmFileList $tmpRoot.hrm -take=* -leave=BPMName,s -nowarning
	eval exec sddsxref $vrmFileList $tmpRoot.vrm -take=* -leave=BPMName,s -nowarning
	eval exec sddsxref $hvrmFileList $tmpRoot.hvrm -take=* -leave=BPMName,s -nowarning
	eval exec sddsxref $vhrmFileList $tmpRoot.vhrm -take=* -leave=BPMName,s -nowarning
    }

    #------ Reformat the way I need:
    if {$coupledMatrix == 2} {
	#------ Here I just blindly change suffix assuming that in the APS lattice file all the BPM are double plane. Filtering is done later anyway.
	lappend deleteFiles $tmpRoot.hrm $tmpRoot.vrm $tmpRoot.vhrm $tmpRoot.hvrm
	#------ Parameter file for direct response:
	exec sddsmakedataset -pipe=out -col=plane,type=string -data=X,Y \
	    | sddsbreak -pipe -rowlimit=1 \
	    | sddsprocess -pipe=in $tmpRoot.params -process=plane,first,Plane -process=plane,first,BpmPlane \
	    -process=plane,first,OrbitPlane 
	lappend deleteFiles $tmpRoot.params
	#------ hrm file:
	exec sddscombine $tmpRoot.hrm $tmpRoot.vhrm -pipe=out \
	    | sddsconvert -pipe -del=col,s \
	    | sddsxref -pipe=in $tmpRoot.params $hrmFile -leave=* -transfer=para,*Plane
	#------ vrm file:
	exec sddscombine $tmpRoot.hvrm $tmpRoot.vrm -pipe=out \
	    | sddsconvert -pipe -del=col,s \
	    | sddsxref -pipe=in $tmpRoot.params $vrmFile -leave=* -transfer=para,*Plane
    }
    exec sddsprocess $hrmFile -redef=para,ElemByElemSR,$elemByElemSR,type=long -redef=para,KickActual,$kickXactual -nowarning
    exec sddsprocess $vrmFile -redef=para,ElemByElemSR,$elemByElemSR,type=long -redef=para,KickActual,$kickYactual -nowarning
    file delete ${hrmFile}~ ${vrmFile}~
    eval file delete $deleteFiles
}

#----------------------------------------------------------------------------------------------------------------------------
proc ReduceOrbitKick {args} {
    APSParseArguments {eleFile divideFactor}
    set tmpRoot /tmp/[APSTmpString]-reduceOrbitKick
    set eleCopyFile $tmpRoot.ele
    if [catch {open $eleFile r} fid0] {return -code error "Error opening eleFile $eleFile: $fid0"}
    if [catch {open $eleCopyFile w} fid1] {return -code error "Error opening eleFile $eleCopyFile: $fid1"}
    while {[gets $fid0 fileLine] >= 0} {
	if {[string first "corrector_tweek" $fileLine] != -1} {
	    set kickString [string range $fileLine [expr [string first "=" $fileLine]+1] [string length $fileLine]]
	    set kickList [split $kickString ","]
	    set kickX [expr [lindex $kickList 0] / $divideFactor]
	    set kickY [expr [lindex $kickList 1] / $divideFactor]
	    set newFileLine [string range $fileLine 0 [string first "=" $fileLine]]
	    append newFileLine "[format %8.2e $kickX], [format %8.2e $kickY],"
	    puts $fid1 $newFileLine
	} else {
	    puts $fid1 $fileLine
	}
    }
    close $fid0
    close $fid1
    file copy -force $tmpRoot.ele $eleFile
    file delete $tmpRoot.ele
}
#----------------------------------------------------------------------------------------------------------------------------
#--- If elegant sees no same name correctors in the response matrix output, it will not put #N at the end of the name
#--- We need it always
proc FixColumnNames {args} {
    APSParseArguments {rootList}
    foreach rootname $rootList {
	set colTargetNameList [join [exec sdds2stream $rootname.use -col=TagName] ]
	set hCorrList [Fit_DeleteElementsFromList -elementList [exec sddsquery -col $rootname.hrm] -deleteElements "BPMName s"]
	set vCorrList [Fit_DeleteElementsFromList -elementList [exec sddsquery -col $rootname.vrm] -deleteElements "BPMName s"]
	set renameStringH ""
	set renameStringV ""
	for {set i 0} {$i < [llength $hCorrList]} {incr i} {
	    set colName [lindex $hCorrList $i]
	    set targetName [lindex $colTargetNameList $i]
	    if {[string compare $colName $targetName] != 0} {
		append renameStringH " -rename=col,$colName=$targetName "
	    }
	}
	for {set i 0} {$i < [llength $vCorrList]} {incr i} {
	    set colName [lindex $vCorrList $i]
	    set targetName [lindex $colTargetNameList [expr $i + [llength $hCorrList]]]
	    if {[string compare $colName $targetName] != 0} {
		append renameStringV " -rename=col,$colName=$targetName "
	    }
	}
	if {[llength $renameStringH] != 0} {
	    eval exec sddsconvert $rootname.hrm $renameStringH -nowarning
	    eval exec sddsconvert $rootname.vhrm $renameStringH -nowarning
	    file delete $rootname.hrm~ $rootname.vhrm~
	}
	if {[llength $renameStringV] != 0} {
 	    eval exec sddsconvert $rootname.vrm $renameStringV -nowarning
	    eval exec sddsconvert $rootname.hvrm $renameStringV -nowarning
	    file delete $rootname.vrm~ $rootname.hvrm~
	}
    }
}

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

proc SingleElegantResponse {args} {
    global bRho verbose usePelegant
    APSParseArguments {accelCodeBinDir lteFile beamlineName latticeParamFileList corrList twiFile elemByElemSR \
			   coupledMatrix fixedOrbitLength useResponseFromOrbits tmpRoot kickX kickY}
    set deleteFiles ""
    set elegantOutFile $tmpRoot.log
    set eleFile $tmpRoot.ele
    set doneFile $tmpRoot.done
    set hrmFile $tmpRoot.hrm
    set vrmFile $tmpRoot.vrm
    set hvrmFile $tmpRoot.hvrm
    set vhrmFile $tmpRoot.vhrm
    lappend deleteFiles $tmpRoot.log $tmpRoot.ele $tmpRoot.done $tmpRoot.hrm $tmpRoot.vrm $tmpRoot.hvrm $tmpRoot.vhrm

    exec sddsmakedataset -pipe=out -col=TagName,type=string -data=[join $corrList ,] \
	| sddsprocess -pipe=in $tmpRoot.use -print=col,ElementParameter,STEERING -def=col,ParameterValue,1
    if [catch {Fit_MakeElementNameFromTagName -input $tmpRoot.use -addElementOccurence 1} result] {
	return -code error "Fit_MakeElementNameFromTagName: $result"
    }
    lappend deleteFiles $tmpRoot.use
    set correctorsFile $tmpRoot.use

    #--- Create elegant file for RM calculations:
    switch -exact $coupledMatrix {
	0 {set correctionMatrixOutput "response \"$hrmFile $vrmFile\" output_at_each_step 1 fixed_length $fixedOrbitLength"}
	1 {set correctionMatrixOutput "response \"$hrmFile $vrmFile\" output_at_each_step 1 fixed_length $fixedOrbitLength coupled 1"}
	2 {set correctionMatrixOutput "response \"$tmpRoot.hrm $tmpRoot.vrm $tmpRoot.vhrm $tmpRoot.hvrm\" output_at_each_step 1 fixed_length $fixedOrbitLength coupled 1 use_response_from_computed_orbits $useResponseFromOrbits"}
    }

    set alterElemList [list "allow_missing_elements 1 allow_missing_parameters 1 name * type * value 0 item STEERING" \
			    "allow_missing_elements 1 allow_missing_parameters 1 name * type * value 0 item HSTEERING" \
			    "allow_missing_elements 1 allow_missing_parameters 1 name * type * value 0 item VSTEERING" \
			    "allow_missing_elements 1 allow_missing_parameters 1 name * type * value 0 item XSTEERING" \
			    "allow_missing_elements 1 allow_missing_parameters 1 name * type * value 0 item YSTEERING" \
			    "name * type RFCA item CHANGE_T value 0 allow_missing_elements 1"]
    if $elemByElemSR {
	set alterElemList [concat [list "name * type *BEN* item SYNCH_RAD value 1" "name * type *BEN* item ISR value 1" \
				       "name * type *BEN* item ISR1PART value 0"] $alterElemList]
    }

    #------ load_parameters2 always goes after alter_elements
    #------ default_order=2 is important when you have non-zero large orbits
    set energy [format %10.4e [expr $bRho * 299792458e-6]]
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -run_setup [list "lattice $lteFile use_beamline $beamlineName semaphore_file $doneFile p_central_mev $energy default_order 2 rootname $tmpRoot"] \
		   -load_parameters [list "filename_list \"[concat $latticeParamFileList]\""] \
		   -alter_elements $alterElemList \
		   -load_parameters2 [list "filename_list \"$correctorsFile\""] \
		   -run_control [list "n_steps 1"] \
		   -correct [list "mode orbit fixed_length $fixedOrbitLength fixed_length_matrix $fixedOrbitLength use_response_from_computed_orbits $useResponseFromOrbits corrector_tweek \"$kickX $kickY\""] \
		   -correction_matrix_output [list "$correctionMatrixOutput"] \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE1: $result"
    }
    lappend deleteFiles $tmpRoot.orb
    
    if $usePelegant {
	set command "/usr/lib64/mpich/bin/mpiexec -np 16 Pelegant $eleFile > $elegantOutFile"
    } else {
	set command "$accelCodeBinDir/elegant $eleFile > $elegantOutFile"
    }
    if $verbose {puts stdout $command}
    return [list $command $elegantOutFile $deleteFiles]
}

#----------------------------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------
#--- Main elegant procedure
#----------------------------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------

proc CalculateElegantRM {args} {
    global tmpDir useElegant4bpmErrors
    APSParseArguments {lteFile beamlineName latticeParamFileList includeTwissCalc elemByElemSR \
			   hrmFile vrmFile hCorrList vCorrList workDir useQsub fixedOrbitLength kickX kickY \
                           directCalc twiFile fitDispersion dispFitWeight dispColumn coupledMatrix closedOrbit \
			   qsubCommand numberOfQsubTasks waitInterval waitTime accelCodeBinDir verbose useDoubleOrbits \
			   twissRefFile twissRefElement malignElement useKickFiles kickFileX kickFileY \
			   bpmGainXFile bpmGainYFile bpmTiltFile bpmCoefFile bpmNoise bpmNoiseSeed dispNoise averageOrbits}

    set tmpRoot $tmpDir/[APSTmpString]-elegRM
    set deleteFiles ""
    set localParamFileList $latticeParamFileList
    #--- Remove BPM errors from parameter files:
    if !$useElegant4bpmErrors {
	for {set i 0} {$i < [llength $latticeParamFileList]} {incr i} {
	    set filename [lindex $latticeParamFileList $i]
	    set columnList [exec sddsquery -col $filename]
	    if {[lsearch -exact $columnList VariableType] != -1} {
		if {[exec sddsprocess $filename -pipe=out -nowarning "-match=col,VariableType=?Bpm,VariableType=bpmTilt,|,VariableType=bpmCoef,|" \
			 | sdds2stream -pipe=in -rows=bare] != 0} {
		    exec sddsprocess $filename $tmpRoot.tmpList -nowarning -match=col,VariableType=xBpm,! -match=col,VariableType=yBpm,! \
			-match=col,VariableType=bpmTilt,! -match=col,VariableType=bpmCoef,! 
		    set localParamFileList [lreplace $localParamFileList $i $i $tmpRoot.tmpList]
		    lappend deleteFiles $tmpRoot.tmpList
		}
	    } elseif {[lsearch -exact $columnList ElementType] != -1} {
		if {[exec sddsprocess $filename -pipe=out -nowarning -match=col,ElementType=*MON* \
			 | sdds2stream -pipe=in -rows=bare] != 0} {
		    exec sddsprocess $filename $tmpRoot.tmpList -nowarning -match=col,ElementType=*MON*,!
		    set localParamFileList [lreplace $localParamFileList $i $i $tmpRoot.tmpList]
		    lappend deleteFiles $tmpRoot.tmpList
		}
	    }
	}
    }

    #--- directCalc = 0 for beta function based RM calculations; 
    #--- directCalc = 1 for my own orbit-based calculations of RM;
    #--- directCalc = 2 for elegant orbit-based calculations of RM;
    if {$directCalc == 1} {
	if [catch {CalculateDirectRM -lteFile $lteFile -beamlineName $beamlineName -latticeParamFileList $localParamFileList \
		       -hrmFile $hrmFile -vrmFile $vrmFile \
		       -hCorrList $hCorrList -vCorrList $vCorrList -workDir $workDir -useQsub $useQsub \
		       -coupledMatrix $coupledMatrix -fitDispersion $fitDispersion -dispColumn $dispColumn \
		       -dispFitWeight $dispFitWeight -twiFile $twiFile -kickX $kickX -kickY $kickY \
		       -qsubCommand $qsubCommand -numberOfQsubTasks $numberOfQsubTasks -verbose $verbose \
		       -waitInterval $waitInterval -waitTime $waitTime \
		       -fixedOrbitLength $fixedOrbitLength -accelCodeBinDir $accelCodeBinDir \
		       -closedOrbit $closedOrbit -useDoubleOrbits $useDoubleOrbits \
		       -useKickFiles $useKickFiles -kickFileX $kickFileX -kickFileY $kickFileY \
		       -twissRefFile $twissRefFile -twissRefElement $twissRefElement -malignElement $malignElement \
		       -bpmNoise $bpmNoise -bpmNoiseSeed $bpmNoiseSeed -averageOrbits $averageOrbits} result] {
            return -code error "CalculateDirectRM: $result"
        }
    } else {
	if {$directCalc == 2} {set useResponseFromOrbits 1} else {set useResponseFromOrbits 0}
        if [catch {NondirectElegantResponse -workDir $workDir -lteFile $lteFile -beamlineName $beamlineName \
		       -latticeParamFileList $localParamFileList -twiFile $twiFile \
                       -hrmFile $hrmFile -vrmFile $vrmFile -hCorrList $hCorrList -vCorrList $vCorrList \
                       -fixedOrbitLength $fixedOrbitLength -coupledMatrix $coupledMatrix \
		       -useResponseFromOrbits $useResponseFromOrbits -verbose $verbose -elemByElemSR $elemByElemSR \
		       -useQsub $useQsub -qsubCommand $qsubCommand -numberOfQsubTasks $numberOfQsubTasks \
		       -waitInterval $waitInterval -waitTime $waitTime -kickX $kickX -kickY $kickY \
		       -accelCodeBinDir $accelCodeBinDir -includeTwissCalc $includeTwissCalc} result] {
            return -code error "NondirectElegantResponse: $result"
        }
    }

    foreach filename [list $hrmFile $vrmFile $twiFile] renameColumns [list 1 1 0] {
	if [catch {AddTagNameColumn -filename $filename -renameColumns $renameColumns} result] {return -code error "AddTagNameColumn: $result"}
    }
    if $fitDispersion {
        if [catch {AddDispersionColumn -hrmFile $hrmFile -vrmFile $vrmFile \
		       -twiFile $twiFile -dispColumn $dispColumn -closedOrbit $closedOrbit \
                       -dispFitWeight $dispFitWeight -coupledMatrix $coupledMatrix} result] {
            return -code error "AddDispersionColumn ($hrmFile -- $vrmFile -- $twiFile): $result"
        }
    }

    if !$useElegant4bpmErrors {
	#--- Apply BPM gain and tilt correction and add BPM noise:
	if {[string length $bpmGainXFile] && [string length $bpmGainYFile] && [string length $bpmTiltFile] && [string length $bpmCoefFile]} {
	    if {[file exists $bpmGainXFile] && [file exists $bpmGainYFile] && [file exists $bpmTiltFile] && [file exists $bpmCoefFile]} {
		#------ Set seed if not set above. Always use 1 as first digit because "date" could return 0 in front, which is invalid integer:
		if {[string length $bpmNoiseSeed] == 0} {set bpmNoiseSeed 1[string range [exec date +%N] 3 [string length [exec date +%N]]]}
		if [catch {Fit_ApplyBPMCorrections -hrmFile $hrmFile -vrmFile $vrmFile \
			       -coupledMatrix $coupledMatrix -bpmGainXFile $bpmGainXFile -bpmGainYFile $bpmGainYFile \
			       -bpmTiltFile $bpmTiltFile -bpmCoefFile $bpmCoefFile} result] {
		    return -code error "Fit_ApplyBPMCorrections: $result"
		}
	    }
	}
	if {$directCalc != 1 && $bpmNoise > 0} {
	    #--- For directCalc=1 bpm noise is added directly to orbits
	    if [catch {AddBpmNoise2RM -hrmFile $hrmFile -vrmFile $vrmFile -bpmNoise $bpmNoise -seed $bpmNoiseSeed -dispNoise $dispNoise} result] {
		return -code error "AddBpmNoise2RM: $result"
	    }
	}
    }
    catch {eval file delete $deleteFiles}
}

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

proc AddTagNameColumn {args} {
    set renameColumns 0
    APSParseArguments {filename renameColumns}
    set columnList [join [exec sddsquery $filename -col] ]
    if {[lsearch -exact $columnList TagName] != -1} {return}

    if {[lsearch -exact $columnList ElementOccurence] != -1} {
	#--- Twiss or clo files have ElementOccurence, RM files do not. We need to triple MONI rows with BPMName#X and BPMName#Y.
	exec sddsprocess $filename $filename.index -redef=col,Index,i_row -redef=col,PageIndex,i_page
	exec sddsprocess $filename.index $filename.bpmx -nowarning -match=col,ElementType=MONI "-print=col,TagName,%s#X#%1d,ElementName,ElementOccurence"
	exec sddsprocess $filename.index $filename.bpmy -nowarning -match=col,ElementType=MONI "-print=col,TagName,%s#Y#%1d,ElementName,ElementOccurence"
	exec sddsprocess $filename.index $filename.bpm  -nowarning -match=col,ElementType=MONI "-print=col,TagName,%s#%1d,ElementName,ElementOccurence"
	exec sddsprocess $filename.index $filename.nonBpms -nowarning "-match=col,ElementType=MONI,!" "-print=col,TagName,%s#%1d,ElementName,ElementOccurence"
	exec sddscombine $filename.nonBpms $filename.bpmx $filename.bpmy $filename.bpm -pipe=out -merge \
	    | sddssort -pipe -col=PageIndex \
	    | sddsbreak -pipe -change=PageIndex \
	    | sddssort -pipe -col=Index \
	    | sddsconvert -pipe=in $filename.tmp -del=col,Index,PageIndex
	file copy -force $filename.tmp $filename
	file delete $filename.nonBpms $filename.bpmx $filename.bpmy $filename.index $filename.bpm
    } else {
	#--- RM files do not have ElementOccurence column:
	#------ elegant can create a mixture of names in rm files when some names have #N in it and others don't.
	#------ This can be in both rows and columns of the rm files.
	#------ To handle it, add #1 at the end, then delete second # symbol and after it.
	if $renameColumns {
	    if [catch {exec sddsconvert $filename -pipe=out "-edit=col,*,e i@#1@ a 2S?@#@ 2d" \
			   | sddsconvert -pipe=in $filename.tmp -rename=col,BPMName#1=BPMName} result] {
		return -code error "Error renaming columns in $filename: $result"
	    }
	} else {
	    file copy -force $filename $filename.tmp
	}
	if [catch {exec sddsprocess $filename.tmp -pipe=out "-edit=col,BPMName1,BPMName,e i@#1@ a 2S?@#@ 2d" "-edit=col,BPMName2,BPMName1,S@#@5d" \
		       "-scan=col,ElementOccurence,BPMName1,%ld,type=long,edit=Z#" -print=col,TagName,%s#%s#%1d,BPMName2,Plane,ElementOccurence \
		       | sddsconvert -pipe=in $filename -del=col,BPMName1,BPMName2,ElementOccurence} result] {
	    return -code error "Error adding TagName column to $filename.tmp: $result"
	}
    }
    file delete $filename.tmp
}

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

proc CalculateElegantBetas {args} {
    global tmpDir bRho
    #------ don't always want to supply reference file
    set deleteFiles ""
    set closedOrbit 1
    set twissRefFile ""
    set twissRefElement ""
    APSParseArguments {lteFile beamlineName latticeParamFileList twiFile workDir accelCodeBinDir elemByElemSR \
			    closedOrbit twissRefFile twissRefElement fixedOrbitLength}
    set tmpRoot $tmpDir/[APSTmpString]-calculateElegantBetas
    set eleFile $tmpRoot.ele
    lappend deleteFiles $eleFile
    set localClosedOrbit $closedOrbit
    if {!$closedOrbit && (![string length $twissRefFile] || ![string length $twissRefElement])} {
	return -code error "Error: twissRefFile and/or twissRefElement are empty for closedOrbit=0 option."
    }
    #------ output_at_each_step is required to calculate twiss parameters on new orbit.
    set energy [format %10.4e [expr $bRho * 299792458e-6]]
    if $closedOrbit {
	set twissOption [list "filename $twiFile matched $closedOrbit output_at_each_step 1"]
	set orbitOutput [exec editstring -stream $twiFile -edit=%@.twi@.clo@]
	set alterElemList [list "name * type RFCA item CHANGE_T value 0 allow_missing_elements 1"]
	if $elemByElemSR {
	    set alterElemList [concat [list "name * type *BEN* item SYNCH_RAD value 1" "name * type *BEN* item ISR value 1" \
					   "name * type *BEN* item ISR1PART value 0"] $alterElemList]
	}
	if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		       -run_setup [list "lattice $lteFile use_beamline $beamlineName default_order 2 p_central_mev $energy"] \
		       -load_parameters [list "filename_list \"$latticeParamFileList\""] \
		       -alter_elements $alterElemList \
		       -run_control [list "n_steps 1"] \
		       -closed_orbit [list "fixed_length $fixedOrbitLength output $orbitOutput"] \
		       -twiss_output $twissOption \
		       -bunched_beam [list "n_particles_per_bunch 1"] \
		   } result] {
	    return -code error "Fit_GenerateElegantFileFromLTE1: $result"
	}
    } else {
	set twissOption [list "filename $twiFile matched $localClosedOrbit reference_file $twissRefFile reference_element $twissRefElement"]
	if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		       -run_setup [list "lattice $lteFile use_beamline $beamlineName default_order 2 p_central_mev $energy"] \
		       -load_parameters [list "filename_list \"$latticeParamFileList\""] \
		       -run_control [list "n_steps 1"] \
		       -twiss_output $twissOption \
		       -bunched_beam [list "n_particles_per_bunch 1"] \
		   } result] {
	    return -code error "Fit_GenerateElegantFileFromLTE1: $result"
	}
    }
    set elegantOutFile $tmpRoot.log
    if [catch {exec $accelCodeBinDir/elegant $eleFile > $elegantOutFile} result] {
        return -code error "Error running elegant with file $eleFile (see $elegantOutFile): $result"
    }
    lappend deleteFiles $elegantOutFile

    #------ Running elegant again to calculate etax starting from 0 (measurement for beamline):
    if !$closedOrbit {
	set betaValues [exec sddsprocess $twissRefFile -pipe=out -match=col,ElementName=$twissRefElement \
			    | sdds2stream -pipe=in -col=betax,betay]
	if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		       -run_setup [list "lattice $lteFile use_beamline $beamlineName default_order 2 p_central_mev $energy"] \
		       -load_parameters [list "filename_list \"$latticeParamFileList\""] \
		       -twiss_output [list "filename $twiFile.tmpDisp matched $localClosedOrbit beta_x [lindex $betaValues 0] beta_y [lindex $betaValues 1] eta_x 0.0"] \
		       -run_control [list "n_steps 1"] \
		       -bunched_beam [list "n_particles_per_bunch 1"] \
		   } result] {
	    return -code error "Fit_GenerateElegantFileFromLTE1: $result"
	}
	set elegantOutFile $tmpRoot.log
	if [catch {exec $accelCodeBinDir/elegant $eleFile > $elegantOutFile} result] {
	    return -code error "Error running elegant with file $eleFile (see $elegantOutFile): $result"
	}
	exec sddsxref $twiFile $twiFile.tmpDisp -take=etax,etay -rename=col,etax=etax_0 -rename=col,etay=etay_0 -nowarning
	lappend deleteFiles ${twiFile}~ $twiFile.tmpDisp
    }
    eval file delete $deleteFiles
}

#----------------------------------------------------------------------------------------------------------------------------
#--- inputFile is multipage file where each page is orbit due to one page of correctorFile with x and y in the same page
#--- outputFile has all x orbits in page 1, and all y orbits in page 2

proc Make2pageOrbitFile { args } {
    global tmpDir
    APSParseArguments {inputFile outputFile correctorFile closedOrbit useDoubleOrbits averageOrbits}
    
    if $averageOrbits {
	if [catch {AverageOrbits -orbitFile $inputFile -averageOrbits $averageOrbits -closedOrbit $closedOrbit} result] {
	    return -code error "AverageOrbits: $result"
	}
    }
    
    set tmpRoot $tmpDir/[APSTmpString]-2pageOrbitFile
    set bpmTypeList [list HMON VMON]
    set planeList [list X Y]
    set bpmPlaneList [list X Y]
    set orbitPlaneList [list X Y]
    if $closedOrbit {set orbitColumnList [list x y]} else {set orbitColumnList [list Cx Cy]}
    set outputList [list $tmpRoot.1 $tmpRoot.2]
    set deleteFiles ""
    if $averageOrbits {
	exec sddsprocess $correctorFile $tmpRoot.corr "-redef=para,Flag,i_page 1 - $averageOrbits mod" -filter=para,Flag,0,0
    } else {
	file copy -force $correctorFile $tmpRoot.corr
    }
    if $useDoubleOrbits {
	set Ncorr [expr [exec sdds2stream $inputFile -npages=bare] / $useDoubleOrbits]
	set fileList ""
	for {set iCorr 0} {$iCorr < $Ncorr} {incr iCorr} {
	    set colName [exec sddsconvert $tmpRoot.corr -pipe=out -keep=[expr $iCorr * $useDoubleOrbits + 1] \
			     | sdds2stream -pipe=in -para=ColumnName]
	    if [catch {exec sddsxref $inputFile $tmpRoot.corr -pipe=out -leave=* -transfer=para,ColumnName \
			   | sddsprocess -pipe -match=para,ColumnName=$colName -match=col,ElementType=*MON* \
			   "-redef=para,NormalizedKick,-1 i_page 1 - 2 * $useDoubleOrbits 1 - / +" \
			   | sddsenvelope -pipe -slope=NormalizedKick,[join $orbitColumnList ,] -copy=TagName \
			   | sddstranspose -pipe -newColumn=TagName -oldColumn=Orbit \
			   | sddsprocess -pipe -print=col,Corrector,$colName \
			   | sddscombine -pipe -merge \
			   | sddsbreak -pipe -rowlimit=1 \
			   | sddstranspose -pipe -newColumn=Corrector -oldColumn=TagName \
			   | sddsconvert -pipe -del=para,* \
			   | sddsprocess -pipe=in $tmpRoot.$colName} result] {
		return -code error "Error processing double orbits ($colName -- $inputFile -- $tmpRoot.corr): $result"
	    }
	    lappend fileList $tmpRoot.$colName
	}
	if {[llength $fileList] == 1} {
	    file copy -force $fileList $tmpRoot
	} else {
	    eval exec sddsxref $fileList $tmpRoot -take=* -leave=TagName
	}
	foreach index [list 1 2] Plane [list X Y] {
	    exec sddsconvert $tmpRoot -pipe=out -keep=$index \
		| sddsprocess -pipe=in $tmpRoot.$index -reprint=para,Plane,$Plane -reprint=para,BpmPlane,$Plane -reprint=para,OrbitPlane,$Plane \
		-redef=para,UseDoubleOrbits,$useDoubleOrbits,type=long
	}
	eval file delete $fileList
	lappend deleteFiles $tmpRoot $tmpRoot.1 $tmpRoot.2 $tmpRoot.corr
    } else {
	foreach bpmType $bpmTypeList plane $planeList bpmPlane $bpmPlaneList orbitPlane $orbitPlaneList \
	    orbitColumn $orbitColumnList output $outputList {
		if [catch {exec sddsprocess $inputFile -pipe=out "-match=col,ElementType=$bpmType,ElementType=MONI,|" \
			       | sddsconvert -pipe -retain=col,TagName,$orbitColumn \
			       | sddstranspose -pipe -newColumn=TagName -oldColumn=Orbit \
			       | sddsxref -pipe $tmpRoot.corr -leave=* -transfer=para,ColumnName \
			       | sddsprocess -pipe -print=col,Corrector,%s,ColumnName \
			       | sddscombine -pipe -merge \
			       | sddstranspose -pipe -newColumn=Corrector -oldColumn=TagName \
			       | sddsconvert -pipe -del=para,* \
			       | sddsprocess -pipe=in $output -reprint=para,Plane,$plane -reprint=para,BpmPlane,$bpmPlane \
			       -reprint=para,OrbitPlane,$orbitPlane -redef=para,UseDoubleOrbits,$useDoubleOrbits,type=long} result] {
		    return -code error "Error processing response ($inputFile -- $tmpRoot.corr): $result"
		}
		lappend deleteFiles $output $tmpRoot.corr
	    }
    }
    exec sddscombine $tmpRoot.1 $tmpRoot.2 $outputFile -overWrite
    eval file delete $deleteFiles
}

#----------------------------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------
#------ Common procedures ---------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------

proc AddDispersionColumn {args} {
    global tmpDir directDispersionCalc dispFile dispDirectParamList
    APSParseArguments { hrmFile vrmFile twiFile dispColumn dispFitWeight coupledMatrix closedOrbit }
    set tmpRoot $tmpDir/[APSTmpString]-addDisp

    set tmpFileList ""
    if $directDispersionCalc {
	set localDispFile $dispFile
	#------ Check that element name in dispDirectParamList is not missing in the lattice:
	array set dispParamArray $dispDirectParamList
	set malinS [exec sddsprocess $twiFile -pipe=out -match=col,TagName=$dispParamArray(malignElementName) -nowarning \
			| sdds2stream -pipe=in -col=s]
	if {[string length $malinS] == 0} {return -code error "====== Error: name $dispParamArray(malignElementName) is missing from $twiFile"}
	if {$malinS != 0} {puts stdout "Warning: dispersion measurement element $dispParamArray(malignElementName) is not in the beginning of the beamline."}
    } else {
	if $closedOrbit {
	    set etaxCol etax; set etayCol etay
	} else {
	    set etaxCol etax_0; set etayCol etay_0
	}
	exec sddsconvert $twiFile $tmpRoot.dx -retain=col,TagName,$etaxCol -rename=col,$etaxCol=$dispColumn
	exec sddsconvert $twiFile $tmpRoot.dy -retain=col,TagName,$etayCol -rename=col,$etayCol=$dispColumn
	lappend tmpFileList $tmpRoot.dx $tmpRoot.dy

	switch -regexp -- $coupledMatrix {
	    0 {
		set localDispFile $tmpRoot.dx
	    }
	    1 {
		set localDispFile $tmpRoot.dy
	    }
	    2 {
		exec sddscombine $tmpRoot.dx $tmpRoot.dy $tmpRoot.d1 -overWrite
		set localDispFile $tmpRoot.d1
		lappend tmpFileList $tmpRoot.d1
	    }
	}
    }
    exec sddsxref $hrmFile $localDispFile -pipe=out -nowarning -take=$dispColumn -match=TagName \
	| sddsprocess -pipe=in $tmpRoot.hrm "-redef=col,$dispColumn,$dispColumn $dispFitWeight *"
    file copy -force $tmpRoot.hrm $hrmFile
    eval file delete $tmpFileList $tmpRoot.hrm
}

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

proc Truncate2pageFiles {args} {
    APSParseArguments {hrmFile vrmFile coupledMatrix}
    if {$coupledMatrix == 2} { return }
    switch -exact -- $coupledMatrix {
        0 { set pageX 1; set pageY 2 }
	1 { set pageX 2; set pageY 1 }
    }
    if [llength $hrmFile] {
	exec sddsconvert $hrmFile -nowarning -keepPage=$pageX
	file delete $hrmFile~
    }
    if [llength $vrmFile] {
	exec sddsconvert $vrmFile -nowarning -keepPage=$pageY
	file delete $vrmFile~
    }
}

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

proc SplitMultipageFile {args} {
    global tmpDir
    APSParseArguments {filename splitTasks}
    set tmpRoot $tmpDir/[APSTmpString]-splitMultipageFile
    set fileList ""
    if [catch {exec sdds2stream $filename -npages=bare} nOrbits] {return -code error $nOrbits}
    #--- Make splitTasks at least equal nOrbits:
    if {$nOrbits < $splitTasks} {set splitTasks $nOrbits}
    #--- No sense in having more than (nOrbits/2+1) tasks unless it equals to nOrbits:
    if {$splitTasks < $nOrbits && $splitTasks > [exec rpnl "$nOrbits 2 / int $nOrbits 2 mod +"]} {set splitTasks [exec rpnl "$nOrbits 2 / int $nOrbits 2 mod +"]}
    #--- Orbits per file:
    set orbitsPerFile [expr $nOrbits / $splitTasks]
    set remainder [exec rpnl "$nOrbits $splitTasks mod"]
    #--- Remainder is spread out over first $remainder pages:
    set sum 0
    for {set i 0} {$i < $splitTasks} {incr i} {
	set nn $orbitsPerFile
	if {$i < $remainder} {incr nn}
	exec sddsprocess $filename $tmpRoot.$i -redef=para,PageNumber,i_page -filter=para,PageNumber,[expr $sum+1],[expr $sum+$nn]
	lappend fileList $tmpRoot.$i
	incr sum $nn
    }
    if {$sum != $nOrbits} {return -code error "Error splitting multi-page file: total number of tasks does not equal to original number of orbits."}
    return $fileList
}

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

proc CalculateDirectRM {args} {

    global tmpDir locoBinDir directDispersionCalc useElegant4bpmErrors additionalErrorSources additionalErrorSourcesSeed dispDirectParamList

    #--- Here for empty variables I need "dummy" because APSParseArguments doesn't work with commandOptions
    set rmDivFactor 1
    set sddsLatticeFile "dummy"
    set nonlin "0"
    set useQsub 0
    set waitInterval 5
    set waitTime 600
    set bRho "dummy"
    set fixedOrbitLength "dummy"
    set verbose 0
    set closedOrbit 1
    set kickFileX dummy
    set kickFileY dummy
    APSParseArguments {lteFile beamlineName latticeParamFileList \
			   sddsLatticeFile hrmFile vrmFile nonlin hCorrList vCorrList workDir useQsub coupledMatrix \
                           fitDispersion dispColumn dispFitWeight twiFile bRho kickX kickY \
			   qsubCommand numberOfQsubTasks waitInterval waitTime fixedOrbitLength \
			   accelCodeBinDir verbose closedOrbit useDoubleOrbits twissRefFile twissRefElement malignElement \
			   useKickFiles kickFileX kickFileY bpmNoise bpmNoiseSeed averageOrbits}

    set refCorr RefTraj
    set tmpRoot $tmpDir/[APSTmpString]-calcDirectRM
    set counter 0
    set planeList [list X Y]
    set rmFileList [list $hrmFile $vrmFile]
    if !$useDoubleOrbits {
	set corListList [list [concat $refCorr $hCorrList] [concat $refCorr $vCorrList]]
    } else {
	set corListList [list $hCorrList $vCorrList]
    }

    set kickListMain [list [expr $kickX / $rmDivFactor] [expr $kickY / $rmDivFactor]]
    set kickFileList [list $kickFileX $kickFileY]

    set toMeters 1.0
    set rmFileRoot rm.ele
    set scriptName $locoBinDir/locoCalculateResponseMatrix
    set dqsWorkDir $workDir/dqs/orbit
    set rootJobName orbit
    set totalPathList ""
    set jobNameList ""
    set doneFileList ""
    set commandList ""
    catch {exec /bin/rm -r $dqsWorkDir}
    if ![file exists $workDir/dqs] { exec mkdir $workDir/dqs }
    if ![file exists $dqsWorkDir]  { exec mkdir $dqsWorkDir }
    foreach plane $planeList rmFile $rmFileList corList $corListList kick $kickListMain kickFile $kickFileList {
	if [llength $corList] {
	    set fileList$plane ""
	    if $useKickFiles {
		if [catch {SplitMultipageFile -filename $kickFile -splitTasks $numberOfQsubTasks} listAfterSplitting] {
		    return -code error "SplitMultipageFile: $listAfterSplitting"
		}
	    } else {
		if [catch {Fit_SplitList -List $corList -splitTasks $numberOfQsubTasks} listAfterSplitting] {
		    return -code error "Fit_SplitList: $listAfterSplitting"
		}
	    }

	    #--- Initiate local values for these variables (they will be redefined later if directDispersionCalc=1):
	    foreach [list useDoubleOrbitsLocal useKickFilesLocal fixedOrbitLengthLocal] [list $useDoubleOrbits $useKickFiles $fixedOrbitLength] {}

	    #--- Parameter file to calculate dispersion directly:
	    if {[string compare "X" $plane] == 0 && $directDispersionCalc} {
		if [catch {MakeDispersionParamFile -paramFile $tmpRoot.disp.param -closedOrbit $closedOrbit} result] {
		    return -code error "MakeDispersionParamFile: $result"
		}
		lappend listAfterSplitting $tmpRoot.disp.param
	    }

	    foreach taskList $listAfterSplitting {
		if {[string compare $taskList $tmpRoot.disp.param] == 0} {
		    foreach [list useDoubleOrbitsLocal useKickFilesLocal fixedOrbitLengthLocal] [list 0 1 0] {}
		}
		incr counter
		set totalPath $dqsWorkDir/$counter
		catch {file delete $totalPath} 
		catch {exec mkdir $totalPath}
		set jobName $rootJobName$counter
		set orbitFile $totalPath/orbit.sdds
		set correctorFile $totalPath/rmTemplate.param
		set rmLatticeFile $totalPath/$rmFileRoot
		if [catch {MakeFileForOrbitCalcs -lteFile $lteFile -beamlineName $beamlineName \
			       -latticeParamFileList $latticeParamFileList \
			       -rmLatticeFile $rmLatticeFile -closedOrbit $closedOrbit \
			       -correctorFile $correctorFile} result] {
		    return -code error "MakeFileForOrbitCalcs: $result"
		}
		set commandOptions " -runMode runCalculateOrbits -latticeFile $rmLatticeFile -orbitFile $orbitFile "
		append commandOptions " -tmpDir $tmpDir -workDir $totalPath -nonlin $nonlin -correctorFile $correctorFile "
		append commandOptions " -corrList \"$taskList\" -kick $kick -sddsLatticeFile $sddsLatticeFile "
		append commandOptions " -bRho $bRho -locoBinDir $locoBinDir -closedOrbit $closedOrbit "
		append commandOptions " -accelCodeBinDir $accelCodeBinDir -fixedOrbitLength $fixedOrbitLengthLocal "
		append commandOptions " -useDoubleOrbits $useDoubleOrbitsLocal -refCorr $refCorr -useKickFiles $useKickFilesLocal "
		append commandOptions " -averageOrbits $averageOrbits -verbose $verbose"
		if {!$useElegant4bpmErrors && $bpmNoise > 0} {
		    append commandOptions " -useElegant4bpmErrors $useElegant4bpmErrors -bpmNoise $bpmNoise -bpmNoiseSeed \"$bpmNoiseSeed\" "
		}
		if [string length $additionalErrorSources] {
		    append commandOptions " -additionalErrorSources \"$additionalErrorSources\" -additionalErrorSourcesSeed $additionalErrorSourcesSeed"
		}
		if [string length $bpmNoiseSeed] {incr bpmNoiseSeed 131}
		if [string length $additionalErrorSourcesSeed] {incr additionalErrorSourcesSeed 71}
		if $useKickFilesLocal {
		    file copy -force $taskList $correctorFile
		    if {[string compare $taskList $tmpRoot.disp.param] != 0} {file delete $taskList}
		}
		set command "$scriptName "
		append command $commandOptions
		lappend totalPathList $totalPath
		lappend jobNameList $jobName
		lappend doneFileList $totalPath/task.completed
		lappend commandList $command
		lappend fileList$plane $orbitFile
		#--- If not useQsub - run jobs one by one right here:
		if !$useQsub {
		    if $verbose {puts stdout $command}
		    if [catch {eval exec $command >@ stdout} result] {
			return -code error "$scriptName: $result"
		    }
		}
	    }
	}
    }

    if $useQsub {
	if [catch {Fit_SubmitJobs -commandList $commandList -jobNameList $jobNameList -doneFileList $doneFileList \
		       -verbose $verbose} returnList] {
	    return -code error "Fit_SubmitJobs: $returnList"
	}
	set pidList [lindex $returnList 0]
	set logFileList [lindex $returnList 1]
	if [catch {Fit_WaitForTasks -pidList $pidList -jobNameList $jobNameList -commandList $commandList \
		       -logFileList $logFileList -doneFileList $doneFileList -waitSec $waitTime \
		       -waitInterval $waitInterval -usePopupWindow 0 \
		       -useQsub $useQsub -abortFile [APSTmpString] \
		       -verbose $verbose} result] {
	    return -code error "Fit_WaitForTasks: $result"
	}
    }

    #------ Separate dispersion related orbits from response related orbits:
    if $directDispersionCalc {
	set dispOrbitFile [lindex $fileListX [expr [llength $fileListX]-1]]
	set paramFile [lindex $totalPathList [expr [llength $fileListX]-1]]/rmTemplate.param
	set fileListX [lreplace $fileListX [expr [llength $fileListX]-1] [expr [llength $fileListX]-1] ]
	if [catch {CalculateDispersionDirectly -orbitFile $dispOrbitFile -paramFile $tmpRoot.disp.param} result] {
	    return -code error "CalculateDispersionDirectly: $result"
	}
	file delete $tmpRoot.disp.param
    }

    #------ Calculate response from orbits:
    foreach fileList [list $fileListX $fileListY] rmFile $rmFileList kick $kickListMain {
	if {[llength $fileList] > 1} {
	    eval exec sddsxref -nowarning $fileList $tmpRoot.xref -leave=*Name 
	} else {
	    file copy -force $fileList $tmpRoot.xref
	}
	#--- Factor toMeters is 1 for elegant (elegant orbit is in m).
	if $useDoubleOrbits {
	    #--- For useDoubleOrbits, the file already contains orbit slopes (normalized to 1, not to $kick)
	    if [catch {exec sddsprocess $tmpRoot.xref $rmFile \
			   "-redef=col,%s,%s $kick / $toMeters /,select=*,exclude=*Name,units=m/rad" \
			   -redef=para,Kick,$kick -redef=para,DirectDispersionCalc,$directDispersionCalc,type=long \
			   -redef=para,AverageOrbits,$averageOrbits,type=long -redef=para,FixedOrbitLength,$fixedOrbitLength,type=long \
			   -print=para,AdditionalErrorSources,$additionalErrorSources \
			   -print=para,DispDirectParamList,$dispDirectParamList} result] {
		return -code error "Error combining orbit files: $result"
	    }
	} else {
	    if [catch {exec sddsconvert $tmpRoot.xref -pipe=out -del=col,TagName \
			   | sddsprocess -pipe -redef=para,Kick,$kick \
			   "-redef=col,%s,%s $refCorr - $kick / $toMeters /,select=*,exclude=$refCorr,units=m/rad" \
			   -redef=para,DirectDispersionCalc,$directDispersionCalc,type=long -redef=para,AverageOrbits,$averageOrbits,type=long \
			   -redef=para,FixedOrbitLength,$fixedOrbitLength,type=long \
			   -print=para,AdditionalErrorSources,$additionalErrorSources \
			   -print=para,DispDirectParamList,$dispDirectParamList \
			   | sddsconvert -pipe -del=col,$refCorr \
			   | sddsxref -pipe=in [lindex $fileList 0] $rmFile -take=TagName} result] {
		return -code error "Error combining orbit files ($tmpRoot.xref): $result"
	    }
	}
	if [catch {Fit_MakeElementNameFromTagName -input $rmFile -elementNameColumn BPMName} result] {
	    return -code error "Fit_MakeElementNameFromTagName: $result"
	}
    }
    file delete $tmpRoot.xref

    foreach path $totalPathList {
	if {$verbose < 2} {
	    catch {eval file delete [glob -nocomplain -- $path/*]}
	    catch {file delete $path}
	}
    }
    if [catch {Truncate2pageFiles -coupledMatrix $coupledMatrix -hrmFile $hrmFile -vrmFile $vrmFile} result] {
        return -code error "Truncate2pageFiles: $result"
    }
}

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

proc MakeFileForOrbitCalcs { args } {
    global bRho
    APSParseArguments { lteFile beamlineName latticeParamFileList rmLatticeFile correctorFile closedOrbit }

    set energy [format %10.4e [expr $bRho * 299792458e-6]]
    if $closedOrbit {
	set latticeOption [list "lattice $lteFile use_beamline $beamlineName p_central_mev $energy"]
	set orbitOption [list "output <ORBITFILE> output_monitors_only 1 fixed_length <FIXEDORBITLENGTH> closed_orbit_accuracy <CLOACCURACY>"]
    } else {
	set latticeOption [list "lattice $lteFile use_beamline $beamlineName centroid <ORBITFILE> p_central_mev $energy"]
	set orbitOption ""
    }
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $rmLatticeFile \
		   -run_setup $latticeOption \
		   -load_parameters [list "filename_list \"$latticeParamFileList\"" \
					 "filename $correctorFile change_defined_values 0"] \
		   -closed_orbit $orbitOption \
		   -run_control [list "n_steps <NSTEPS>"] \
		   -bunched_beam [list "n_particles_per_bunch 1"] \
	       } result] {
	return -code error "Fit_GenerateElegantFileFromLTE: $result"
    }
}

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

proc MakeSddsLattice { args } {
    global tmpDir bRho
    set tmpRoot $tmpDir/[APSTmpString]-makeSddsLattice
    APSParseArguments { lteFile beamlineName latticeParamFileList accelCodeBinDir \
			    workDir sddsLatticeFile }

    set eleFile $tmpRoot.ele
    set twiFile $tmpRoot.twi
    set paramFile $tmpRoot.param
    set energy [format %10.4e [expr $bRho * 299792458e-6]]
    if [catch {Fit_GenerateElegantFileFromLTE1 -eleOutputFile $eleFile \
		   -run_setup [list "lattice $lteFile use_beamline $beamlineName parameters $paramFile p_central_mev $energy"] \
		   -load_parameters [list "filename_list \"$latticeParamFileList\""] \
		   -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 elegantOutFile $tmpRoot.log
    file delete $elegantOutFile
    if [catch {exec $accelCodeBinDir/elegant $eleFile > $elegantOutFile} result] {
	return -code error "Error running elegant with file \"$eleFile\" (see $elegantOutFile): $result"
    }
    if [catch {exec sddsprocess $paramFile -pipe=out -match=col,ElementParameter=L \
		   | sddsxref -pipe=in $twiFile $sddsLatticeFile -take=s -match=ElementName -nowarning} result] {
	return -code error "Error taking column s from $twiFile to $paramFile: $result"
    }
    if [catch {AddTagNameColumn -filename $sddsLatticeFile} result] {return -code error "AddTagNameColumn: $result"}
    file delete $elegantOutFile $twiFile $paramFile $eleFile
}

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

proc CalculateOrbits { args } {
    global tmpDir useElegant4bpmErrors additionalErrorSources
    APSParseArguments {latticeFile orbitFile workDir tmpDir nonlin correctorFile corrList kick verbose \
			    sddsLatticeFile bRho accelCodeBinDir fixedOrbitLength closedOrbit refCorr useDoubleOrbits useKickFiles averageOrbits}
    set tmpRoot $tmpDir/[APSTmpString]-orbit

    if !$useKickFiles {
	#------ Make kick list:
	#------ +kick, -kick for double orbits; +kick for single orbits but with 0 for ref corrector
	if $useDoubleOrbits {
	    foreach {corrKickList kickList corrList1} [list "" "" ""] {}
	    for {set i 0} {$i < $useDoubleOrbits} {incr i} {lappend corrKickList [expr -1.0*$kick + 2*$i*$kick/($useDoubleOrbits-1)]}
	    foreach corr $corrList {
		set kickList [concat $kickList $corrKickList]
		for {set i 0} {$i < $useDoubleOrbits} {incr i} {lappend corrList1 $corr}
	    }
	    set corrList $corrList1
	} else {
	    set kickList ""
	    foreach corr $corrList {
		if {[string compare $corr $refCorr] == 0} {lappend kickList 0} else {lappend kickList $kick}
	    }
	}
	if [catch {MakeKickFile -correctorFile $correctorFile -corrList $corrList -kickList $kickList \
		       -useDoubleOrbits $useDoubleOrbits} result] {
	    return -code error "MakeKickFile: $result"
	}
    } else {
	if $useDoubleOrbits {
	    #--- Remove refCorr page
	    exec sddsprocess $correctorFile -nowarning "-match=para,ColumnName=$refCorr,!"
	    file delete ${correctorFile}~
	    #--- Add pages to handle useDoubleOrbits 
	    if [catch {ExpandKickFilePages -correctorFile $correctorFile -useDoubleOrbits $useDoubleOrbits} result] {
		return -code error "ExpandKickFilePages: $result"
	    }
	}
    }

    #--- Multiply (copy with different noise seeds) every page $averageOrbits times
    if $averageOrbits {
	if [catch {MultiplyKickFilePages -correctorFile $correctorFile -averageOrbits $averageOrbits} result] {
	    return -code error "MultiplyKickFilePages: $result"
	}
    }

    #--- Add additional errors to the kick file:
    if {[string length $additionalErrorSources] != 0} {
	if [catch {AddAdditionalErrorSources -additionalErrorSources $additionalErrorSources -correctorFile $correctorFile} result] {
	    return -code error "AddAdditionalErrorSources: $result"
	}
    }

    set Nsteps [exec sdds2stream $correctorFile -npages=bare]

    set elegantOutFile $tmpRoot.out
    set cloAccuracy 1e-7
    set cloAccuracyMax 2e-6
    set orbitError 1
    while {$orbitError == 1} {
	if {$cloAccuracy > $cloAccuracyMax} {
	    return -code error "Error:: elegant reports not converging orbits even when increasing closed_orbit_accuracy to $cloAccuracyMax"
	}
	catch {file delete $elegantOutFile}
	set command "$accelCodeBinDir/elegant $latticeFile -macro=NSTEPS=$Nsteps,ORBITFILE=$orbitFile.eleg,FIXEDORBITLENGTH=$fixedOrbitLength,CLOACCURACY=$cloAccuracy"
	if $verbose {puts stdout $command}
	if [catch {eval exec $command > $elegantOutFile} result] {
	    return -code error "Error running elegant with file $latticeFile (look in $elegantOutFile file): $result"
	} else {
	    #--- Number of pages in the clo file could be less that Nsteps, if some orbits didn't converge (elegant does not return error 
	    #--- for non-converging orbits).
	    #--- Repeat with increased closed_orbit_accuracy:
	    set NstepsOut [exec sdds2stream $orbitFile.eleg -npages=bare]
	    if {$NstepsOut != $Nsteps} {set cloAccuracy [expr $cloAccuracy * 2]} else {set orbitError 0}
	}
    }
    file delete $elegantOutFile

    #------ Add BPM noise:
    if !$useElegant4bpmErrors {
	if [catch {AddBpmNoise2Orbit -orbitFile $orbitFile.eleg} result] {return -code error "AddBpmNoise2Orbit: $result"}
    }

    #------ Check if any particle was lost -- not done.

    if [catch {AddTagNameColumn -filename $orbitFile.eleg} result] {return -code error "AddTagNameColumn: $result"}

    #--- Transform from multi-page orbit file where X and Y are on the same page to 2-page file with x and Y on different pages:
    if [catch {Make2pageOrbitFile -inputFile $orbitFile.eleg -outputFile $orbitFile \
		   -correctorFile $correctorFile -closedOrbit $closedOrbit \
		   -useDoubleOrbits $useDoubleOrbits -averageOrbits $averageOrbits} result] {
	return -code error "Make2pageOrbitFile: $result"
    }
}

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

#----------------------------------------------------------------------------------------------------------------------------
# Average every $averageOrbits pages
proc AverageOrbits {args} {
    global tmpDir
    APSParseArguments {orbitFile averageOrbits closedOrbit}
    set nPagesTotal [exec sdds2stream $orbitFile -npages=bare]
    set nKicks [expr $nPagesTotal / $averageOrbits]
    if $closedOrbit {set orbitColumnList [list x y]} else {set orbitColumnList [list Cx Cy]}
    set tmpRoot $tmpDir/[APSTmpString]-averOrbits
    set fileList ""
    for {set i 0} {$i < $nKicks} {incr i} {
	exec sddsprocess $orbitFile -pipe=out "-redef=para,PageIndex,i_page 1 -" \
	    -filter=para,PageIndex,[expr $i*$averageOrbits],[expr ($i+1)*$averageOrbits - 1] \
	    | sddsenvelope -pipe=in $tmpRoot.$i -mean=[join $orbitColumnList ,] -copy=TagName,ElementType
	lappend fileList $tmpRoot.$i
    }
    eval exec sddscombine $fileList -pipe=out \
	| sddsconvert -pipe=in $orbitFile \"-edit=col,*Mean,%@Mean@@\"
}
#----------------------------------------------------------------------------------------------------------------------------
proc CalculateDispersionDirectly {args} {
    global tmpDir dispDirectParamList dispColumn dispFile
    APSParseArguments {orbitFile paramFile}
    set tmpRoot $tmpDir/[APSTmpString]-dirDisp
    exec sddscombine $paramFile -pipe=out -merge \
	| sddsprocess -pipe=in $tmpRoot.dp -match=col,ElementParameter=DP
    foreach Plane [list X Y] {
	exec sddsprocess $orbitFile -pipe=out -match=para,Plane=$Plane \
	    | sddsbreak -pipe -rowlimit=1 \
	    | sddsprocess -pipe -process=TagName,first,TagName \
	    | sddstranspose -pipe -root=Orbit \
	    | sddsxref -pipe $tmpRoot.dp -take=ParameterValue -reuse=page \
	    | sddsprocess -pipe -process=Orbit,slope,Slope,functionof=ParameterValue \
	    | sddscollapse -pipe \
	    | sddsconvert -pipe -retain=col,TagName,Slope -rename=col,Slope=$dispColumn \
	    | sddsprocess -pipe=in $tmpRoot.$Plane -redef=col,$dispColumn,$dispColumn,units=m
    }
    array set dispParamArray $dispDirectParamList
    exec sddscombine $tmpRoot.X $tmpRoot.Y -pipe=out \
	| sddsprocess -pipe=in $dispFile -redef=para,DP0,$dispParamArray(dp0) -redef=para,DP1,$dispParamArray(dp1) \
	-redef=para,Steps,$dispParamArray(Steps),type=long
    file delete $tmpRoot.X $tmpRoot.Y $tmpRoot.dp
}

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

proc MakeDispersionParamFile {args} {
    global dispDirectParamList
    APSParseArguments {paramFile closedOrbit}
    array set dispParamArray $dispDirectParamList
    set parameterValueList ""
    set elementNameList ""
    set elementParameterList ""
    set rowLimit 1
    for {set i 0} {$i < $dispParamArray(Steps)} {incr i} {
	set dp [expr $dispParamArray(dp0) + ($dispParamArray(dp1) - $dispParamArray(dp0))/($dispParamArray(Steps)-1)*$i]
	lappend parameterValueList $dp
	lappend elementNameList $dispParamArray(malignElementName)
	lappend elementParameterList DP
	if $closedOrbit {
	    lappend parameterValueList 0
	    lappend elementNameList $dispParamArray(malignElementName)
	    lappend elementParameterList ON_PASS
	    set rowLimit 2
	}
    }
    exec sddsmakedataset -pipe=out -col=ParameterValue,type=double -data=[join $parameterValueList ,] \
	-col=TagName,type=string -data=[join $elementNameList ,] -col=ElementParameter,type=string -data=[join $elementParameterList ,] \
	| sddsprocess -pipe -print=col,ParameterMode,differential -redef=para,UseDoubleOrbits,0,type=long \
	-redef=col,Index,i_row,type=long -print=col,ColumnName,DISP%02d,Index \
	| sddsbreak -pipe -rowlimit=$rowLimit \
	| sddsprocess -pipe -process=ColumnName,first,ColumnName \
	| sddsconvert -pipe=in $paramFile -del=col,Index,ColumnName
    if [catch {Fit_MakeElementNameFromTagName -input $paramFile -addElementOccurence 1} result] {
	return -code error "Fit_MakeElementNameFromTagName : $result"
    }
}

#----------------------------------------------------------------------------------------------------------------------------
proc AddBpmNoise2Orbit {args} {
    global tmpDir bpmNoise bpmNoiseSeed
    APSParseArguments {orbitFile}
    set gLimit 3
    set tmpRoot $tmpDir/[APSTmpString]-bpmNoise
    if {$bpmNoise > 0} {
	if ![string length $bpmNoiseSeed] {set bpmNoiseSeed [exec rpnl "rnd 1000 * 1 + int"]}
	set colList [join [exec sddsquery $orbitFile -col] ]
	if {[lsearch $colList Cx] != -1 && [lsearch $colList Cy] != -1} {
	    foreach {colX colY} {Cx Cy} {}
	} else {
	    foreach {colX colY} {x y} {}
	}
	exec sddsprocess $orbitFile $tmpRoot "-rpnexpression=$bpmNoiseSeed srnd" \
	    "-redef=col,$colX,$colX $gLimit grndl $bpmNoise * +" "-redef=col,$colY,$colY $gLimit grndl $bpmNoise * +"
	incr bpmNoiseSeed 131
	file copy -force $tmpRoot $orbitFile
	file delete $tmpRoot
    }
}
#----------------------------------------------------------------------------------------------------------------------------
# Every page needs to be copied useDoubleOrbits times with varied strength of all kicks in the page to make a linear scan:
proc ExpandKickFilePages {args} {
    global tmpDir
    APSParseArguments {correctorFile useDoubleOrbits}
    set tmpRoot $tmpDir/[APSTmpString]-expPages
    set nKicks [exec sdds2stream $correctorFile -npages=bare]
    set deleteFiles ""
    for {set i 0} {$i < $nKicks} {incr i} {
	exec sddsconvert $correctorFile $tmpRoot.$i -keep=[expr $i+1]
	lappend deleteFiles $tmpRoot.$i
	for {set j 0} {$j < $useDoubleOrbits} {incr j} {
	    set factor [expr -1.0 + 2*$j/($useDoubleOrbits-1)]
	    exec sddsprocess $tmpRoot.$i $tmpRoot.$i.$j "-redef=col,ParameterValue,ParameterValue $factor *" -redef=para,OrbitScan,$j,type=long
	    lappend fileList $tmpRoot.$i.$j
	    lappend deleteFiles $tmpRoot.$i.$j
	}
    }
    eval exec sddscombine $fileList -pipe=out \
	| sddsconvert -pipe -retain=para,ColumnName,OrbitScan \
	| sddsprocess -pipe=in $correctorFile -redef=para,UseDoubleOrbits,$useDoubleOrbits,type=long
    eval file delete $deleteFiles
}
#----------------------------------------------------------------------------------------------------------------------------
# Every page needs to be copied averageOrbits times with different noise seeds
proc MultiplyKickFilePages {args} {
    global tmpDir
    APSParseArguments {correctorFile averageOrbits}
    set tmpRoot $tmpDir/[APSTmpString]-cpPages
    exec sddsprocess $correctorFile $tmpRoot -def=para,PageIndex,i_page
    set fileList ""; for {set i 0} {$i < $averageOrbits} {incr i} {lappend fileList $tmpRoot}
    eval exec sddscombine $fileList -pipe=out \
	| sddssort -pipe -para=PageIndex \
	| sddsprocess -pipe -redef=para,AverageOrbits,$averageOrbits,type=long \
	| sddsconvert -pipe=in $tmpRoot.out -del=para,PageIndex,Filename
    file copy -force $tmpRoot.out $correctorFile
    file delete $tmpRoot $tmpRoot.out
}
#----------------------------------------------------------------------------------------------------------------------------
# additionalErrorSources should have form "<TagName> <ElementParameter> <NoiseType> <NoiseAmplitude> ..."
# NoiseType: 0 for uniform with +-NoiseAmplitude; 1 for Gaussian with sigma NoiseAmplitude
# For example: "MALIN#1 DP 1 1e-4 B2C9ES2C#1 KICK 0 1e-4 B2C9ES2C#1 KICK 0 1e-4 S39IS1#1 FSE 1 1e-4 S39IS2#1 FSE 1 1e-4"
proc AddAdditionalErrorSources {args} {
    global tmpDir additionalErrorSourcesSeed
    APSParseArguments {additionalErrorSources correctorFile}
    set tmpRoot $tmpDir/[APSTmpString]-addErr
    set nOrbits [exec sdds2stream $correctorFile -npages=bare]
    set gLimit 2
    set sourceParams 4
    foreach {nameList paramNameList noiseTypeList noiseAmplList} {"" "" "" ""} {}
    set N [expr [llength $additionalErrorSources] / $sourceParams]
    for {set i 0} {$i < $N} {incr i} {
	lappend nameList [lindex $additionalErrorSources [expr $i * $sourceParams]]
	lappend paramNameList [lindex $additionalErrorSources [expr $i * $sourceParams + 1]]
	lappend noiseTypeList [lindex $additionalErrorSources [expr $i * $sourceParams + 2]]
	lappend noiseAmplList [lindex $additionalErrorSources [expr $i * $sourceParams + 3]]
    }
    exec sddsmakedataset $tmpRoot.0 \
	-col=TagName,type=string -data=[join $nameList ,] \
	-col=NoiseType,type=long -data=[join $noiseTypeList ,] \
	-col=NoiseAmpl,type=double -data=[join $noiseAmplList ,] \
	-col=ElementParameter,type=string -data=[join $paramNameList ,] 
    set fileList ""; for {set i 0} {$i < $nOrbits} {incr i} {lappend fileList $tmpRoot.0}
    eval exec sddscombine $fileList -pipe=out \
	| sddsprocess -pipe \"-rpnexpression=$additionalErrorSourcesSeed srnd\" \
	\"-redef=col,ParameterValue,NoiseType 0 == ? rnd 2 * 1 - NoiseAmpl * : $gLimit grndl NoiseAmpl * \$\" \
	-print=col,ParameterMode,differential -redef=col,PageNumber,i_page \"-redef=col,PageNumber2,i_page 2 *\" -redef=col,RowNumber,i_row \
	| sddsconvert -pipe=in $tmpRoot.1 -retain=col,TagName,ElementParameter,ParameterMode,ParameterValue,PageNumber*,RowNumber -del=para,*
    exec sddsprocess $correctorFile -pipe=out -redef=col,PageNumber,i_page "-redef=col,PageNumber2,i_page 2 * 1 -" -redef=col,RowNumber,i_row \
	| sddsconvert -pipe=in $tmpRoot.2 -retain=col,TagName,ElementParameter,ParameterMode,ParameterValue,PageNumber*,RowNumber -del=para,*
	
    exec sddscombine $tmpRoot.2 $tmpRoot.1 -pipe=out -merge \
	| sddssort -pipe -col=PageNumber2 -col=PageNumber -col=RowNumber \
	| sddsbreak -pipe -change=PageNumber \
	| sddsxref -pipe $correctorFile -leave=* -transfer=para,ColumnName,UseDoubleOrbits,AverageOrbits \
	| sddsconvert -pipe=in $tmpRoot.3 -del=col,PageNumber*,RowNumber
    file copy -force $tmpRoot.3 $correctorFile
    file delete $tmpRoot.0 $tmpRoot.1 $tmpRoot.2 $tmpRoot.3
    if [catch {Fit_MakeElementNameFromTagName -input $correctorFile -addElementOccurence 1} result] {
	return -code error "Fit_MakeElementNameFromTagName : $result"
    }
}

#----------------------------------------------------------------------------------------------------------------------------
proc MakeKickFile { args } {
    APSParseArguments { correctorFile corrList kickList useDoubleOrbits}
    exec sddsmakedataset -pipe=out \
	-col=TagName,type=string -data=[join $corrList ,] \
	-col=ParameterValue,type=double -data=[join $kickList ,] \
	| sddsprocess -pipe -print=col,ElementParameter,KICK -print=col,ParameterMode,differential \
	| sddsbreak -pipe -rowLimit=1 \
	| sddsprocess -pipe=in $correctorFile -process=TagName,first,ColumnName -redef=para,UseDoubleOrbits,$useDoubleOrbits,type=long
    if [catch {Fit_MakeElementNameFromTagName -input $correctorFile -addElementOccurence 1} result] {
	return -code error "Fit_MakeElementNameFromTagName : $result"
    }
}

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

proc AddBpmNoise2RM {args} {
    global tmpDir kickX kickY useDoubleOrbits dispColumn
    APSParseArguments {hrmFile vrmFile bpmNoise seed dispNoise}
    set tmpRoot $tmpDir/[APSTmpString]-addBpmNoise
    if ![string length $seed] {set seed [exec rpnl "rnd 1000 * int"]}
    set gLimit 2
    foreach ext [list hrm vrm] kick [list $kickX $kickY] {
	#------ Scale bpm noise:
	#------ Output file of response calculations is (X1-X2)/kick (/2 for useDoubleOrbits)
	if $useDoubleOrbits {
	    set bpmNoiseRM [expr $bpmNoise * 1.414 / $kick / 2]
	} else {
	    set bpmNoiseRM [expr $bpmNoise * 1.414 / $kick]
	}

	#--- Command below is built this way so that if there is column "s" in the rmFile, it is preserved in the end
	set colList [join [exec sddsquery [set ${ext}File] -col] ]
	if {[lsearch $colList $dispColumn] != -1 } {
	    if [catch {exec sddsprocess [set ${ext}File] -pipe=out "-rpnexpression=$seed srnd" "-redef=col,%s,$gLimit grndl $bpmNoiseRM *,select=*,exclude=*Name" \
			   "-redef=col,$dispColumn,$dispColumn $bpmNoiseRM / $dispNoise *" \
			   | sddsmatrixop -pipe=in $tmpRoot.rm -push=[set ${ext}File] -add} result] {
		return -code error "Error adding BPM noise ([set ${ext}File]): $result"
	    }
	} else {
	    if [catch {exec sddsprocess [set ${ext}File] -pipe=out "-rpnexpression=$seed srnd" "-redef=col,%s,$gLimit grndl $bpmNoiseRM *,select=*,exclude=*Name" \
			   | sddsmatrixop -pipe=in $tmpRoot.rm -push=[set ${ext}File] -add} result] {
		return -code error "Error adding BPM noise ([set ${ext}File]): $result"
	    }
	}
	exec sddsconvert [set ${ext}File] -pipe=out -retain=col,*Name,s \
	    | sddsxref -pipe=in $tmpRoot.rm $tmpRoot.rm1 -take=* -leave=*Name,s
	file copy -force $tmpRoot.rm1 [set ${ext}File]
	file delete $tmpRoot.rm $tmpRoot.rm1
    }
}

#----------------------------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------
#--- Main program body ------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------------------------------------------

set dispFile $workDir/dispersion.sdds

switch -exact $runMode {
    getSDDSLattice {
	if [catch {MakeSddsLattice -lteFile $lteFile -beamlineName $beamlineName -latticeParamFileList $latticeParamFileList \
		       -workDir $workDir -sddsLatticeFile $sddsLatticeFile -accelCodeBinDir $accelCodeBinDir} result] {
	    puts stderr "MakeSddsLattice: $result \n Calling argument (LOCO_BINDIR: $env(LOCO_BINDIR)): $rememberArgs"; exit
	}
    }
    getBetaFunctions {
	if [catch {CalculateElegantBetas -lteFile $lteFile -beamlineName $beamlineName \
		       -latticeParamFileList $latticeParamFileList -twiFile $twiFile -workDir $workDir \
		       -accelCodeBinDir $accelCodeBinDir -closedOrbit $closedOrbit \
		       -twissRefFile $twissRefFile -twissRefElement $twissRefElement \
		       -fixedOrbitLength $fixedOrbitLength -elemByElemSR $elemByElemSR} result] {
	    puts stderr "CalculateElegantBetas: $result \n Calling argument (LOCO_BINDIR: $env(LOCO_BINDIR)): $rememberArgs"; exit
	}
    }
    getResponseMatrix {
	if {$useDoubleOrbits == 1} {set useDoubleOrbits 2}
	if [string length $twiFile] {
	    #------ If not closed orbit, dispersion has to be calculated starting from zero. It requires special attention.
	    #------ Therefore use separate procedure. If not, twiss file is calculated as part of RM for speed.
	    if {!$closedOrbit || $directCalc == 1} {
		if [catch {CalculateElegantBetas -lteFile $lteFile -beamlineName $beamlineName \
			       -latticeParamFileList $latticeParamFileList -twiFile $twiFile -workDir $workDir \
			       -accelCodeBinDir $accelCodeBinDir -closedOrbit $closedOrbit \
			       -twissRefFile $twissRefFile -twissRefElement $twissRefElement \
			       -fixedOrbitLength $fixedOrbitLength -elemByElemSR $elemByElemSR} result] {
		    puts stderr "CalculateElegantBetas: $result \n Calling argument (LOCO_BINDIR: $env(LOCO_BINDIR)): $rememberArgs"; exit
		}
		set includeTwissCalc 0
	    } else {
		set includeTwissCalc 1
	    }
	}
	if [catch {CalculateElegantRM -lteFile $lteFile -beamlineName $beamlineName \
		       -latticeParamFileList $latticeParamFileList -hrmFile $hrmFile -vrmFile $vrmFile \
		       -hCorrList $hCorrList -vCorrList $vCorrList -directCalc $directCalc \
		       -fixedOrbitLength $fixedOrbitLength -kickX $kickX -kickY $kickY \
		       -workDir $workDir -twiFile $twiFile -coupledMatrix $coupledMatrix \
		       -fitDispersion $fitDispersion -dispFitWeight $dispFitWeight -dispColumn $dispColumn \
		       -useQsub $useQsub -numberOfQsubTasks $numberOfQsubTasks \
		       -waitInterval $waitInterval -waitTime $waitTime -qsubCommand $qsubCommand \
		       -accelCodeBinDir $accelCodeBinDir -closedOrbit $closedOrbit \
		       -twissRefFile $twissRefFile -twissRefElement $twissRefElement -malignElement $malignElement \
		       -useKickFiles $useKickFiles -kickFileX $kickFileX -kickFileY $kickFileY \
		       -includeTwissCalc $includeTwissCalc -elemByElemSR $elemByElemSR \
		       -useDoubleOrbits $useDoubleOrbits -verbose $verbose \
		       -bpmGainXFile $bpmGainXFile -bpmGainYFile $bpmGainYFile \
		       -bpmTiltFile $bpmTiltFile -bpmCoefFile $bpmCoefFile \
		       -bpmNoise $bpmNoise -bpmNoiseSeed $bpmNoiseSeed -dispNoise $dispNoise \
		       -averageOrbits $averageOrbits} result] { 
	    puts stderr "CalculateElegantRM: $result \n Calling argument (LOCO_BINDIR: $env(LOCO_BINDIR)): $rememberArgs"; exit
	}
    
	#--- Transfering tunes from twiss file to response files:
	if [string length $twiFile] {
	    if [llength $hCorrList] {
		if [catch {exec sddsxref $hrmFile $twiFile $hrmFile.tmp -nowarning -transfer=para,nux,nuy -leave=*} result] {
		    return -code error "Error (hrmFile $hrmFile): $result"
		}
		file copy -force $hrmFile.tmp $hrmFile
		file delete $hrmFile.tmp
	    }
	    if [llength $vCorrList] {
		exec sddsxref $vrmFile $twiFile $vrmFile.tmp -nowarning -transfer=para,nux,nuy -leave=*
		file copy -force $vrmFile.tmp $vrmFile
		file delete $vrmFile.tmp
	    }
	}
	if [string length $doneFile] {exec touch $doneFile}
    }
    runCalculateOrbits {
	if [catch {CalculateOrbits -latticeFile $latticeFile -orbitFile $orbitFile -verbose $verbose \
		       -workDir $workDir -tmpDir $tmpDir -nonlin $nonlin -correctorFile $correctorFile \
		       -corrList $corrList -kick $kick -sddsLatticeFile $sddsLatticeFile \
		       -bRho $bRho -accelCodeBinDir $accelCodeBinDir -closedOrbit $closedOrbit -refCorr $refCorr \
		       -fixedOrbitLength $fixedOrbitLength -useDoubleOrbits $useDoubleOrbits -useKickFiles $useKickFiles \
		       -averageOrbits $averageOrbits} result] {
	    return -code error "CalculateOrbits: $result"
	}
	exec echo END > $workDir/task.completed
    }
}
exit
