]]\n\
[-reverseOrder] [-localFit]\n\n\
textOutput Make the output file a text file instead of an SDDS file.\n\
totalCharge Specify the total charge in Coulombs.\n\
chargeParameter Specify the name of a parameter in the input file that gives\n\
the total charge in Coulombs.\n\
wavelength Specify the wavelength of light in meters.\n\
slices Specify the number of slices to cut the beam into.\n\
steer Force the x, x', y, and y' centroids for the whole beam to zero.\n\
Slices may still have nonzero centroids.\n\
removePTail Removes the momentum tails from the beam. deltaLimit is the maximum\n\
|p-|/
that will be kept. If 'fit' is given, then a linear fit to\n\
(t, p) is done and removal is based on residuals from this fit.\n\
If 'beamOutput' is given, the filtered beam is written to the named\n\
file for review.\n\
reverseOrder By default, the data for the head of the beam comes first. This\n\
option causes elegant to put the data for the tail of the beam \n\
first.\n\
localFit If given, then for each slice the program performs a local linear\n\
fit of momentum vs time, and subtracts it from the momentum data.\n\
This removes a contribution to the energy spread, but may not be\n\
appropriate if slices are longer than the cooperation length.\n\n\
Program by Robert Soliday. This is version 2, August 2008, M. Borland. ("__DATE__")\n";
#define PTAILS_DELTALIMIT 0x0001U
#define PTAILS_FITMODE 0x0002U
#define PTAILS_OUTPUT 0X0004U
long removeMomentumTails(double *x, double *xp, double *y,
double *yp, double *s, double *p, double *t,
long rows, double limit, unsigned long flags);
void removeLocalFit(double *p, double *s, short *selected, double pAverage, long rows);
/* ********** */
int main(int argc, char **argv)
{
SDDS_DATASET SDDS_input, SDDS_output, SDDS_pTails;
FILE *fileID=NULL;
SCANNED_ARG *s_arg;
long i, j, i_arg;
char *input, *output;
unsigned long pipeFlags=0;
long noWarnings=0, tmpfile_used=0;
long retval, origRows, rows, sddsOutput=1, steer=0, firstOne, sliceOffset=0;
long reverseOrder=0;
char *chargeParameter;
unsigned long removePTailsFlags = 0;
double pTailsDeltaLimit, term;
char *pTailsOutputFile = NULL;
long localFit = 0;
short *selected;
char column[7][15] = {"x", "xp", "y", "yp", "t", "p", "particleID"};
long columnIndex[7];
double *tValues, *sValues, *xValues, *xpValues, *yValues, *ypValues, *pValues;
double sAverage, pAverage, pStDev;
double xAverage, xpAverage, xRMS, x2Average, xp2Average, xxpAverage, xEmittance;
double yAverage, ypAverage, yRMS, y2Average, yp2Average, yypAverage, yEmittance;
double sMin=1, sMax=-1, s1, s2;
/*
double pMin, pMax;
*/
double alphax, alphay;
double wavelength=0.0001;
double totalCharge=0;
double current;
long slices=4;
long tmp=0;
long N=0;
input = output = NULL;
sValues = NULL;
chargeParameter = NULL;
selected = NULL;
SDDS_RegisterProgramName(argv[0]);
argc = scanargs(&s_arg, argc, argv);
if (argc<3)
bomb(NULL, USAGE);
/* parse the command line */
for (i_arg=1; i_arg=1) {
rows = origRows = SDDS_RowCount(&SDDS_input);
if (rows < 1) {
fprintf(stderr, "error: no rows found for page %ld\n", retval);
exit(1);
}
/* get all the needed column data */
if (!(tValues = SDDS_GetNumericColumn(&SDDS_input, "t", SDDS_DOUBLE)) ||
!(xValues = SDDS_GetNumericColumn(&SDDS_input, "x", SDDS_DOUBLE)) ||
!(xpValues = SDDS_GetNumericColumn(&SDDS_input, "xp", SDDS_DOUBLE)) ||
!(yValues = SDDS_GetNumericColumn(&SDDS_input, "y", SDDS_DOUBLE)) ||
!(ypValues = SDDS_GetNumericColumn(&SDDS_input, "yp", SDDS_DOUBLE)) ||
!(pValues = SDDS_GetNumericColumn(&SDDS_input, "p", SDDS_DOUBLE))) {
SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
exit(1);
}
if (chargeParameter &&
!SDDS_GetParameterAsDouble(&SDDS_input, chargeParameter, &totalCharge))
SDDS_Bomb("unable to read charge from file");
if ((!(tValues)) || (!(xValues)) || (!(xpValues)) || (!(yValues)) || (!(ypValues)) || (!(pValues))) {
fprintf(stderr, "error: invalid data for page %ld\n", retval);
exit(1);
}
selected = realloc(selected, sizeof(*selected)*rows);
sValues = realloc(sValues, sizeof(*sValues)*rows);
sAverage = xAverage = xpAverage = yAverage = ypAverage = 0;
for (i=0;i sValues[j])) {
selected[j] = 1;
N++;
xAverage += xValues[j];
xpAverage += xpValues[j];
yAverage += yValues[j];
ypAverage += ypValues[j];
pAverage += pValues[j];
}
}
if (N>2) {
current = N * 2.99792458e8 * (totalCharge/(origRows * wavelength));
xAverage /= N;
yAverage /= N;
xpAverage /= N;
ypAverage /= N;
pAverage /= N;
/*
pMin = pMax = pAverage;
*/
pStDev = 0;
if (localFit)
removeLocalFit(pValues, sValues, selected, pAverage, rows);
for (j=0;j pMax)
pMax = pValues[j];
*/
}
}
pStDev = sqrt(pStDev/(N-1.0));
x2Average /= N;
y2Average /= N;
xp2Average /= N;
yp2Average /= N;
xxpAverage /= N;
yypAverage /= N;
xRMS = sqrt(x2Average);
yRMS = sqrt(y2Average);
if ((term = (x2Average * xp2Average) - sqr(xxpAverage))>0)
xEmittance = sqrt(term)*pAverage;
else
xEmittance = 0;
if ((term = (y2Average * yp2Average) - sqr(yypAverage))>0)
yEmittance = sqrt(term)*pAverage;
else
yEmittance = 0;
alphax = -xxpAverage/(xEmittance/pAverage);
alphay = -yypAverage/(yEmittance/pAverage);
} else {
pAverage = pStDev = xEmittance = yEmittance = xRMS = yRMS = xAverage =
yAverage = xpAverage = ypAverage = alphax = alphay = current = 0;
}
if (sddsOutput) {
if (SDDS_SetRowValues(&SDDS_output, SDDS_SET_BY_NAME|SDDS_PASS_BY_VALUE, i+sliceOffset,
"s", s1,
"t", -s1/c_mks,
"gamma", pAverage,
"dgamma", pStDev,
/* This was commented out to be consistant with GENESIS
where dgamma is the rms value of gamma
"dgamma", pMax - pMin,
"gammaStDev", pStDev,
*/
"xemit", xEmittance,
"yemit", yEmittance,
"xrms", xRMS,
"yrms", yRMS,
"xavg", xAverage,
"yavg", yAverage,
"pxavg", xpAverage,
"pyavg", ypAverage,
"alphax", alphax,
"alphay", alphay,
"current", current,
"wakez", 0.0,
"N", N,
NULL) != 1) {
SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
exit(1);
}
} else {
fprintf(fileID, "\n % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E % .6E 0.000000E+00", \
s1, pAverage, pStDev, xEmittance, yEmittance, \
xRMS, yRMS, xAverage, yAverage, xpAverage, ypAverage, alphax, alphay, current);
}
}
free(tValues);
free(xValues);
free(xpValues);
free(yValues);
free(ypValues);
free(pValues);
}
free(sValues);
/* close the output file */
if (sddsOutput) {
if (!SDDS_WritePage(&SDDS_output)) {
SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
exit(1);
}
if (!SDDS_Terminate(&SDDS_output)) {
SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
exit(1);
}
} else {
fclose(fileID);
}
/* close the input file */
if (!SDDS_Terminate(&SDDS_input)) {
SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
exit(1);
}
if (removePTailsFlags&PTAILS_OUTPUT &&
!SDDS_Terminate(&SDDS_pTails)) {
SDDS_PrintErrors(stderr, SDDS_VERBOSE_PrintErrors);
exit(1);
}
return(0);
}
long removeMomentumTails(double *x, double *xp, double *y,
double *yp, double *s, double *p, double *t,
long rows, double limit, unsigned long flags)
/* t is equivalent to s. It is filtered here for convenience in output
* to an elegant-type file
*/
{
double *delta, pAve;
long ip, jp;
#if defined(DEBUG)
static FILE *fp = NULL;
#endif
if (!rows)
return rows;
if (!(delta = malloc(sizeof(*delta)*rows)))
SDDS_Bomb("memory allocation failure");
for (ip=pAve=0; iplimit) {
if (jp>ip) {
x[ip] = x[jp];
xp[ip] = xp[jp];
y[ip] = y[jp];
yp[ip] = yp[jp];
s[ip] = s[jp];
p[ip] = p[jp];
t[ip] = t[jp];
delta[ip] = delta[jp];
jp --;
ip --;
}
rows --;
}
}
free(delta);
#if defined(DEBUG)
if (!fp) {
fp = fopen("e2g.sdds", "w");
fprintf(fp, "SDDS1\n&column name=x type=double &end\n");
fprintf(fp, "&column name=xp type=double &end\n");
fprintf(fp, "&column name=y type=double &end\n");
fprintf(fp, "&column name=yp type=double &end\n");
fprintf(fp, "&column name=s type=double &end\n");
fprintf(fp, "&column name=p type=double &end\n");
fprintf(fp, "&data mode=ascii &end\n");
fprintf(fp, "%ld\n", rows);
}
for (ip=0; ip