AnalysisFramework/limits/run2011limits.C

#include "roostats_cl95.C"

int getLimits( double lumi, double elumi, double bkg, double ebkg, double elm4, double eelm4, double elm8, double eelm8, int obs  ) {
//  gROOT->LoadMacro("roostats_cl95.C+");

  int niter = 100;
  
  double evt         = roostats_cl95(1.0,       0.001,  1.0, 0.001, bkg, ebkg, obs, kFALSE, 1);
  LimitResult expevt = roostats_clm( 1.0,       0.001,  1.0, 0.001, bkg, ebkg, niter, 1);

  double      lm4    = roostats_cl95(lumi, lumi*elumi, elm4, eelm4, bkg, ebkg, obs, kFALSE, 1);
  LimitResult exp4   = roostats_clm( lumi, lumi*elumi, elm4, eelm4, bkg, ebkg, niter, 1 );

  double      lm8    = roostats_cl95(lumi, lumi*elumi, elm8, eelm8, bkg, ebkg, obs, kFALSE, 1);
  LimitResult exp8   = roostats_clm( lumi, lumi*elumi, elm8, eelm8, bkg, ebkg, niter, 1 );

  
  std::cout << "UL on #(events) = " << evt << ", expected = " << expevt.GetExpectedLimit();
  std::cout <<      "+" << expevt.GetOneSigmaHighRange()-expevt.GetExpectedLimit();
  std::cout <<      "-" << expevt.GetExpectedLimit()-expevt.GetOneSigmaLowRange()  <<  std::endl;
  std::cout << "UL on LM4 SxBRxA= " << lm4 << ", expected = " << exp4.GetExpectedLimit();
  std::cout <<      "+" << exp4.GetOneSigmaHighRange()-exp4.GetExpectedLimit();
  std::cout <<      "-" << exp4.GetExpectedLimit()-exp4.GetOneSigmaLowRange() << std::endl; 
  std::cout << "UL on LM8 SxBRxA= " << lm8 << ", expected = " << exp8.GetExpectedLimit(); 
  std::cout <<      "+" << exp8.GetOneSigmaHighRange()-exp8.GetExpectedLimit();
  std::cout <<      "-" << exp8.GetExpectedLimit()-exp8.GetOneSigmaLowRange() << std::endl;

  return 0;
}

//_____________________________________________________________________________________________
int run2011_3jets_50(double elumi = 0.04 ) {

double lumi = 191.;
double bkg = 24;
double ebkg = sqrt(6*6+1.4*1.4+2.4*2.4);
int obs = 20;

double elm4 = 0.919;
double eelm4s = sqrt(4*4+2.*2.+1.0*1.0+2.*2.+1.4*1.4)/100.;
double eelm4t = (0.7*0.7)/100.;
double elm8 = 0.885;
double eelm8s = sqrt(4*4+2.*2.+1.4*1.4+2.*2.+1.5*1.5)/100.;
double eelm8t = (0.9*0.9)/100.;

std::cout << "LM4 syst error: " << eelm4s << std::endl;
std::cout << "LM8 syst error: " << eelm8s << std::endl;
  
double eelm4 = sqrt(eelm4s*eelm4s+eelm4t*eelm4t+0.01*elm4*elm4); // Add 10% theo. unc.
double eelm8 = sqrt(eelm8s*eelm8s+eelm8t*eelm8t+0.01*elm8*elm8);

getLimits( lumi, elumi, bkg, ebkg, elm4, eelm4, elm8, eelm8, obs );

return 0;

}

//_____________________________________________________________________________________________
int run2011_3jets_100(double elumi = 0.04) {

double lumi = 191.;
double bkg = 8;
double ebkg = sqrt(4*4+0.1*0.1+0.4*0.4);
int obs = 6;

double elm4 = 0.904;
double eelm4s = sqrt(4*4+2.*2.+1.0*1.0+4.*4.+4.3*4.3)/100.;
double eelm4t = (0.9*0.9)/100.;
double elm8 = 0.853;
double eelm8s = sqrt(4*4+2.*2.+2.0*2.0+4.*4.+4.5*4.5)/100.;
double eelm8t = (1.1*1.1)/100.;

std::cout << "LM4 syst erro: " << eelm4s << std::endl;
std::cout << "LM8 syst erro: " << eelm8s << std::endl;
  
double eelm4 = sqrt(eelm4s*eelm4s+eelm4t*eelm4t+0.01*elm4*elm4); // Add 10% theo. unc.
double eelm8 = sqrt(eelm8s*eelm8s+eelm8t*eelm8t+0.01*elm8*elm8);

getLimits( lumi, elumi, bkg, ebkg, elm4, eelm4, elm8, eelm8, obs );

return 0;

}

