ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C

Comparing UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C (file contents):
Revision 1.22 by buchmann, Fri Jun 14 15:52:03 2013 UTC vs.
Revision 1.26 by buchmann, Fri Jun 21 12:44:46 2013 UTC

#	Line 32 | Line 32
32		#include "EdgeModules/RooSUSYTPdf.cxx"
33		#include "EdgeModules/RooSUSYBkgPdf.cxx"
34
35	+	#include "md5/md5.h"
36	+	#include "md5/md5.C"
37
38		using namespace std;
39		using namespace PlottingSetup;
#	Line 76 | Line 78 | namespace EdgeFitter {
78		TCut cut;
79
80		RooDataSet* AllData;
81	<	RooDataSet* SFSample;
81	>	RooDataSet* eeSample;
82	>	RooDataSet* mmSample;
83		RooDataSet* OFSample;
84
85		bool MarcoDebug=true;
#	Line 92 | Line 95 | namespace EdgeFitter {
95		bool AllowTriangle=true;
96		}
97
95	–	TLatex* WriteAppleStudel() {
96	–	string sel="Apple Strudel Preliminary";
97	–	TLatex *sele = new TLatex(0.97,0.135,sel.c_str());
98	–	sele->SetNDC(true);
99	–	sele->SetTextColor(TColor::GetColor("#848484"));
100	–	sele->SetTextFont(42);
101	–	sele->SetTextAlign(32);
102	–	sele->SetTextSize(0.03);
103	–	sele->SetTextAngle(270);
104	–	return sele;
105	–	}
106	–
98		TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
99		float massLM[1];
100		massLM[0]=mass;
#	Line 283 | Line 274 | float EdgeFitter::calcExclusion(RooWorks
274
275		}
276
277	+	string GetSkimName(int isample) {
278	+	return removefunnystring(allsamples.collection[isample].filename);
279	+	}
280	+
281		TTree* SkimTree(int isample) {
282	<	string TreeName = removefunnystring(allsamples.collection[isample].filename);
282	>	string TreeName = GetSkimName(isample);
283		cout << "About to skim " << TreeName << " for sample id " << isample << endl;
284	<	TTree* newTree = new TTree("NanoTree",TreeName.c_str());
284	>	TTree* newTree = new TTree(TreeName.c_str(),TreeName.c_str());
285
286		float mll,edgeWeight;
287		int id1,id2;
#	Line 339 | Line 334 | void EdgeFitter::InitializeVariables(flo
334		}
335
336		TTree* MergeTrees(vector<TTree*> trees) {
337	+	assert(trees.size()>0);
338		TTree * newtree = (TTree*)trees[0]->CloneTree(0);
339		newtree->SetTitle("FullTree");
340		newtree->SetName("FullTree");
#	Line 368 | Line 364 | TTree* MergeTrees(vector<TTree*> trees)
364
365		void EdgeFitter::PrepareDatasets(int is_data) {
366		write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
367	<	TFile *tempout = new TFile("tempout.root","RECREATE");
367	>
368	>	ensure_directory_exists(PlottingSetup::cbafbasedir+"/EdgeCache/");
369	>
370	>	stringstream FileName;
371	>	FileName << PlottingSetup::cbafbasedir << "/EdgeCache/" << md5((const char*)cut) << ".root";
372	>
373	>
374	>	TFile *fEdgeCache;
375	>	if(PlottingSetup::do_CleanCache) fEdgeCache = new TFile(FileName.str().c_str(),"RECREATE");
376	>	else fEdgeCache = new TFile(FileName.str().c_str(),"UPDATE");
377	>
378	>	if(fEdgeCache->GetNkeys()==0) {
379	>	ofstream CacheOverview;
380	>	CacheOverview.open((PlottingSetup::cbafbasedir+"/EdgeCache/CacheOverview").c_str(),ios::app);
381	>	CacheOverview << md5((const char*)cut) << ";" << (const char*) cut << endl;
382	>	CacheOverview.close();
383	>	}
384	>
385		vector<TTree*> SkimmedTrees;
386		TTree *SkimmedTree[(int)allsamples.collection.size()];
387		for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
#	Line 378 | Line 391 | void EdgeFitter::PrepareDatasets(int is_
391		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
392		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
393		if(EdgeFitter::MarcoDebug) dout << "Considering : " << allsamples.collection[isample].samplename << endl;
394	<	SkimmedTree[isample] = SkimTree(isample);
395	<	tempout->cd();
396	<	SkimmedTree[isample]->Write();
394	>
395	>	string SkimName = GetSkimName(isample);
396	>	SkimmedTree[isample] = (TTree*)fEdgeCache->Get(SkimName.c_str());
397	>	if(!SkimmedTree[isample]) {
398	>	dout << "Need to generate tree for " << allsamples.collection[isample].filename << endl;
399	>	SkimmedTree[isample] = SkimTree(isample);
400	>	fEdgeCache->cd();
401	>	SkimmedTree[isample]->Write();
402	>	} else {
403	>	dout << "Loaded tree for " << allsamples.collection[isample].filename << " from cache file " << endl;
404	>	}
405	>
406		SkimmedTrees.push_back(SkimmedTree[isample]);
407		}
408
#	Line 397 | Line 419 | void EdgeFitter::PrepareDatasets(int is_
419		string title="CMS Data";
420		if(is_data!=1) title="CMS SIMULATION";
421		RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
422	<	tempout->cd();
423	<	completetree->Write();
424	<	tempout->Close();
425	<
426	<	EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
422	>	fEdgeCache->Close();
423	>	dout << "Edge cache closed." << endl;
424	>
425	>	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==0");
426	>	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==1");
427		EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
428		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
429
430	<	SFSample->SetName("SFSample");
430	>	eeSample->SetName("eeSample");
431	>	mmSample->SetName("mmSample");
432		OFSample->SetName("OFSample");
433		AllData->SetName("AllData");
434
435		if(EdgeFitter::MarcoDebug) {
436		dout << "Number of (weighted) events in full sample = " << AllData->sumEntries() << " (raw events : " << AllData->numEntries() << ")" <<  endl;
437	<	dout << "Number of (weighted) events in eemm sample = " << SFSample->sumEntries() << " (raw events : " << SFSample->numEntries() << ")" <<  endl;
437	>	dout << "Number of (weighted) events in ee sample = " << eeSample->sumEntries() << " (raw events : " << eeSample->numEntries() << ")" <<  endl;
438	>	dout << "Number of (weighted) events in mm sample = " << mmSample->sumEntries() << " (raw events : " << mmSample->numEntries() << ")" <<  endl;
439		dout << "Number of (weighted) events in em sample = " << OFSample->sumEntries() << " (raw events : " << OFSample->numEntries() << ")" <<  endl;
440		}
441
442		}
443
444		void EdgeFitter::DrawDatasetContent(int is_data) {
445	<	TCanvas* cSFdata = new TCanvas("cSFdata","cSFdata") ;
445	>	TCanvas* ceedata = new TCanvas("ceedata","ceedata") ;
446
447		RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
448	<	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
448	>	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("ee sample")) ;
449		frame1->GetXaxis()->CenterTitle(1);
450		frame1->GetYaxis()->CenterTitle(1);
451	<	SFSample->plotOn(frame1,RooFit::Name("SFdata")) ;
451	>	eeSample->plotOn(frame1,RooFit::Name("eedata")) ;
452
453	<	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
453	>	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("mm sample")) ;
454		frame2->GetXaxis()->CenterTitle(1);
455		frame2->GetYaxis()->CenterTitle(1);
456	<	OFSample->plotOn(frame2,RooFit::Name("OFdata")) ;
456	>	mmSample->plotOn(frame2,RooFit::Name("mmdata")) ;
457
458	<	//  THStack sf_ref = allsamples.DrawStack("sf_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2"),is_data, luminosity);
459	<	//  THStack of_ref = allsamples.DrawStack("of_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1!=id2"),is_data, luminosity);
458	>	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
459	>	frame3->GetXaxis()->CenterTitle(1);
460	>	frame3->GetYaxis()->CenterTitle(1);
461	>	OFSample->plotOn(frame3,RooFit::Name("OFdata")) ;
462
463	<	TH1F *sf_ref = allsamples.Draw("sf_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2"),is_data, luminosity);
463	>	TH1F *ee_ref = allsamples.Draw("ee_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2&&id1==0"),is_data, luminosity);
464	>	TH1F *mm_ref = allsamples.Draw("mm_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1==id2&&id1==1"),is_data, luminosity);
465		TH1F *of_ref = allsamples.Draw("of_ref","mll",int((mllmax-mllmin)/5.0),mllmin,mllmax, "m_{ll} [GeV]", "events", cut&&TCut("id1!=id2"),is_data, luminosity);
466
467	<	TH1F *sf_fit = (TH1F*)SFSample->createHistogram("sf_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
467	>	ee_ref->SetFillColor(TColor::GetColor("#3ADF00"));
468	>	mm_ref->SetFillColor(TColor::GetColor("#3ADF00"));
469	>	of_ref->SetFillColor(TColor::GetColor("#3ADF00"));
470	>
471	>	TH1F *ee_fit = (TH1F*)eeSample->createHistogram("ee_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
472	>	TH1F *mm_fit = (TH1F*)mmSample->createHistogram("mm_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
473		TH1F *of_fit = (TH1F*)OFSample->createHistogram("of_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
474
475	<	TLegend *leg = allsamples.allbglegend();
476	<	leg->AddEntry(sf_fit,"RooDataSet content","p");
477	<	leg->SetX1(0.58);
475	>	TLegend *leg = make_legend();
476	>	leg->AddEntry(ee_fit,"RooDataSet content","p");
477	>	leg->AddEntry(ee_ref,"CBAF cross-check","f");
478	>	leg->SetX1(0.5);
479	>	leg->SetY1(0.75);
480
481	<	cSFdata->cd() ;
481	>	ceedata->cd() ;
482		gPad->SetLeftMargin(0.15);
483	+	frame1->SetMaximum(1.2*frame1->GetMaximum());
484		frame1->GetYaxis()->SetTitleOffset(1.4);
485		frame1->Draw();
486	<	sf_ref->Draw("histo,same");
487	<	sf_fit->Draw("e1,same");
486	>	ee_ref->Draw("histo,same");
487	>	ee_fit->Draw("e1,same");
488		leg->Draw("same");
489		if(is_data==data) DrawPrelim();
490		else DrawPrelim(PlottingSetup::luminosity,true);
491	<	CompleteSave(cSFdata,"Edge/SF_NoFit");
491	>	CompleteSave(ceedata,"Edge/PreFit_ee");
492
493	<	TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
459	<	cOFdata->cd() ;
493	>	ceedata->cd() ;
494		gPad->SetLeftMargin(0.15);
495	<	frame2->SetMaximum(frame1->GetMaximum());
495	>	frame2->SetMaximum(1.2*frame2->GetMaximum());
496		frame2->GetYaxis()->SetTitleOffset(1.4);
497		frame2->Draw();
498	+	mm_ref->Draw("histo,same");
499	+	mm_fit->Draw("e1,same");
500	+	leg->Draw("same");
501	+	if(is_data==data) DrawPrelim();
502	+	else DrawPrelim(PlottingSetup::luminosity,true);
503	+	CompleteSave(ceedata,"Edge/PreFit_mm");
504	+
505	+	TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
506	+	cOFdata->cd() ;
507	+	gPad->SetLeftMargin(0.15);
508	+	frame3->SetMaximum(1.2*frame3->GetMaximum());
509	+	frame3->GetYaxis()->SetTitleOffset(1.4);
510	+	frame3->Draw();
511		of_ref->Draw("histo,same");
512		of_fit->Draw("e1,same");
513		leg->Draw("same");
514		if(is_data==data) DrawPrelim();
515		else DrawPrelim(PlottingSetup::luminosity,true);
516	<	CompleteSave(cOFdata,"Edge/OF_NoFit");
516	>	CompleteSave(cOFdata,"Edge/PreFit_OF");
517
518	<	delete sf_fit;
518	>	delete ee_fit;
519	>	delete mm_fit;
520		delete of_fit;
521
522	<	delete sf_ref;
522	>	delete ee_ref;
523	>	delete mm_ref;
524		delete of_ref;
476	–
525		}
526
527		string WriteWithError(float central, float error, int digits) {
#	Line 527 | Line 575 | void EdgeFitter::DoFit(int is_data, floa
575		RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
576		RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
577		RooCategory sample("sample","sample") ;
578	<	sample.defineType("SF");
579	<	//sample.defineType("mm");
578	>	sample.defineType("ee");
579	>	sample.defineType("mm");
580		sample.defineType("OF");
581
582	<	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
583	<	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
582	>	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("mm",*mmSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
583	>
584	>	//------------------------------------------------------------------------------------------------------------------------------------
585	>	//    _____ _               __
586	>	//   / ____| |             /_ |
587	>	//  | (___ | |_ ___ _ __    | |
588	>	//   \___ \| __/ _ \ '_ \   | |
589	>	//   ____) | ||  __/ |_) |  | |
590	>	//  |_____/ \__\___| .__/   |_|
591	>	//                 | |
592	>	//                 |_|
593	>	//
594	>	// Fit OF to get good initial values for flavor-symmetric background
595	>
596
597
598		//First we make a fit to opposite flavor
599	<	RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 1.5*combData.sumEntries());
600	<	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
599	>	RooRealVar fttbarOF("fttbarOF", "fttbarOF", OFSample->sumEntries(), 0, combData.sumEntries());
600	>	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.77, 0.01, 4.0);
601		RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
602	<	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
602	>	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.0272, 0.001, 1.0);
603		RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
604		RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
605		RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
#	Line 568 | Line 628 | void EdgeFitter::DoFit(int is_data, floa
628		if(EdgeFitter::FixedMEdge>=0) CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str()+"__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
629		else CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
630		delete pof;
631	+
632	+	TCanvas* pof2 = new TCanvas("pof2","pof2") ;
633	+	pof2->cd();
634	+	gPad->SetLeftMargin(0.15);
635	+	frameO->GetYaxis()->SetTitleOffset(1.4);
636	+	frameO->Draw();
637	+	if(is_data==data) DrawPrelim();
638	+	else DrawPrelim(PlottingSetup::luminosity,true);
639	+	stringstream OFFitParameterSummary;
640	+	OFFitParameterSummary << "#splitline{ftbbarOF = " <<  WriteWithError(fttbarOF.getVal(),fttbarOF.getError(),3)  << "}{";
641	+	OFFitParameterSummary << "#splitline{par1ttbarOF = " <<  WriteWithError(par1ttbarOF.getVal(),par1ttbarOF.getError(),3)  << "}{";
642	+	OFFitParameterSummary << "#splitline{par2ttbarOF = " <<  WriteWithError(par2ttbarOF.getVal(),par2ttbarOF.getError(),3)  << "}{";
643	+	OFFitParameterSummary << "#splitline{par3ttbarOF = " <<  WriteWithError(par3ttbarOF.getVal(),par3ttbarOF.getError(),3)  << "}{";
644	+	OFFitParameterSummary << "par4ttbarOF = " <<  WriteWithError(par4ttbarOF.getVal(),par4ttbarOF.getError(),3);
645	+	OFFitParameterSummary << "}}}}";
646	+	TLatex *OFFitinfobox = new TLatex(0.57,0.75,OFFitParameterSummary.str().c_str());
647	+	OFFitinfobox->SetNDC();
648	+	OFFitinfobox->SetTextSize(0.03);
649	+	OFFitinfobox->Draw();
650	+
651	+
652	+	if(EdgeFitter::FixedMEdge>=0) CompleteSave(pof2,"Edge/OF__OFFitonly_"+prefix.str()+"__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge)+"__INFO",false,false);
653	+	else CompleteSave(pof2,"Edge/OF__OFFitonly_"+prefix.str()+"__INFO",false,false);
654	+
655	+	delete pof2;
656
657		if(resultOF->covQual()!=3) {
658		write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
#	Line 579 | Line 664 | void EdgeFitter::DoFit(int is_data, floa
664		write_info(__FUNCTION__,"OF fit converged");
665		}
666
667	+
668	+	//    _____ _               ___
669	+	//   / ____| |             |__ \
670	+	//  | (___ | |_ ___ _ __      ) |
671	+	//   \___ \| __/ _ \ '_ \    / /
672	+	//   ____) | ||  __/ |_) |  / /_
673	+	//  |_____/ \__\___| .__/  |____|
674	+	//                 | |
675	+	//                 |_|
676	+	//
677	+	// This is where the actual simultaneous fitting happens
678	+
679	+
680	+	// Step 2a: set up edge position (if fixed), and maximum Z yield
681		float StartingMedge=70;
682		if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
584	–
683
684		RooDataSet *ZDataSet = (RooDataSet*)EdgeFitter::AllData->reduce("id1==id2 && abs(mll-91.2)<20");
685
#	Line 590 | Line 688 | void EdgeFitter::DoFit(int is_data, floa
688		delete ZDataSet;
689
690
691	<	// Now same flavor
594	<	RooRealVar fzSF("fzSF", "fzSF", 39, 39, maxZ);
595	<	RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
596	<	//RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
597	<	RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
598	<	RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
599	<	widthzSF.setConstant(1);
600	<
601	<	/*RooRealVar fttbarSFx("fttbarSFx","fttbarSFx",0.2,1.8);
602	<	RooRealVar mRsfof("mRsfof","mRsfof",1.02);
603	<	RooRealVar wRsfof("wRsfof","wRsfof",0.07);
604	<
605	<	RooGaussian RSFOF("RSFOF","RSFOF",fttbarSFx,mRsfof,wRsfof);
606	<	RooProdPdf fttbarSF("fttbarSF","fttbarOF*RSFOF",RooArgList(RSFOF,fttbarOF));*/
607	<	RooRealVar fttbarSF("fttbarSF", "fttbarSF", fttbarOF.getVal(), 0.2*fttbarOF.getVal(), 1.5*fttbarOF.getVal());
691	>	// Step 2b: set up the models!
692
693	<	RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", 1.02*par1ttbarOF.getVal(), (1.02-0.07)*par1ttbarOF.getVal(), (1.02+0.07)*par1ttbarOF.getVal());
694	<
695	<	RooRealVar fsignalSF("fsignalSF", "fsignalSF", 0, 0, 300);
696	<	RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
697	<	RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
698	<	RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
699	<
700	<	RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
701	<
702	<	if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
693	>	RooRealVar eefrac("eefrac","eefrac",eeSample->sumEntries()/(eeSample->sumEntries()+mmSample->sumEntries()),0.1,1.0);
694	>	RooFormulaVar mmfrac("mmfrac","1-eefrac",RooArgSet(eefrac));
695	>
696	>	RooRealVar rOneOverEMFrac("rOneOverEMFrac","rOneOverEMFrac",1.02,1.02-0.07,1.02+0.07);
697	>	RooRealVar rmean("rmean","rmean",1.02);
698	>	RooRealVar rwidth("rwidth","rwidth",0.07);
699	>	RooGaussian OneOverEMFrac("OneOverEMFrac","OneOverEMFrac",rOneOverEMFrac,rmean,rwidth);
700	>
701	>
702	>	// FIRST: Drell-Yan: different parameters for resolution of DY (and of course yield), but same mean and so on
703	>	RooRealVar widthz("widthz", "widthz", 2.94);
704	>	widthz.setConstant(1);
705	>	RooRealVar meanz("meanz", "meanz", 91.1876, 89, 93);
706	>	RooRealVar fzSF("fzSF","fzSF",39,0,maxZ);
707	>	RooFormulaVar fzee("fzee","eefrac*fzSF",RooArgSet(eefrac,fzSF));
708	>	RooFormulaVar fzmm("fzmm","mmfrac*fzSF",RooArgSet(mmfrac,fzSF));
709	>	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0.5, 5);
710	>	RooRealVar sigmazmm("sigmazmm", "sigmazmm", 5, 0.5, 5);
711	>
712	>	RooVoigtian zee("zee", "zee", mll, meanz, widthz, sigmazee);
713	>	RooVoigtian zmm("zmm", "zmm", mll, meanz, widthz, sigmazmm);
714	>	RooAddPdf zSF("zSF","zSF",RooArgList(zee,zmm));
715
620	–	RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
621	–	RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
622	–
623	–	RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
716
717	+	//SECOND: Flavor-Symmetry contribution. Only relative fraction is different, all other parameters are the same;
718	+	RooFormulaVar fttbaree("fttbaree","eefrac*fttbarOF*OneOverEMFrac",RooArgSet(eefrac,fttbarOF,OneOverEMFrac));
719	+	RooFormulaVar fttbarmm("fttbarmm","mmfrac*fttbarOF*OneOverEMFrac",RooArgSet(mmfrac,fttbarOF,OneOverEMFrac));
720	+
721	+	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
722	+	RooSUSYBkgPdf ttbarmm("ttbarmm","ttbarmm", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
723	+	RooAddPdf ttbarSF("ttbarSF","ttbarSF", RooArgList(ttbaree,ttbarmm), RooArgList(fttbaree,fttbarmm));
724	+
725	+
726	+	//for signal only the resolution and absolute fraction are different for ee / mm
727	+	RooRealVar fsignal("fsignal", "fsignal", 0, 0, 300);
728	+	RooFormulaVar fsignalee("fsignalee","eefrac*fsignal",RooArgSet(eefrac,fsignal));
729	+	RooFormulaVar fsignalmm("fsignalmm","mmfrac*fsignal",RooArgSet(mmfrac,fsignal));
730	+
731	+	RooRealVar par1signal("par1signal", "par1signal", 45, 20, 100);
732	+	RooRealVar par3signal("par3signal", "par3signal", StartingMedge, 0, 300);
733	+	if(EdgeFitter::FixedMEdge>0) par3signal.setConstant();
734	+
735	+	RooSUSYTPdf signalee("signalee","signalee", mll , par1signal, sigmazee, par3signal);
736	+	RooSUSYTPdf signalmm("signalmm","signalmm", mll , par1signal, sigmazmm, par3signal);
737	+	RooAddPdf signalSF("signalSF","signalSF",RooArgList(signalee,signalmm),RooArgList(fsignalee,fsignalmm));
738
739		if(!EdgeFitter::AllowTriangle) {
740	<	fsignalSF.setVal(0.0); // kill off the signal if we don't want the triangle
741	<	fsignalSF.setConstant(1);
742	<	par1signalSF.setConstant(1);
743	<	par2signalSF.setConstant(1);
631	<	par3signalSF.setConstant(1);
740	>	fsignal.setVal(0.0); // kill off the signal if we don't want the triangle
741	>	fsignal.setConstant(1);
742	>	par1signal.setConstant(1);
743	>	par3signal.setConstant(1);
744		}
745
746	+
747	+	//COMPLETE MODEL:
748	+
749	+	RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
750	+	RooAddPdf model_mm("model_mm","model_mm", RooArgList(zmm, ttbarmm, signalmm), RooArgList(fzmm, fttbarmm, fsignalmm));
751	+	//  RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
752	+	RooAddPdf model_SF("model_SF","model_SF",RooArgList(model_ee,model_mm));
753	+
754		RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
755	<	simPdf.addPdf(model_SF,"SF") ;
755	>	simPdf.addPdf(model_ee,"ee") ;
756	>	simPdf.addPdf(model_mm,"mm") ;
757		simPdf.addPdf(model_OF,"OF") ;
758
638	–
759		RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
760
761		if(result->covQual()!=3) {
#	Line 652 | Line 772 | void EdgeFitter::DoFit(int is_data, floa
772		result->Print();
773		dout << "============================= < /Full results > =============================" << endl;
774
775	+
776
777		RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
778		frame1->GetXaxis()->CenterTitle(1);
779		frame1->GetYaxis()->CenterTitle(1);
780	<	combData.plotOn(frame1,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::SF")) ;
781	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
782	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
783	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("DYSFonly"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
784	<	if(EdgeFitter::AllowTriangle) simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("SignalSFonly"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
785	<
786	<	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame1->chiSquare("FullFit", "SFdata", 3);
787	<	else EdgeFitter::FixedMEdgeChi2_H0=frame1->chiSquare("FullFit", "SFdata", 3);
788	<
789	<	dout << "Result   : " << endl;
790	<	if(EdgeFitter::AllowTriangle) dout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
791	<	//  dout << "f ttbar  : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
792	<	dout << "f ttbar  : " << fttbarSF.getVal() << " +/- NO ERROR CUZ ITS A PDF "<< endl;
793	<	dout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
794	<	dout << "f z SF   : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
795	<	if(EdgeFitter::AllowTriangle) dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
675	<	else dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
780	>	combData.plotOn(frame1,RooFit::Name("eedata"),RooFit::Cut("sample==sample::ee")) ;
781	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("FullFitEE"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
782	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("TTbarEEonly"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
783	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("DYEEonly"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
784	>	if(EdgeFitter::AllowTriangle) simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("SignalEEonly"),RooFit::Components("signalee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
785	>
786	>
787	>	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("mm sample")) ;
788	>	frame2->GetXaxis()->CenterTitle(1);
789	>	frame2->GetYaxis()->CenterTitle(1);
790	>	combData.plotOn(frame2,RooFit::Name("mmdata"),RooFit::Cut("sample==sample::mm")) ;
791	>	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("FullFitMM"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
792	>	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("TTbarMMonly"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
793	>	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("DYMMonly"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
794	>	if(EdgeFitter::AllowTriangle) simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("SignalMMonly"),RooFit::Components("signalmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
795	>
796
677	–	// The same plot for the cointrol sample slice
797		RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
798		frame3->GetXaxis()->CenterTitle(1);
799		frame3->GetYaxis()->CenterTitle(1);
800	<	frame3->SetMaximum(frame1->GetMaximum());
801	<	combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
802	<	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
803	<	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
804	<
800	>	combData.plotOn(frame3,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
801	>	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("FullFitOF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
802	>	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
803	>	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("DYOFonly"),RooFit::Components("zOF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
804	>
805	>
806	>	RooPlot* frame4 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
807	>	frame4->GetXaxis()->CenterTitle(1);
808	>	frame4->GetYaxis()->CenterTitle(1);
809	>	combData.plotOn(frame4,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::ee||sample==sample::mm")) ;
810	>	simPdf.plotOn(frame4,RooFit::Slice(sample,"ee"),RooFit::Name("FullFitSFee"),RooFit::ProjWData(sample,combData),RooFit::LineColor(kRed),RooFit::LineStyle(kDashed));
811	>	simPdf.plotOn(frame4,RooFit::Slice(sample,"ee"),RooFit::Slice(sample,"mm"),RooFit::Name("FullFitSF"),RooFit::ProjWData(sample,combData),RooFit::LineColor(kRed));
812	>	simPdf.plotOn(frame4,RooArgSet(RooFit::Slice(sample,"ee"),RooFit::Slice(sample,"mm")),RooFit::Name("FullFitSF"),RooFit::ProjWData(sample,combData),RooFit::LineColor(kGreen),RooFit::LineStyle(kDotted)));
813	>	simPdf.plotOn(frame4,RooFit::Slice(sample,"mm"),RooFit::Name("FullFitSFmm"),RooFit::ProjWData(sample,combData),RooFit::LineColor(kBlue),RooFit::LineStyle(kDotted));
814	>	simPdf.plotOn(frame4,RooFit::Components(RooArgSet(ttbaree,ttbarmm)),RooFit::LineStyle(kDotted));
815	>	simPdf.plotOn(frame4,RooFit::Components(RooArgSet(zee,zmm)),RooFit::LineStyle(kDotted));
816	>	//  model_SF.plotOn(frame4);
817	>	// ttbarSF.plotOn(frame4,RooFit::LineStyle(kDashed));
818	>	// zSF.plotOn(frame4,RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
819	>	if(EdgeFitter::AllowTriangle) signalSF.plotOn(frame4,RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
820	>
821
822	+	/// ********************************************************************************************************************************
823		TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
824		c->cd() ;
825		gPad->SetLeftMargin(0.15);
#	Line 691 | Line 827 | void EdgeFitter::DoFit(int is_data, floa
827		frame1->Draw();
828		if(is_data==data) DrawPrelim();
829		else DrawPrelim(PlottingSetup::luminosity,true);
830	<	stringstream infotext;
831	<	infotext << "#splitline{Fit results (" << EdgeFitter::Mode << ">" << jzb_cut << "): }{#splitline{";
832	<	infotext << "N(Data) = " << EdgeFitter::SFSample->sumEntries() << "}{#splitline{";
833	<	infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
834	<	//infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
835	<	infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),0,3) << "}{#splitline{";
836	<	write_warning(any2string(__LINE__),"Don't have the error yet, need to complete this");
837	<	if(EdgeFitter::AllowTriangle) {
838	<	infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
839	<	infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
840	<	} else infotext << "}{}}}}}";
705	<
706	<	TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
707	<	infobox->SetNDC();
708	<	infobox->SetTextSize(0.03);
709	<	infobox->Draw();
710	<	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
711	<	else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
712	<	delete c;
830	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_ee__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
831	>	else CompleteSave(c,"Edge/"+prefix.str()+"_ee",false,false);
832	>
833	>	c->cd() ;
834	>	gPad->SetLeftMargin(0.15);
835	>	frame2->GetYaxis()->SetTitleOffset(1.4);
836	>	frame2->Draw();
837	>	if(is_data==data) DrawPrelim();
838	>	else DrawPrelim(PlottingSetup::luminosity,true);
839	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_mm__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
840	>	else CompleteSave(c,"Edge/"+prefix.str()+"_mm",false,false);
841
842	<	TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
715	<	e->cd();
842	>	c->cd() ;
843		gPad->SetLeftMargin(0.15);
844		frame3->GetYaxis()->SetTitleOffset(1.4);
845		frame3->Draw();
846		if(is_data==data) DrawPrelim();
847		else DrawPrelim(PlottingSetup::luminosity,true);
848	<	if(EdgeFitter::FixedMEdge>=0) CompleteSave(e,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
849	<	else CompleteSave(e,"Edge/"+prefix.str()+"_OF",false,false);
850	<	delete e;
851	<
725	<
726	<
727	<
728	<	/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
729	<	f->cd();
848	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
849	>	else CompleteSave(c,"Edge/"+prefix.str()+"_OF",false,false);
850	>
851	>	c->cd() ;
852		gPad->SetLeftMargin(0.15);
853		frame4->GetYaxis()->SetTitleOffset(1.4);
854		frame4->Draw();
855		if(is_data==data) DrawPrelim();
856		else DrawPrelim(PlottingSetup::luminosity,true);
857	<	CompleteSave(f,"Edge/"+prefix.str()+"_SF");
858	<	delete f;*/
859	<
857	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
858	>	else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
859	>
860	>	delete c;
861	>	/// ********************************************************************************************************************************
862
739	–	/*
740	–	float maxZ=200;
741	–	RooWorkspace* wspace = new RooWorkspace();
742	–	stringstream mllvar;
743	–	mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
744	–	wspace->factory(mllvar.str().c_str());
745	–	wspace->var("mll")->setBins(30);
746	–	wspace->factory("nSig[1.,0.,100.]");
747	–	wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
748	–	wspace->factory("rME[1.12,1.05,1.19]");
749	–	wspace->factory("effUncert[1.]");
750	–	EdgeFitter::prepareLimits(wspace, true);
751	–	*/
863
864	<	write_warning(__FUNCTION__," A lot missing here to calculate limits");
864	>	//  if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame2->chiSquare("FullFit", "EEdata", 3);
865	>	//  else EdgeFitter::FixedMEdgeChi2_H0=frame2->chiSquare("FullFit", "EEdata", 3);
866	>
867	>
868
869		}
870
#	Line 770 | Line 884 | void EdgeFitter::DoEdgeFit(string mcjzb,
884		RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
885
886
773	–
887		EdgeFitter::AllowTriangle=false;
888		EdgeFitter::DoFit(is_data, jzb_cut);
889

Diff Legend

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