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.3 by buchmann, Wed Jul 20 08:52:17 2011 UTC vs.
Revision 1.17 by buchmann, Thu Sep 1 16:23:46 2011 UTC

#	Line 1 | Line 1
1		#include <iostream>
2		#include <vector>
3		#include <sys/stat.h>
4	+	#include <fstream>
5
6		#include <TCut.h>
7		#include <TROOT.h>
#	Line 167 | Line 168 | ratio_binning.push_back(80);
168
169		}
170
171	<	void calculate_upper_limits(string mcjzb, string datajzb) {
172	<	write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
173	<	//  write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
174	<	//  rediscover_the_top(mcjzb,datajzb);
175	<	/*
176	<	TCanvas *c3 = new TCanvas("c3","c3");
177	<	c3->SetLogy(1);
178	<	vector<float> binning;
179	<	//binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
180	<	binning.push_back(50);
181	<	binning.push_back(100);
182	<	binning.push_back(150);
183	<	binning.push_back(200);
184	<	binning.push_back(500);
185	<	TH1F *datapredictiona = allsamples.Draw("datapredictiona",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity);
186	<	TH1F *datapredictionb = allsamples.Draw("datapredictionb",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
187	<	TH1F *datapredictionc = allsamples.Draw("datapredictionc",    datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
188	<	TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
189	<	dataprediction->Add(datapredictionb,-1);
190	<	dataprediction->Add(datapredictionc);
191	<	TH1F *puresignal     = allsamples.Draw("puresignal",        mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
192	<	TH1F *signalpred     = allsamples.Draw("signalpred",    "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
193	<	TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
194	<	TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
195	<	TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
196	<	signalpred->Add(signalpredlo,-1);
197	<	signalpred->Add(signalpredro);
198	<	puresignal->Add(signalpred,-1);//subtracting signal contamination
199	<	ofstream myfile;
200	<	myfile.open ("ShapeFit_log.txt");
201	<	establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
202	<	myfile.close();
203	<	*/
203	<	}
204	<
205	<	vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, bool doobserved=false) {
206	<	float sigma95=0.0,sigma95A=0.0;
207	<	dout << "Now calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << 1<< ") " << endl;
208	<	sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, 1);
209	<	if(doobserved) {
210	<	dout << "Now calling : CL95A(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << 1<< ") " << endl;
211	<	sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], 1);
212	<	}
171	>	vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, string plotfilename, bool doobserved) {
172	>	float sigma95=-9.9,sigma95A=-9.9;
173	>	int nuisancemodel=0;
174	>	/*
175	>	USAGE OF ROOSTATS_CL95
176	>	" Double_t             limit = roostats_cl95(ilum, slum, eff, seff, bck, sbck, n, gauss = false, nuisanceModel, method, plotFileName, seed); \n"
177	>	" LimitResult expected_limit = roostats_clm(ilum, slum, eff, seff, bck, sbck, ntoys, nuisanceModel, method, seed); \n"
178	>	" Double_t     average_limit = roostats_cla(ilum, slum, eff, seff, bck, sbck, nuisanceModel, method, seed); \n"
179	>	"                                                                     \n"
180	>	"
181	>	" Double_t obs_limit = limit.GetObservedLimit();                      \n"
182	>	" Double_t exp_limit = limit.GetExpectedLimit();                      \n"
183	>	" Double_t exp_up    = limit.GetOneSigmaHighRange();                  \n"
184	>	" Double_t exp_down  = limit.GetOneSigmaLowRange();                   \n"
185	>	" Double_t exp_2up   = limit.GetTwoSigmaHighRange();                  \n"
186	>	" Double_t exp_2down = limit.GetTwoSigmaLowRange();                   \n"
187	>	*/
188	>	if(mceff<=0) {
189	>	write_warning(__FUNCTION__,"Cannot compute upper limit in this configuration as the efficiency is negative:");
190	>	dout << "mc efficiency=" << mceff << " +/- " << mcefferr;
191	>	vector<float> sigmas;
192	>	sigmas.push_back(-1);
193	>	sigmas.push_back(-1);
194	>	return sigmas;
195	>	} else {
196	>	int nlimittoysused=1;
197	>	if(doobserved) nlimittoysused=nlimittoys;
198	>	/*  dout << "Now calling : CL95(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << Nobs[ibin] << "," << false << "," << nuisancemodel<< ") " << endl;
199	>	sigma95 = CL95(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], Nobs[ibin], false, nuisancemodel);
200	>	*/
201	>	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;
202	>	LimitResult limit = roostats_limit(luminosity,lumiuncert*luminosity,mceff,mcefferr,Npred[ibin],Nprederr[ibin],nlimittoysused,false,nuisancemodel,limitmethod,plotfilename,1);
203	>
204		vector<float> sigmas;
205	+	sigmas.push_back(limit.GetExpectedLimit());//expected
206	+	sigmas.push_back(limit.GetObservedLimit());//observed
207	+	//up to here for backward compatibility
208	+	sigmas.push_back(limit.GetOneSigmaHighRange());//expected, up
209	+	sigmas.push_back(limit.GetTwoSigmaHighRange());//expected, 2 up
210	+	sigmas.push_back(limit.GetOneSigmaLowRange());//expected, down
211	+	sigmas.push_back(limit.GetTwoSigmaLowRange());//expected, 2 down
212	+	//  if(doobserved) {
213	+	//    dout << "Now calling : CLA(" << luminosity << "," <<  lumiuncert*luminosity << "," << mceff << "," << mcefferr << "," << Npred[ibin] << "," << Nprederr[ibin] << "," << nuisancemodel<< ") " << endl;
214	+	//    sigma95A = CLA(luminosity, lumiuncert*luminosity, mceff, mcefferr, Npred[ibin], Nprederr[ibin], nuisancemodel);
215	+	//  }
216	+
217	+	/*  vector<float> sigmas;
218		sigmas.push_back(sigma95);
219	<	sigmas.push_back(sigma95A);
219	>	sigmas.push_back(sigma95A);*/
220		return sigmas;
221	+
222	+
223	+	}
224	+	write_warning(__FUNCTION__,"STILL MISSING SIGMAS, LIMITS, EVERYTHING ...");
225		}
226
227		void compute_upper_limits_from_counting_experiment(vector<vector<float> > uncertainties,vector<float> jzbcuts, string mcjzb, bool doobserved) {
#	Line 234 | Line 242 | void compute_upper_limits_from_counting_
242		float staterr=uncertainties[isample*jzbcuts.size()+ibin][2];
243		float systerr=uncertainties[isample*jzbcuts.size()+ibin][3];
244		float toterr =uncertainties[isample*jzbcuts.size()+ibin][4];
245	<	float observed,null,result;
238	<	fill_result_histos(observed, null,null,null,null,null,null,null,mcjzb,JZBcutat,(int)5,result,(signalsamples.FindSample(signalsamples.collection[isample].filename)),signalsamples);
239	<	observed-=result;//this is the actual excess we see!
240	<	float expected=observed/luminosity;
245	>	float observed,observederr,null,result;
246
247	+	//      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);
248	+	//      observed-=result;//this is the actual excess we see!
249	+	//      float expected=observed/luminosity;
250	+	string plotfilename=(string)(TString(signalsamples.collection[isample].samplename)+TString("___JZB_geq_")+TString(any2string(JZBcutat))+TString(".png"));
251		dout << "Sample: " << signalsamples.collection[isample].samplename << ", JZB>"<<JZBcutat<< " : " << mceff << " +/- " << staterr << " (stat) +/- " << systerr << " (syst) --> toterr = " << toterr << endl;
252	<	vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,doobserved);
252	>	vector<float> sigmas = compute_one_upper_limit(mceff,toterr,ibin,mcjzb,plotfilename,doobserved);
253
254		if(doobserved) {
255	<	rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
255	>	//      rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(expected)+")");
256	>	rows.push_back(any2string(sigmas[0])+";"+any2string(sigmas[1])+";"+"("+any2string(signalsamples.collection[isample].xs)+")");
257		vrows.push_back(sigmas[0]);
258		vrows.push_back(sigmas[1]);
259	<	vrows.push_back(expected);
259	>	//      vrows.push_back(expected);
260	>	vrows.push_back(signalsamples.collection[isample].xs);
261		}
262		else {
263	<	rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")");
263	>	//      rows.push_back(any2string(sigmas[0])+"("+any2string(expected)+")");
264	>	rows.push_back(any2string(sigmas[0]));
265		vrows.push_back(sigmas[0]);
266	<	vrows.push_back(expected);
266	>	vrows.push_back(signalsamples.collection[isample].xs);
267	>	//      vrows.push_back(expected);
268		}
269		}//end of bin loop
270		limits.push_back(rows);
271		vlimits.push_back(vrows);
272		}//end of sample loop
273	<	dout << endl << endl << "PAS table 3: " << endl << endl;
273	>	dout << endl << endl << endl << "_________________________________________________________________________________________________" << endl << endl;
274	>	dout << endl << endl << "PAS table 3:   (notation: limit [95%CL])" << endl << endl;
275		dout << "\t";
276		for (int irow=0;irow<jzbcuts.size();irow++) {
277		dout << jzbcuts[irow] << "\t";
#	Line 265 | Line 279 | void compute_upper_limits_from_counting_
279		dout << endl;
280		for(int irow=0;irow<limits.size();irow++) {
281		for(int ientry=0;ientry<limits[irow].size();ientry++) {
282	<	dout << limits[irow][ientry] << "\t";
282	>	if (limits[irow][ientry]>0) dout << limits[irow][ientry] << "\t";
283	>	else dout << " (N/A) \t";
284		}
285		dout << endl;
286		}
287
288		if(!doobserved) {
289	<	dout << endl << endl << "LIMITS: " << endl;
290	<	dout << "\t";
289	>	dout << endl << endl << "LIMITS: (Tex)" << endl;
290	>	tout << "\\begin{table}[hbtp]" << endl;
291	>	tout << "\\renewcommand{\arraystretch}{1.3}" << endl;
292	>	tout << "\\begin{center}" << endl;
293	>	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;
294	>	tout << "" << endl;
295	>	tout << "\\begin{tabular}{ | l | ";
296	>	for (int irow=0;irow<jzbcuts.size();irow++) tout << " l |";
297	>	tout << "} " << endl << " \\hline " << endl << "& \t ";
298		for (int irow=0;irow<jzbcuts.size();irow++) {
299	<	dout << jzbcuts[irow] << "\t";
299	>	tout << "JZB $>$ " << jzbcuts[irow] << " GeV & \t ";
300		}
301	<	dout << endl;
301	>	tout << " \\\\ \\hline " << endl;
302		for(int irow=0;irow<limits.size();irow++) {
303	<	dout << limits[irow][0] << "\t";
303	>	tout << limits[irow][0] << " \t";
304		for(int ientry=0;ientry<jzbcuts.size();ientry++) {
305	<	dout << Round(vlimits[irow][2*ientry] / vlimits[irow][2*ientry+1],3)<< "\t";
305	>	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";
306	>	else tout << " & ( N / A ) \t";
307	>	//      dout << Round(vlimits[irow][2*ientry],3) << " / " << Round(vlimits[irow][2*ientry+1],3)<< "\t";
308		}
309	<	dout << endl;
309	>	tout << " \\\\ \\hline " << endl;
310		}
311	+	tout << "\\end{tabular}" << endl;
312	+	tout << "      \\end{tabular}"<< endl;
313	+	tout << "\\end{center}"<< endl;
314	+	tout << "\\end{table} "<< endl;
315	+
316		}//do observed
317
318		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;
319	<	dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t Prediction [pb]" << endl;
319	>	dout << "Scenario \t Efficiency [%] \t Upper limits [pb] \t \\sigma [pb]" << endl;
320		for(int icut=0;icut<jzbcuts.size();icut++) {
321	<	dout << "Region with JZB>" << jzbcuts[icut] << endl;
321	>	dout << "Region with JZB>" << jzbcuts[icut] << (ConsiderSignalContaminationForLimits?"  (accounting for signal contamination)":"  (not accounting for signal contamination)") << endl;
322		for(int isample=0;isample<signalsamples.collection.size();isample++) {
323	<	dout << limits[icut][0] << "\t" << Round(100*uncertainties[isample*jzbcuts.size()+icut][1],1) << "+/-" << Round(100*uncertainties[isample*jzbcuts.size()+icut][2],1) << " (stat) +/- " << Round(100*uncertainties[isample*jzbcuts.size()+icut][3],1) << " (syst) \t" << Round((vlimits[isample][2*icut]),3) << "\t" << Round(vlimits[isample][2*icut+1],3) << endl;
323	>	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;
324		}
325		dout << endl;
326		}
327	+
328	+	write_warning(__FUNCTION__,"Still need to update the script");
329		}
330
300	–	void susy_scan_axis_labeling(TH2F *histo) {
301	–	histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
302	–	histo->GetXaxis()->CenterTitle();
303	–	histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
304	–	histo->GetYaxis()->CenterTitle();
305	–	}
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	<	float arrbinning[binning.size()];
337	<	for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
338	<	TH1F *puredata   = allsamples.Draw("puredata",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
339	<	puredata->SetMarkerSize(DataMarkerSize);
340	<	TH1F *allbgs   = allsamples.Draw("allbgs",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
341	<	TH1F *allbgsb   = allsamples.Draw("allbgsb",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
342	<	TH1F *allbgsc   = allsamples.Draw("allbgsc",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
343	<	allbgs->Add(allbgsb,-1);
344	<	allbgs->Add(allbgsc);
345	<	int ndata=puredata->Integral();
346	<	ofstream myfile;
347	<	myfile.open ("susyscan_log.txt");
348	<	TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
349	<	TTree *suevents = (TTree*)susyscanfile->Get("events");
350	<	TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
351	<	TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
352	<	TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
353	<	TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
354	<
355	<	susy_scan_axis_labeling(exclusionmap);
356	<	susy_scan_axis_labeling(exclusionmap1s);
357	<	susy_scan_axis_labeling(exclusionmap2s);
358	<	susy_scan_axis_labeling(exclusionmap3s);
359	<
360	<	Int_t MyPalette[100];
361	<	Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
362	<	Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
363	<	Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
364	<	Double_t stop[] = {0., .25, .50, .75, 1.0};
365	<	Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
366	<	for (int i=0;i<100;i++) MyPalette[i] = FI+i;
367	<
368	<	gStyle->SetPalette(100, MyPalette);
344	<
345	<	for(int m23=50;m23<500;m23+=25) {
346	<	for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
347	<	{
348	<	c3->cd();
349	<	stringstream drawcondition;
350	<	drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
351	<	TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
352	<	TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
353	<	stringstream drawvar,drawvar2;
354	<	drawvar<<mcjzb<<">>puresignal";
355	<	drawvar2<<"-"<<mcjzb<<">>puresignall";
356	<	suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
357	<	suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
358	<	if(puresignal->Integral()<60) {
359	<	delete puresignal;
360	<	continue;
361	<	}
362	<	puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
363	<	puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
364	<	stringstream saveas;
365	<	saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
366	<	dout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
367	<	//        TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
368	<	//        TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
369	<	//        TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
370	<	//        signalpred->Add(signalpredlo,-1);
371	<	//        signalpred->Add(signalpredro);
372	<	//        puresignal->Add(signalpred,-1);//subtracting signal contamination
373	<	//---------------------
374	<	//      dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
375	<	//    TH1F *puresignal = allsamples.Draw("puresignal",mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
376	<	vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
377	<	if(results.size()==0) {
378	<	delete puresignal;
379	<	continue;
380	<	}
381	<	exclusionmap->Fill(m23,m0,results[0]);
382	<	exclusionmap1s->Fill(m23,m0,results[1]);
383	<	exclusionmap2s->Fill(m23,m0,results[2]);
384	<	exclusionmap3s->Fill(m23,m0,results[3]);
385	<	delete puresignal;
386	<	dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
387	<	}
388	<	}//end of model scan for loop
332	>
333	>	/********************************************************************** new : Limits using SHAPES ***********************************
334	>
335	>
336	>	SSSSSSSSSSSSSSS hhhhhhh
337	>	SS:::::::::::::::Sh:::::h
338	>	S:::::SSSSSS::::::Sh:::::h
339	>	S:::::S     SSSSSSSh:::::h
340	>	S:::::S             h::::h hhhhh         aaaaaaaaaaaaa  ppppp   ppppppppp       eeeeeeeeeeee        ssssssssss
341	>	S:::::S             h::::hh:::::hhh      a::::::::::::a p::::ppp:::::::::p    ee::::::::::::ee    ss::::::::::s
342	>	S::::SSSS          h::::::::::::::hh    aaaaaaaaa:::::ap:::::::::::::::::p  e::::::eeeee:::::eess:::::::::::::s
343	>	SS::::::SSSSS     h:::::::hhh::::::h            a::::app::::::ppppp::::::pe::::::e     e:::::es::::::ssss:::::s
344	>	SSS::::::::SS   h::::::h   h::::::h    aaaaaaa:::::a p:::::p     p:::::pe:::::::eeeee::::::e s:::::s  ssssss
345	>	SSSSSS::::S  h:::::h     h:::::h  aa::::::::::::a p:::::p     p:::::pe:::::::::::::::::e    s::::::s
346	>	S:::::S h:::::h     h:::::h a::::aaaa::::::a p:::::p     p:::::pe::::::eeeeeeeeeee        s::::::s
347	>	S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::p    p::::::pe:::::::e           ssssss   s:::::s
348	>	SSSSSSS     S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::ppppp:::::::pe::::::::e          s:::::ssss::::::s
349	>	S::::::SSSSSS:::::S h:::::h     h:::::ha:::::aaaa::::::a p::::::::::::::::p  e::::::::eeeeeeee  s::::::::::::::s
350	>	S:::::::::::::::SS  h:::::h     h:::::h a::::::::::aa:::ap::::::::::::::pp    ee:::::::::::::e   s:::::::::::ss
351	>	SSSSSSSSSSSSSSS    hhhhhhh     hhhhhhh  aaaaaaaaaa  aaaap::::::pppppppp        eeeeeeeeeeeeee    sssssssssss
352	>	p:::::p
353	>	p:::::p
354	>	p:::::::p
355	>	p:::::::p
356	>	p:::::::p
357	>	ppppppppp
358	>
359	>
360	>	*********************************************************************** new : Limits using SHAPES ***********************************/
361	>
362	>
363	>	void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
364	>	dout << "Creatig shape templates ... ";
365	>	if(identifier!="") dout << "for systematic called "<<identifier;
366	>	dout << endl;
367	>	int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
368	>	if(dataormc!=data) write_warning(__FUNCTION__,"WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
369
370	<	dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
371	<	c3->cd();
372	<	exclusionmap->Draw("CONTZ");
373	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
374	<	exclusionmap->Draw("COLZ");
375	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
376	<
377	<	exclusionmap1s->Draw("CONTZ");
378	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
379	<	exclusionmap1s->Draw("COLZ");
380	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
381	<
382	<	exclusionmap2s->Draw("CONTZ");
383	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
384	<	exclusionmap2s->Draw("COLZ");
385	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
386	<
387	<	exclusionmap3s->Draw("CONTZ");
388	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
389	<	exclusionmap3s->Draw("COLZ");
390	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
391	<
392	<	TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
393	<	exclusionmap->Write();
394	<	exclusionmap1s->Write();
395	<	exclusionmap2s->Write();
396	<	exclusionmap3s->Write();
397	<	exclusion_limits->Close();
398	<	susyscanfile->Close();
370	>	TCut limitnJetcut;
371	>	if(JES==noJES) limitnJetcut=cutnJets;
372	>	else {
373	>	if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
374	>	if(JES==JESup) limitnJetcut=cutnJetsJESup;
375	>	}
376	>	TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
377	>	TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
378	>	TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
379	>	TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
380	>
381	>	TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
382	>	TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
383	>	TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
384	>	TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
385	>
386	>	TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
387	>	TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
388	>	TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
389	>	TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
390	>
391	>	TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
392	>	TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
393	>	TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
394	>	TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
395	>
396	>	string obsname="data_obs";
397	>	string predname="background";
398	>	string signalname="signal";
399	>	if(identifier!="") {
400	>	obsname=("data_"+identifier);
401	>	predname=("background_"+identifier);
402	>	signalname="signal_"+identifier;
403	>	}
404
405	<	myfile.close();
405	>	TH1F *obs = (TH1F*)ZOSSFP->Clone();
406	>	obs->SetName(obsname.c_str());
407	>	obs->Write();
408	>	TH1F *pred = (TH1F*)ZOSSFN->Clone();
409	>	pred->Add(ZOSOFP,1.0/3);
410	>	pred->Add(ZOSOFN,-1.0/3);
411	>	pred->Add(SBOSSFP,1.0/3);
412	>	pred->Add(SBOSSFN,-1.0/3);
413	>	pred->Add(SBOSOFP,1.0/3);
414	>	pred->Add(SBOSOFN,-1.0/3);
415	>	pred->SetName(predname.c_str());
416	>	pred->Write();
417	>
418	>	//  TH1F *Lobs = (TH1F*)LZOSSFP->Clone();
419	>	//  TH1F *Lpred = (TH1F*)LZOSSFN->Clone();
420	>
421	>	TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
422	>	TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
423	>	Lobs->Add(LZOSSFP);
424	>	Lpred->Add(LZOSSFN);
425	>	Lpred->Add(LZOSOFP,1.0/3);
426	>	Lpred->Add(LZOSOFN,-1.0/3);
427	>	Lpred->Add(LSBOSSFP,1.0/3);
428	>	Lpred->Add(LSBOSSFN,-1.0/3);
429	>	Lpred->Add(LSBOSOFP,1.0/3);
430	>	Lpred->Add(LSBOSOFN,-1.0/3);
431	>	TH1F *signal = (TH1F*)Lobs->Clone();
432	>	signal->Add(Lpred,-1);
433	>	signal->SetName(signalname.c_str());
434	>	signal->Write();
435	>
436	>	delete Lobs;
437	>	delete Lpred;
438	>
439	>	delete ZOSSFP;
440	>	delete ZOSOFP;
441	>	delete ZOSSFN;
442	>	delete ZOSOFN;
443	>
444	>	delete SBOSSFP;
445	>	delete SBOSOFP;
446	>	delete SBOSSFN;
447	>	delete SBOSOFN;
448	>
449	>	delete LZOSSFP;
450	>	delete LZOSOFP;
451	>	delete LZOSSFN;
452	>	delete LZOSOFN;
453	>
454	>	delete LSBOSSFP;
455	>	delete LSBOSOFP;
456	>	delete LSBOSSFN;
457	>	delete LSBOSOFN;
458	>
459		}
460
461	+	void prepare_datacard(TFile *f) {
462	+	TH1F *dataob = (TH1F*)f->Get("data_obs");
463	+	TH1F *signal = (TH1F*)f->Get("signal");
464	+	TH1F *background = (TH1F*)f->Get("background");
465	+
466	+	ofstream datacard;
467	+	ensure_directory_exists(get_directory()+"/limits");
468	+	datacard.open ((get_directory()+"/limits/susydatacard.txt").c_str());
469	+	datacard << "Writing this to a file.\n";
470	+	datacard << "imax 1\n";
471	+	datacard << "jmax 1\n";
472	+	datacard << "kmax *\n";
473	+	datacard << "---------------\n";
474	+	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
475	+	datacard << "---------------\n";
476	+	datacard << "bin 1\n";
477	+	datacard << "observation "<<dataob->Integral()<<"\n";
478	+	datacard << "------------------------------\n";
479	+	datacard << "bin             1          1\n";
480	+	datacard << "process         signal     background\n";
481	+	datacard << "process         0          1\n";
482	+	datacard << "rate            "<<signal->Integral()<<"         "<<background->Integral()<<"\n";
483	+	datacard << "--------------------------------\n";
484	+	datacard << "lumi     lnN    1.10       1.0\n";
485	+	datacard << "bgnorm   lnN    1.00       1.4  uncertainty on our prediction (40%)\n";
486	+	datacard << "JES    shape    1          1    uncertainty on background shape and normalization\n";
487	+	datacard << "peak   shape    1          1    uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
488	+	datacard << "#                                so divide the unit gaussian by 2 before doing the interpolation\n";
489	+	datacard.close();
490	+	}
491
492
493	<
493	>	void prepare_limits(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbbins) {
494	>	ensure_directory_exists(get_directory()+"/limits");
495	>	TFile *limfile = new TFile((get_directory()+"/limits/limitfile.root").c_str(),"RECREATE");
496	>	TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
497	>	limit_shapes_for_systematic_effect(limfile,"",mcjzb,datajzb,noJES,jzbbins,limcan);
498	>	limit_shapes_for_systematic_effect(limfile,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan);
499	>	limit_shapes_for_systematic_effect(limfile,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan);
500	>	limit_shapes_for_systematic_effect(limfile,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan);
501	>	limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
502	>
503	>	prepare_datacard(limfile);
504	>	limfile->Close();
505	>	write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
506	>
507	>	}

Diff Legend

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