# This program demonstrates the computation of 95 % C.L. limits.
# Statistical errors are included in the treatment.  
#
# Signal hypothesis is excluded at the 95% CL if CLs = 0.05
# and at more than the 95% CL if CLs < 0.05, assuming that signal is present
#
# Authors:  Christophe.Delaere@cern.ch on 21.08.02 and Sergei Chekanov

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


# just convinient print
def printCL(mess,confidence ):
  print "\n------- "+mess+" -----------"
  print "CLs    : " ,confidence.getCLs()
  print "CLb    : " ,confidence.getCLb()
  print "CLsb    : " ,confidence.getCLsb()
  print "expected     : " ,confidence.getExpectedCLs_b()
  print "expected     : " ,confidence.getExpectedCLb_b()
  print "expected     : " ,confidence.getExpectedCLb_b()
  print "Prob for background fluctuation", 1-confidence.getCLb()


c1 = HPlot("Canvas")
c1.setGTitle("Statistical significance of a signal")
c1.visible()
c1.setRange(-4,4,0.0,300)
c1.setNameX("Variable")
c1.setNameY("Events")

# set
background = H1D("Background",30,-4.0,4.0)
background.setColor(Color.green)
background.setFill(1)
background.setFillColor(Color.green)
background.setErrAll(0)

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

data     = H1D("Data",30,-4.0,4.0)
data.setColor(Color.black)
data.setStyle("p")

r=Random()
for i in range(1000):
      background.fill(r.nextGaussian(),1.0)
      data.fill(r.nextGaussian(),1.0)

for i in range(200):
      signal.fill(1+0.2*r.nextGaussian(),1.0)

for i in range(70):
      data.fill(1+0.2*r.nextGaussian(),1.0)


sigback=background.oper(signal,"Signal+Background","+")
sigback.setErrAll(0)

c1.cd(1,1)
c1.draw(signal)
c1.draw(background)
c1.draw(sigback)
c1.draw(data)

# add systematic errors from 2 sources
# err_names=["source1","source2"]
# err_b=[0.05, 0.0] # error on background is 5% (source1) and 0% (source2)
# err_s=[0.0, 0.01] # error on signal is  0% (source1) and 1% (source2)

print "Wait.. Calculating.."
# assume first that we do not have any systematics, only statistical uncertanties
d=DataSource()
d.addChannel(signal,background,data)
climit=CLimits(d, 100000)
printCL("No systematics", climit.getLimit())

print "Wait.. Calculating.."
# assume data source has +/-10% systematics 
d=DataSource()
d.addChannel(signal,background,data,0,0,0.1,"source")
d.addChannel(signal,background,data,0,0,-0.1,"source")
climit=CLimits(d, 100000)
printCL("Data with +/-5% systematics", climit.getLimit())

print "Wait.. Calculating.."
# assume data source has +/-2% systematics and background has 10% systematics 
d=DataSource()
d.addChannel(signal,background,data,0,0.1,0.02,"source")
d.addChannel(signal,background,data,0,-0.1,-0.02,"source")
climit=CLimits(d, 100000)
printCL("Data with +/-5% and background with +/-4% systematics", climit.getLimit())




# export to some image (png,eps,pdf,jpeg...)
# c1.export(Editor.DocMasterName()+".png");
# edit the image
# IEditor(Editor.DocMasterName()+".png");