# using multipass calculations. We calculate the upper and lower
# value for 90% confidence levels.
# 
# Based on Computer Physics Communications 149 (2002) 97, R.Barlow.

from jhpro.stat.limit import *
from java.awt import Color
from jhplot import *

# build a singleton
c1 = HPlot("Canvas")
c1.visible()
c1.setAutoRange()
c1.setGTitle("Limit calculations")
c1.setNameX("Values")
c1.setNameY("Probability")

upper=P1D("Upper values")
lower=P1D("Lower values")

# build calculations. No need to define limit here
sensitivity=1
sigma=0;
limitguess=0;
MCevents=10000;
s=StatConfidence(sensitivity,sigma,limitguess,MCevents);

# build an experimental data
# - 200 observed events
# - 100 expected background
# - 20  error on expected background 
exp1=ExpData(200,100,20)
s.addData(exp1)

# second experiment with 150 observed events
# exp2=ExpData(150,140,20)
# s.addData(exp2)


# iterate ove different guesses
limitGuess=0
print "-> Calculating lower limit for 90% confidence level using Cousins+Highland .."
for i in range(20,150):
    s.setLimitGuess(i)
    s.run(0)
    lower.add(i,s.getProbabilityCH())
    if (s.getProbabilityCH()>0.89 and s.getProbabilityCH()<0.91):
       limitGuess=i
print "Lower Value for 90% limit =", limitGuess  

print "-> Calculating upper Limit for 90% confidence level using Cousins+Highland:"
limitGuess=0
for i in range(20,150):
    s.setLimitGuess(i)
    s.run(1)
    upper.add(i,s.getProbabilityCH())
    if (s.getProbabilityCH()>0.09 and s.getProbabilityCH()<0.11):
       limitGuess=i
 
lower.setColor(Color.blue)     
c1.draw(lower)
c1.draw(upper)

print "Upper Value for 90% limit =", limitGuess