Project

General

Profile

fast_ms_emit.py

Python Source Code - Chip Edstrom, 02/14/2018 05:33 PM

 
1
import matplotlib.pyplot as plt
2
from math import pi, isnan
3
from time import localtime, strftime, time as now
4
from numpy import arange, exp, abs, delete, float, min, max, array, inf, random, std, median
5
from scipy.optimize import leastsq as myFit
6
from scipy.integrate import quad as integrate
7
from subprocess import call
8
import tarfile
9
import commands
10
import sys
11
import os
12

    
13
t0 = now()
14
filedate = str( t0 )
15

    
16
# EXPERIMENT SETUP PARAMS
17
myW = 40e-6                                        # Slit Width [m] (d in paper)
18
myS = 400e-6                                        # Slit Spacing (between slits) [m] (w in paper)
19
myX107 = 9221.9455                                # Survey distance along beamline to x107 [mm] (Slit possition)
20
myX108 = 9676.6585                                # Survey distance along beamline to x108 [mm]
21
myX109 = 10058.638                                # Survey distance along beamline to x109 [mm]
22
myX111 = 10759.953                                # Survey distance along beamline to x111 [mm]
23
myX115 = 14000.000                                # Distance along beamline to x115 [mm] (Slit possition) !!! not from survey and adds distance through chicane !!!
24
myX120 = 16600.485                                # Survey distance along beamline to x120 [mm]
25
myX124 = 19147.663                                # Survey distance along beamline to x124 [mm]
26
l0 = 107
27
l1 = 111
28
z0 = myX107
29
z1 = myX111
30
myL = (z1 - z0)/1000 # 0.445                        # Drift length between slit face and screen [m]
31
x_res = 9.3e-6                                        # X-Axis Resolution [m/px], (0.008 / 366) = 2.8758e-5 for initial sample
32

    
33
# CAMERA PARAMS
34
myCam = 4
35
myPlane = "X"
36
myAver = 5
37
MINIMUM_SLITS = 3
38

    
39
# TEST PARAMS
40
DEBUG = False
41
READ_LIVE = True
42
RANDOMIZE_NOISE = False
43
mySample = "/export/home1/edstrom/acl/emit_calc/sample.dat"
44
test_noise = 0.01
45
test_level = 0.0
46

    
47
# Paths
48
OUTPUT_DIRECTORY = "/usr/local/userb/fast/sync/"# Output Directory to be used
49
ARCHIVE_DIRECTORY = "/usr/local/cbs_files/cns_write/acl/data_files/fast_ms_emit/" # Archiving directory for tars
50
myHistogram = ARCHIVE_DIRECTORY + "fast_ms_emit.hist"
51
ACL = "/usr/local/xcons/bin/acl"
52
myHistFetch = "/usr/local/cbs_files/sequencer/acl/fast_gethist.acl"
53

    
54
# PROCESS PARAMS
55
smooth_range = 40                                # Smoothing Range for Baseline Fit
56
pick_smooth = 4                                        # Smoothing factor for picking 
57
min_peak_width = 20                                # Minimum peak width for picking
58
min_threshold =  0.002                                # Threshold over the standard background fit to start peak detection (0.1 = 10%)
59
n_pts = 0                                        # Total Number of Points
60

    
61
# CONTROL FLAGS
62
SUBTRACT_BACKGROUND = 1                                # 0 = no subtraction; 1 = fit subtraction; 2 = median subtraction
63
CENTER_ON_ENVELOPE = True                        # Center all data around the envelope (rather than the central/highest peak)
64
VERBOSE = True                                        # Extra verbose output
65
FRAME_OUTPUT = False                                # Only outputs the detected gaussians plus 10% (at most)
66
FULL_FIT = True                                        # Initially fit for the full ensemble of gaussians
67
FULL_OUTPUT = True                                # Enable output to bd-www
68

    
69
y_in,x_in,y,x = [],[],[],[]                        # Histogram functions, each with size of n_pts
70
a = []                                                # Full fit parameter array, will be (n_bumps + 1) x 4 for the baseline and all other gaussians
71
p_full = []                                        # A 1-D list of the a[] parameterization for the fit routine
72
                                                # Individual sets of gaussian fit parameters:
73
                                                #        0 = Amplitude
74
                                                #        1 = Centroid
75
                                                #        2 = Sigma
76
                                                #        3 = Linear offset
77
bump= []                                        # Individual bump information
78
i,j,k = 0,0,0                                        # Loop iterators
79
myArg = []                                        # Command Line Argumments
80

    
81
########################################################### COMMAND LINE ARGUMENTS
82
print "Welcome to the FAST Multi-slit emittance scanner."
83

    
84
if (len(sys.argv) > 1):
85
        myArg = sys.argv[1].split(',')                        # Read in data
86
nArgs = len(myArg)
87
if (nArgs < 1):
88
        print "No command-line arguments passed..."
89
        print "> Syntax: python fast_ms_emit.py arg1=val1,arg2=val2..."
90
        print "> Valid arguments: <FILENAME>,cam=(1-4),plane(X/Y),aver=(int),"
91
        print ">          z0=(107/115),z1=(108/109/111/115/120/124),res=(float),"
92
        print ">          thresh=(float),minwidth=(int),bgsubtract=(0/1/2),"
93
        print ">          min_peaks=(int),plot=(0/1)"
94
        print
95
        print "Proceeding with default arugements..."
96
else:
97
        print "Parsing command-line arguments..."
98
i=0
99
while (i < nArgs):
100
        myArg[i] = myArg[i].split('=')
101
        if(myArg[i][0] == 'cam'):
102
                if((int(myArg[i][1]) > 0) and (int(myArg[i][1]) < 5)):
103
                        print "> Using DC" + str(int(myArg[i][1])) + " for histogram information."
104
                        myCam=int(myArg[i][1])
105
                else:
106
                        print "> !! Invalid camera selection: [" + myArg[i][1] + "]"
107
        if(myArg[i][0] == 'res'):
108
                if(float(myArg[i][1]) > 0):
109
                        print "> Using " + str(float(myArg[i][1])) + " as the camera resolution."
110
                        x_res=float(myArg[i][1])
111
                else:
112
                        print "> !! Invalid resolution: [" + myArg[i][1] + "]"
113
        if(myArg[i][0] == 'z0'):
114
                if (int(myArg[i][1])  == 107):
115
                        l0 = 107
116
                        z0 = myX107
117
                if (int(myArg[i][1])  == 115):
118
                        l0 = 115
119
                        z0 = myX115
120
                print "> Slits are presumed to be inserted at x" + str(l0) + "."
121
        if(myArg[i][0] == 'z1'):
122
                if (int(myArg[i][1])  == 108):
123
                        l1 = 108
124
                        z1 = myX108
125
                if (int(myArg[i][1])  == 109):
126
                        l1 = 109
127
                        z1 = myX109
128
                if (int(myArg[i][1])  == 111):
129
                        l1 = 111
130
                        z1 = myX111
131
                if (int(myArg[i][1])  == 115):
132
                        l1 = 115
133
                        z1 = myX115
134
                if (int(myArg[i][1])  == 120):
135
                        l1 = 120
136
                        z1 = myX120
137
                if (int(myArg[i][1])  == 124):
138
                        l1 = 124
139
                        z1 = myX124
140
                print "> YAG screen and camera are presumed to be active at x" + str(l0) + "."
141
        if(myArg[i][0] == 'min_peaks'):
142
                if (int(myArg[i][1])  > 2 ):
143
                        print "> Using " + str(int(myArg[i][1])) + " as the minimum number of peaks to accept for fitting."
144
                        MINIMUM_SLITS=int(myArg[i][1])
145
                else:
146
                        print "> !! Invalid peak number: [" + myArg[i][1] + "]"
147
        if(myArg[i][0] == 'plot'):
148
                if (int(myArg[i][1]) == 0):
149
                        print "> Not displaying process plots."
150
                        DEBUG = False
151
                else:
152
                        print "> Displaying process plots... close each plot to continue."
153
                        DEBUG = True
154
        if(myArg[i][0] == 'bgsubtract'):
155
                if (int(myArg[i][1]) == 1):
156
                        print "> Subtracting out the fit background."
157
                        SUBTRACT_BACKGROUND = 1
158
                elif (int(myArg[i][1]) == 2):
159
                        print "> Subtracting out the median value (no fit to the background)"
160
                        SUBTRACT_BACKGROUND = 2
161
                else:
162
                        print "> Not subtracting out the measured background."
163
                        SUBTRACT_BACKGROUND = 0
164
        if(myArg[i][0] == 'thresh'):
165
                if(float(myArg[i][1]) > 0):
166
                        print "> Using " + str(float(myArg[i][1])) + " as the miminum fit threshold (for peak determination)."
167
                        min_threshold=float(myArg[i][1])
168
                else:
169
                        print "> !! Invalid minimum fit threshold: [" + myArg[i][1] + "]"
170
        if(myArg[i][0] == 'minwidth'):
171
                if(int(myArg[i][1]) > 0):
172
                        print "> Using " + str(int(myArg[i][1])) + " as the miminum peak width [px] (for peak determination)."
173
                        min_peak_width=int(myArg[i][1])
174
                else:
175
                        print "> !! Invalid minimum peak width in px: [" + myArg[i][1] + "]"        
176
        if(myArg[i][0] == 'aver'):
177
                if(int(myArg[i][1]) > -1):
178
                        print "> Averaging " + str(int(myArg[i][1])) + " histogram readbacks for measurement."
179
                        myAver=int(myArg[i][1])
180
                else:
181
                        print "> !! Invalid averaging request: [" + myArg[i][1] + "]"
182
        if(myArg[i][0] == 'plane'):
183
                if ((myArg[i][1] == 'X') or (myArg[i][1] == 'x') or (myArg[i][1] == 'Y') or (myArg[i][1] == 'y')):
184
                        print "> " + myArg[i][1] + "-axis histogram to be used for collection."
185
                        myPlane = myArg[i][1]
186
                else:
187
                        print "> !! Invalid plane selection: [" + myArg[i][1] + "]"
188
        if(myArg[i][0].find(".tar") > 0):
189
                READ_LIVE = False
190
                print "> tar source: " + myArg[i][0]
191
                tar = tarfile.open(myArg[i][0], 'r')
192
                tar.extractall()
193
                tar.close()
194
                mySample = "ms_emit.hist"
195
        if(myArg[i][0].find(".dat") > 0) or (myArg[i][0].find(".hist") > 0):
196
                READ_LIVE = False
197
                print "> data source: " + myArg[i][0]
198
                mySample = myArg[i][0]
199
        i+=1
200
myL = (z1 - z0)/1000
201
print
202
#quit()
203
########################################################### UTILITY FUNCITONS
204
def make_float_array(a_in):
205
        a_out = [0.0] * len(a_in)
206
        i=0
207
        while ( i < len(a_in) ):
208
                a_out[i] = float(a_in[i])
209
                i+=1
210
        return a_out
211

    
212
def make_random_array(a_pts):
213
        a_out = [0.0] * len(a_in)
214
        for i in a_out:
215
                i = random.random()
216
        return a_out
217

    
218
def smooth_array(a_in, sf):
219
        a_out = [0.0] * len(a_in)
220
        i = 0
221
        while ( i < len(a_in)):
222
                j=i-sf
223
                k=i+sf
224
                if (j < 0): j=0
225
                if (k > len(a_in)): k = len(a_in)
226
                a_out[i] = (sum(a_in[j:k])/(k-j))
227
                i+=1
228
        return a_out
229
        
230
def form_array(a_in, n_d):
231
        a_out = []
232
        i=0
233
        while ( i < (len(a_in) / n_d)):
234
                a_part = []
235
                j=0
236
                while ( j < n_d ):
237
                        a_part.append(a_in[(i*n_d)+j])
238
                        j+=1
239
                a_out.append(a_part)
240
                i+=1
241
        return a_out        
242

    
243
def flatten_array(a_in):
244
        a_out = []
245
        for i in a_in:
246
                for j in i: a_out.append(j)
247
        return a_out
248

    
249
def plot_arrays(Y,X):
250
        print sdt() + " -- Close the current plot to continue..."
251
        plt.figure()
252
        i=0
253
        if (len(X) < len(Y)):
254
                while(i < len(Y)):
255
                        plt.plot(X[0],Y[i])
256
                        i=i+1
257
        else:
258
                while(i < len(Y)):
259
                        plt.plot(X[i],Y[i])
260
                        i=i+1
261
        plt.get_current_fig_manager().window.set_position(1)
262
        plt.show()
263
        return 0
264

    
265
def find_array_index(a_in, s):
266
        b_in = arange(0, len(a_in))
267
        i=0
268
        while (i < len(a_in)):
269
                if ( a_in[i] == s ): return b_in[i]
270
                i=i+1
271
        return None
272

    
273
def find_FWHM(a_in, min_pts, thresh):
274
        a_max = find_array_index(a_in, max(a_in))
275
        i = a_max
276
        j = a_max
277
        while ( i > 0 ):
278
                if ( a_in[i] < ((a_in[a_max] + thresh)/2)): break
279
                i=i-1
280
        while ( j < len(a_in) ):
281
                if ( a_in[j] < ((a_in[a_max] + thresh)/2)): break
282
                j=j+1
