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.6 by buchmann, Fri Jul 22 10:15:34 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 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	–	}
274	–
275	–	void scan_susy_space(string mcjzb, string datajzb) {
276	–	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	–	}
390	–
391	–
269
270
271		//********************************************************************** new : Limits using SHAPES ***********************************
#	Line 398 | Line 275 | void limit_shapes_for_systematic_effect(
275		if(identifier!="") dout << "for systematic called "<<identifier;
276		dout << endl;
277		int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
278	<	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!");
278	>	if(dataormc!=data) write_warning(__FUNCTION__,"WATCH OUT! Not using data for limits!!!! this is ok for tests, but not ok for anything official!");
279
280		TCut limitnJetcut;
281		if(JES==noJES) limitnJetcut=cutnJets;
#	Line 498 | Line 375 | void prepare_datacard(TFile *f) {
375
376		ofstream datacard;
377		ensure_directory_exists(get_directory()+"/limits");
378	<	datacard.open ((get_directory()+"/limits/susylm4datacard.txt").c_str());
378	>	datacard.open ((get_directory()+"/limits/susydatacard.txt").c_str());
379		datacard << "Writing this to a file.\n";
380		datacard << "imax 1\n";
381		datacard << "jmax 1\n";
#	Line 534 | Line 411 | void prepare_limits(string mcjzb, string
411		limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
412
413		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!");
414		limfile->Close();
415	+	write_info("prepare_limits","limitfile.root and datacard.txt have been generated. You can now use them to calculate limits!");
416
417		}

Diff Legend

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