Computation of 95 % CL limits with statistical errors
Source code name: "stat_limits_discovery.py"
Programming language: Python
Topic: Statistics/Limits
DMelt Version 1. Last modified: 06/04/2015. License: Pro
https://datamelt.org/code/cache/stat_limits_discovery_7860.py
To run this script using the DMelt IDE, copy the above URL link to the menu [File]→[Read script from URL] of the DMelt IDE.


 
"""
1-CLb gives a probability that a background-only fluctuates and gives 
expected signal. If you will see 0 value - increase the number of MC simu

for 5-sigma discovery, 1-CLb should be less than  2.87*10-7
for 3-sigma discovery, 1-CLb should be less than 1.3 *10-3

Authors:  Christophe.Delaere@cern.ch (C++ example) and Sergei Chekanov (Java)
"""

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

c1 = HPlot("Canvas",600,400)
c1.setGTitle("Prob of background fluctuation")
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.draw(signal)
c1.draw(background)
c1.draw(sigback)
c1.draw(data)


print "Wait .. Calculating .."
datasource = DataSource(signal, background,data)
climit = CLimits(datasource,100000) 
confidence = climit.getLimit()
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()
# probability that background only fluctuates to this peak is
print "Prob for background fluctuation", 1-confidence.getCLb()


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


You see the box below because you did not login.