283
        if ((j-i) < min_pts):
284
                i=0
285
                j=len(a_in)
286
        return [i, j]
287

    
288
def find_peaks(a_in, min_gap, thresh):
289
        a_out=[]
290
        peak_det = False
291
        a_pts = len(a_in)
292
        a_mod = smooth_array(a_in, (pick_smooth))
293
        i, j, k, n, myPeakIndex, myCent, myPeak, maxPeak = 0, 0, 0, 0, 0, 0, 0.0, 0.0
294
        while (i < len(a_mod)):
295
                if ( a_mod[i] > thresh) and (not peak_det):
296
                        if ( DEBUG ):
297
                                print "Begin: " + str(i) + ": " + str(i*x_res) + ", " + str(a_mod[i])
298
                        peak_det = True
299
                        j = i
300
                        i += (min_gap/2)
301
                        if (i >= a_pts): i = a_pts - 1
302
                if ( (a_mod[i] < thresh) or (i == a_pts - 1) ) and (peak_det):
303
                        if ( DEBUG ):
304
                                print "End:   " + str(i) + ":" + str(i*x_res) + ", " + str(a_mod[i])
305
                                print
306
                        peak_det = False
307
                        k = i
308
                        i += (min_gap/2)
309
                        a_eval = a_in[j:k]
310
                        myPeak = max(a_eval)
311
                        myCent = (find_array_index(a_eval, max(a_eval))+j)
312
                        a_out.append( [ myPeak, myCent ] )
313
                        if (myPeak > maxPeak): 
314
                                maxPeak, myPeakIndex = myPeak, n
315
                        n+= 1
316
                i+=1
317
        i=0
318
        while (i < len(a_out)):        # FIND PEAK WINDOWS
319
                myCent = a_out[i][1]
320
                if (len(a_out) <= 1):        # Catch if only a single peak is detected above
321
                        return [0],0 
322
                if (i == 0):
323
                        myWid = (a_out[i+1][1] - myCent)
324
                else:
325
                        myWid = (myCent - a_out[i-1][1])
326
                j = myCent - (myWid/2)
327
                if (j < 0): j=0
328
                k = myCent + (myWid/2)
329
                if (k >= len(a_in)): k = (len(a_in)-1)
330
                #                    Full Width        x0    x1
331
                a_out[i] = a_out[i] + [myWid] + [j] + [k]
332
                i=i+1
333
        return a_out, myPeakIndex
334
        
335
def snow(): return strftime("%H:%M:%S %m-%d-%Y", localtime())
336
def snowfile(): return strftime("%Y-%m-%d_%H-%M-%S", localtime())
337
def sfiletime(): return strftime("%H:%M:%S %m-%d-%Y", localtime(myFileStat.st_ctime))
338
def sdt(): return ('{0: >6.2f}'.format(now() - t0))
339
def f2s(a_in,b):
340
        a_out = '{0: .' + str(b) + 'f}'
341
        return a_out.format(a_in)
342
def e2s(a_in,b):
343
        a_out = '{0: .' + str(b) + 'e}'
344
        return a_out.format(a_in)
345
def i2s(a_in,b):
346
        a_out = '{0: ' + str(b) + 'd}'
347
        return a_out.format(a_in)
348

    
349
def get_histogram(c, myPlane, nAver): # c = camera = 1-4, myPlane = "X" or "Y", nAver >= 1
350
        myACLString = ACL + " " + myHistFetch + " N:DC" + str(c) + myPlane + "H " + str(nAver) + " > " + myHistogram
351
        os.system(myACLString)
352
        return
353

    
354
def set_output(b):
355
        myACLString = ACL + " \"enable settings;"
356
        i = 0
357
        while (i < len(b)):
358
                myACLString += "set N:MSEMIT[" + str(i) + "] "
359
                if (not isnan(b[i])):
360
                        myACLString += str(float(b[i])) + ";"
361
                else:
362
                        myACLString += "-999;"
363
                i+=1
364
        myACLString += "\""
365
        result=commands.getstatusoutput(myACLString)
366
        return result
367

    
368
########################################################### FIT FUNCTIONS
369
def peval_full(q, b):
370
        # This is a full evaluation of the gaussians without a baseline... this ignores the baseline fit in a[n_bumps]
371
        ptot = [0.0]*len(q)
372
        i=0
373
        while (i < n_bumps):
374
                if ( len(b) > n_bumps ):
375
                        p_test = [b[i*4],b[(i*4)+1],b[(i*4)+2],b[(i*4)+3]]         # For if b is a [1 x (4*n_bumps)] array
376
                else:
377
                        p_test = b[i]                                                 # For if b is a [4 x n_bumps] array
378
                ptot = ptot + peval(q, p_test) # Single gaussian test function
379
                i=i+1
380
        return ptot
381

    
382
def pevalb_full(q, b):
383
        # This is a full evaluation of the gaussians including the baseline (in a[n_bumps])
384
        ptot = [0.0]*len(q)
385
        i=0
386
        while (i < n_bumps):
387
                if ( len(b) > n_bumps ):
388
                        p_test = [b[i*4],b[(i*4)+1],b[(i*4)+2],b[(i*4)+3]]        # For if b is a [1 x (4*n_bumps)] array
389
                else:
390
                        p_test = b[i]                                                 # For if b is a [4 x n_bumps] array
391
                ptot = ptot + peval(q, p_test) # Single gaussian test function
392
                i=i+1
393
        if ( len(b) > n_bumps ):
394
                p_test = [b[i*4],b[(i*4)+1],b[(i*4)+2],b[(i*4)+3]]
395
        else:
396
                p_test = b[i]
397
        ptot = ptot + pevalb(q, p_test) # Baseline test function
398
        return ptot
399
        
400
def pevalb(q, b): return ((b[0]*exp(-0.5*(((q-b[1])/b[2])**2))) + b[3]) # A gaussian on a linear offset (to be used for the baseline)
401
def peval(q, b): return (b[0]*exp(-0.5*(((q-b[1])/b[2])**2))) # A gaussian with no linear offset (to be used for each of the superimposed gaussians)
402
def pfunc(b, r, q): return (r - peval(q, b))
403
def pfuncb(b, r, q): return (r - pevalb(q, b))
404
def pfuncf(b, r, q): return (r - peval_full(q, b))
405
def pfuncfb(b, r, q): return (r - pevalb_full(q, b))
406
def fit(function, p_ls, y_ls, x_ls):
407
        errfunc = lambda p_ls, x_lsl, y_lsl: function(p_ls, y_lsl, x_lsl)
408
        pfit, pcov, infodict, errmsg, success = myFit(errfunc, p_ls, args=(x_ls, y_ls), full_output=1, epsfcn=0.0001)
