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

if {![info exists env(OAG_TOP_DIR)]} { set env(OAG_TOP_DIR) /usr/local }
set auto_path [linsert $auto_path 0  $env(OAG_TOP_DIR)/oag/apps/lib/$env(HOST_ARCH)]

APSStandardSetup

# Script to repeatedly double the number of particles in a distribution.
# Algorithm:
# 1. For each doubling, insert a new particle "near" every pair of existing particles in time.
#    The particle has a new t value, but the same (x, xp, y, yp, p) as one of the original particles.
# 2. Bin the beam according to t into a large number of bins.  Randomize the assignment of p values
#    across particles in the same bin.
#
# usage: doubleDist6 -input <name> -output <name> -doublings <number> -nt <bins>
# -doublings <number> : give the number of times to double the number of particles.
#                       -doublings 3 will multiply the number of particles by 2^3
# -nt <bins>          : give the number of time bins to use.  Particle's total momenta
#                       are randomly swapped within the same time bin to reduce the
#                       number of particles with (x, xp, y, yp, p) all identical
#
# M. Borland, 2009.

set usage {usage: doubleDist6 -input <name> -output <name> -doublings <number> -nt <bins>}
set args $argv
set output ""
set input ""
set nt 1000
set doublings 1
if {[APSStrictParseArguments {input output doublings nt}] \
      || ![string length $input] || ![string length $output] || $doublings<1 || $nt<100} {
    return -code error "$usage"
}

if ![file exists $input] {
    return -code error "$input: not found"
}
if [file exists $output] {
    return -code error "$output: already exists"
}

for {set d 0} {$d<$doublings} {incr d} {
    puts stderr "Doubling pass $d"
    exec sddscombine $input $input -pipe=out \
      | sddsprocess -pipe -redefine=column,page,i_page,type=short \
      | sddscombine -pipe -merge \
      | sddssort -pipe -column=t \
      | sddsprocess -pipe=in $output.$d \
      "-redefine=parameter,n_rows2,n_rows 2 -,type=long" \
      "-redefine=column,doAve,i_row n_rows2 < page 2 == && ? 1 : 0 $ ,type=short" \
      "-redefine=column,t,doAve 1 == ? t i_row 1 + &t \[ + 2 / t i_row 1 + &t \[ - rnd 2 * 1 - * + : t $ " 
    set input $output.$d
    set last $output.$d
    lappend fileList $output.$d
}

file rename $last $output.tm
eval file delete $fileList

puts stderr "Binning by time, randomizing p"
exec sddsprocess $output.tm -pipe=out \
    -process=t,min,tMin -process=t,spread,tSpread \
    "-define=column,it,t tMin - tSpread / $nt *,type=long" \
    | sddssort -pipe -column=it \
    | sddsbreak -pipe -change=it \
    | sddsprocess -pipe -define=column,RN,rnd -process=p,spread,pSpread \
    "-redefine=column,p,rnd 2 * 1 - pSpread * 30 / p +" \
    | sddssort -pipe -column=RN \
    | sddscombine -pipe=in -merge $output.p -retain=column,p -overwrite

puts stderr "Inserting new p values"
exec sddsxref $output.p $output.tm $output -take=*

file delete $output.p $output.tm $output.tb



