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


17

myW = 40e6

18

myS = 400e6

19

myX107 = 9221.9455

20

myX108 = 9676.6585

21

myX109 = 10058.638

22

myX111 = 10759.953

23

myX115 = 14000.000

24

myX120 = 16600.485

25

myX124 = 19147.663

26

l0 = 107

27

l1 = 111

28

z0 = myX107

29

z1 = myX111

30

myL = (z1  z0)/1000

31

x_res = 9.3e6

32


33


34

myCam = 4

35

myPlane = "X"

36

myAver = 5

37

MINIMUM_SLITS = 3

38


39


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


48

OUTPUT_DIRECTORY = "/usr/local/userb/fast/sync/"

49

ARCHIVE_DIRECTORY = "/usr/local/cbs_files/cns_write/acl/data_files/fast_ms_emit/"

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


55

smooth_range = 40

56

pick_smooth = 4

57

min_peak_width = 20

58

min_threshold = 0.002

59

n_pts = 0

60


61


62

SUBTRACT_BACKGROUND = 1

63

CENTER_ON_ENVELOPE = True

64

VERBOSE = True

65

FRAME_OUTPUT = False

66

FULL_FIT = True

67

FULL_OUTPUT = True

68


69

y_in,x_in,y,x = [],[],[],[]

70

a = []

71

p_full = []

72


73


74


75


76


77

bump= []

78

i,j,k = 0,0,0

79

myArg = []

80


81


82

print "Welcome to the FAST Multislit emittance scanner."

83


84

if (len(sys.argv) > 1):

85

myArg = sys.argv[1].split(',')

86

nArgs = len(myArg)

87

if (nArgs < 1):

88

print "No commandline arguments passed..."

89

print "> Syntax: python fast_ms_emit.py arg1=val1,arg2=val2..."

90

print "> Valid arguments: <FILENAME>,cam=(14),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 commandline 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


203


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=isf

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])/(kj))

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=i1

280

while ( j < len(a_in) ):

281

if ( a_in[j] < ((a_in[a_max] + thresh)/2)): break

282

j=j+1

283

if ((ji) < 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)):

319

myCent = a_out[i][1]

320

if (len(a_out) <= 1):

321

return [0],0

322

if (i == 0):

323

myWid = (a_out[i+1][1]  myCent)

324

else:

325

myWid = (myCent  a_out[i1][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


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):

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


369

def peval_full(q, b):

370


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]]

376

else:

377

p_test = b[i]

378

ptot = ptot + peval(q, p_test)

379

i=i+1

380

return ptot

381


382

def pevalb_full(q, b):

383


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]]

389

else:

390

p_test = b[i]

391

ptot = ptot + peval(q, p_test)

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)

398

return ptot

399


400

def pevalb(q, b): return ((b[0]*exp(0.5*(((qb[1])/b[2])**2))) + b[3])

401

def peval(q, b): return (b[0]*exp(0.5*(((qb[1])/b[2])**2)))

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


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]

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


426


427


428


429


430


431


432


433


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))

437


438

def eval_ellipse(b, n, q_lim):

439

a_out = []

440

b_out = []

441


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=i1

458

a_out.append(a_out[0])

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):

469


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):

473


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


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)

483

myFile = open(myHistogram,"r")

484

else:

485

print sdt() + "  Reading data from " + mySample

486

myFileStat = os.stat(mySample)

487

myFile = open(mySample,"r")

488

y_in = myFile.read().split('\n')

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)

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 xaxis information in units of [px]"

523

else:

524

print " > Inconsistent or nonexistant xaxis 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


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


542

i+=1

543


544


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


557

ps,pse,psc,pscs = fit(pfuncb, ps, y_s, x)

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):

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


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


595


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):

602


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])

608

full_chisq += (p_chisq / n_bumps)

609

i+=1

610


611


612

lin_offset = median(y)

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)

619

pb,pbe,pbc,pbcs = fit(pfuncb, pb, y_residual, x)

620

y_base = pevalb(x, pb)

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):

641

a[i][3] = 0.0

642

i+=1

643


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)

648

a = form_array(p_full, 4)

649

print sdt() + "  Full Fit Chi^2 = " + e2s(full_chisq,4)

650


651


652


653

print sdt() + "  Eliminating Bad Gaussians... "

654


655

if ((not SUBTRACT_BACKGROUND) and (abs(a[n_bumps][2]) < 1e4)):

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)