409
        chisq = ((infodict['fvec'])**2).sum()/(len(y_ls) - len(p_ls))
410
                # chisq is the same as (errfunc(pfit, x_ls, y_ls)**2).sum()/(len(y_ls)-len(p_ls))
411
        if (len(y_ls) > len(p_ls)) and pcov is not None:
412
                pcov = pcov * chisq
413
        else:
414
                pcov = inf
415
        pfit_out = [0.0]*len(pfit)
416
        perr_out = [0.0]*len(pfit)
417
        i=0
418
        while (i < len(pfit)):
419
                pfit_out[i] = pfit[i] # format the output array into a standard list
420
                try: perr_out[i] = (abs(pcov[i][i])**0.5)
421
                except: perr_out[i] = 0.0
422
                i=i+1
423
        return [pfit_out, perr_out, pcov, chisq]
424

    
425
########################################################### ELLIPSE EVALUATION FUNCTIONS
426
# The phase ellipse is based on twiss parameters: Alpha, Beta, Gamma, Epsilon
427
# (Epsilon = emittance)
428
# x` = (B +/- sqrt(B^2 - 4AC)) / 2A
429
#        A = Beta
430
#        B = -2 * Alpha * x
431
#        C = Gamma * x^2 - Epsilon
432
# x limits @ ( B^2 - (4 * A * C) ) = 0
433
#         -> x = +-((myBeta * myEmit)/(myBeta * myGamma - myAlpha^2))^0.5
434
def ellipse_limits(b):
435
        q = [abs(((b[1] * b[3])/(b[1] * b[2] - b[1]**2))**0.5),  abs(-1*((b[1] * b[3])/(b[1] * b[2] - b[1]**2))**0.5)]
436
        return (1.25*max(q)) # A little padding
437
        
438
def eval_ellipse(b, n, q_lim): # but I'll just use the q_lim=secMom[0]
439
        a_out = []
440
        b_out = []
441
        #q_lim = ellipse_limits(b)
442
        if (not isnan(q_lim)):
443
                q = arange(-1*q_lim,q_lim,((2*q_lim)/n))
444
                f = eval_ellipse_pos(q, b)
445
                g = eval_ellipse_neg(q, b)
446
                i=0
447
                while (i < len(f)):
448
                        if (not isnan(f[i])):
449
                                a_out.append(q[i])
450
                                b_out.append(f[i])
451
                        i=i+1
452
                i = len(g) - 1
453
                while (i > 0):
454
                        if (not isnan(g[i])):
455
                                a_out.append(q[i])
456
                                b_out.append(g[i])
457
                        i=i-1
458
                a_out.append(a_out[0])        # Close the ellipse
459
                b_out.append(b_out[0])        
460
        else:
461
                print "     -- WARNING: Cannot create ellipse with the given parameters:"
462
                print "           p_ellipse: " + str(b)
463
                print "           secMom[1]: " + str(q_lim)
464
                a_out.append(0.0)
465
                b_out.append(0.0)
466
        return a_out, b_out
467
        
468
def eval_ellipse_pos(q, b):        # Positive half of the sqrt
469
        #      ( (     B      + (    B^2        -  4 *  A *           C              )^0.5 )/(2 *  A  ))
470
        return ( ( (-2*b[0]*q) + ( (-2*b[0]*q)**2 - (4 * b[1]* ((b[2]*(q**2)) - b[3]) ) )**0.5)/(2 * b[1]))
471

    
472
def eval_ellipse_neg(q, b):        # Negative half of the sqrt
473
        #      ( (     B      - (    B^2        -  4 *  A *           C              )^0.5 )/(2 *  A  ))
474
        return ( ( (-2*b[0]*q) - ( (-2*b[0]*q)**2 - (4 * b[1]* ((b[2]*(q**2)) - b[3]) ) )**0.5)/(2 * b[1]))        
475

    
476
###########################################################
477

    
478
########################################################### READ IN DATA AND FORMAT IT
479
if (READ_LIVE):
480
        print  sdt() + " -- Reading data from the DC" + str(myCam) + " " + myPlane + "-Axis Histogram (" + str(myAver) + "x)" 
481
        print "  (" + myHistogram + ")"
482
        get_histogram(myCam,myPlane,myAver) # e.g. Camera #3, Y-Axis, 1 averages...
483
        myFile = open(myHistogram,"r")
484
else: 
485
        print  sdt() + " -- Reading data from " + mySample
486
        myFileStat = os.stat(mySample)
487
        myFile = open(mySample,"r")                        # Get File info
488
y_in = myFile.read().split('\n')                        # Read in data
489
myFile.close()
490
for myLine in y_in:
491
        myPart = myLine.split('\t')
492
        if(len(myPart) > 4):
493
                try:
494
                        y.append(float(myPart[4]))
495
                        x.append(float(myPart[2]))
496
                except:
497
                        pass
498
        else:
499
                try:
500
                        y.append(float(myPart[0]))
501
                except:
502
                        pass
503
n_pts = len(y)
504
if(n_pts == 0):
505
        print
506
        print "!! Failed to parse input data file !!"
507
        quit()
508
if (True):
509
        i = n_pts - 1
510
        while i > 0:
511
                if y[i] == 0:
512
                        y = delete(y,i)                        # Remove all trailing 0's
513
                        n_pts-=1
514
                        i-=1
515
                else:
516
                        i=0
517
print  sdt() + " -- Number of points in dataset: " + str(n_pts)
518
if(n_pts < 100):
519
        print "    Too few points! (Parsed=" + str(n_pts) + " / File=" + str(n_pts_ini) + ")"
520
        quit()
521
if(len(x) == n_pts):
522
        print "       -> Using x-axis information in units of [px]"
523
else:
524
        print "       -> Inconsistent or non-existant x-axis information available"
525
        x = arange(0.0,n_pts,1.0)
526
        x *= x_res
527
        
528
print "  Slits @ x" + str(l0) + ", YAG Screen @ x" + str(l1) + ", Total distance: " + str(myL) + "m."
529
print
530

    
531
if (DEBUG): plot_arrays([y],[x])
532

    
533
########################################################### ADD NOISE & NOISE FLOOR TO TEST DATA
534
if (not READ_LIVE) and (RANDOMIZE_NOISE):
535
        test_noise = test_noise * random.random()
536
        test_level = test_level * random.random()
537
        i=0
538
        print " -- Adding random noise (" + test_noise + ") and level (" + test_level + ") for testing..."
539
        while(i < n_pts):
540
                y[i] += (((random.random()-0.5)*2*test_noise) + test_level)
541
                #if ( y[i] < 0.0 ): y[i] = 0.0
542
                i+=1
543

    
544
########################################################### FIND PEAKS
545
if (DEBUG):
546
        i = 0
