Commit 781a41d0 authored by Zolnierczuk, Piotr's avatar Zolnierczuk, Piotr
Browse files

echo error estimate

parent 06031135
Loading
Loading
Loading
Loading

misc/get_beam_power.py

0 → 100755
+72 −0
Original line number Diff line number Diff line
#!/usr/bin/env python

"""
Web 
"""
import os.path
import time

import urllib.request, urllib.parse, urllib.error
import html.parser

class BeamStatusParser(html.parser.HTMLParser):
    results = {}
    current_tag = None

    def handle_starttag(self, tag, attrs):
        if tag == 'span':
            for attr in attrs:
                if attr[0] == 'id':
                    key = attr[1]
                    if key.startswith("val_"):
                        self.current_tag = key[4:]
                    else:
                        self.current_tag = key

    def handle_endtag(self, tag):
        self.current_tag = None

    def handle_data(self, data):
        if not self.current_tag:
            return
        values = data.strip().split()
        reslist = []
        succ = False
        for val in values:
            try:
                val = eval(val)
                succ = True
            except NameError as ValueError:
                pass
            finally:
                reslist.append(val)
        if not succ:
            reslist = " ".join(reslist)
        self.results[self.current_tag] = reslist


def get_beam_power(uri):
    fd = urllib.request.urlopen(uri)
    parser = BeamStatusParser()
    parser.feed(fd.read().decode("utf-8"))
    #power = parser.results.get('beam_power', [0.0, 'kW'])
    # TODO handle units
    return parser.results

def main():
    import argparse
    default_uri  = 'http://status.sns.ornl.gov/getRoundBox.jsp?which=SNSBeamInfo'

    parser = argparse.ArgumentParser(description='beam_status')
    parser.set_defaults(uri=default_uri)
    parser.add_argument('--uri', '-u',  dest='uri',  help='uri to upload the info   (default=%(default)s)')

    args = parser.parse_args()

    res = get_beam_power(args.uri)
    for key in res:
	    print("%-24s: %s" % (key, res[key]))

if __name__ == "__main__":
    main()

misc/get_proposal.py

0 → 100755
+34 −0
Original line number Diff line number Diff line
#!/usr/bin/env python

"""

"""
import os.path
import time

import urllib.request, urllib.error, urllib.parse, ssl
import json

IPTS_CALIB=4531

url1=r'http://snsapp1.sns.ornl.gov/xprod_ro/proposal_production/express/getScheduledIPTSByInstrId/BL-15'
url2=r'https://snsapp1.sns.ornl.gov/xprod_ro/proposal_production/express/getTitleTeamMemberLocalContactByInstrIdPrpslId/BL-15/%d'


def get_json_uri(uri):
	req  = urllib.request.Request(uri)
	cont = ssl.SSLContext(ssl.PROTOCOL_TLSv1)
	resp = urllib.request.urlopen(req, context=cont)
	return json.loads(resp.read())

	

proposals = get_json_uri(url1)
print((list(proposals.keys())))
for item in proposals['items']:
	ipts = item.get('ipts', -1)
	print(ipts)
	if ipts>0: # and ipts!=IPTS_CALIB:
		pdata = get_json_uri(url2 % ipts)
		print(pdata)
+28 −0
Original line number Diff line number Diff line
#!/usr/bin/env python
"""

"""
import sys
import os.path
import time

import urllib.request, urllib.error, urllib.parse, ssl
import json
import pprint

url1=r'http://snsapp1.sns.ornl.gov/xprod_ro/proposal_production/express/getScheduledIPTSByInstrId/BL-15'
url2=r'https://snsapp1.sns.ornl.gov/xprod_ro/proposal_production/express/getTitleTeamMemberLocalContactByInstrIdPrpslId/BL-15/%d'

def get_json_uri(uri):
	req  = urllib.request.Request(uri)
	cont = ssl.SSLContext(ssl.PROTOCOL_TLSv1)
	resp = urllib.request.urlopen(req, context=cont)
	return json.loads(resp.read())


ipts = int(sys.argv[1])
#print(ipts)
if ipts>0: # and ipts!=IPTS_CALIB:
    pdata = get_json_uri(url2 % ipts)
    pprint.pprint(pdata.get('items'))