660


661

i=0

662

while (i < n_bumps):

663


664

if (a[i][0] < 0.0 or ((abs(a[i][2]) > 1e3) or (abs(a[i][2]) < 1e7))):

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


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]

676

result = integrate(peval, min(x), max(x), p_eval)

677

I.append(result[0])

678

i+=1

679


680


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):

686

m = (ix_cent)

687

x_env.append((m*myS) + a[x_cent][1])

688

y_env.append(a[i][0])

689


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)

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)

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


711


712

result = integrate(pevalb, min(x), max(x), args=(a[n_bumps]))

713

I.append(result[0])

714


715

result = integrate(peval, min(x), max(x), args=(p_env))

716

I.append(result[0])

717


718


719

print sdt() + "  Centering on the Envelope Centroid..."

720


721

if (CENTER_ON_ENVELOPE):

722

trans_offset = p_env[1]

723

i=n_bumps

724

i_offset = 0

725

while (i >= 0):

726

if ( (a[i][1]  trans_offset) > 0 ): i_offset = i

727

i=1

728

if (i_offset > 0):

729

frac_offset = (a[i_offset][1]  trans_offset) / (a[i_offset][1]  a[i_offset1][1])

730

else:

731

CENTER_ON_ENVELOPE = False

732

trans_offset = a[x_cent][1]

733

frac_offset = 0.0

734

else:

735

trans_offset = a[x_cent][1]

736


737


738

x = trans_offset

739


740


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

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


756

i=0

757

if CENTER_ON_ENVELOPE:

758

myOffset = (i_offset  frac_offset)

759

else:

760

myOffset = x_cent

761


762

myDiv = []

763

secMom = [0.0] * 3

764

myIntSum = 0.0

765

while ( i < n_bumps ):

766

if (a[i][0] > 0.0):

767

m = (imyOffset)

768

myIntSum = myIntSum + I[i]

769

myDiv.append((a[i][1]  (m * myS)) / (myL))

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 = (imyOffset)

777

secMom[0] = secMom[0] + (I[i] * (m*myS)**2)

778

secMom[1] = secMom[1] + (I[i] * ((myDiv[i])**2 + ((a[i][2])**2  ((myW**2)/12.0))/(myL**2)))

779

secMom[2] = secMom[2] + (I[i] * (m*myS) * myDiv[i])

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

788

myEmitum = myEmit * 1000.0 * 1000.0

789

mydEmitum = mydEmit * myEmitum

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


799


800


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


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!! > Nonphysical parameters found!"

823

x_ellipse = [0.0] * n_pts

824

y_ellipse = [0.0] * n_pts

825


826


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[k1],3)

861

if (k > n_bumps):

862

print

863

print " Background fit:"

864

else:

865

myOut += " " + e2s(myDiv[k1],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;fontsize: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

Emittance (UnNorm) : <b>" + e2s(myEmitum,3) + " μm</b>\n\

880

Dist to Screen (L) : " + e2s(myL,3) + " m\n\

881

Slit Width (w) : " + e2s(myW,3) + " m\n\

882

Slit Spacing (s) : " + e2s(myS,3) + " m\n\

883

Histogram Amplitude : " + f2s(max(y),2) + " ArbU\n\

884

Histogram Width : " + str(n_pts) + " px\n\

885

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

<x^2> : " + e2s(secMom[0],3) + " m<sup>2</sup>\n\

891

<x`^2> : " + e2s(secMom[1],3) + " rad<sup>2</sup>\n\

892

<xx`> : " + e2s(secMom[2],3) + " m*rad\n\

893

<b>Checks:</b>\n\

894

Full χ<sup>2 </sup> : " + e2s(full_chisq,3) + "\n\

895

σ<sub>Env</sub><sup>2</sup> : " + e2s((p_env[2]**2),3) + " m<sup>2</sup>\n\

896

<x^2> / σ<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

α : " + f2s(myAlpha,3) + "\n\

902

β : " + f2s(myBeta,3) + "\n\

903

γ : " + f2s(myGamma,3) + "\n"

904


905

myOut += " <b>Background Subtraction:</b>\n "

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 > " + 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> </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> </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 wwwbd"

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://wwwbd.fnal.gov/cgisrf/fast.sync.pl", shell=True)

964

else:

965

print sdt() + "  Not Syncing to wwwbd"

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()