547
        while (( i < 10 ) and ( i < n_pts)):
548
                if (i == 0):
549
                        print "       y: [" + i2s(i,3) + "] " + e2s(y[i],4)
550
                else:
551
                        print "          [" + i2s(i,3) + "] " + e2s(y[i],4)
552
                i+=1
553
        if (i < n_pts): print "            ..."
554
y_s = smooth_array(y,smooth_range)
555
ps = [(max(y_s) - min(y_s)), ((max(x) - min(x))/2), ((max(x) - min(x))/20), min(y_s)]
556
#if (DEBUG): plot_arrays([y_s,y],[x])
557
ps,pse,psc,pscs = fit(pfuncb, ps, y_s, x)        # Fit to the smoothed data to find a baseline fit
558
lin_baseline = ps[3]
559
print sdt() + " -- Linear Baseline   = " + f2s(lin_baseline,4)
560
noise_level = abs(min(y - lin_baseline))
561
print sdt() + " -- Noise Level       = " + e2s(noise_level,4)
562
pick_threshold = (lin_baseline) + ((max(y) - min(y))*min_threshold) + (noise_level/2)
563
print sdt() + " -- Peak Pick Thresh  = " + e2s(pick_threshold,4)
564

    
565
if (DEBUG): # Check smooth fit to linear baseline
566
        y_s_eval = pevalb(x,ps)
567
        plot_arrays([y,y_s,y_s_eval],[x])
568

    
569
detect_threshold = 1.0 - ((noise_level/max(y))/0.35)
570
if (detect_threshold < 0.0): detect_threshold = 0.0
571
print sdt() + " -- Peak Detection  = " + f2s(detect_threshold,2)
572
if (detect_threshold > 1):
573
        print sdt() + " -- No Discernible Peaks For Fitting!! Quitting..."
574
        quit()
575

    
576
bump, x_cent = find_peaks(y, min_peak_width, pick_threshold)
577
n_bumps = len(bump)
578
if (n_bumps < MINIMUM_SLITS):
579
        print sdt() + " -- Detected only " + str(n_bumps) + " peak(s)! Quitting..."
580
        quit()
581
print sdt() + " -- Number of peaks found: " + str(n_bumps)
582
print sdt() + " -- Center index: " + str(x_cent)
583

    
584
# Find the bump data above the FWHM and the FWFM data for each
585
x_fwhm,y_fwhm,x_part,y_part=[],[],[],[]
586
i = 0
587
while (i < n_bumps):
588
        y_part.append(y[bump[i][3]:bump[i][4]])
589
        x_part.append(x[bump[i][3]:bump[i][4]])
590
        i+=1
591

    
592
if (DEBUG): plot_arrays(y_part,x_part)
593

    
594
########################################################### PERFORM INITIAL (SINGLE-GAUSSIAN) FIT(S)
595
# Set initial guesses for fitting
596
print sdt() + " -- Fitting single gaussians"
597
a =([[0.0, 0.0, 1.0, 0.0]]) * (n_bumps + 1)
598
ae=([[0.0, 0.0, 0.0, 0.0]]) * (n_bumps + 1)
599
full_chisq = 0.0
600
i=0
601
while (i < n_bumps): # Find Initial Fit Guesses for each bump
602
        #           Amp                                Cent                Sigma                 Lin Off
603
        a[i] = ([bump[i][0], x[bump[i][1]], (bump[i][2] * x_res)/8, lin_baseline])
604
        i+=1
605
i=0
606
while(i < n_bumps):
607
        a[i], ae[i], p_conv, p_chisq = fit(pfuncb, a[i], y_part[i], x_part[i])        # Perform Fit (including baseline)
608
        full_chisq += (p_chisq / n_bumps)
609
        i+=1
610

    
611
########################################################### FIND A BASELINE FIT / MEDIAN
612
lin_offset = median(y) # Also used as a fit parameter for the envelope if no background subtraction is performed
613
if (SUBTRACT_BACKGROUND == 1):
614
        print sdt() + " -- Fitting and subtracting the baseline"
615
        y_residual = (y - pevalb_full(x, flatten_array(a)))
616
        ys_res = smooth_array(y_residual,smooth_range)
617
        pb=([(max(ys_res)-min(ys_res))/2, (max(x)-min(x))/2, (max(x)-min(x))*10, lin_baseline])
618
        pb,pbe,pbc,pbcs = fit(pfuncb, pb, ys_res, x)        # Perform smoothed fit (including baseline)
619
        pb,pbe,pbc,pbcs = fit(pfuncb, pb, y_residual, x)# Perform baseline fit (including baseline)
620
        y_base = pevalb(x, pb) # Adjust these functions as needed to fit to the background
621
        if (DEBUG):
622
                plot_arrays([y_residual,y_base],[x])
623
        a[n_bumps] = pb
624
        y -= y_base
625
        y_s -= y_base
626
elif (SUBTRACT_BACKGROUND == 2):
627
        print sdt() + " -- Subtracting the median"
628
        a[n_bumps] = [0,0,0,lin_offset]
629
        y_base = []
630
        i = 0
631
        while(i < n_pts):
632
                y_base.append(a[n_bumps][3])
633
                i += 1
634
        if (DEBUG):
635
                plot_arrays([y,y_base],[x])
636
        y -= a[n_bumps][3]
637
        y_s -= a[n_bumps][3]
638

    
639
i=0
640
while(i <= n_bumps): # Remove the baseline offset on all individual gaussians...
641
        a[i][3] = 0.0
642
        i+=1
643
########################################################### PERFORM FULL FIT (IF DESIRED)
644

    
645
if FULL_FIT:
646
        print sdt() + " -- Performing Full Fit"
647
        p_full, pe_full, full_cov, full_chisq = fit(pfuncfb, flatten_array(a), y, x)        # Perform Fit
648
        a = form_array(p_full, 4)
649
        print sdt() + " -- Full Fit Chi^2 = " + e2s(full_chisq,4)
650

    
651

    
652
########################################################### RE-FIT/ELIMINATE BAD GAUSSIANS
653
print sdt() + " -- Eliminating Bad Gaussians... "
654
# Sigma sanity check on background
655
if ((not SUBTRACT_BACKGROUND) and (abs(a[n_bumps][2]) < 1e-4)):
656
        print "      --> Bad Background Fit  (Params: " + str(a[i]) + ")"
657
        a[n_bumps][0] = 0.0
658
        a[n_bumps][3] = lin_baseline
659
        p_full, pe_full, full_cov, full_chisq = fit(pfuncfb, flatten_array(a), y, x)        # Perform Fit
660

    
661
i=0
662
while (i < n_bumps):
663
        # Amplitude and Sigma sanity checks on individual gaussians
