ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/LimitCalculation.C

Comparing UserCode/cbrown/AnalysisFramework/Plotting/Modules/LimitCalculation.C (file contents):
Revision 1.15 by buchmann, Thu Sep 1 11:01:42 2011 UTC vs.
Revision 1.21 by buchmann, Wed Sep 14 12:02:01 2011 UTC

#	Line 15 | Line 15
15		#include <TF1.h>
16		#include <TSQLResult.h>
17		#include <TProfile.h>
18	+	#include <TSystem.h>
19	+	#include "LimitDroplet.C"
20
21		//#include "TTbar_stuff.C"
22		using namespace std;
#	Line 170 | Line 172 | ratio_binning.push_back(80);
172
173		vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doobserved) {
174		float sigma95=-9.9,sigma95A=-9.9;
173	–	int nuisancemodel=0;
175		/*
176		USAGE OF ROOSTATS_CL95
177		" Double_t             limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
#	Line 194 | Line 195 | USAGE OF ROOSTATS_CL95
195		return sigmas;
196		} else {
197		int nlimittoysused=1;
198	<	if(doobserved) nlimittoysused=nlimittoys;
199	<	/*  dout << "Now calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
200	<	sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel);
201	<	*/
202	<	dout << "Now calling : roostats_limit(" << luminosity << "," << lumiuncert*luminosity << ","<<mceff <<","<<mcefferr<<","<<Npred[ibin]<<","<<Nprederr[ibin] << ",n=" << nlimittoysused << ",gauss=" << false << ", nuisanceModel="<<nuisancemodel<<",method="<<limitmethod<<",plotfilename="<<plotfilename<<",seed=1) " << endl;
203	<	LimitResult limit = roostats_limit(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,1);
204	<
198	>
199	>	///------------------------------------------ < NEW > ----------------------------------------------------------
200	>
201	>	int secondssince1970=time(NULL);
202	>	stringstream repname;
203	>	repname << PlottingSetup::cbafbasedir << "/exchange/report_" << secondssince1970 << "_"<<plotfilename<< "__"<< ".txt";
204	>
205	>	/* - report filename [1]
206	>	- luminosity [2]
207	>	- lumi uncert [3]
208	>	- MC efficiency [4]
209	>	- MC efficiency error [5]
210	>	- Npred [6]
211	>	- Nprederr [7]
212	>	- Nobs [8]
213	>	- JZB cut [9]
214	>	- plot name  [10]*/
215	>
216	>	dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
217	>
218	>	stringstream command;
219	>	command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/LimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << plotfilename;
220	>	dout << command.str() << endl;
221	>
222	>	int retval = 256;
223	>	int attempts=0;
224	>	while(!(retval==0||attempts>=5)) {//try up to 5 times
225	>	attempts++;
226	>	dout << "Starting limit calculation (LimitCapsule) now : Attempt " << attempts << endl;
227	>	retval=gSystem->Exec(command.str().c_str());
228	>	}
229	>
230	>	LimitDroplet limres;
231	>	limres.readDroplet(repname.str());
232	>	dout << limres << endl;
233	>	remove(repname.str().c_str());
234	>	sigma95=limres.expected;
235	>
236	>
237	>	///------------------------------------------ < /NEW > ----------------------------------------------------------
238	>	//if(doobserved) nlimittoysused=nlimittoys;
239	>	//    nlimittoysused=nlimittoys;
240	>	//  dout << "Now calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
241	>	//  sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel);
242	>	//  sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel);
243	>
244	>	/*  dout << "Now calling : roostats_cl95(" << luminosity << "," << lumiuncert*luminosity << ","<<mceff <<","<<mcefferr<<","<<Npred[ibin]<<","<<Nprederr[ibin] << ",n=" << nlimittoysused << ",gauss=" << false << ",nuisanceModel="<<nuisancemodel<<",method="<<limitmethod<<",plotfilename="<<plotfilename<<",seed=0) " << endl;
245	>	/*  dout << "Now calling : roostats_limit(" << luminosity << "," << lumiuncert*luminosity << ","<<mceff <<","<<mcefferr<<","<<Npred[ibin]<<","<<Nprederr[ibin] << ",n=" << nlimittoysused << ",gauss=" << false << ", nuisanceModel="<<nuisancemodel<<",method="<<limitmethod<<",plotfilename="<<plotfilename<<",seed=1) " << endl;
246	>	LimitResult limit = roostats_limit(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,0);
247	>	dout << "Now interpreting and saving results ... " << endl;
248		vector<float> sigmas;
249		sigmas.push_back(limit.GetExpectedLimit());//expected
250		sigmas.push_back(limit.GetObservedLimit());//observed
#	Line 209 | Line 253 | USAGE OF ROOSTATS_CL95
253		sigmas.push_back(limit.GetTwoSigmaHighRange());//expected, 2 up
254		sigmas.push_back(limit.GetOneSigmaLowRange());//expected, down
255		sigmas.push_back(limit.GetTwoSigmaLowRange());//expected, 2 down
256	<	//  if(doobserved) {
256	>	*/
257	>	//  float limit = roostats_cl95(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,0);
258	>	if(doobserved) {
259	>	write_warning(__FUNCTION__,"OBSERVED LIMITS HAVE BEEN SWITCHED OFF TEMPORARILY");
260		//    dout << "Now calling : CLA(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << nuisancemodel<< ") " << endl;
261		//    sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], nuisancemodel);
262	<	//  }
263	<
264	<	/*  vector<float> sigmas;
262	>	}
263	>	//  vector<float> sigmas;
264	>	//  sigmas.push_back(limit);
265	>	vector<float> sigmas;
266		sigmas.push_back(sigma95);
267	<	sigmas.push_back(sigma95A);*/
267	>	sigmas.push_back(sigma95A);
268		return sigmas;
269
270
#	Line 247 | Line 295 | void compute_upper_limits_from_counting_
295		//      fill_result_histos(observed,observederr, null,null,null,null,null,null,null,mcjzb,JZBcutat,14000,(int)5,result,(signalsamples.FindSample(signalsamples.collection[isample].filename)),signalsamples);
296		//      observed-=result;//this is the actual excess we see!
297		//      float expected=observed/luminosity;
298	<	string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat)));
298	>	string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
299		dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
300		vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doobserved);
301
#	Line 288 | Line 336 | void compute_upper_limits_from_counting_
336		if(!doobserved) {
337		dout << endl << endl << "LIMITS: (Tex)" << endl;
338		tout << "\\begin{table}[hbtp]" << endl;
339	<	tout << "\\renewcommand{\arraystretch}{1.3}" << endl;
339	>	tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
340		tout << "\\begin{center}" << endl;
341		tout << "\\caption{Observed upper limits on the cross section of different LM benchmark points " << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << "}\\label{tab:lmresults}" << endl;
342		tout << "" << endl;

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines