# This program demonstrates the computation of 95 % C.L. limits.
# Statistical errors are included in the treatment.
# As systematic, we will consider 2 data sources.  
#
# 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:  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 "--> Signal hypothesis is excluded at level (%) ", (1-confidence.getCLs())*100.


c1 = HPlot("Canvas")
c1.setGTitle("Computation of 95 % C.L. limits")
c1.setRange(-4,4,0.0,120)
c1.visible()
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")

# second data ses (after some systematic variation) 
data1     = H1D("Systematic1",30,-4.0,4.0)
data1.setColor(Color.red)
data1.setStyle("p")

# third data set (after some systematic variation)
data2     = H1D("Systematic2",30,-4.0,4.0)
data2.setColor(Color.green)
data2.setStyle("p")



r=Random()
for i in range(25000):
      background.fill(r.nextGaussian(),0.02)
      signal.fill(1+0.2*r.nextGaussian(),0.001)
for i in range(500):
      data.fill(r.nextGaussian())    
for i in range(480):
      data1.fill(1.05*r.nextGaussian())
for i in range(530):
      data2.fill(0.95*r.nextGaussian())


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)
c1.draw(data1)
c1.draw(data2)

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


# assume two data sources with systematic variations 
print "Wait.."
d=DataSource()
d.addChannel(signal,background,data,0,0,0,"data")
d.addChannel(signal,background,data1,0,0,0,"systematics1")
d.addChannel(signal,background,data2,0,0,0,"systematics1")
climit=CLimits(d, 100000)
printCL("Data with with systematics", climit.getLimit())



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