Page MenuHomeHEPForge

tools.py
No OneTemporary

tools.py

#!/usr/bin/env python
from numpy import random, zeros,array,log, exp, ceil,array
import os, sys, yoda
import pylab
def smearsingle(A,nsamples,default):
return array([default if p.y <=0 or p.yErrAvg<=0 else random.normal(p.y, p.yErrAvg, nsamples) for p in A.points ])
def smearYODA(BG,S,nsamples):
z=zeros(nsamples)
return smearsingle(BG,nsamples,z),smearsingle(BG,nsamples,z),smearsingle(S,nsamples,z)
factmap={}
def LLR(s,b1,b2,d):
global factmap
fact=0.
if int(ceil(d)) not in factmap.keys():
factmap[int(ceil(d))]=sum([log(i) for i in xrange(2, int(ceil(d)))])
fact=factmap[int(ceil(d))]
sb = -1.*(s+b1)+d*log(s+b1) -fact
a1 = -1.*b1 +d*log(b1)-fact
a2 = -1.*b2 +d*log(b2)-fact
return -2.*(sb-a2), -2*(a1-a2)
OBS=sys.argv[1]
BG=yoda.readYODA(sys.argv[2])[OBS].mkScatter()
SIG=yoda.readYODA(sys.argv[3])[OBS].mkScatter()
D=yoda.readYODA(sys.argv[4])["/REF%s"%OBS].mkScatter()
LLROBS=0
for num, p in enumerate(D.points):
if p.y>0 and BG.points[num].y>0 and (BG.points[num].y + SIG.points[num].y)>0:
LLROBS += LLR(SIG.points[num].y, BG.points[num].y,BG.points[num].y, p.y)[0]
aLLR,clspb,clspb2,clb,N,i=[],0.,0.,0.,0,1000
f_BG,f_BG1,f_SIG = smearYODA(BG,SIG, 1000)
while N<10000:
if i==1000:
i=0
f_BG,f_BG1,f_SIG = smearYODA(BG,SIG, 1000)
llr,llrb,llrb2 = 0.,0.,0.
for num, p in enumerate(D.points):
if p.y>0 and f_SIG[num][i]>0 and f_BG[num][i]>0 and f_BG1[num][i]>0:
t1,t2 =LLR(f_SIG[num][i], f_BG[num][i], f_BG1[num][i], p.y)
llr += t1
llrb += t2
i+=1
N+=1
aLLR.append([llr, llrb ])
if (llr>LLROBS):
clspb+=1.
if (llrb>LLROBS):
clb+=1.
print "clspb",clspb/N, " clb",clb/N
X=array(aLLR)
pylab.hist(X[:,1],bins=100, label="BG only",histtype="step")
pylab.hist(X[:,0],bins=100, label="S+B",histtype="step")
pylab.axvline(LLROBS, label="LLROBS")
pylab.legend()
pylab.show()

File Metadata

Mime Type
text/x-python
Expires
Sat, Dec 21, 1:36 PM (19 h, 36 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4019017
Default Alt Text
tools.py (1 KB)

Event Timeline