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.1 by buchmann, Tue Jul 19 08:15:45 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>
17	+	#include <fstream>
18
19		#include <TCut.h>
20		#include <TROOT.h>
#	Line 14 | Line 28
28		#include <TF1.h>
29		#include <TSQLResult.h>
30		#include <TProfile.h>
31	+	#include <TSystem.h>
32	+	#include "LimitDroplet.C"
33
34		//#include "TTbar_stuff.C"
35		using namespace std;
#	Line 22 | Line 38 | using namespace PlottingSetup;
38
39
40		void rediscover_the_top(string mcjzb, string datajzb) {
41	<	cout << "Hi! today we are going to (try to) rediscover the top!" << endl;
41	>	dout << "Hi! today we are going to (try to) rediscover the top!" << endl;
42		TCanvas *c3 = new TCanvas("c3","c3");
43		c3->SetLogy(1);
44		vector<float> binning;
#	Line 62 | Line 78 | void rediscover_the_top(string mcjzb, st
78		TH1F *datapredictiont = allsamples.Draw("datapredictiont",    "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
79		TH1F *datapredictiono = allsamples.Draw("datapredictiono",    "-"+datajzb, nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
80		datapredictiont->Add(datapredictiono,-1);
81	<	cout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl;
81	>	dout << "Second way of doing this !!!! Analytical shape to the left :-D" << endl;
82		vector<TF1*> functions = do_cb_fit_to_plot(datapredictiont,10);
83		datapredictiont->SetMarkerColor(kRed);
84		datapredictiont->SetLineColor(kRed);
#	Line 141 | Line 157 | ratio_binning.push_back(80);
157		ratio->Divide(datapredictiontb);
158
159		for (int i=0;i<=ratio_binning.size();i++) {
160	<	cout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
160	>	dout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
161		}
162
163		//  ratio->Divide(datapredictiontb);
#	Line 167 | Line 183 | ratio_binning.push_back(80);
183
184		}
185
186	<	void calculate_upper_limits(string mcjzb, string datajzb) {
187	<	write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
188	<	//  write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
189	<	//  rediscover_the_top(mcjzb,datajzb);
190	<	/*
191	<	TCanvas *c3 = new TCanvas("c3","c3");
192	<	c3->SetLogy(1);
193	<	vector<float> binning;
194	<	//binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
195	<	binning.push_back(50);
196	<	binning.push_back(100);
197	<	binning.push_back(150);
198	<	binning.push_back(200);
199	<	binning.push_back(500);
200	<	TH1F *datapredictiona = allsamples.Draw("datapredictiona",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity);
201	<	TH1F *datapredictionb = allsamples.Draw("datapredictionb",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
202	<	TH1F *datapredictionc = allsamples.Draw("datapredictionc",    datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
203	<	TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
204	<	dataprediction->Add(datapredictionb,-1);
205	<	dataprediction->Add(datapredictionc);
206	<	TH1F *puresignal     = allsamples.Draw("puresignal",        mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
207	<	TH1F *signalpred     = allsamples.Draw("signalpred",    "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
208	<	TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
209	<	TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
210	<	TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
211	<	signalpred->Add(signalpredlo,-1);
212	<	signalpred->Add(signalpredro);
213	<	puresignal->Add(signalpred,-1);//subtracting signal contamination
214	<	ofstream myfile;
215	<	myfile.open ("ShapeFit_log.txt");
216	<	establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
217	<	myfile.close();
218	<	*/
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
190	>	" Double_t             limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
191	>	" LimitResult expected_limit = roostats_clm(ilum, slum, eff, seff, bck, sbck, ntoys, nuisanceModel, method, seed); \n"
192	>	" Double_t     average_limit = roostats_cla(ilum, slum, eff, seff, bck, sbck, nuisanceModel, method, seed); \n"
193	>	"                                                                     \n"
194	>	"
195	>	" Double_t obs_limit = limit.GetObservedLimit();                      \n"
196	>	" Double_t exp_limit = limit.GetExpectedLimit();                      \n"
197	>	" Double_t exp_up    = limit.GetOneSigmaHighRange();                  \n"
198	>	" Double_t exp_down  = limit.GetOneSigmaLowRange();                   \n"
199	>	" Double_t exp_2up   = limit.GetTwoSigmaHighRange();                  \n"
200	>	" Double_t exp_2down = limit.GetTwoSigmaLowRange();                   \n"
201	>	*/
202	>	if(mceff<=0) {
203	>	write_warning(__FUNCTION__,"Cannot compute upper limit in this configuration as the efficiency is negative:");
204	>	dout << "mc efficiency=" << mceff << " +/- " << mcefferr;
205	>	vector<float> sigmas;
206	>	sigmas.push_back(-1);
207	>	sigmas.push_back(-1);
208	>	return sigmas;
209	>	} else {
210	>	int nlimittoysused=1;
211	>
212	>	///------------------------------------------ < NEW > ----------------------------------------------------------
213	>
214	>	int secondssince1970=time(NULL);
215	>	stringstream repname;
216	>	repname << PlottingSetup::cbafbasedir << "/exchange/report_" << secondssince1970 << "_"<<plotfilename<< "__"<< ".txt";
217	>
218	>	/* - report filename [1]
219	>	- luminosity [2]
220	>	- lumi uncert [3]
221	>	- MC efficiency [4]
222	>	- MC efficiency error [5]
223	>	- Npred [6]
224	>	- Nprederr [7]
225	>	- Nobs [8]
226	>	- JZB cut [9]
227	>	- plot name  [10]*/
228	>
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	>	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>=3)) {//try up to 3 times
240	>	attempts++;
241	>	dout << "Starting limit calculation (TimedLimitCapsule) now : Attempt " << attempts << endl;
242	>	retval=gSystem->Exec(command.str().c_str());
243	>	}
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.observed;
257	>
258	>
259	>	///------------------------------------------ < /NEW > ----------------------------------------------------------
260	>	vector<float> sigmas;
261	>	sigmas.push_back(sigma95);
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	>	}//end of mc efficiency is ok
274		}
275
276	<	void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts) {
277	<	cout << "Doing counting experiment ... " << endl;
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;
280	>
281
282		for(int isample=0;isample<signalsamples.collection.size();isample++) {
283	<	cout << "Considering sample " << signalsamples.collection[isample].samplename << endl;
283	>	vector<string> rows;
284	>	vector<float> vrows;
285	>	dout << "Considering sample " << signalsamples.collection[isample].samplename << endl;
286	>	rows.push_back(signalsamples.collection[isample].samplename);
287		for(int ibin=0;ibin<jzbcuts.size();ibin++) {
288	<	cout << "   Considering bin " << jzbcuts[ibin] << endl;
289	<	for(int isyst=0;isyst<uncertainties[isample*jzbcuts.size()+ibin].size();isyst++) {
290	<	if(isyst==0) {
291	<	if(uncertainties[isample*jzbcuts.size()+ibin][isyst]!=jzbcuts[ibin]) cout << "WATCH OUT THERE SEEMS TO BE A PROBLEM - THE TEST BIN DOESN'T CORRESPOND TO YOUR BIN!" << endl;
292	<	continue;
293	<	}
294	<	cout << "     Considering syst " << isyst << " : " << uncertainties[isample*jzbcuts.size()+ibin][isyst] << endl;
295	<	}//end of syst loop
288	>	dout << "_________________________________________________________________________________" << endl;
289	>	float JZBcutat=uncertainties[isample*jzbcuts.size()+ibin][0];
290	>	float mceff=uncertainties[isample*jzbcuts.size()+ibin][1];
291	>	float staterr=uncertainties[isample*jzbcuts.size()+ibin][2];
292	>	float systerr=uncertainties[isample*jzbcuts.size()+ibin][3];
293	>	float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
294	>	float observed,observederr,null,result;
295	>
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,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) {
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]);
311	>	vrows.push_back(signalsamples.collection[isample].xs);
312	>	}
313	>	else {
314	>	rows.push_back(any2string(sigmas[0]));
315	>	vrows.push_back(sigmas[0]);
316	>	vrows.push_back(signalsamples.collection[isample].xs);
317	>	}
318		}//end of bin loop
319	+	limits.push_back(rows);
320	+	vlimits.push_back(vrows);
321		}//end of sample loop
221	–	}
322
323	<	void susy_scan_axis_labeling(TH2F *histo) {
324	<	histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
325	<	histo->GetXaxis()->CenterTitle();
326	<	histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
327	<	histo->GetYaxis()->CenterTitle();
328	<	}
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	<	void scan_susy_space(string mcjzb, string datajzb) {
333	<	TCanvas *c3 = new TCanvas("c3","c3");
334	<	vector<float> binning;
335	<	binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
336	<	/*
337	<	binning.push_back(50);
338	<	binning.push_back(75);
339	<	binning.push_back(100);
340	<	binning.push_back(150);
341	<	binning.push_back(200);
342	<	binning.push_back(500);
241	<	*/
242	<	float arrbinning[binning.size()];
243	<	for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
244	<	TH1F *puredata   = allsamples.Draw("puredata",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
245	<	puredata->SetMarkerSize(DataMarkerSize);
246	<	TH1F *allbgs   = allsamples.Draw("allbgs",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
247	<	TH1F *allbgsb   = allsamples.Draw("allbgsb",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
248	<	TH1F *allbgsc   = allsamples.Draw("allbgsc",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
249	<	allbgs->Add(allbgsb,-1);
250	<	allbgs->Add(allbgsc);
251	<	int ndata=puredata->Integral();
252	<	ofstream myfile;
253	<	myfile.open ("susyscan_log.txt");
254	<	TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
255	<	TTree *suevents = (TTree*)susyscanfile->Get("events");
256	<	TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
257	<	TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
258	<	TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
259	<	TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
260	<
261	<	susy_scan_axis_labeling(exclusionmap);
262	<	susy_scan_axis_labeling(exclusionmap1s);
263	<	susy_scan_axis_labeling(exclusionmap2s);
264	<	susy_scan_axis_labeling(exclusionmap3s);
265	<
266	<	Int_t MyPalette[100];
267	<	Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
268	<	Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
269	<	Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
270	<	Double_t stop[] = {0., .25, .50, .75, 1.0};
271	<	Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
272	<	for (int i=0;i<100;i++) MyPalette[i] = FI+i;
273	<
274	<	gStyle->SetPalette(100, MyPalette);
275	<
276	<	for(int m23=50;m23<500;m23+=25) {
277	<	for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
278	<	{
279	<	c3->cd();
280	<	stringstream drawcondition;
281	<	drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
282	<	TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
283	<	TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
284	<	stringstream drawvar,drawvar2;
285	<	drawvar<<mcjzb<<">>puresignal";
286	<	drawvar2<<"-"<<mcjzb<<">>puresignall";
287	<	suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
288	<	suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
289	<	if(puresignal->Integral()<60) {
290	<	delete puresignal;
291	<	continue;
292	<	}
293	<	puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
294	<	puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
295	<	stringstream saveas;
296	<	saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
297	<	cout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
298	<	//        TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
299	<	//        TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
300	<	//        TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
301	<	//        signalpred->Add(signalpredlo,-1);
302	<	//        signalpred->Add(signalpredro);
303	<	//        puresignal->Add(signalpred,-1);//subtracting signal contamination
304	<	//---------------------
305	<	//      cout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
306	<	//    TH1F *puresignal = allsamples.Draw("puresignal",mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
307	<	vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
308	<	if(results.size()==0) {
309	<	delete puresignal;
310	<	continue;
311	<	}
312	<	exclusionmap->Fill(m23,m0,results[0]);
313	<	exclusionmap1s->Fill(m23,m0,results[1]);
314	<	exclusionmap2s->Fill(m23,m0,results[2]);
315	<	exclusionmap3s->Fill(m23,m0,results[3]);
316	<	delete puresignal;
317	<	cout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
332	>	dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
333	>	dout << endl << endl << "PAS table 3:   (notation: limit [95%CL])" << endl << endl;
334	>	dout << "\t";
335	>	for (int irow=0;irow<jzbcuts.size();irow++) {
336	>	dout << jzbcuts[irow] << "\t";
337	>	}
338	>	dout << endl;
339	>	for(int irow=0;irow<limits.size();irow++) {
340	>	for(int ientry=0;ientry<limits[irow].size();ientry++) {
341	>	if (limits[irow][ientry]>0) dout << limits[irow][ientry] << "\t";
342	>	else dout << " (N/A) \t";
343		}
344	<	}//end of model scan for loop
344	>	dout << endl;
345	>	}
346
347	<	cout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
348	<	c3->cd();
349	<	exclusionmap->Draw("CONTZ");
350	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
351	<	exclusionmap->Draw("COLZ");
352	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
353	<
354	<	exclusionmap1s->Draw("CONTZ");
355	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
356	<	exclusionmap1s->Draw("COLZ");
357	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
358	<
359	<	exclusionmap2s->Draw("CONTZ");
360	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
361	<	exclusionmap2s->Draw("COLZ");
362	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
363	<
364	<	exclusionmap3s->Draw("CONTZ");
365	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
366	<	exclusionmap3s->Draw("COLZ");
367	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
368	<
369	<	TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
370	<	exclusionmap->Write();
371	<	exclusionmap1s->Write();
372	<	exclusionmap2s->Write();
373	<	exclusionmap3s->Write();
374	<	exclusion_limits->Close();
375	<	susyscanfile->Close();
347	>	if(!doexpected) {
348	>	dout << endl << endl << "LIMITS: (Tex)" << endl;
349	>	tout << "\\begin{table}[hbtp]" << endl;
350	>	tout << "\\renewcommand{\\arraystretch}{1.3}" << endl;
351	>	tout << "\\begin{center}" << endl;
352	>	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;
353	>	tout << "" << endl;
354	>	tout << "\\begin{tabular}{ | l | ";
355	>	for (int irow=0;irow<jzbcuts.size();irow++) tout << " l |";
356	>	tout << "} " << endl << " \\hline " << endl << "& \t ";
357	>	for (int irow=0;irow<jzbcuts.size();irow++) {
358	>	tout << "JZB $>$ " << jzbcuts[irow] << " GeV & \t ";
359	>	}
360	>	tout << " \\\\ \\hline " << endl;
361	>	for(int irow=0;irow<limits.size();irow++) {
362	>	tout << limits[irow][0] << " \t";
363	>	for(int ientry=0;ientry<jzbcuts.size();ientry++) {
364	>	if(vlimits[irow][2*ientry]>0) tout << " & " << Round(vlimits[irow][2*ientry],2) << " \t (" << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "x \\sigma ) \t";
365	>	else tout << " & ( N / A ) \t";
366	>	//      dout << Round(vlimits[irow][2*ientry],3) << " / " << Round(vlimits[irow][2*ientry+1],3)<< "\t";
367	>	}
368	>	tout << " \\\\ \\hline " << endl;
369	>	}
370	>	tout << "\\end{tabular}" << endl;
371	>	tout << "      \\end{tabular}"<< endl;
372	>	tout << "\\end{center}"<< endl;
373	>	tout << "\\end{table} "<< endl;
374	>
375	>	}//do observed
376	>
377	>	dout << endl << endl << "Final selection efficiencies with total statistical and systematic errors, and corresponding observed and expected upper limits (UL) on ($\\sigma\\times$  BR $\\times$ acceptance) for the LM4 and LM8 scenarios, in the different regions. The last column contains the predicted ($\\sigma \\times $BR$\\times$ acceptance) at NLO obtained from Monte Carlo simulation." << endl;
378	>	dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t \\sigma [pb]" << endl;
379	>	for(int icut=0;icut<jzbcuts.size();icut++) {
380	>	dout << "Region with JZB>" << jzbcuts[icut] << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << endl;
381	>	for(int isample=0;isample<signalsamples.collection.size();isample++) {
382	>	dout << limits[isample][0] << "\t" << Round(100*uncertainties[isample*jzbcuts.size()+icut][1],3) << "+/-" << Round(100*uncertainties[isample*jzbcuts.size()+icut][2],3) << " (stat) +/- " << Round(100*uncertainties[isample*jzbcuts.size()+icut][3],3) << " (syst) \t" << Round((vlimits[isample][2*icut]),3) << "\t" << Round(vlimits[isample][2*icut+1],3) << endl;
383	>	}
384	>	dout << endl;
385	>	}
386	>	allresults.Print();
387	>
388
351	–	myfile.close();
389		}
390
391
392	+
393	+	/********************************************************************** new : Limits using SHAPES ***********************************
394	+
395	+
396	+	SSSSSSSSSSSSSSS hhhhhhh
397	+	SS:::::::::::::::Sh:::::h
398	+	S:::::SSSSSS::::::Sh:::::h
399	+	S:::::S     SSSSSSSh:::::h
400	+	S:::::S             h::::h hhhhh         aaaaaaaaaaaaa  ppppp   ppppppppp       eeeeeeeeeeee        ssssssssss
401	+	S:::::S             h::::hh:::::hhh      a::::::::::::a p::::ppp:::::::::p    ee::::::::::::ee    ss::::::::::s
402	+	S::::SSSS          h::::::::::::::hh    aaaaaaaaa:::::ap:::::::::::::::::p  e::::::eeeee:::::eess:::::::::::::s
403	+	SS::::::SSSSS     h:::::::hhh::::::h            a::::app::::::ppppp::::::pe::::::e     e:::::es::::::ssss:::::s
404	+	SSS::::::::SS   h::::::h   h::::::h    aaaaaaa:::::a p:::::p     p:::::pe:::::::eeeee::::::e s:::::s  ssssss
405	+	SSSSSS::::S  h:::::h     h:::::h  aa::::::::::::a p:::::p     p:::::pe:::::::::::::::::e    s::::::s
406	+	S:::::S h:::::h     h:::::h a::::aaaa::::::a p:::::p     p:::::pe::::::eeeeeeeeeee        s::::::s
407	+	S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::p    p::::::pe:::::::e           ssssss   s:::::s
408	+	SSSSSSS     S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::ppppp:::::::pe::::::::e          s:::::ssss::::::s
409	+	S::::::SSSSSS:::::S h:::::h     h:::::ha:::::aaaa::::::a p::::::::::::::::p  e::::::::eeeeeeee  s::::::::::::::s
410	+	S:::::::::::::::SS  h:::::h     h:::::h a::::::::::aa:::ap::::::::::::::pp    ee:::::::::::::e   s:::::::::::ss
411	+	SSSSSSSSSSSSSSS    hhhhhhh     hhhhhhh  aaaaaaaaaa  aaaap::::::pppppppp        eeeeeeeeeeeeee    sssssssssss
412	+	p:::::p
413	+	p:::::p
414	+	p:::::::p
415	+	p:::::::p
416	+	p:::::::p
417	+	ppppppppp
418	+
419	+
420	+	*********************************************************************** new : Limits using SHAPES ***********************************/
421	+
422	+
423	+	void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
424	+	dout << "Creatig shape templates ... ";
425	+	if(identifier!="") dout << "for systematic called "<<identifier;
426	+	dout << endl;
427	+	int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
428	+	if(dataormc!=data) write_warning(__FUNCTION__,"WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
429	+
430	+	TCut limitnJetcut;
431	+	if(JES==noJES) limitnJetcut=cutnJets;
432	+	else {
433	+	if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
434	+	if(JES==JESup) limitnJetcut=cutnJetsJESup;
435	+	}
436	+	TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
437	+	TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
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;
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;
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";
471	+	string signalname="signal";
472	+	if(identifier!="") {
473	+	obsname=("data_"+identifier);
474	+	predname=("background_"+identifier);
475	+	signalname="signal_"+identifier;
476	+	}
477	+
478	+	TH1F *obs = (TH1F*)ZOSSFP->Clone("observation");
479	+	obs->SetName(obsname.c_str());
480	+	obs->Write();
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	+
498	+	//  TH1F *Lobs = (TH1F*)LZOSSFP->Clone();
499	+	//  TH1F *Lpred = (TH1F*)LZOSSFN->Clone();
500	+
501	+	TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
502	+	TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
503	+	Lobs->Add(LZOSSFP);
504	+	Lpred->Add(LZOSSFN);
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());
521	+	signal->Write();
522	+
523	+	delete Lobs;
524	+	delete Lpred;
525	+
526	+	delete ZOSSFP;
527	+	delete ZOSOFP;
528	+	delete ZOSSFN;
529	+	delete ZOSOFN;
530	+
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	+	if(PlottingSetup::RestrictToMassPeak) {
544	+	delete LSBOSSFP;
545	+	delete LSBOSOFP;
546	+	delete LSBOSSFN;
547	+	delete LSBOSOFN;
548	+	}
549
550	+	}
551	+
552	+	void prepare_datacard(TFile *f) {
553	+	TH1F *dataob = (TH1F*)f->Get("data_obs");
554	+	TH1F *signal = (TH1F*)f->Get("signal");
555	+	TH1F *background = (TH1F*)f->Get("background");
556
557	+	ofstream datacard;
558	+	ensure_directory_exists(get_directory()+"/limits");
559	+	datacard.open ((get_directory()+"/limits/susydatacard.txt").c_str());
560	+	datacard << "Writing this to a file.\n";
561	+	datacard << "imax 1\n";
562	+	datacard << "jmax 1\n";
563	+	datacard << "kmax *\n";
564	+	datacard << "---------------\n";
565	+	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
566	+	datacard << "---------------\n";
567	+	datacard << "bin 1\n";
568	+	datacard << "observation "<<dataob->Integral()<<"\n";
569	+	datacard << "------------------------------\n";
570	+	datacard << "bin             1          1\n";
571	+	datacard << "process         signal     background\n";
572	+	datacard << "process         0          1\n";
573	+	datacard << "rate            "<<signal->Integral()<<"         "<<background->Integral()<<"\n";
574	+	datacard << "--------------------------------\n";
575	+	datacard << "lumi     lnN    1.10       1.0\n";
576	+	datacard << "bgnorm   lnN    1.00       1.4  uncertainty on our prediction (40%)\n";
577	+	datacard << "JES    shape    1          1    uncertainty on background shape and normalization\n";
578	+	datacard << "peak   shape    1          1    uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
579	+	datacard << "#                                so divide the unit gaussian by 2 before doing the interpolation\n";
580	+	datacard.close();
581	+	}
582	+
583	+
584	+	void prepare_limits(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbbins) {
585	+	ensure_directory_exists(get_directory()+"/limits");
586	+	TFile *limfile = new TFile((get_directory()+"/limits/limitfile.root").c_str(),"RECREATE");
587	+	TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
588	+	limit_shapes_for_systematic_effect(limfile,"",mcjzb,datajzb,noJES,jzbbins,limcan);
589	+	limit_shapes_for_systematic_effect(limfile,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan);
590	+	limit_shapes_for_systematic_effect(limfile,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan);
591	+	limit_shapes_for_systematic_effect(limfile,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan);
592	+	limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
593	+
594	+	prepare_datacard(limfile);
595	+	limfile->Close();
596	+	write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
597	+
598	+	}

Diff Legend

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