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.25 by buchmann, Mon Sep 26 15:48:20 2011 UTC vs.
Revision 1.31 by buchmann, Thu Nov 24 08:19:17 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 doexpected) {
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/NewLimitCapsule.exec " << repname.str() << " " << luminosity << " " << luminosity*lumiuncert << " " << mceff << " " << mcefferr << " " << Npred[ibin] << " " << Nprederr[ibin] << " " << Nobs[ibin] << " " << -1 << " " << PlottingSetup::basedirectory << "/" << plotfilename << " " << doexpected;
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;
#	Line 251 | Line 273 | USAGE OF ROOSTATS_CL95
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 doexpected) {
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 271 | Line 293 | void compute_upper_limits_from_counting_
293		float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
294		float observed,observederr,null,result;
295
274	–	//      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);
275	–	//      observed-=result;//this is the actual excess we see!
276	–	//      float expected=observed/luminosity;
296		string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
297		dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
298	<	vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected);
299	<
298	>	vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doexpected,flipped);
299	>
300	>	tripple LibraryUpperLimits;
301	>	LibraryUpperLimits.name=signalsamples.collection[isample].samplename;
302	>	LibraryUpperLimits.first=mceff*signalsamples.collection[isample].xs * PlottingSetup::luminosity;
303	>	LibraryUpperLimits.second=staterr*signalsamples.collection[isample].xs * PlottingSetup::luminosity;
304	>	int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
305	>	(allresults.predictions[resultindex]).SignalYield.push_back(LibraryUpperLimits);
306	>
307		if(doexpected) {
282	–	//      rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
308		rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
309		vrows.push_back(sigmas[0]);
310		vrows.push_back(sigmas[1]);
286	–	//      vrows.push_back(expected);
311		vrows.push_back(signalsamples.collection[isample].xs);
312		}
313		else {
290	–	//      rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")");
314		rows.push_back(any2string(sigmas[0]));
315		vrows.push_back(sigmas[0]);
316		vrows.push_back(signalsamples.collection[isample].xs);
294	–	//      vrows.push_back(expected);
317		}
318		}//end of bin loop
319		limits.push_back(rows);
320		vlimits.push_back(vrows);
321		}//end of sample loop
322	+
323	+	dout << endl << endl << "_______________________________________________________________________________________" << endl;
324	+	dout << "Going to store upper limit on event yield in result library: " << endl;
325	+	for(int ibin=0;ibin<jzbcuts.size();ibin++) {
326	+	int resultindex=PlottingSetup::allresults.Find(jzbcuts[ibin]);
327	+	vector<float> Normsigmas = compute_one_upper_limit(1.0,0.0, resultindex, mcjzb, "UPPERLIMIT", false, 0);
328	+	(allresults.predictions[resultindex]).UpperLimit=Normsigmas[0]*PlottingSetup::luminosity;
329	+	}
330	+	dout << endl << "_______________________________________________________________________________________" << endl;
331	+
332		dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
333		dout << endl << endl << "PAS table 3:   (notation: limit [95%CL])" << endl << endl;
334		dout << "\t";
#	Line 351 | Line 383 | void compute_upper_limits_from_counting_
383		}
384		dout << endl;
385		}
386	+	allresults.Print();
387	+
388	+
389		}
390
391
#	Line 403 | Line 438 | void limit_shapes_for_systematic_effect(
438		TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
439		TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
440
441	<	TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
442	<	TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
443	<	TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
444	<	TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
441	>	TH1F *SBOSSFP;
442	>	TH1F *SBOSOFP;
443	>	TH1F *SBOSSFN;
444	>	TH1F *SBOSOFN;
445
446		TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
447		TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
448		TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
449		TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
450
451	<	TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
452	<	TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
453	<	TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
454	<	TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
451	>	TH1F *LSBOSSFP;
452	>	TH1F *LSBOSOFP;
453	>	TH1F *LSBOSSFN;
454	>	TH1F *LSBOSOFN;
455	>
456	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
457	>	if(PlottingSetup::RestrictToMassPeak) {
458	>	SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
459	>	SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
460	>	SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
461	>	SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
462	>
463	>	LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
464	>	LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
465	>	LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
466	>	LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
467	>	}
468
469		string obsname="data_obs";
470		string predname="background";
#	Line 427 | Line 475 | void limit_shapes_for_systematic_effect(
475		signalname="signal_"+identifier;
476		}
477
478	<	TH1F *obs = (TH1F*)ZOSSFP->Clone();
478	>	TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
479		obs->SetName(obsname.c_str());
480		obs->Write();
481	<	TH1F *pred = (TH1F*)ZOSSFN->Clone();
482	<	pred->Add(ZOSOFP,1.0/3);
483	<	pred->Add(ZOSOFN,-1.0/3);
484	<	pred->Add(SBOSSFP,1.0/3);
485	<	pred->Add(SBOSSFN,-1.0/3);
486	<	pred->Add(SBOSOFP,1.0/3);
487	<	pred->Add(SBOSOFN,-1.0/3);
481	>	TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
482	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
483	>	if(PlottingSetup::RestrictToMassPeak) {
484	>	pred->Add(ZOSOFP,1.0/3);
485	>	pred->Add(ZOSOFN,-1.0/3);
486	>	pred->Add(SBOSSFP,1.0/3);
487	>	pred->Add(SBOSSFN,-1.0/3);
488	>	pred->Add(SBOSOFP,1.0/3);
489	>	pred->Add(SBOSOFN,-1.0/3);
490	>	} else {
491	>	pred->Add(ZOSOFP,1.0);
492	>	pred->Add(ZOSOFN,-1.0);
493	>	}
494	>
495		pred->SetName(predname.c_str());
496		pred->Write();
497
#	Line 447 | Line 502 | void limit_shapes_for_systematic_effect(
502		TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
503		Lobs->Add(LZOSSFP);
504		Lpred->Add(LZOSSFN);
505	<	Lpred->Add(LZOSOFP,1.0/3);
506	<	Lpred->Add(LZOSOFN,-1.0/3);
507	<	Lpred->Add(LSBOSSFP,1.0/3);
508	<	Lpred->Add(LSBOSSFN,-1.0/3);
509	<	Lpred->Add(LSBOSOFP,1.0/3);
510	<	Lpred->Add(LSBOSOFN,-1.0/3);
505	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
506	>	if(PlottingSetup::RestrictToMassPeak) {
507	>	Lpred->Add(LZOSOFP,1.0/3);
508	>	Lpred->Add(LZOSOFN,-1.0/3);
509	>	Lpred->Add(LSBOSSFP,1.0/3);
510	>	Lpred->Add(LSBOSSFN,-1.0/3);
511	>	Lpred->Add(LSBOSOFP,1.0/3);
512	>	Lpred->Add(LSBOSOFN,-1.0/3);
513	>	} else {
514	>	Lpred->Add(LZOSOFP,1.0);
515	>	Lpred->Add(LZOSOFN,-1.0);
516	>	}
517	>
518		TH1F *signal = (TH1F*)Lobs->Clone();
519		signal->Add(Lpred,-1);
520		signal->SetName(signalname.c_str());
#	Line 466 | Line 528 | void limit_shapes_for_systematic_effect(
528		delete ZOSSFN;
529		delete ZOSOFN;
530
531	<	delete SBOSSFP;
532	<	delete SBOSOFP;
533	<	delete SBOSSFN;
534	<	delete SBOSOFN;
531	>	if(PlottingSetup::RestrictToMassPeak) {
532	>	delete SBOSSFP;
533	>	delete SBOSOFP;
534	>	delete SBOSSFN;
535	>	delete SBOSOFN;
536	>	}
537
538		delete LZOSSFP;
539		delete LZOSOFP;
540		delete LZOSSFN;
541		delete LZOSOFN;
542
543	<	delete LSBOSSFP;
544	<	delete LSBOSOFP;
545	<	delete LSBOSSFN;
546	<	delete LSBOSOFN;
543	>	if(PlottingSetup::RestrictToMassPeak) {
544	>	delete LSBOSSFP;
545	>	delete LSBOSOFP;
546	>	delete LSBOSSFN;
547	>	delete LSBOSOFN;
548	>	}
549
550		}
551

Diff Legend

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