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.4 by buchmann, Wed Jul 20 12:26:11 2011 UTC vs.
Revision 1.8 by buchmann, Wed Jul 27 12:16:15 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
170	–	void calculate_upper_limits(string mcjzb, string datajzb) {
171	–	write_warning("calculate_upper_limits","Upper limit calculation temporarily deactivated");
172	–	//  write_warning("calculate_upper_limits","Calculation of SUSY upper limits has been temporarily suspended in favor of top discovery");
173	–	//  rediscover_the_top(mcjzb,datajzb);
174	–	/*
175	–	TCanvas *c3 = new TCanvas("c3","c3");
176	–	c3->SetLogy(1);
177	–	vector<float> binning;
178	–	//binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
179	–	binning.push_back(50);
180	–	binning.push_back(100);
181	–	binning.push_back(150);
182	–	binning.push_back(200);
183	–	binning.push_back(500);
184	–	TH1F *datapredictiona = allsamples.Draw("datapredictiona",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity);
185	–	TH1F *datapredictionb = allsamples.Draw("datapredictionb",    "-"+datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
186	–	TH1F *datapredictionc = allsamples.Draw("datapredictionc",    datajzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity);
187	–	TH1F *dataprediction = (TH1F*)datapredictiona->Clone();
188	–	dataprediction->Add(datapredictionb,-1);
189	–	dataprediction->Add(datapredictionc);
190	–	TH1F *puresignal     = allsamples.Draw("puresignal",        mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
191	–	TH1F *signalpred     = allsamples.Draw("signalpred",    "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
192	–	TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
193	–	TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
194	–	TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
195	–	signalpred->Add(signalpredlo,-1);
196	–	signalpred->Add(signalpredro);
197	–	puresignal->Add(signalpred,-1);//subtracting signal contamination
198	–	ofstream myfile;
199	–	myfile.open ("ShapeFit_log.txt");
200	–	establish_upper_limits(puredata,dataprediction,puresignal,"LM4",myfile);
201	–	myfile.close();
202	–	*/
203	–	}
204	–
171		vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, bool doobserved=false) {
172		float sigma95=0.0,sigma95A=0.0;
173		int nuisancemodel=1;
#	Line 292 | Line 258 | void compute_upper_limits_from_counting_
258		for(int icut=0;icut<jzbcuts.size();icut++) {
259		dout << "Region with JZB>" << jzbcuts[icut] << endl;
260		for(int isample=0;isample<signalsamples.collection.size();isample++) {
261	<	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;
261	>	dout << limits[isample][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;
262		}
263		dout << endl;
264		}
#	Line 300 | Line 266 | void compute_upper_limits_from_counting_
266		write_warning("compute_upper_limits_from_counting_experiment","Still need to update the script");
267		}
268
303	–	void susy_scan_axis_labeling(TH2F *histo) {
304	–	histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
305	–	histo->GetXaxis()->CenterTitle();
306	–	histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
307	–	histo->GetYaxis()->CenterTitle();
308	–	}
269
270	<	void scan_susy_space(string mcjzb, string datajzb) {
271	<	TCanvas *c3 = new TCanvas("c3","c3");
272	<	vector<float> binning;
273	<	binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
274	<	float arrbinning[binning.size()];
275	<	for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
276	<	TH1F *puredata   = allsamples.Draw("puredata",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
277	<	puredata->SetMarkerSize(DataMarkerSize);
278	<	TH1F *allbgs   = allsamples.Draw("allbgs",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
279	<	TH1F *allbgsb   = allsamples.Draw("allbgsb",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
280	<	TH1F *allbgsc   = allsamples.Draw("allbgsc",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
281	<	allbgs->Add(allbgsb,-1);
282	<	allbgs->Add(allbgsc);
283	<	int ndata=puredata->Integral();
284	<	ofstream myfile;
285	<	myfile.open ("susyscan_log.txt");
286	<	TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
287	<	TTree *suevents = (TTree*)susyscanfile->Get("events");
288	<	TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
289	<	TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
290	<	TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
291	<	TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
292	<
293	<	susy_scan_axis_labeling(exclusionmap);
294	<	susy_scan_axis_labeling(exclusionmap1s);
295	<	susy_scan_axis_labeling(exclusionmap2s);
296	<	susy_scan_axis_labeling(exclusionmap3s);
297	<
298	<	Int_t MyPalette[100];
299	<	Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
300	<	Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
301	<	Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
302	<	Double_t stop[] = {0., .25, .50, .75, 1.0};
303	<	Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
304	<	for (int i=0;i<100;i++) MyPalette[i] = FI+i;
305	<
306	<	gStyle->SetPalette(100, MyPalette);
347	<
348	<	for(int m23=50;m23<500;m23+=25) {
349	<	for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
350	<	{
351	<	c3->cd();
352	<	stringstream drawcondition;
353	<	drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
354	<	TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
355	<	TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
356	<	stringstream drawvar,drawvar2;
357	<	drawvar<<mcjzb<<">>puresignal";
358	<	drawvar2<<"-"<<mcjzb<<">>puresignall";
359	<	suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
360	<	suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
361	<	if(puresignal->Integral()<60) {
362	<	delete puresignal;
363	<	continue;
364	<	}
365	<	puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
366	<	puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
367	<	stringstream saveas;
368	<	saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
369	<	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;
370	<	//        TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
371	<	//        TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
372	<	//        TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
373	<	//        signalpred->Add(signalpredlo,-1);
374	<	//        signalpred->Add(signalpredro);
375	<	//        puresignal->Add(signalpred,-1);//subtracting signal contamination
376	<	//---------------------
377	<	//      dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
378	<	//    TH1F *puresignal = allsamples.Draw("puresignal",mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
379	<	vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
380	<	if(results.size()==0) {
381	<	delete puresignal;
382	<	continue;
383	<	}
384	<	exclusionmap->Fill(m23,m0,results[0]);
385	<	exclusionmap1s->Fill(m23,m0,results[1]);
386	<	exclusionmap2s->Fill(m23,m0,results[2]);
387	<	exclusionmap3s->Fill(m23,m0,results[3]);
388	<	delete puresignal;
389	<	dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
390	<	}
391	<	}//end of model scan for loop
270	>
271	>	/********************************************************************** new : Limits using SHAPES ***********************************
272	>
273	>
274	>	SSSSSSSSSSSSSSS hhhhhhh
275	>	SS:::::::::::::::Sh:::::h
276	>	S:::::SSSSSS::::::Sh:::::h
277	>	S:::::S     SSSSSSSh:::::h
278	>	S:::::S             h::::h hhhhh         aaaaaaaaaaaaa  ppppp   ppppppppp       eeeeeeeeeeee        ssssssssss
279	>	S:::::S             h::::hh:::::hhh      a::::::::::::a p::::ppp:::::::::p    ee::::::::::::ee    ss::::::::::s
280	>	S::::SSSS          h::::::::::::::hh    aaaaaaaaa:::::ap:::::::::::::::::p  e::::::eeeee:::::eess:::::::::::::s
281	>	SS::::::SSSSS     h:::::::hhh::::::h            a::::app::::::ppppp::::::pe::::::e     e:::::es::::::ssss:::::s
282	>	SSS::::::::SS   h::::::h   h::::::h    aaaaaaa:::::a p:::::p     p:::::pe:::::::eeeee::::::e s:::::s  ssssss
283	>	SSSSSS::::S  h:::::h     h:::::h  aa::::::::::::a p:::::p     p:::::pe:::::::::::::::::e    s::::::s
284	>	S:::::S h:::::h     h:::::h a::::aaaa::::::a p:::::p     p:::::pe::::::eeeeeeeeeee        s::::::s
285	>	S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::p    p::::::pe:::::::e           ssssss   s:::::s
286	>	SSSSSSS     S:::::S h:::::h     h:::::ha::::a    a:::::a p:::::ppppp:::::::pe::::::::e          s:::::ssss::::::s
287	>	S::::::SSSSSS:::::S h:::::h     h:::::ha:::::aaaa::::::a p::::::::::::::::p  e::::::::eeeeeeee  s::::::::::::::s
288	>	S:::::::::::::::SS  h:::::h     h:::::h a::::::::::aa:::ap::::::::::::::pp    ee:::::::::::::e   s:::::::::::ss
289	>	SSSSSSSSSSSSSSS    hhhhhhh     hhhhhhh  aaaaaaaaaa  aaaap::::::pppppppp        eeeeeeeeeeeeee    sssssssssss
290	>	p:::::p
291	>	p:::::p
292	>	p:::::::p
293	>	p:::::::p
294	>	p:::::::p
295	>	ppppppppp
296	>
297	>
298	>	*********************************************************************** new : Limits using SHAPES ***********************************/
299	>
300	>
301	>	void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
302	>	dout << "Creatig shape templates ... ";
303	>	if(identifier!="") dout << "for systematic called "<<identifier;
304	>	dout << endl;
305	>	int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
306	>	if(dataormc!=data) write_warning(__FUNCTION__,"WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
307
308	<	dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
309	<	c3->cd();
310	<	exclusionmap->Draw("CONTZ");
311	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
312	<	exclusionmap->Draw("COLZ");
313	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
314	<
315	<	exclusionmap1s->Draw("CONTZ");
316	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
317	<	exclusionmap1s->Draw("COLZ");
318	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
319	<
320	<	exclusionmap2s->Draw("CONTZ");
321	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
322	<	exclusionmap2s->Draw("COLZ");
323	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
324	<
325	<	exclusionmap3s->Draw("CONTZ");
326	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
327	<	exclusionmap3s->Draw("COLZ");
328	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
329	<
330	<	TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
331	<	exclusionmap->Write();
332	<	exclusionmap1s->Write();
333	<	exclusionmap2s->Write();
334	<	exclusionmap3s->Write();
335	<	exclusion_limits->Close();
336	<	susyscanfile->Close();
308	>	TCut limitnJetcut;
309	>	if(JES==noJES) limitnJetcut=cutnJets;
310	>	else {
311	>	if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
312	>	if(JES==JESup) limitnJetcut=cutnJetsJESup;
313	>	}
314	>	TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
315	>	TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
316	>	TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
317	>	TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
318	>
319	>	TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
320	>	TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
321	>	TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
322	>	TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
323	>
324	>	TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
325	>	TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
326	>	TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
327	>	TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
328	>
329	>	TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
330	>	TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
331	>	TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
332	>	TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
333	>
334	>	string obsname="data_obs";
335	>	string predname="background";
336	>	string signalname="signal";
337	>	if(identifier!="") {
338	>	obsname=("data_"+identifier);
339	>	predname=("background_"+identifier);
340	>	signalname="signal_"+identifier;
341	>	}
342
343	<	myfile.close();
343	>	TH1F *obs = (TH1F*)ZOSSFP->Clone();
344	>	obs->SetName(obsname.c_str());
345	>	obs->Write();
346	>	TH1F *pred = (TH1F*)ZOSSFN->Clone();
347	>	pred->Add(ZOSOFP,1.0/3);
348	>	pred->Add(ZOSOFN,-1.0/3);
349	>	pred->Add(SBOSSFP,1.0/3);
350	>	pred->Add(SBOSSFN,-1.0/3);
351	>	pred->Add(SBOSOFP,1.0/3);
352	>	pred->Add(SBOSOFN,-1.0/3);
353	>	pred->SetName(predname.c_str());
354	>	pred->Write();
355	>
356	>	//  TH1F *Lobs = (TH1F*)LZOSSFP->Clone();
357	>	//  TH1F *Lpred = (TH1F*)LZOSSFN->Clone();
358	>
359	>	TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
360	>	TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
361	>	Lobs->Add(LZOSSFP);
362	>	Lpred->Add(LZOSSFN);
363	>	Lpred->Add(LZOSOFP,1.0/3);
364	>	Lpred->Add(LZOSOFN,-1.0/3);
365	>	Lpred->Add(LSBOSSFP,1.0/3);
366	>	Lpred->Add(LSBOSSFN,-1.0/3);
367	>	Lpred->Add(LSBOSOFP,1.0/3);
368	>	Lpred->Add(LSBOSOFN,-1.0/3);
369	>	TH1F *signal = (TH1F*)Lobs->Clone();
370	>	signal->Add(Lpred,-1);
371	>	signal->SetName(signalname.c_str());
372	>	signal->Write();
373	>
374	>	delete Lobs;
375	>	delete Lpred;
376	>
377	>	delete ZOSSFP;
378	>	delete ZOSOFP;
379	>	delete ZOSSFN;
380	>	delete ZOSOFN;
381	>
382	>	delete SBOSSFP;
383	>	delete SBOSOFP;
384	>	delete SBOSSFN;
385	>	delete SBOSOFN;
386	>
387	>	delete LZOSSFP;
388	>	delete LZOSOFP;
389	>	delete LZOSSFN;
390	>	delete LZOSOFN;
391	>
392	>	delete LSBOSSFP;
393	>	delete LSBOSOFP;
394	>	delete LSBOSSFN;
395	>	delete LSBOSOFN;
396	>
397		}
398
399	+	void prepare_datacard(TFile *f) {
400	+	TH1F *dataob = (TH1F*)f->Get("data_obs");
401	+	TH1F *signal = (TH1F*)f->Get("signal");
402	+	TH1F *background = (TH1F*)f->Get("background");
403	+
404	+	ofstream datacard;
405	+	ensure_directory_exists(get_directory()+"/limits");
406	+	datacard.open ((get_directory()+"/limits/susydatacard.txt").c_str());
407	+	datacard << "Writing this to a file.\n";
408	+	datacard << "imax 1\n";
409	+	datacard << "jmax 1\n";
410	+	datacard << "kmax *\n";
411	+	datacard << "---------------\n";
412	+	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
413	+	datacard << "---------------\n";
414	+	datacard << "bin 1\n";
415	+	datacard << "observation "<<dataob->Integral()<<"\n";
416	+	datacard << "------------------------------\n";
417	+	datacard << "bin             1          1\n";
418	+	datacard << "process         signal     background\n";
419	+	datacard << "process         0          1\n";
420	+	datacard << "rate            "<<signal->Integral()<<"         "<<background->Integral()<<"\n";
421	+	datacard << "--------------------------------\n";
422	+	datacard << "lumi     lnN    1.10       1.0\n";
423	+	datacard << "bgnorm   lnN    1.00       1.4  uncertainty on our prediction (40%)\n";
424	+	datacard << "JES    shape    1          1    uncertainty on background shape and normalization\n";
425	+	datacard << "peak   shape    1          1    uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
426	+	datacard << "#                                so divide the unit gaussian by 2 before doing the interpolation\n";
427	+	datacard.close();
428	+	}
429
430
431	<
431	>	void prepare_limits(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbbins) {
432	>	ensure_directory_exists(get_directory()+"/limits");
433	>	TFile *limfile = new TFile((get_directory()+"/limits/limitfile.root").c_str(),"RECREATE");
434	>	TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
435	>	limit_shapes_for_systematic_effect(limfile,"",mcjzb,datajzb,noJES,jzbbins,limcan);
436	>	limit_shapes_for_systematic_effect(limfile,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan);
437	>	limit_shapes_for_systematic_effect(limfile,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan);
438	>	limit_shapes_for_systematic_effect(limfile,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan);
439	>	limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
440	>
441	>	prepare_datacard(limfile);
442	>	limfile->Close();
443	>	write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
444	>
445	>	}

Diff Legend

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