664
        if (a[i][0] < 0.0 or ((abs(a[i][2]) > 1e-3) or (abs(a[i][2]) < 1e-7))):
665
                print "    " + str(i+1) + " --> Bad Fit  (Params: " + str(a[i]) + ")... ignoring in calculations"
666
                a[i][0] = 0.0
667
                a[i][2] = 1.0
668
        i+=1
669

    
670
########################################################### INTEGRATE UNDER BEAMLET GAUSSIANS
671
i=0
672
I=[]
673
p_eval=[0.0, 0.0, 0.0, 0.0]
674
while (i < n_bumps):
675
        p_eval = [(a[i][0] - a[i][3]), a[i][1], a[i][2], 0.0] # Is this fair to do if we don't do the background subtraction? Otherwise we just want p_eval = a[i]
676
        result = integrate(peval, min(x), max(x), p_eval)
677
        I.append(result[0])        # Integration for each gaussian
678
        i+=1
679

    
680
########################################################### FIT TO ENVELOPE AT THE SLIT MASK
681
print sdt() + " -- Fitting to the envelope at the slit mask"
682
x_env,y_env=[],[]
683
i=0
684
while (i<n_bumps):
685
        if (a[i][0] > 0.0 and I[i] > 0.0): # Ignore bad gaussians marked above
686
                m = (i-x_cent)
687
                x_env.append((m*myS) + a[x_cent][1])        # At the slit mask all points will by definition be at intervals of myS
688
                y_env.append(a[i][0])                         # Use max amplitude for each fit beamlet gaussian
689
                #y_env.append(max(y)*I[i]/max(I))            # Use a scaled integral for each beamlet gaussian
690
        i+=1
691
if(SUBTRACT_BACKGROUND > 0):
692
        p_env = [(max(y_env)), (max(x_env) + min(x_env))/2, max(x_env)/2, 0.0]
693
        p_env, pe_env, env_cov, env_chisq = fit(pfunc, p_env, y_env, x_env)                # Perform Fit
694
        y_env_eval = peval(x, p_env)
695
else:
696
        p_env = [(max(y_env)), (max(x_env) + min(x_env))/2, max(x_env)/2, lin_offset]
697
        p_env, pe_env, env_cov, env_chisq = fit(pfuncb, p_env, y_env, x_env)                # Perform Fit with baseline
698
        y_env_eval = pevalb(x, p_env)
699
        i = 0
700
        while (i < n_pts):
701
                if (y_env_eval[i] < 0.0):
702
                        y_env_eval[i] = 0.0
703
                i += 1
704
print sdt() + " -- Envelope Fit Chi^2 = " + e2s(env_chisq,4)
705
print sdt() + " -- Envelope Centroid = " + e2s(p_env[1],4)
706

    
707
if (DEBUG): 
708
        plot_arrays([y,y_env_eval,y_env],[x,x,x_env])
709

    
710
########################################################### INTEGRATE UNDER BACKGROUND/ENV GAUSSIANS
711

    
712
result = integrate(pevalb, min(x), max(x), args=(a[n_bumps]))
713
I.append(result[0])        # Background integration in I[n_bumps]
714

    
715
result = integrate(peval, min(x), max(x), args=(p_env))        # Evaluate the envelope with no baseline
716
I.append(result[0])         # Envelope integration in I[n_bumps + 1]
717

    
718
########################################################### CENTER ON ENVELOPE CENTROID
719
print sdt() + " -- Centering on the Envelope Centroid..."
720

    
721
if (CENTER_ON_ENVELOPE):
722
        trans_offset = p_env[1]                # Remove the offset between original x and envelope centroid
723
        i=n_bumps
724
        i_offset = 0
725
        while (i >= 0):
726
                if ( (a[i][1] - trans_offset) > 0 ): i_offset = i        # Find the peak just above the envelope centroid... aka a good solid peak
727
                i-=1
728
        if (i_offset > 0):
729
                frac_offset = (a[i_offset][1] - trans_offset) / (a[i_offset][1] - a[i_offset-1][1])
730
        else:
731
                CENTER_ON_ENVELOPE = False
732
                trans_offset = a[x_cent][1]        # Just center on the peak after all...
733
                frac_offset = 0.0
734
else:
735
        trans_offset = a[x_cent][1]        # Remove the offset between original x and highest peak
736

    
737
# Translate x-axis components
738
x -= trans_offset
739

    
740
# Translate all fit parameters
741
x_min, x_max = x[len(x)-1],0.0
742
i=0
743
while (i < n_bumps):
744
        a[i][1] -= trans_offset
745
        if ( a[i][1] > x_max ): x_max = a[i][1]
746
        if ( a[i][1] < x_min ): x_min = a[i][1]
747
        i+=1
748
a[n_bumps][1] -= trans_offset # Shift Background as well
749
p_env[1] -= trans_offset
750

    
751
print "      -- Centering on Envelope: " + str(CENTER_ON_ENVELOPE)
752
print "      -- Positive Center Index: " + str(i_offset)
753
print "      -- Fractional Offset: " + e2s(frac_offset,4)
754

    
755
########################################################### FIND EMITTANCE PARAMETERS
756
i=0
757
if CENTER_ON_ENVELOPE:
758
        myOffset = (i_offset - frac_offset)        # Center about Envelope
759
else:
760
        myOffset = x_cent                         # Center about highest peak
761

    
762
myDiv = []
763
secMom = [0.0] * 3
764
myIntSum = 0.0
765
while ( i < n_bumps ):
766
        if (a[i][0] > 0.0):        # So long as there's a peak, include it in the calculation...
767
                m = (i-myOffset)
768
                myIntSum = myIntSum + I[i]                                        # Sum Integrals
769
                myDiv.append((a[i][1] - (m * myS)) / (myL))                        # Calculate Divergence
770
        else:
771
                myDiv.append(0.0)
772
        i+=1
773
i=0
774
while ( i < n_bumps ):
775
        if (a[i][0] > 0.0):
776
                m = (i-myOffset)
777
                secMom[0] = secMom[0] + (I[i] * (m*myS)**2)                                                                # <x^2>
778
                secMom[1] = secMom[1] + (I[i] * ((myDiv[i])**2 + ((a[i][2])**2 - ((myW**2)/12.0))/(myL**2)))                # <x'^2>
779
                secMom[2] = secMom[2] + (I[i] * (m*myS) * myDiv[i])                                                        # <xx'>
780
        i+=1
781
i=0
782
while ( i < len(secMom)):
783
        secMom[i] = (secMom[i] / (myIntSum))
784
        i+=1
