SusyScan/Limits/Limit.cc

//User
#include "cls.h"
#include "ConfigFile.h"
#include "table.h"

//system
#include <string>
#include <iostream>
#include <cmath>
#include <cstdlib>

//ROOT
#include "TH1.h"

using namespace std;


void Limit(int argc, char *argv[])
{
  if (argc<1){ 
    std::cerr<<"Error: Expected '"<<argv[0]<<" <limit config>'!"<<std::endl;
    exit(1);
  }
  
  for (int file=1; file<argc; ++file){   
    std::cout << "opening: "<<argv[file]<<std::endl;
    ConfigFile config(argv[file]);

    int    ndata = config.read<int>("data");
    double nback  = config.read<double>("background");
    double bkgUncert = config.read<double>("background.uncertainty");
    double signal  = config.read<double>("signal.LO");
    double sigUncert = config.read<double>("signal.LO.uncertainty");
    double sigSignal = config.read<double>("signal.LO.signalregion.IsoMuon");
    double sigSigTau  = config.read<double>("signal.LO.signalregion.Tau");
    double signalNLO  = config.read<double>("signal.NLO");
    double sigUncertNLO = config.read<double>("signal.NLO.uncertainty");
    double sigSignalNLO = config.read<double>("signal.NLO.signalregion.IsoMuon");
    double sigSigTauNLO = config.read<double>("signal.NLO.signalregion.Tau");
    double bkgControl = config.read<double>("background.controlregion.IsoMuon");
    double bkgSignal  = config.read<double>("background.signalregion.IsoMuon");
  
    double xsec = config.read<double>("Xsection");
    double kfactor = config.read<double>("signal.kFactor");
    
    //bkgUncert += sigSignal * sigUncert/signal;
    
    //config.add("data", ndata);
    //config.add("background", nback);
    //config.add("background.uncertainty", bkgUncert);

    char * n = new char[32];
    sprintf(n,"bkgd%d",file); TH1 * bgd  = new TH1F(n, "",1,0,1);
    sprintf(n,"bkgd_NoSig%d",file); TH1 * bgd_NoSig  = new TH1F(n, "",1,0,1);
    sprintf(n,"data%d",file); TH1 * data = new TH1F(n, "",1,0,1);
    sprintf(n,"sig%d", file); TH1 * sig  = new TH1F(n,"",1,0,1);
    data->SetBinContent(1,ndata);
    data->SetBinError(1,sqrt(ndata));
    bgd_NoSig->SetBinContent(1,nback);
    bgd_NoSig->SetBinError(1,bkgUncert);
    double bLO=(nback - sigSignal -sigSigTau<0?0:nback - sigSignal -sigSigTau);
    bgd->SetBinContent(1,bLO);
    bgd->SetBinError(1,bkgUncert);
    sig->SetBinContent(1,signal);
    sig->SetBinError(1,sigUncert);

    cls limit_sys("cls_limit.ps",
          sig,  // signal with total syst. and stat. errors
          bgd,  // with total syst. and stat. errors
          data,
          bgd_NoSig); 
    limit_sys.SetStat(true); //Consider Bin-Errors (containing the full syst. uncertainties)
    limit_sys.SetNMC(50000); //number of pseudo-experiments for each CL calculation
    //limit_sys.Draw();

    sprintf(n,"bkgdNLO%d",file); TH1 * bgdNLO  = new TH1F(n, "",1,0,1);
    sprintf(n,"bkgdNoSigNLO%d",file); TH1 * bgdNoSigNLO  = new TH1F(n, "",1,0,1);
    sprintf(n,"dataNLO%d",file); TH1 * dataNLO = new TH1F(n, "",1,0,1);
    sprintf(n,"sigNLO%d", file); TH1 * sigNLO  = new TH1F(n,"",1,0,1);
    
    double bNLO=(nback - sigSignalNLO -sigSigTauNLO<0?0:nback - sigSignalNLO - sigSigTauNLO);

    dataNLO->SetBinContent(1,ndata);
    dataNLO->SetBinError(1,sqrt(ndata));
    bgdNLO->SetBinContent(1,bNLO);
    bgdNLO->SetBinError(1,bkgUncert);
    bgdNoSigNLO->SetBinContent(1,nback);
    bgdNoSigNLO->SetBinError(1,bkgUncert);
    sigNLO->SetBinContent(1,signalNLO);
    sigNLO->SetBinError(1,sigUncertNLO);

    cls limit_sys_NLO("cls_limit_NLO.ps",
          sigNLO,  // signal with total syst. and stat. errors
          bgdNLO,  // with total syst. and stat. errors
          dataNLO,
          bgdNoSigNLO); 
    limit_sys_NLO.SetStat(true); //Consider Bin-Errors (containing the full syst. uncertainties)
    limit_sys_NLO.SetNMC(50000); //number of pseudo-experiments for each CL calculation
    //limit_sys.Draw();

    try {
      //limit_sys.WriteResult( &config, xsec );
      limit_sys.WriteConfidence( &config, "LO." );
      limit_sys_NLO.WriteConfidence( &config, "NLO." );
    }
    catch(exception& e){
      cout << "Exception catched: " << e.what();
    }

    //write stuff:  
    time_t rawtime;
    struct tm * timeinfo;
    time ( &rawtime );
    timeinfo = localtime ( &rawtime );
    ofstream ofile;
    ofile.open (argv[file]);
    if (ofile.good())
      std::cout << "writing '"<<argv[file]<<"'"<<std::endl;
    else if ( (ofile.rdstate() & ifstream::failbit ) != 0 )
      std::cerr << "ERROR opening '"<<argv[file]<<"'! Does the directory exist?"<<std::endl;
    
    ofile << config.getComment() << asctime (timeinfo) 
          << config.getComment()<< "\n"
          << config;

    std::cout<<"d:"<<ndata<<", b:"<<bNLO<<", s:"<<signalNLO<<std::endl;

  }
}