//___________________________________________________________________________________________
int main(void) {

  run2011_3jets_50();
  run2011_3jets_100();
//  run2011_3jets_50(0.1);
//  run2011_3jets_100(0.1);
 
  return 0;
}
Revision:	1.1
Committed:	Fri Jun 24 07:48:47 2011 UTC (13 years, 10 months ago) by fronga
Content type:	text/plain
Branch:	MAIN
CVS Tags:	cbaf_4_98ifb_paper, cbaf_4p7ifb, Honeypot, cbaf_2p1ifb, HEAD
Log Message:	Adding macros used to compute limits.
#	User	Rev	Content
1	fronga	1.1	#include "roostats_cl95.C"
2
3			int getLimits( double lumi, double elumi, double bkg, double ebkg, double elm4, double eelm4, double elm8, double eelm8, int obs ) {
4			// gROOT->LoadMacro("roostats_cl95.C+");
5
6			int niter = 100;
7
8			double evt = roostats_cl95(1.0, 0.001, 1.0, 0.001, bkg, ebkg, obs, kFALSE, 1);
9			LimitResult expevt = roostats_clm( 1.0, 0.001, 1.0, 0.001, bkg, ebkg, niter, 1);
10
11			double lm4 = roostats_cl95(lumi, lumi*elumi, elm4, eelm4, bkg, ebkg, obs, kFALSE, 1);
12			LimitResult exp4 = roostats_clm( lumi, lumi*elumi, elm4, eelm4, bkg, ebkg, niter, 1 );
13
14			double lm8 = roostats_cl95(lumi, lumi*elumi, elm8, eelm8, bkg, ebkg, obs, kFALSE, 1);
15			LimitResult exp8 = roostats_clm( lumi, lumi*elumi, elm8, eelm8, bkg, ebkg, niter, 1 );
16
17
18			std::cout << "UL on #(events) = " << evt << ", expected = " << expevt.GetExpectedLimit();
19			std::cout << "+" << expevt.GetOneSigmaHighRange()-expevt.GetExpectedLimit();
20			std::cout << "-" << expevt.GetExpectedLimit()-expevt.GetOneSigmaLowRange() << std::endl;
21			std::cout << "UL on LM4 SxBRxA= " << lm4 << ", expected = " << exp4.GetExpectedLimit();
22			std::cout << "+" << exp4.GetOneSigmaHighRange()-exp4.GetExpectedLimit();
23			std::cout << "-" << exp4.GetExpectedLimit()-exp4.GetOneSigmaLowRange() << std::endl;
24			std::cout << "UL on LM8 SxBRxA= " << lm8 << ", expected = " << exp8.GetExpectedLimit();
25			std::cout << "+" << exp8.GetOneSigmaHighRange()-exp8.GetExpectedLimit();
26			std::cout << "-" << exp8.GetExpectedLimit()-exp8.GetOneSigmaLowRange() << std::endl;
27
28			return 0;
29			}
30
31			//_____________________________________________________________________________________________
32			int run2011_3jets_50(double elumi = 0.04 ) {
33
34			double lumi = 191.;
35			double bkg = 24;
36			double ebkg = sqrt(66+1.41.4+2.4*2.4);
37			int obs = 20;
38
39			double elm4 = 0.919;
40			double eelm4s = sqrt(44+2.2.+1.01.0+2.2.+1.4*1.4)/100.;
41			double eelm4t = (0.7*0.7)/100.;
42			double elm8 = 0.885;
43			double eelm8s = sqrt(44+2.2.+1.41.4+2.2.+1.5*1.5)/100.;
44			double eelm8t = (0.9*0.9)/100.;
45
46			std::cout << "LM4 syst error: " << eelm4s << std::endl;
47			std::cout << "LM8 syst error: " << eelm8s << std::endl;
48
49			double eelm4 = sqrt(eelm4seelm4s+eelm4teelm4t+0.01elm4elm4); // Add 10% theo. unc.
50			double eelm8 = sqrt(eelm8seelm8s+eelm8teelm8t+0.01elm8elm8);
51
52			getLimits( lumi, elumi, bkg, ebkg, elm4, eelm4, elm8, eelm8, obs );
53
54			return 0;
55
56			}
57
58			//_____________________________________________________________________________________________
59			int run2011_3jets_100(double elumi = 0.04) {
60
61			double lumi = 191.;
62			double bkg = 8;
63			double ebkg = sqrt(44+0.10.1+0.4*0.4);
64			int obs = 6;
65
66			double elm4 = 0.904;
67			double eelm4s = sqrt(44+2.2.+1.01.0+4.4.+4.3*4.3)/100.;
68			double eelm4t = (0.9*0.9)/100.;
69			double elm8 = 0.853;
70			double eelm8s = sqrt(44+2.2.+2.02.0+4.4.+4.5*4.5)/100.;
71			double eelm8t = (1.1*1.1)/100.;
72
73			std::cout << "LM4 syst erro: " << eelm4s << std::endl;
74			std::cout << "LM8 syst erro: " << eelm8s << std::endl;
75
76			double eelm4 = sqrt(eelm4seelm4s+eelm4teelm4t+0.01elm4elm4); // Add 10% theo. unc.
77			double eelm8 = sqrt(eelm8seelm8s+eelm8teelm8t+0.01elm8elm8);
78
79			getLimits( lumi, elumi, bkg, ebkg, elm4, eelm4, elm8, eelm8, obs );
80
81			return 0;
82
83			}
84
85			//___________________________________________________________________________________________
86			int main(void) {
87
88			run2011_3jets_50();
89			run2011_3jets_100();
90			// run2011_3jets_50(0.1);
91			// run2011_3jets_100(0.1);
92
93			return 0;
94			}