``` 1 # Statistics. Computation of 95 % C.L. limits using stat. and systematical errors. 2 # This program demonstrates the computation of 95 % C.L. limits. 3 # Statistical errors are included in the treatment. 4 # As systematic, we will consider 2 data sources. 5 # 6 # Signal hypothesis is excluded at the 95% CL if CLs = 0.05 7 # and at more than the 95% CL if CLs < 0.05, assuming that signal is present 8 # 9 # Authors: Sergei Chekanov 10 11 from java.awt import Color,Font 12 from java.util import Random 13 from jhplot import * 14 from jhpro.stat import * 15 16 17 # just convinient print 18 def printCL(mess,confidence ): 19 print "\n------- "+mess+" -----------" 20 print "CLs : " ,confidence.getCLs() 21 print "CLb : " ,confidence.getCLb() 22 print "CLsb : " ,confidence.getCLsb() 23 print "expected : " ,confidence.getExpectedCLs_b() 24 print "expected : " ,confidence.getExpectedCLb_b() 25 print "expected : " ,confidence.getExpectedCLb_b() 26 print "--> Signal hypothesis is excluded at level (%) ", (1-confidence.getCLs())*100. 27 28 29 c1 = HPlot("Canvas") 30 c1.setGTitle("Computation of 95 % C.L. limits") 31 c1.setRange(-4,4,0.0,120) 32 c1.visible() 33 c1.setNameX("Variable") 34 c1.setNameY("Events") 35 36 37 # set 38 background = H1D("Background",30,-4.0,4.0) 39 background.setColor(Color.green) 40 background.setFill(1) 41 background.setFillColor(Color.green) 42 background.setErrAll(0) 43 44 signal = H1D("Signal",30,-4.0,4.0) 45 signal.setFill(1) 46 signal.setFillColor(Color.red) 47 signal.setColor(Color.red) 48 49 data = H1D("Data",30,-4.0,4.0) 50 data.setColor(Color.black) 51 data.setStyle("p") 52 53 # second data ses (after some systematic variation) 54 data1 = H1D("Systematic1",30,-4.0,4.0) 55 data1.setColor(Color.red) 56 data1.setStyle("p") 57 58 # third data set (after some systematic variation) 59 data2 = H1D("Systematic2",30,-4.0,4.0) 60 data2.setColor(Color.green) 61 data2.setStyle("p") 62 63 64 65 r=Random() 66 for i in range(25000): 67 background.fill(r.nextGaussian(),0.02) 68 signal.fill(1+0.2*r.nextGaussian(),0.001) 69 for i in range(500): 70 data.fill(r.nextGaussian()) 71 for i in range(480): 72 data1.fill(1.05*r.nextGaussian()) 73 for i in range(530): 74 data2.fill(0.95*r.nextGaussian()) 75 76 77 sigback=background.oper(signal,"Signal+Background","+") 78 sigback.setErrAll(0) 79 80 c1.cd(1,1) 81 c1.draw(signal) 82 c1.draw(background) 83 c1.draw(sigback) 84 c1.draw(data) 85 c1.draw(data1) 86 c1.draw(data2) 87 88 # assume first that we do not have any systematics, only statistical uncertanties 89 print "Wait.." 90 d=DataSource() 91 d.addChannel(signal,background,data) 92 climit=CLimits(d, 100000) 93 printCL("No systematics", climit.getLimit()) 94 95 96 # assume two data sources with systematic variations 97 print "Wait.." 98 d=DataSource() 99 d.addChannel(signal,background,data,0,0,0,"data") 100 d.addChannel(signal,background,data1,0,0,0,"systematics1") 101 d.addChannel(signal,background,data2,0,0,0,"systematics1") 102 climit=CLimits(d, 100000) 103 printCL("Data with with systematics", climit.getLimit()) 104 105 106 107 # export to some image (png,eps,pdf,jpeg...) 108 # c1.export(Editor.DocMasterName()+".png"); 109 # edit the image 110 # IEditor(Editor.DocMasterName()+".png"); ### © jHepWork. S.Chekanov ### ```