785
myEmitSq = (secMom[0]*secMom[1]) - (secMom[2]*secMom[2])
786
myEmit = (myEmitSq)**(0.5)
787
mydEmit = (p_env[2] * secMom[1]) / myEmit        # Gives differential as a fraction of the Emittance
788
myEmitum = myEmit * 1000.0 * 1000.0                # x(mm/m) & x'(mrad/rad) = um
789
mydEmitum = mydEmit * myEmitum                        # Statistical Error
790
myAlpha = -1.0 * (secMom[2]/myEmit)
791
myBeta = secMom[0]/myEmit
792
myGamma = (1.0 + (myAlpha**2))/myBeta
793

    
794
print sdt() + " -- Emitance = " + f2s(myEmitum,4) + " +/- " + e2s(mydEmitum,4) + " um (" + e2s(myEmit,4) + " m*Rad)"
795
print "          Emit = sqrt(" + e2s(secMom[0]*secMom[1],4) + " -" + e2s((secMom[2]*secMom[2]),4) + " ) m*Rad"
796

    
797

    
798
########################################################### CONTOUR PLOT
799

    
800
########################################################### FRAME OUTPUT AND GENERATE CURVES
801

    
802
xi_min = 0
803
xi_max = (len(x) - 1)
804
if (FRAME_OUTPUT):
805
        x_margin = (x_max - x_min) * 0.15
806
        x_min -= x_margin
807
        x_max += x_margin
808
        while (x[xi_min] < x_min): xi_min += 1
809
        while (x[xi_max] > x_max): xi_max -= 1
810
x_out = x[xi_min:xi_max]
811
y_out = y[xi_min:xi_max]
812
y_fit = pevalb_full(x_out, flatten_array(a))
813
y_env = pevalb(x_out,p_env)
814
y_s_out = y_s[xi_min:xi_max]
815

    
816
########################################################### CALCULATE PHASE ELLIPSE
817

    
818
if (not isnan(myEmit)):
819
        p_ellipse = [myAlpha, myBeta, myGamma, myEmit]
820
        x_ellipse, y_ellipse = eval_ellipse(p_ellipse, len(x), (secMom[0]**(0.5)))
821
else:
822
        print " WARNING!! -> Non-physical parameters found!"
823
        x_ellipse = [0.0] * n_pts
824
        y_ellipse = [0.0] * n_pts
825

    
826
########################################################### OUTPUT DATA
827
if (VERBOSE):
828
        print
829
        print " -- SUMMARY --"
830
        print "    Second Moments:"
831
        print "         <x^2>  = " + e2s(secMom[0],4)
832
        print "         <x'^2> = " + e2s(secMom[1],4)
833
        print "         <xx'>  = " + e2s(secMom[2],4)
834
        print
835
        print "    Alpha : " + f2s(myAlpha,4)
836
        print "    Beta  : " + f2s(myBeta,4)
837
        print "    Gamma : " + f2s(myGamma,4)
838
        print
839
        print "    Phase Ellipse: "
840
        print "          Max x  : " + e2s(max(abs(x_ellipse)), 4)
841
        print "          Max x' : " + e2s(max(abs(y_ellipse)), 4)
842
        print
843
        print " -- Bump Detection"
844
        print " Bump  YMax      X(YMax)      Width        Xi       Xf"
845
        k=1
846
        for i in bump:
847
                myOut = i2s(k,3)
848
                for j in i:
849
                        myOut += " " + e2s(j,3)
850
                print myOut
851
                k+=1
852
        print
853
        print " -- Final Fit Params & Div Calc"
854
        print " Bump   Amp        Cent        Sig      LinOff   Integral  Divergence"
855
        k=1
856
        for i in a:
857
                myOut = i2s(k,3)
858
                for j in i:
859
                        myOut += " " + e2s(j,3)
860
                myOut += " " + e2s(I[k-1],3)
861
                if (k > n_bumps):
862
                        print
863
                        print "   Background fit:"
864
                else:
865
                        myOut += " " + e2s(myDiv[k-1],3)
866
                print myOut
867
                k+=1
868
        print
869

    
870
if (FULL_OUTPUT):
871
        myOut = ""
872
        print sdt() + " -- Saving Summary to " + OUTPUT_DIRECTORY + "ms_emit.htm"
873
        myOut += "<html><body><div><span style='float:left;font-size:20px;'>Multislit Emittance Summary</span><pre style='float:right;'>" + snow() +"</pre><br>"
874
        if (not READ_LIVE):
875
                myOut += "<pre>File: " + mySample + " (" + sfiletime() + ")</pre>"
876
        myOut += "<table width=100%><tr>"
877
        myOut += "<td valign=top><div><pre>\n\
878
        <b>Emittance & System Params:</b>\n\
879
        &nbsp;&nbsp;&nbsp;&nbsp;Emittance (Un-Norm) : <b>" + e2s(myEmitum,3) + " &mu;m</b>\n\
880
        &nbsp;&nbsp;&nbsp;&nbsp;Dist to Screen (L)  : " + e2s(myL,3) + " m\n\
881
        &nbsp;&nbsp;&nbsp;&nbsp;Slit Width (w)      : " + e2s(myW,3) + " m\n\
882
        &nbsp;&nbsp;&nbsp;&nbsp;Slit Spacing (s)    : " + e2s(myS,3) + " m\n\
883
        &nbsp;&nbsp;&nbsp;&nbsp;Histogram Amplitude : " + f2s(max(y),2) + " ArbU\n\
884
        &nbsp;&nbsp;&nbsp;&nbsp;Histogram Width     :  " + str(n_pts) + " px\n\
885
        &nbsp;&nbsp;&nbsp;&nbsp;Camera Resolution   : " + e2s(x_res,3) + " m/px\n\
886
        </pre></div></td>"
887
        
888
        myOut += "<td valign=top><div><pre>\n\
889
        <b>Second Moments:</b>\n\
890
        &nbsp;&nbsp;&nbsp;&nbsp;&lt;x^2&gt;  : " + e2s(secMom[0],3) + " m<sup>2</sup>\n\
891
        &nbsp;&nbsp;&nbsp;&nbsp;&lt;x`^2&gt; : " + e2s(secMom[1],3) + " rad<sup>2</sup>\n\
892
        &nbsp;&nbsp;&nbsp;&nbsp;&lt;xx`&gt;  : " + e2s(secMom[2],3) + " m*rad\n\
893
        <b>Checks:</b>\n\
894
        &nbsp;&nbsp;&nbsp;&nbsp;Full &chi;<sup>2 </sup> : " + e2s(full_chisq,3) + "\n\
895
        &nbsp;&nbsp;&nbsp;&nbsp;&sigma;<sub>Env</sub><sup>2</sup>    : " + e2s((p_env[2]**2),3) + " m<sup>2</sup>\n\
896
        &nbsp;&nbsp;&nbsp;&nbsp;&lt;x^2&gt / &sigma;<sub>Env</sub><sup>2</sup>  : " + f2s((secMom[0]/(p_env[2]**2)),3) + "\
