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.24 by buchmann, Thu Sep 15 13:54:15 2011 UTC vs.
Revision 1.27 by buchmann, Mon Oct 24 15:26:14 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 216 | Line 229 | USAGE OF ROOSTATS_CL95
229		dout << "Calling limit capsule instead of calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
230
231		stringstream command;
232	<	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;
232	>	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;
233		dout << command.str() << endl;
234
235		int retval = 256;
236		int attempts=0;
237	<	while(!(retval==0||attempts>=5)) {//try up to 5 times
237	>	while(!(retval==0||attempts>=3)) {//try up to 3 times
238		attempts++;
239		dout << "Starting limit calculation (LimitCapsule) now : Attempt " << attempts << endl;
240		retval=gSystem->Exec(command.str().c_str());
#	Line 235 | Line 248 | USAGE OF ROOSTATS_CL95
248
249
250		///------------------------------------------ < /NEW > ----------------------------------------------------------
238	–	//if(doexpected) 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;
251		vector<float> sigmas;
252	<	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);
252	>	sigmas.push_back(sigma95);
253		if(doexpected) {
254	<	write_warning(__FUNCTION__,"EXPECTED LIMITS HAVE BEEN SWITCHED OFF TEMPORARILY");
255	<	//    dout << "Now calling : CLA(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << nuisancemodel<< ") " << endl;
256	<	//    sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], nuisancemodel);
254	>	sigmas.push_back(limres.expected);
255	>	sigmas.push_back(limres.upper68);
256	>	sigmas.push_back(limres.lower68);
257	>	sigmas.push_back(limres.upper95);
258	>	sigmas.push_back(limres.lower95);
259		}
260	<	//  vector<float> sigmas;
264	<	//  sigmas.push_back(limit);
265	<	vector<float> sigmas;
266	<	sigmas.push_back(sigma95);
267	<	sigmas.push_back(sigma95A);
260	>
261		return sigmas;
262
263
264	<	}
272	<	write_warning(__FUNCTION__,"STILL MISSING SIGMAS, LIMITS, EVERYTHING ...");
264	>	}//end of mc efficiency is ok
265		}
266
267		void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doexpected) {
#	Line 424 | Line 416 | void limit_shapes_for_systematic_effect(
416		TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
417		TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
418
419	<	TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
420	<	TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
421	<	TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
422	<	TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
419	>	TH1F *SBOSSFP;
420	>	TH1F *SBOSOFP;
421	>	TH1F *SBOSSFN;
422	>	TH1F *SBOSOFN;
423
424		TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
425		TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
426		TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
427		TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
428
429	<	TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
430	<	TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
431	<	TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
432	<	TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
429	>	TH1F *LSBOSSFP;
430	>	TH1F *LSBOSOFP;
431	>	TH1F *LSBOSSFN;
432	>	TH1F *LSBOSOFN;
433	>
434	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
435	>	if(PlottingSetup::RestrictToMassPeak) {
436	>	SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
437	>	SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
438	>	SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
439	>	SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
440	>
441	>	LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
442	>	LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
443	>	LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
444	>	LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
445	>	}
446
447		string obsname="data_obs";
448		string predname="background";
#	Line 448 | Line 453 | void limit_shapes_for_systematic_effect(
453		signalname="signal_"+identifier;
454		}
455
456	<	TH1F *obs = (TH1F*)ZOSSFP->Clone();
456	>	TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
457		obs->SetName(obsname.c_str());
458		obs->Write();
459	<	TH1F *pred = (TH1F*)ZOSSFN->Clone();
460	<	pred->Add(ZOSOFP,1.0/3);
461	<	pred->Add(ZOSOFN,-1.0/3);
462	<	pred->Add(SBOSSFP,1.0/3);
463	<	pred->Add(SBOSSFN,-1.0/3);
464	<	pred->Add(SBOSOFP,1.0/3);
465	<	pred->Add(SBOSOFN,-1.0/3);
459	>	TH1F *pred = (TH1F*)ZOSSFN->Clone("prediction");
460	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
461	>	if(PlottingSetup::RestrictToMassPeak) {
462	>	pred->Add(ZOSOFP,1.0/3);
463	>	pred->Add(ZOSOFN,-1.0/3);
464	>	pred->Add(SBOSSFP,1.0/3);
465	>	pred->Add(SBOSSFN,-1.0/3);
466	>	pred->Add(SBOSOFP,1.0/3);
467	>	pred->Add(SBOSOFN,-1.0/3);
468	>	} else {
469	>	pred->Add(ZOSOFP,1.0);
470	>	pred->Add(ZOSOFN,-1.0);
471	>	}
472	>
473		pred->SetName(predname.c_str());
474		pred->Write();
475
#	Line 468 | Line 480 | void limit_shapes_for_systematic_effect(
480		TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
481		Lobs->Add(LZOSSFP);
482		Lpred->Add(LZOSSFN);
483	<	Lpred->Add(LZOSOFP,1.0/3);
484	<	Lpred->Add(LZOSOFN,-1.0/3);
485	<	Lpred->Add(LSBOSSFP,1.0/3);
486	<	Lpred->Add(LSBOSSFN,-1.0/3);
487	<	Lpred->Add(LSBOSOFP,1.0/3);
488	<	Lpred->Add(LSBOSOFN,-1.0/3);
483	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
484	>	if(PlottingSetup::RestrictToMassPeak) {
485	>	Lpred->Add(LZOSOFP,1.0/3);
486	>	Lpred->Add(LZOSOFN,-1.0/3);
487	>	Lpred->Add(LSBOSSFP,1.0/3);
488	>	Lpred->Add(LSBOSSFN,-1.0/3);
489	>	Lpred->Add(LSBOSOFP,1.0/3);
490	>	Lpred->Add(LSBOSOFN,-1.0/3);
491	>	} else {
492	>	Lpred->Add(LZOSOFP,1.0);
493	>	Lpred->Add(LZOSOFN,-1.0);
494	>	}
495	>
496		TH1F *signal = (TH1F*)Lobs->Clone();
497		signal->Add(Lpred,-1);
498		signal->SetName(signalname.c_str());
#	Line 487 | Line 506 | void limit_shapes_for_systematic_effect(
506		delete ZOSSFN;
507		delete ZOSOFN;
508
509	<	delete SBOSSFP;
510	<	delete SBOSOFP;
511	<	delete SBOSSFN;
512	<	delete SBOSOFN;
509	>	if(PlottingSetup::RestrictToMassPeak) {
510	>	delete SBOSSFP;
511	>	delete SBOSOFP;
512	>	delete SBOSSFN;
513	>	delete SBOSOFN;
514	>	}
515
516		delete LZOSSFP;
517		delete LZOSOFP;
518		delete LZOSSFN;
519		delete LZOSOFN;
520
521	<	delete LSBOSSFP;
522	<	delete LSBOSOFP;
523	<	delete LSBOSSFN;
524	<	delete LSBOSOFN;
521	>	if(PlottingSetup::RestrictToMassPeak) {
522	>	delete LSBOSSFP;
523	>	delete LSBOSOFP;
524	>	delete LSBOSSFN;
525	>	delete LSBOSOFN;
526	>	}
527
528		}
529

Diff Legend

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