from java.awt import Color
from java.util import Random 
from jhplot import *
from jhpro.stat import *

data=H1D("Fake data",30,0.0,4.0)
data.setColor(Color.black); data.setStyle("p")

backg=H1D("Backg",30,0.0,4.0)
backg.setColor(Color.green); backg.setFill(1); backg.setErrAll(0)
backg.setFillColor(Color.green)

signal= H1D("Signal",30,0.0,4.0)
signal.setFill(1); signal.setFillColor(Color.red)
signal.setColor(Color.red)

r=Random(2L) # make random numbers reproducible
Max=10000    # events in empirical data  
strength=10  # initial signal strength at start of the scan 

# fill pseudo-data and some established theory (-10 events, for example) 
for i in xrange(Max):
    data.fill(abs(r.nextGaussian()))
    backg.fill(abs(r.nextGaussian()))  #  predictions have lower normalization

for ev in range(100):  # scan through the parameters of new theory (events) 
   strength=strength+5 # increment Nr of signal events by 10 
   signal.clear()
   for i in xrange(strength): 
       signal.fill(abs(2+0.3*r.nextGaussian()))
   events=signal.sumAllBinHeights()
   datasource = DataSource(signal, backg, data)
   print "Calculating ..  Nr of signal events:",events
   climit = CLimits(datasource,50000)
   conf = climit.getLimit()
   valCLs=conf.getCLs(); print "  CLs = ",valCLs 
   if  valCLs == float("inf"): continue 
   exc=int((1-valCLs)*100) 
   print "  Signal is excluded at level ",exc,"%"  
   del climit
   if (valCLs<0.05):
       print "Done! Nr ",events," events needed to exclude signal with ",exc,"% CL" 
       break

c1 = HPlot()  # show the results
c1.visible(); c1.setRange(0,3,0.0,1500)
sigback=backg.oper(signal,"Signal+Backg","+")
sigback.setColor(Color.blue)
sigback.setErrAll(0)
c1.draw([backg,signal,sigback,data])
c1.export("limit.pdf")