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.5 by buchmann, Wed Jul 20 14:34:31 2011 UTC

#	Line 167 | Line 167 | ratio_binning.push_back(80);
167
168		}
169
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	–
170		vector<float> compute_one_upper_limit(float mceff,float mcefferr, int ibin, string mcjzb, bool doobserved=false) {
171		float sigma95=0.0,sigma95A=0.0;
172		int nuisancemodel=1;
#	Line 426 | Line 391 | void scan_susy_space(string mcjzb, strin
391
392
393
394	+	//********************************************************************** new : Limits using SHAPES ***********************************
395	+
396	+	void limit_shapes_for_systematic_effect(TFile *limfile, string identifier, string mcjzb, string datajzb, int JES,vector<float> binning, TCanvas *limcan) {
397	+	dout << "Creatig shape templates ... ";
398	+	if(identifier!="") dout << "for systematic called "<<identifier;
399	+	dout << endl;
400	+	int dataormc=mcwithsignal;//this is only for tests - for real life you want dataormc=data !!!
401	+	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!");
402	+
403	+	TCut limitnJetcut;
404	+	if(JES==noJES) limitnJetcut=cutnJets;
405	+	else {
406	+	if(JES==JESdown) limitnJetcut=cutnJetsJESdown;
407	+	if(JES==JESup) limitnJetcut=cutnJetsJESup;
408	+	}
409	+	TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
410	+	TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
411	+	TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,dataormc,luminosity);
412	+	TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,dataormc,luminosity);
413	+
414	+	TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
415	+	TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
416	+	TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
417	+	TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,dataormc,luminosity);
418	+
419	+	TH1F *LZOSSFP = allsamples.Draw("LZOSSFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
420	+	TH1F *LZOSOFP = allsamples.Draw("LZOSOFP",mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
421	+	TH1F *LZOSSFN = allsamples.Draw("LZOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSSF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
422	+	TH1F *LZOSOFN = allsamples.Draw("LZOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutmass&&cutOSOF&&limitnJetcut&&basiccut,mc,luminosity,allsamples.FindSample("LM4"));
423	+
424	+	TH1F *LSBOSSFP = allsamples.Draw("LSBOSSFP",mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
425	+	TH1F *LSBOSOFP = allsamples.Draw("LSBOSOFP",mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
426	+	TH1F *LSBOSSFN = allsamples.Draw("LSBOSSFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSSF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
427	+	TH1F *LSBOSOFN = allsamples.Draw("LSBOSOFN","-"+mcjzb,binning, "JZB4limits", "events",cutOSOF&&limitnJetcut&&basiccut&&sidebandcut,mc,luminosity,allsamples.FindSample("LM4"));
428	+
429	+	string obsname="data_obs";
430	+	string predname="background";
431	+	string signalname="signal";
432	+	if(identifier!="") {
433	+	obsname=("data_"+identifier);
434	+	predname=("background_"+identifier);
435	+	signalname="signal_"+identifier;
436	+	}
437	+
438	+	TH1F *obs = (TH1F*)ZOSSFP->Clone();
439	+	obs->SetName(obsname.c_str());
440	+	obs->Write();
441	+	TH1F *pred = (TH1F*)ZOSSFN->Clone();
442	+	pred->Add(ZOSOFP,1.0/3);
443	+	pred->Add(ZOSOFN,-1.0/3);
444	+	pred->Add(SBOSSFP,1.0/3);
445	+	pred->Add(SBOSSFN,-1.0/3);
446	+	pred->Add(SBOSOFP,1.0/3);
447	+	pred->Add(SBOSOFN,-1.0/3);
448	+	pred->SetName(predname.c_str());
449	+	pred->Write();
450	+
451	+	//  TH1F *Lobs = (TH1F*)LZOSSFP->Clone();
452	+	//  TH1F *Lpred = (TH1F*)LZOSSFN->Clone();
453	+
454	+	TH1F *Lobs = new TH1F("Lobs","Lobs",binning.size()-1,&binning[0]);
455	+	TH1F *Lpred = new TH1F("Lpred","Lpred",binning.size()-1,&binning[0]);
456	+	Lobs->Add(LZOSSFP);
457	+	Lpred->Add(LZOSSFN);
458	+	Lpred->Add(LZOSOFP,1.0/3);
459	+	Lpred->Add(LZOSOFN,-1.0/3);
460	+	Lpred->Add(LSBOSSFP,1.0/3);
461	+	Lpred->Add(LSBOSSFN,-1.0/3);
462	+	Lpred->Add(LSBOSOFP,1.0/3);
463	+	Lpred->Add(LSBOSOFN,-1.0/3);
464	+	TH1F *signal = (TH1F*)Lobs->Clone();
465	+	signal->Add(Lpred,-1);
466	+	signal->SetName(signalname.c_str());
467	+	signal->Write();
468	+
469	+	delete Lobs;
470	+	delete Lpred;
471	+
472	+	delete ZOSSFP;
473	+	delete ZOSOFP;
474	+	delete ZOSSFN;
475	+	delete ZOSOFN;
476	+
477	+	delete SBOSSFP;
478	+	delete SBOSOFP;
479	+	delete SBOSSFN;
480	+	delete SBOSOFN;
481	+
482	+	delete LZOSSFP;
483	+	delete LZOSOFP;
484	+	delete LZOSSFN;
485	+	delete LZOSOFN;
486	+
487	+	delete LSBOSSFP;
488	+	delete LSBOSOFP;
489	+	delete LSBOSSFN;
490	+	delete LSBOSOFN;
491	+
492	+	}
493	+
494	+	void prepare_datacard(TFile *f) {
495	+	TH1F *dataob = (TH1F*)f->Get("data_obs");
496	+	TH1F *signal = (TH1F*)f->Get("signal");
497	+	TH1F *background = (TH1F*)f->Get("background");
498	+
499	+	ofstream datacard;
500	+	ensure_directory_exists(get_directory()+"/limits");
501	+	datacard.open ((get_directory()+"/limits/susylm4datacard.txt").c_str());
502	+	datacard << "Writing this to a file.\n";
503	+	datacard << "imax 1\n";
504	+	datacard << "jmax 1\n";
505	+	datacard << "kmax *\n";
506	+	datacard << "---------------\n";
507	+	datacard << "shapes * * limitfile.root $PROCESS $PROCESS_$SYSTEMATIC\n";
508	+	datacard << "---------------\n";
509	+	datacard << "bin 1\n";
510	+	datacard << "observation "<<dataob->Integral()<<"\n";
511	+	datacard << "------------------------------\n";
512	+	datacard << "bin             1          1\n";
513	+	datacard << "process         signal     background\n";
514	+	datacard << "process         0          1\n";
515	+	datacard << "rate            "<<signal->Integral()<<"         "<<background->Integral()<<"\n";
516	+	datacard << "--------------------------------\n";
517	+	datacard << "lumi     lnN    1.10       1.0\n";
518	+	datacard << "bgnorm   lnN    1.00       1.4  uncertainty on our prediction (40%)\n";
519	+	datacard << "JES    shape    1          1    uncertainty on background shape and normalization\n";
520	+	datacard << "peak   shape    1          1    uncertainty on signal resolution. Assume the histogram is a 2 sigma shift, \n";
521	+	datacard << "#                                so divide the unit gaussian by 2 before doing the interpolation\n";
522	+	datacard.close();
523	+	}
524	+
525	+
526	+	void prepare_limits(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbbins) {
527	+	ensure_directory_exists(get_directory()+"/limits");
528	+	TFile *limfile = new TFile((get_directory()+"/limits/limitfile.root").c_str(),"RECREATE");
529	+	TCanvas *limcan = new TCanvas("limcan","Canvas for calculating limits");
530	+	limit_shapes_for_systematic_effect(limfile,"",mcjzb,datajzb,noJES,jzbbins,limcan);
531	+	limit_shapes_for_systematic_effect(limfile,"peakUp",newjzbexpression(mcjzb,jzbpeakerrormc),newjzbexpression(datajzb,jzbpeakerrordata),noJES,jzbbins,limcan);
532	+	limit_shapes_for_systematic_effect(limfile,"peakDown",newjzbexpression(mcjzb,-jzbpeakerrormc),newjzbexpression(datajzb,-jzbpeakerrordata),noJES,jzbbins,limcan);
533	+	limit_shapes_for_systematic_effect(limfile,"JESUp",mcjzb,datajzb,JESup,jzbbins,limcan);
534	+	limit_shapes_for_systematic_effect(limfile,"JESDown",mcjzb,datajzb,JESdown,jzbbins,limcan);
535	+
536	+	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!");
539	+	limfile->Close();
540	+
541	+	}

Diff Legend

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