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.21 by buchmann, Wed Sep 14 12:02:01 2011 UTC vs.
Revision 1.30 by buchmann, Wed Nov 16 14:50:33 2011 UTC

#	Line 1 | Line 1
1	+	/****
2	+
3	+	Off peak status (RestrictToMassPeak) :
4	+
5	+	x  Necessary adaptations identified
6	+	x  Started working on necessary adaptations
7	+	x  Necessary adaptations implemented
8	+	x  Necessary adaptations tested
9	+
10	+	DONE!
11	+
12	+
13	+	****/
14		#include <iostream>
15		#include <vector>
16		#include <sys/stat.h>
#	Line 170 | Line 183 | ratio_binning.push_back(80);
183
184		}
185
186	<	vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doobserved) {
186	>	vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doexpected, int flipped) {
187		float sigma95=-9.9,sigma95A=-9.9;
188		/*
189		USAGE OF ROOSTATS_CL95
#	Line 213 | Line 226 | USAGE OF ROOSTATS_CL95
226		- JZB cut [9]
227		- plot name  [10]*/
228
229	<	dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
229	>	if(flipped==0) dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
230	>	if(flipped>0) dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << flippedNpred[ibin] << "," << flippedNprederr[ibin] << "," << flippedNobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
231
232		stringstream command;
233	<	command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/LimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << plotfilename;
233	>	if(flipped==0) command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected;
234	>	if(flipped>0) command << PlottingSetup::cbafbasedir << "/DistributedModelCalculations/Limits/TimedLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << flippedNpred[ibin] << " " << flippedNprederr[ibin] << " " << flippedNobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected;
235		dout << command.str() << endl;
236
237		int retval = 256;
238		int attempts=0;
239	<	while(!(retval==0||attempts>=5)) {//try up to 5 times
239	>	while(!(retval==0||attempts>=3)) {//try up to 3 times
240		attempts++;
241	<	dout << "Starting limit calculation (LimitCapsule) now : Attempt " << attempts << endl;
241	>	dout << "Starting limit calculation (TimedLimitCapsule) now : Attempt " << attempts << endl;
242		retval=gSystem->Exec(command.str().c_str());
243		}
244	<
244	>	char hostname[1023];
245	>	gethostname(hostname,1023);
246	>	if((!((Contains(hostname,"t3ui")||Contains(hostname,"t3wn"))))&&retval==256) {
247	>	//running via CRAB and encountered the same problem too often: place a problem file to mark this problem!
248	>	stringstream markproblem;
249	>	markproblem << "touch " << PlottingSetup::basedirectory << "/exchange/problemswhilesettinglimits.txt";
250	>	gSystem->Exec(markproblem.str().c_str());
251	>	}
252		LimitDroplet limres;
253		limres.readDroplet(repname.str());
254		dout << limres << endl;
255		remove(repname.str().c_str());
256	<	sigma95=limres.expected;
256	>	sigma95=limres.observed;
257
258
259		///------------------------------------------ < /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
251	–	//up to here for backward compatibility
252	–	sigmas.push_back(limit.GetOneSigmaHighRange());//expected, up
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	–	*/
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	–	//  vector<float> sigmas;
264	–	//  sigmas.push_back(limit);
260		vector<float> sigmas;
261		sigmas.push_back(sigma95);
262	<	sigmas.push_back(sigma95A);
262	>	if(doexpected) {
263	>	sigmas.push_back(limres.expected);
264	>	sigmas.push_back(limres.upper68);
265	>	sigmas.push_back(limres.lower68);
266	>	sigmas.push_back(limres.upper95);
267	>	sigmas.push_back(limres.lower95);
268	>	}
269	>
270		return sigmas;
271
272
273	<	}
272	<	write_warning(__FUNCTION__,"STILL MISSING SIGMAS, LIMITS, EVERYTHING ...");
273	>	}//end of mc efficiency is ok
274		}
275
276	<	void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) {
276	>	void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected, int flipped) {
277		dout << "Doing counting experiment ... " << endl;
278		vector<vector<string> > limits;
279		vector<vector<float> > vlimits;
#	Line 297 | Line 298 | void compute_upper_limits_from_counting_
298		//      float expected=observed/luminosity;
299		string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
300		dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
301	<	vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doobserved);
301	>	vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected,flipped);
302
303	<	if(doobserved) {
303	>	if(doexpected) {
304		//      rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
305		rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
306		vrows.push_back(sigmas[0]);
#	Line 333 | Line 334 | void compute_upper_limits_from_counting_
334		dout << endl;
335		}
336
337	<	if(!doobserved) {
337	>	if(!doexpected) {
338		dout << endl << endl << "LIMITS: (Tex)" << endl;
339		tout << "\\begin{table}[hbtp]" << endl;
340		tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
#	Line 372 | Line 373 | void compute_upper_limits_from_counting_
373		}
374		dout << endl;
375		}
375	–
376	–	write_warning(__FUNCTION__,"Still need to update the script");
376		}
377
378
#	Line 426 | Line 425 | void limit_shapes_for_systematic_effect(
425		TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
426		TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
427
428	<	TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
429	<	TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
430	<	TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
431	<	TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
428	>	TH1F *SBOSSFP;
429	>	TH1F *SBOSOFP;
430	>	TH1F *SBOSSFN;
431	>	TH1F *SBOSOFN;
432
433		TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
434		TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
435		TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
436		TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
437
438	<	TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
439	<	TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
440	<	TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
441	<	TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
438	>	TH1F *LSBOSSFP;
439	>	TH1F *LSBOSOFP;
440	>	TH1F *LSBOSSFN;
441	>	TH1F *LSBOSOFN;
442	>
443	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
444	>	if(PlottingSetup::RestrictToMassPeak) {
445	>	SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
446	>	SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
447	>	SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
448	>	SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
449	>
450	>	LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
451	>	LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
452	>	LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
453	>	LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
454	>	}
455
456		string obsname="data_obs";
457		string predname="background";
#	Line 450 | Line 462 | void limit_shapes_for_systematic_effect(
462		signalname="signal_"+identifier;
463		}
464
465	<	TH1F *obs = (TH1F*)ZOSSFP->Clone();
465	>	TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
466		obs->SetName(obsname.c_str());
467		obs->Write();
468	<	TH1F *pred = (TH1F*)ZOSSFN->Clone();
469	<	pred->Add(ZOSOFP,1.0/3);
470	<	pred->Add(ZOSOFN,-1.0/3);
471	<	pred->Add(SBOSSFP,1.0/3);
472	<	pred->Add(SBOSSFN,-1.0/3);
473	<	pred->Add(SBOSOFP,1.0/3);
474	<	pred->Add(SBOSOFN,-1.0/3);
468	>	TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
469	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
470	>	if(PlottingSetup::RestrictToMassPeak) {
471	>	pred->Add(ZOSOFP,1.0/3);
472	>	pred->Add(ZOSOFN,-1.0/3);
473	>	pred->Add(SBOSSFP,1.0/3);
474	>	pred->Add(SBOSSFN,-1.0/3);
475	>	pred->Add(SBOSOFP,1.0/3);
476	>	pred->Add(SBOSOFN,-1.0/3);
477	>	} else {
478	>	pred->Add(ZOSOFP,1.0);
479	>	pred->Add(ZOSOFN,-1.0);
480	>	}
481	>
482		pred->SetName(predname.c_str());
483		pred->Write();
484
#	Line 470 | Line 489 | void limit_shapes_for_systematic_effect(
489		TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
490		Lobs->Add(LZOSSFP);
491		Lpred->Add(LZOSSFN);
492	<	Lpred->Add(LZOSOFP,1.0/3);
493	<	Lpred->Add(LZOSOFN,-1.0/3);
494	<	Lpred->Add(LSBOSSFP,1.0/3);
495	<	Lpred->Add(LSBOSSFN,-1.0/3);
496	<	Lpred->Add(LSBOSOFP,1.0/3);
497	<	Lpred->Add(LSBOSOFN,-1.0/3);
492	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
493	>	if(PlottingSetup::RestrictToMassPeak) {
494	>	Lpred->Add(LZOSOFP,1.0/3);
495	>	Lpred->Add(LZOSOFN,-1.0/3);
496	>	Lpred->Add(LSBOSSFP,1.0/3);
497	>	Lpred->Add(LSBOSSFN,-1.0/3);
498	>	Lpred->Add(LSBOSOFP,1.0/3);
499	>	Lpred->Add(LSBOSOFN,-1.0/3);
500	>	} else {
501	>	Lpred->Add(LZOSOFP,1.0);
502	>	Lpred->Add(LZOSOFN,-1.0);
503	>	}
504	>
505		TH1F *signal = (TH1F*)Lobs->Clone();
506		signal->Add(Lpred,-1);
507		signal->SetName(signalname.c_str());
#	Line 489 | Line 515 | void limit_shapes_for_systematic_effect(
515		delete ZOSSFN;
516		delete ZOSOFN;
517
518	<	delete SBOSSFP;
519	<	delete SBOSOFP;
520	<	delete SBOSSFN;
521	<	delete SBOSOFN;
518	>	if(PlottingSetup::RestrictToMassPeak) {
519	>	delete SBOSSFP;
520	>	delete SBOSOFP;
521	>	delete SBOSSFN;
522	>	delete SBOSOFN;
523	>	}
524
525		delete LZOSSFP;
526		delete LZOSOFP;
527		delete LZOSSFN;
528		delete LZOSOFN;
529
530	<	delete LSBOSSFP;
531	<	delete LSBOSOFP;
532	<	delete LSBOSSFN;
533	<	delete LSBOSOFN;
530	>	if(PlottingSetup::RestrictToMassPeak) {
531	>	delete LSBOSSFP;
532	>	delete LSBOSOFP;
533	>	delete LSBOSSFN;
534	>	delete LSBOSOFN;
535	>	}
536
537		}
538

Diff Legend

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