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.5 by buchmann, Wed Jul 20 14:34:31 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 257 | 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 265 | 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
268	–	void susy_scan_axis_labeling(TH2F *histo) {
269	–	histo->GetXaxis()->SetTitle("#Chi_{2}^{0}-LSP");
270	–	histo->GetXaxis()->CenterTitle();
271	–	histo->GetYaxis()->SetTitle("m_{#tilde{q}}");
272	–	histo->GetYaxis()->CenterTitle();
273	–	}
269
270	<	void scan_susy_space(string mcjzb, string datajzb) {
271	<	TCanvas *c3 = new TCanvas("c3","c3");
277	<	vector<float> binning;
278	<	binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
279	<	float arrbinning[binning.size()];
280	<	for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
281	<	TH1F *puredata   = allsamples.Draw("puredata",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
282	<	puredata->SetMarkerSize(DataMarkerSize);
283	<	TH1F *allbgs   = allsamples.Draw("allbgs",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
284	<	TH1F *allbgsb   = allsamples.Draw("allbgsb",  "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
285	<	TH1F *allbgsc   = allsamples.Draw("allbgsc",  datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
286	<	allbgs->Add(allbgsb,-1);
287	<	allbgs->Add(allbgsc);
288	<	int ndata=puredata->Integral();
289	<	ofstream myfile;
290	<	myfile.open ("susyscan_log.txt");
291	<	TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
292	<	TTree *suevents = (TTree*)susyscanfile->Get("events");
293	<	TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
294	<	TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
295	<	TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
296	<	TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
297	<
298	<	susy_scan_axis_labeling(exclusionmap);
299	<	susy_scan_axis_labeling(exclusionmap1s);
300	<	susy_scan_axis_labeling(exclusionmap2s);
301	<	susy_scan_axis_labeling(exclusionmap3s);
302	<
303	<	Int_t MyPalette[100];
304	<	Double_t r[]    = {0., 0.0, 1.0, 1.0, 1.0};
305	<	Double_t g[]    = {0., 0.0, 0.0, 1.0, 1.0};
306	<	Double_t b[]    = {0., 1.0, 0.0, 0.0, 1.0};
307	<	Double_t stop[] = {0., .25, .50, .75, 1.0};
308	<	Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
309	<	for (int i=0;i<100;i++) MyPalette[i] = FI+i;
310	<
311	<	gStyle->SetPalette(100, MyPalette);
312	<
313	<	for(int m23=50;m23<500;m23+=25) {
314	<	for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
315	<	{
316	<	c3->cd();
317	<	stringstream drawcondition;
318	<	drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
319	<	TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
320	<	TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
321	<	stringstream drawvar,drawvar2;
322	<	drawvar<<mcjzb<<">>puresignal";
323	<	drawvar2<<"-"<<mcjzb<<">>puresignall";
324	<	suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
325	<	suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
326	<	if(puresignal->Integral()<60) {
327	<	delete puresignal;
328	<	continue;
329	<	}
330	<	puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
331	<	puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
332	<	stringstream saveas;
333	<	saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
334	<	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;
335	<	//        TH1F *signalpredlo   = allsamples.Draw("signalpredlo",  "-"+mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
336	<	//        TH1F *signalpredro   = allsamples.Draw("signalpredro",      mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
337	<	//        TH1F *puredata       = allsamples.Draw("puredata",          datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
338	<	//        signalpred->Add(signalpredlo,-1);
339	<	//        signalpred->Add(signalpredro);
340	<	//        puresignal->Add(signalpred,-1);//subtracting signal contamination
341	<	//---------------------
342	<	//      dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
343	<	//    TH1F *puresignal = allsamples.Draw("puresignal",mcjzb,  binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc,  luminosity,allsamples.FindSample("LM4"));
344	<	vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
345	<	if(results.size()==0) {
346	<	delete puresignal;
347	<	continue;
348	<	}
349	<	exclusionmap->Fill(m23,m0,results[0]);
350	<	exclusionmap1s->Fill(m23,m0,results[1]);
351	<	exclusionmap2s->Fill(m23,m0,results[2]);
352	<	exclusionmap3s->Fill(m23,m0,results[3]);
353	<	delete puresignal;
354	<	dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
355	<	}
356	<	}//end of model scan for loop
357	<
358	<	dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
359	<	c3->cd();
360	<	exclusionmap->Draw("CONTZ");
361	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
362	<	exclusionmap->Draw("COLZ");
363	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
364	<
365	<	exclusionmap1s->Draw("CONTZ");
366	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
367	<	exclusionmap1s->Draw("COLZ");
368	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
369	<
370	<	exclusionmap2s->Draw("CONTZ");
371	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
372	<	exclusionmap2s->Draw("COLZ");
373	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
374	<
375	<	exclusionmap3s->Draw("CONTZ");
376	<	CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
377	<	exclusionmap3s->Draw("COLZ");
378	<	CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
379	<
380	<	TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
381	<	exclusionmap->Write();
382	<	exclusionmap1s->Write();
383	<	exclusionmap2s->Write();
384	<	exclusionmap3s->Write();
385	<	exclusion_limits->Close();
386	<	susyscanfile->Close();
387	<
388	<	myfile.close();
389	<	}
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
393	–
394	–	//********************************************************************** new : Limits using SHAPES ***********************************
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("limit_shapes_for_systematic_effect","WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
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		TCut limitnJetcut;
309		if(JES==noJES) limitnJetcut=cutnJets;
#	Line 498 | Line 403 | void prepare_datacard(TFile *f) {
403
404		ofstream datacard;
405		ensure_directory_exists(get_directory()+"/limits");
406	<	datacard.open ((get_directory()+"/limits/susylm4datacard.txt").c_str());
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";
#	Line 534 | Line 439 | void prepare_limits(string mcjzb, string
439		limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
440
441		prepare_datacard(limfile);
537	–
538	–	write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
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