+9 −6
Original line number Diff line number Diff line
@@ -62,20 +62,23 @@ def linear_fit(x, y, yerr, chionly=False):

def error_estimate(sfun, dj, dj0, counts, ecounts):
    """estimate the error using chi^2+1 technique"""
    from numpy.polynomial.polynomial import Polynomial
    dx = np.abs(dj[1]-dj[0])/2 # FIXME dx hardcoded (kindof)
    # three points for a parabola
    x  = np.asarray((dj0 - dx, dj0, dj0 + dx))
    y  = np.asarray([linear_fit(sfun(dj-_x), counts, ecounts, chionly=True) for _x in x ])
    a  = np.polyfit(x,y,deg=2)
    # y = a[0]*x^2 + a[1]*x + a[2]
    # we need to know x for y_min+1.0
    # x0 = -b/2a
    # dx ~ sqrt(delta)/2a
    a[2]  = a[2] - 1.0 - y[1]
    if a[0]==0:
        return np.nan
    # x_min = -b/2a - parabola vertex (in our case, the chi2 minimum)
    # we need to know x for y_min+1
    a[2]  = a[2] - 1.0 - np.polyval(a, -a[1]/(2*a[0]))
    # dx    = sqrt(delta)/2a
    delta = a[1]**2-4*a[0]*a[2]
    if delta<0:
        return np.nan
    return  0.5*np.sqrt(delta)/a[0]
    return np.sqrt(delta)/a[0]/2



+15 −12
Original line number Diff line number Diff line
@@ -194,20 +194,23 @@ def plot_example(filename, plt, dj0=None, signum=1):
    ax2.errorbar(dj/MICRO, counts, yerr=ecounts, fmt='ko')
    ax2.plot(djf/MICRO, efit,  'r-')
    if init_dj0 is not None:
        ax2.axvline(init_dj0/MICRO, color='black', ls=':') # phase input
    #ax2.axvline((dj0-edj0)/MICRO , color='red' , ls='--' ) # phase fit
    ax2.axvline(dj0/MICRO , color='red' , lw=2 ) # phase fit
    #ax2.axvline((dj0+edj0)/MICRO , color='red' , ls='--' ) # phase fit
    ax2.axvline(dj0_pt/MICRO , color='magenta' , lw=1 ) # phase fit
    ax2.axhline(res['average'][0], color='green')  # average
    ax2.axhline(res['average'][0]-res['amplitude'][0], color='green', ls='--') # average - ampl
    ax2.axhline(res['average'][0]+res['amplitude'][0], color='green', ls='--') # average + ampl
        ax2.axvline(init_dj0/MICRO, color='black', lw=1.0, ls=':') # phase input
    ax2.axvline(dj0_pt/MICRO , color='magenta' , lw=0.5 )    # phase fit 4pt
    #
    ax2.axvline((dj0-edj0)/MICRO , color='red' , lw=0.5, ls='--' ) #
    ax2.axvline(dj0/MICRO ,        color='red' , lw=1.0, ls='-'  ) # phase fit
    ax2.axvline((dj0+edj0)/MICRO , color='red' , lw=0.5, ls='--' ) #
    #
    ax2.axhline(res['average'][0], color='green', lw=1.0)  # average
    ax2.axhline(res['average'][0]-res['amplitude'][0], color='green', lw=1.0, ls='--') # average - ampl
    ax2.axhline(res['average'][0]+res['amplitude'][0], color='green', lw=1.0, ls='--') # average + ampl
    ax2.set_ylabel('counts')
    ax2.grid(True)

    ax2a = ax2.twinx()
    ax2a.set_ylabel(r'$\chi^2$', color='b')
    ax2a.plot(djf/MICRO, chif, 'b-')
    ax2a.plot(dj /MICRO, chi2, 'b.')
    ax2a.set_ylabel(r'$\chi^2$', color='b', )
    ax2a.plot(djf/MICRO, chif, 'b-', lw=1.0, ls='--')
    ax2a.plot(dj /MICRO, chi2, 'bo', ms=2.0)



@@ -226,7 +229,7 @@ if __name__ == "__main__":
    if len(sys.argv)>2:
        signum = np.sign(float(sys.argv[2]))
    else:
        signum = 1.0
        signum = 0.0
    plt.figure()
    plot_example(filename, plt, dj0=0, signum=signum)
    plt.show()