897
        </pre></div></td>"
898
        
899
        myOut += "<td valign=top><div><pre>\n\
900
        <b>Twiss Parameters:</b>\n\
901
        &nbsp;&nbsp;&nbsp;&nbsp;&alpha; : " + f2s(myAlpha,3) + "\n\
902
        &nbsp;&nbsp;&nbsp;&nbsp;&beta; : " + f2s(myBeta,3) + "\n\
903
        &nbsp;&nbsp;&nbsp;&nbsp;&gamma; : " + f2s(myGamma,3) + "\n"
904
        
905
        myOut += "        <b>Background Subtraction:</b>\n        &nbsp;&nbsp;&nbsp;&nbsp;"
906
        if (SUBTRACT_BACKGROUND == 1):
907
                myOut += "Gaussian Fit"
908
        elif (SUBTRACT_BACKGROUND == 2):
909
                myOut += "Median"
910
        else:
911
                myOut += "None"
912
        myOut += "\n"
913
        if (FULL_FIT): 
914
                myOut += "        <b>Full Final Fit:</b>\n        &nbsp;&nbsp;&nbsp;&nbsp;--> " + str(len(p_full)) + " Parameters"
915
        else:
916
                myOut += "        <b>Single Gaussian Fits Only</b>"
917

    
918
        myOut += "</pre></div></td>\
919
        </tr></table>"
920
        
921
        myOut += "<table width=100%><tr><td width=50%><img width=100% src='./ms_emit.png'></td>\
922
        <td><img width=100% src='./ms_phase.png'></td></tr><tr><td colspan=2><center>\
923
        <div><pre>" + str(n_bumps) + " Gaussians Detected (Threshold = " + f2s(pick_threshold,3) + " | Noise = " + f2s(noise_level,3) + ")</pre><hr><pre><table width=100%>\
924
        <tr><td>Gaussian</td><td>Amplitude</td><td>Centroid [m]</td><td>Sigma [m]</td><td>LinOff</td><td>Integral [m]</td><td>Div [Rad]</td></tr>"
925
        i=0
926
        while(i < n_bumps):
927
                myOut += "<tr><td>" + str(i+1) + "</td><td>" + e2s(a[i][0],3) + "</td><td>" + e2s(a[i][1],3) + "</td><td>" + e2s(abs(a[i][2]),3) + "</td><td>--</td><td>" + e2s(I[i],3) + "</td><td>" + e2s(myDiv[i],3) + "</td></tr>"
928
                i+=1
929
        i = n_bumps
930
        if(SUBTRACT_BACKGROUND == 1):
931
                myOut += "<tr><td>Baseline</td><td>" + e2s(a[i][0],3) + "</td><td>" + e2s(a[i][1],3) + "</td><td>" + e2s(abs(a[i][2]),3) + "</td><td>" + e2s(a[i][3],3) + "</td><td>" + e2s(I[i],3) + "</td><td>&nbsp;</td></tr>"
932
        myOut += "<tr><td>Proj Env</td><td>" + e2s(p_env[0],3) + "</td><td>" + e2s(p_env[1],3) + "</td><td>" + e2s(abs(p_env[2]),3) + "</td><td>" + e2s(lin_baseline,3) + "</td><td>" + e2s(I[i+1],3) + "</td><td>&nbsp;</td></tr>\
933
        </table>"
934
        if(SUBTRACT_BACKGROUND == 2):
935
                myOut += "<span style='float:left'>Median Baseline (subtracted): " + e2s(lin_offset,3) + "</span><br>"
936
        myOut += "<hr></td></tr></table></div></body></html>"
937
        myFile = open(OUTPUT_DIRECTORY+"ms_emit.htm","w")
938
        myFile.write(myOut)
939
        myFile.close()
940
        
941
        myOut = ""
942
        print sdt() + " -- Saving Plot Data to " + OUTPUT_DIRECTORY + "ms_emit.dat"
943
        myOut = "x        x\'        x [px]        x [m]        y        y_fit        y_envelope\n"
944
        i=0
945
        j=0
946
        while ((i < len(x_ellipse)) or (i < len(x_out))):
947
                if ( i < len(x_ellipse) ):
948
                        myOut += str(x_ellipse[i]) + "\t" + str(y_ellipse[i])
949
                else:
950
                        myOut += "--\t--"
951
                if ( i < len(x_out) ):
952
                        myOut += "\t" + str(x_out[i]) + "\t" + str(y_out[i]) + "\t" + str(y_fit[i]) + "\t" + str(y_env_eval[i]) + "\t" + str(y_s_out[i])
953
                myOut += "\n"
954
                i+=1
955
        myFile = open(OUTPUT_DIRECTORY+"ms_emit.dat","w")
956
        myFile.write(myOut)
957
        myFile.close()
958
        
959
        print sdt() + " -- Generating plots"
960
        call("/usr/bin/gnuplot " + OUTPUT_DIRECTORY + "ms_emit.gnuPlot", shell=True)
961
        
962
        print sdt() + " -- Syncing data to www-bd"
963
        call("/usr/bin/wget -O /export/home1/edstrom/acl/emit_calc/ms_emit.sync -o /export/home1/edstrom/acl/emit_calc/ms_emit.wget http://www-bd.fnal.gov/cgi-srf/fast.sync.pl", shell=True)
964
else:
965
        print sdt() + " -- Not Syncing to www-bd"
966

    
967
if(READ_LIVE):
968
        print sdt() + " -- Output to ACNET"
969
        tnow = (now() - t0)
970
        result = set_output([myEmitum, myAlpha, myBeta, myGamma, n_bumps, noise_level, pick_threshold, full_chisq, tnow])
971
        tar = tarfile.open(ARCHIVE_DIRECTORY + "ms_emit." + filedate + ".tar", "w")
972
        tar.add(OUTPUT_DIRECTORY + "ms_emit.htm","ms_emit.htm")
973
        tar.add(OUTPUT_DIRECTORY + "ms_emit.png","ms_emit.png")
974
        tar.add(OUTPUT_DIRECTORY + "ms_phase.png","ms_phase.png")
975
        tar.add(ARCHIVE_DIRECTORY + "fast_ms_emit.hist","ms_emit.hist")
976
        tar.close()
977

    
978
print sdt() + " -- DONE!"
979

    
980
if (DEBUG):
981
        y_fit = pevalb_full(x, flatten_array(a))
982
        print sdt() + " -- Plotting Summary..."
983
        plot_arrays([y_env_eval,y_s,y,y_fit],[x])
984
        print sdt() + " -- Plotting phase ellipse..."
985
        plot_arrays([y_ellipse],[x_ellipse])
986
quit()