int main(int argc, char *argv[])
{
  Limit(argc,argv);
}
Revision:	1.6
Committed:	Sun Mar 20 12:17:58 2011 UTC (14 years, 1 month ago) by auterman
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.5:	+17 -5 lines
Log Message:	Modification of CLs: Adding functunality to integrate until a given Test-statistic value, i.e. until a value determined by a no-signal hypotheses for expected observed
#	User	Rev	Content
1	auterman	1.1	//User
2			#include "cls.h"
3			#include "ConfigFile.h"
4			#include "table.h"
5
6			//system
7			#include <string>
8			#include <iostream>
9			#include <cmath>
10	auterman	1.2	#include <cstdlib>
11	auterman	1.1
12			//ROOT
13			#include "TH1.h"
14
15			using namespace std;
16
17
18			void Limit(int argc, char *argv[])
19			{
20			if (argc<1){
21			std::cerr<<"Error: Expected '"<<argv[0]<<" <limit config>'!"<<std::endl;
22			exit(1);
23			}
24
25			for (int file=1; file<argc; ++file){
26			std::cout << "opening: "<<argv[file]<<std::endl;
27			ConfigFile config(argv[file]);
28	auterman	1.2
29	auterman	1.1	int ndata = config.read<int>("data");
30			double nback = config.read<double>("background");
31			double bkgUncert = config.read<double>("background.uncertainty");
32	auterman	1.5	double signal = config.read<double>("signal.LO");
33			double sigUncert = config.read<double>("signal.LO.uncertainty");
34			double sigSignal = config.read<double>("signal.LO.signalregion.IsoMuon");
35			double sigSigTau = config.read<double>("signal.LO.signalregion.Tau");
36			double signalNLO = config.read<double>("signal.NLO");
37			double sigUncertNLO = config.read<double>("signal.NLO.uncertainty");
38			double sigSignalNLO = config.read<double>("signal.NLO.signalregion.IsoMuon");
39			double sigSigTauNLO = config.read<double>("signal.NLO.signalregion.Tau");
40	auterman	1.4	double bkgControl = config.read<double>("background.controlregion.IsoMuon");
41			double bkgSignal = config.read<double>("background.signalregion.IsoMuon");
42	auterman	1.5
43			double xsec = config.read<double>("Xsection");
44			double kfactor = config.read<double>("signal.kFactor");
45	auterman	1.4
46			//bkgUncert += sigSignal * sigUncert/signal;
47
48	auterman	1.2	//config.add("data", ndata);
49			//config.add("background", nback);
50			//config.add("background.uncertainty", bkgUncert);
51
52	auterman	1.1	char * n = new char[32];
53			sprintf(n,"bkgd%d",file); TH1 * bgd = new TH1F(n, "",1,0,1);
54	auterman	1.6	sprintf(n,"bkgd_NoSig%d",file); TH1 * bgd_NoSig = new TH1F(n, "",1,0,1);
55	auterman	1.1	sprintf(n,"data%d",file); TH1 * data = new TH1F(n, "",1,0,1);
56			sprintf(n,"sig%d", file); TH1 * sig = new TH1F(n,"",1,0,1);
57			data->SetBinContent(1,ndata);
58			data->SetBinError(1,sqrt(ndata));
59	auterman	1.6	bgd_NoSig->SetBinContent(1,nback);
60			bgd_NoSig->SetBinError(1,bkgUncert);
61	auterman	1.5	double bLO=(nback - sigSignal -sigSigTau<0?0:nback - sigSignal -sigSigTau);
62	auterman	1.6	bgd->SetBinContent(1,bLO);
63	auterman	1.1	bgd->SetBinError(1,bkgUncert);
64			sig->SetBinContent(1,signal);
65			sig->SetBinError(1,sigUncert);
66
67			cls limit_sys("cls_limit.ps",
68			sig, // signal with total syst. and stat. errors
69			bgd, // with total syst. and stat. errors
70	auterman	1.6	data,
71			bgd_NoSig);
72	auterman	1.1	limit_sys.SetStat(true); //Consider Bin-Errors (containing the full syst. uncertainties)
73			limit_sys.SetNMC(50000); //number of pseudo-experiments for each CL calculation
74			//limit_sys.Draw();
75
76	auterman	1.5	sprintf(n,"bkgdNLO%d",file); TH1 * bgdNLO = new TH1F(n, "",1,0,1);
77	auterman	1.6	sprintf(n,"bkgdNoSigNLO%d",file); TH1 * bgdNoSigNLO = new TH1F(n, "",1,0,1);
78	auterman	1.5	sprintf(n,"dataNLO%d",file); TH1 * dataNLO = new TH1F(n, "",1,0,1);
79			sprintf(n,"sigNLO%d", file); TH1 * sigNLO = new TH1F(n,"",1,0,1);
80	auterman	1.6
81			double bNLO=(nback - sigSignalNLO -sigSigTauNLO<0?0:nback - sigSignalNLO - sigSigTauNLO);
82
83	auterman	1.5	dataNLO->SetBinContent(1,ndata);
84			dataNLO->SetBinError(1,sqrt(ndata));
85			bgdNLO->SetBinContent(1,bNLO);
86			bgdNLO->SetBinError(1,bkgUncert);
87	auterman	1.6	bgdNoSigNLO->SetBinContent(1,nback);
88			bgdNoSigNLO->SetBinError(1,bkgUncert);
89	auterman	1.5	sigNLO->SetBinContent(1,signalNLO);
90			sigNLO->SetBinError(1,sigUncertNLO);
91
92			cls limit_sys_NLO("cls_limit_NLO.ps",
93			sigNLO, // signal with total syst. and stat. errors
94			bgdNLO, // with total syst. and stat. errors
95	auterman	1.6	dataNLO,
96			bgdNoSigNLO);
97	auterman	1.5	limit_sys_NLO.SetStat(true); //Consider Bin-Errors (containing the full syst. uncertainties)
98			limit_sys_NLO.SetNMC(50000); //number of pseudo-experiments for each CL calculation
99			//limit_sys.Draw();
100
101	auterman	1.1	try {
102	auterman	1.5	//limit_sys.WriteResult( &config, xsec );
103			limit_sys.WriteConfidence( &config, "LO." );
104			limit_sys_NLO.WriteConfidence( &config, "NLO." );
105	auterman	1.1	}
106			catch(exception& e){
107			cout << "Exception catched: " << e.what();
108			}
109
110			//write stuff:
111			time_t rawtime;
112			struct tm * timeinfo;
113			time ( &rawtime );
114			timeinfo = localtime ( &rawtime );
115			ofstream ofile;
116			ofile.open (argv[file]);
117	auterman	1.3	if (ofile.good())
118			std::cout << "writing '"<<argv[file]<<"'"<<std::endl;
119			else if ( (ofile.rdstate() & ifstream::failbit ) != 0 )
120			std::cerr << "ERROR opening '"<<argv[file]<<"'! Does the directory exist?"<<std::endl;
121
122	auterman	1.1	ofile << config.getComment() << asctime (timeinfo)
123			<< config.getComment()<< "\n"
124			<< config;
125	auterman	1.6
126			std::cout<<"d:"<<ndata<<", b:"<<bNLO<<", s:"<<signalNLO<<std::endl;
127
128	auterman	1.1	}
129			}
130
131
132			int main(int argc, char *argv[])
133			{
134			Limit(argc,argv);
135			}