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.5 by buchmann, Wed Jun 27 08:53:23 2012 UTC vs.
Revision 1.23 by buchmann, Thu Jun 20 13:16:43 2013 UTC

#	Line 1 | Line 1
1		#include <iostream>
2
3	+
4	+
5	+	#include <TVirtualIndex.h>
6	+
7		#include <RooRealVar.h>
8		#include <RooArgSet.h>
9		#include <RooDataSet.h>
#	Line 7 | Line 11
11		#include <RooCategory.h>
12
13		#include <RooPlot.h>
14	+	#include <RooGaussian.h>
15	+	#include <RooConstVar.h>
16		#include <RooSimultaneous.h>
17		#include <RooAddPdf.h>
18		#include <RooFitResult.h>
#	Line 23 | Line 29
29		#include "RooStats/HypoTestInverterOriginal.h"
30
31		//#include "ParametrizedEdge.C"
32	<	#include "/shome/pablom/RooFit/Pdfs/RooSUSYTPdf.cxx"
33	<	#include "/shome/pablom/RooFit/Pdfs/RooSUSYBkgPdf.cxx"
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 51 | Line 59 | namespace EdgeFitter {
59
60		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree*);
61		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree*);
62	<	void getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
62	>	void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
63		void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
64		void PrepareDatasets(int);
65	+	void DrawDatasetContent(int);
66		void DoFit(int is_data, float jzb_cut);
67		string RandomStorageFile();
68		Yield Get_Z_estimate(float,int);
69		Yield Get_T_estimate(float,int);
70	<	float calcExclusion(RooWorkspace *ws, RooDataSet *data = NULL);
71	<	void prepareLimits(RooWorkspace *ws);
63	<	vector<RooDataSet*> generateToys(RooWorkspace *ws, int nToys);
70	>	float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
71	>	vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
72		void prepareLimits(RooWorkspace *ws, bool ComputeBands);
73		TGraph* prepareLM(float mass, float nEv);
74
#	Line 70 | Line 78 | namespace EdgeFitter {
78		TCut cut;
79
80		RooDataSet* AllData;
81	<	RooDataSet* eeSample;
82	<	RooDataSet* mmSample;
83	<	RooDataSet* emSample;
81	>	RooDataSet* SFSample;
82	>	RooDataSet* OFSample;
83	>
84	>	bool MarcoDebug=true;
85	>
86	>	float FixedMEdge=-1;
87	>	float FixedMEdgeChi2_H1=-1;
88	>	float FixedMEdgeChi2_H0=-1;
89
90	<	bool MarcoDebug;
90	>	bool RejectPointIfNoConvergence=false;
91	>
92	>	string Mode="UndefinedMode";
93	>
94	>	bool AllowTriangle=true;
95		}
96
97		TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
#	Line 93 | Line 110 | TGraph* EdgeFitter::prepareLM(float mass
110		return lm;
111		}
112
113	<	vector<RooDataSet*> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
113	>	vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
114	>	ws->ls();
115		ws->var("nSig")->setVal(0.);
116		ws->var("nSig")->setConstant(true);
117		RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
118	<	vector<RooDataSet*> theToys;
118	>	vector<RooDataSet> theToys;
119
120		RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
121		mcEE.generate(nToys,44,true);
#	Line 106 | Line 124 | vector<RooDataSet*> EdgeFitter::generate
124		RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
125		mcOSOF.generate(nToys,44,true);
126
127	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
127	>	RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
128		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
129		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
130		RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
#	Line 120 | Line 138 | vector<RooDataSet*> EdgeFitter::generate
138		stringstream toyname;
139		toyname << "theToy_" << i;
140		write_warning(__FUNCTION__,"Problem while adding toys");
141	<	//    RooDataSet *toyData = RooDataSet(toyname.str(),toyname.str(),observables,RooFit::Index(ws->cat("cat")),RooFit::Import("OSOF",*toyOSOF),RooFit::Import("EE",*toyEE),RooFit::Import("MM",*toyMM));
142	<	//    theToys.push_back(toyData);
141	>	RooDataSet toyData = RooDataSet(toyname.str().c_str(),toyname.str().c_str(),observables,RooFit::Index(const_cast<RooCategory&>(*ws->cat("cat"))),RooFit::Import("OSOF",*toyOSOF),RooFit::Import("EE",*toyEE),RooFit::Import("MM",*toyMM));
142	>	theToys.push_back(toyData);
143		}
144		ws->var("nSig")->setVal(17.0);
145		ws->var("nSig")->setConstant(false);
146		return theToys;
147		}
148
149	<	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
149	>	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
150		TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
151		vector<float> theLimits;
152	<	for(int itoy=0;itoy<theToys.size();itoy++) {
153	<	float theLimit = calcExclusion(ws,theToys[itoy]);
152	>	for(int itoy=0;itoy<(int)theToys.size();itoy++) {
153	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
154		if(theLimit > 0 ) gauLimit->Fill(theLimit);
155		}
156		const Int_t nQ = 4;
#	Line 145 | Line 163 | void EdgeFitter::getMedianLimit(RooWorks
163		sigmaDown = abs(yQ[1]-median);
164		sigmaUp = abs(yQ[2]-median);
165		twoSigmaUp = abs(yQ[3]-median);
166	<	cout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
166	>	dout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
167		}
168
169		void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
170		if(ComputeBands) {
171	<	vector<RooDataSet*> theToys = EdgeFitter::generateToys(ws,50);
171	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
172		vector<float> medVals;
173		vector<float> medLimits;
174		vector<float> sigmaLimitsDown;
#	Line 179 | Line 197 | void EdgeFitter::prepareLimits(RooWorksp
197		theVals.push_back((float)i);
198		ws->var("m0")->setVal((float)i);
199		ws->var("m0")->setConstant(true);
200	<	theLimits.push_back(calcExclusion(ws));
200	>	//      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
201	>	write_error(__FUNCTION__,"Error while casting roo data set");
202		}
203
204	<	for(int i=0;i<theLimits.size();i++) {
204	>	for(int i=0;i<(int)theLimits.size();i++) {
205		if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
206	<	cout << i << " : " << theVals[i] << endl;
206	>	dout << i << " : " << theVals[i] << endl;
207		theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
208		write_warning(__FUNCTION__,"Need to store limits");
209		}
#	Line 194 | Line 213 | void EdgeFitter::prepareLimits(RooWorksp
213		}
214
215
216	<	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet *data) {
216	>	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
217	>	int numberOfToys=50;
218		RooRealVar mu("mu","nSig",0,10000,"");
219		RooArgSet poi = RooArgSet(mu);
220		RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
#	Line 203 | Line 223 | float EdgeFitter::calcExclusion(RooWorks
223		model->SetWorkspace(*ws);
224		model->SetPdf("combModel");
225		model->SetParametersOfInterest(poi);
226	<	if(!data) data = (RooDataSet*)ws->data("data_obs");
226	>	//  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
227
228	<	RooStats::ProfileLikelihoodCalculator plc(*data, *model);
228	>	RooStats::ProfileLikelihoodCalculator plc(data, *model);
229		plc.SetNullParameters(*nullParams);
230		plc.SetTestSize(0.05);
231	+
232		RooStats::LikelihoodInterval* interval = plc.GetInterval();
233		RooStats::HypoTestResult *htr = plc.GetHypoTest();
234		double theLimit = interval->UpperLimit( mu );
235	<	cout << "Significance " << htr->Significance() << endl;
235	>	//  double significance = htr->Significance();
236
237		ws->defineSet("obs","nB");
238		ws->defineSet("poi","nSig");
#	Line 235 | Line 256 | float EdgeFitter::calcExclusion(RooWorks
256		slrts.SetAltParameters(*sb_model.GetSnapshot());
257		RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
258
259	<	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
259	>	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
260		hc.SetTestStatistic(2);
261	<	hc.SetNumberOfToys(50);
261	>	hc.SetNumberOfToys(numberOfToys);
262
263		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
264		hcInv.SetTestSize(0.05);
#	Line 246 | Line 267 | float EdgeFitter::calcExclusion(RooWorks
267		RooStats::HypoTestInverterResult* hcInt = hcInv.GetInterval();
268		float ulError = hcInt->UpperLimitEstimatedError();
269		theLimit = hcInt->UpperLimit();
270	<	cout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
270	>	dout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
271
272		return theLimit/PlottingSetup::luminosity;
273
274		}
275
276	+	string GetSkimName(int isample) {
277	+	return removefunnystring(allsamples.collection[isample].filename);
278	+	}
279	+
280		TTree* SkimTree(int isample) {
281	<	TTree* newTree = allsamples.collection[isample].events->CloneTree(0);
281	>	string TreeName = GetSkimName(isample);
282	>	cout << "About to skim " << TreeName << " for sample id " << isample << endl;
283	>	TTree* newTree = new TTree(TreeName.c_str(),TreeName.c_str());
284	>
285	>	float mll,edgeWeight;
286	>	int id1,id2;
287	>
288	>	newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
289	>	newTree->Branch("mll",&mll,"mll/F");
290	>	newTree->Branch("id1",&id1,"id1/I");
291	>	newTree->Branch("id2",&id2,"id2/I");
292	>	newTree->Branch("isample",&isample,"isample/I");
293	>
294		float xsweight=1.0;
295		if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
296		if(EdgeFitter::MarcoDebug) {
297	<	cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
298	<	cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
297	>	dout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
298	>	dout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
299		}
300	+
301	+	float tmll;
302	+	int tid1,tid2;
303	+	allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
304	+	allsamples.collection[isample].events->SetBranchAddress("id1",&tid1);
305	+	allsamples.collection[isample].events->SetBranchAddress("id2",&tid2);
306	+
307		TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
308	+	TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
309	+
310		float wgt=1.0;
265	–	allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
311		for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
312		allsamples.collection[isample].events->LoadTree(entry);
313		if (select->EvalInstance()) {
314		allsamples.collection[isample].events->GetEntry(entry);
315	<	wgt=wgt*xsweight;
315	>	mll=tmll;
316	>	id1=tid1;
317	>	id2=tid2;
318	>	wgt=Weight->EvalInstance();
319	>	edgeWeight=wgt*xsweight;
320		newTree->Fill();
321		}
322		}
323
324	<	if(EdgeFitter::MarcoDebug) cout << "     Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
324	>	if(EdgeFitter::MarcoDebug) dout << "   Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
325		return newTree;
326		}
327
#	Line 282 | Line 331 | void EdgeFitter::InitializeVariables(flo
331		jzbmax=_jzbmax;
332		cut=_cut;
333		}
334	<
334	>
335	>	TTree* MergeTrees(vector<TTree*> trees) {
336	>	assert(trees.size()>0);
337	>	TTree * newtree = (TTree*)trees[0]->CloneTree(0);
338	>	newtree->SetTitle("FullTree");
339	>	newtree->SetName("FullTree");
340	>	trees[0]->GetListOfClones()->Remove(newtree);
341	>	trees[0]->ResetBranchAddresses();
342	>	newtree->ResetBranchAddresses();
343	>
344	>	for(int itree=0;itree<trees.size();itree++) {
345	>	newtree->CopyAddresses(trees[itree]);
346	>	Long64_t nentries = trees[itree]->GetEntries();
347	>	for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
348	>	trees[itree]->GetEntry(iEntry);
349	>	newtree->Fill();
350	>	}
351	>	trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
352	>	if (newtree->GetTreeIndex()) {
353	>	newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
354	>	}
355	>	if (newtree && newtree->GetTreeIndex()) {
356	>	newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
357	>	}
358	>	}
359	>	return newtree;
360	>	}
361	>
362	>
363	>
364		void EdgeFitter::PrepareDatasets(int is_data) {
287	–	TTree *completetree;
365		write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
366	<	bool hashit=0;
367	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
366	>
367	>	ensure_directory_exists(PlottingSetup::cbafbasedir+"/EdgeCache/");
368	>
369	>	stringstream FileName;
370	>	FileName << PlottingSetup::cbafbasedir << "/EdgeCache/" << md5((const char*)cut) << ".root";
371	>
372	>	TFile *fEdgeCache = new TFile(FileName.str().c_str(),"UPDATE");
373	>	if(fEdgeCache->GetNkeys()==0) {
374	>	ofstream CacheOverview;
375	>	CacheOverview.open((PlottingSetup::cbafbasedir+"/EdgeCache/CacheOverview").c_str(),ios::app);
376	>	CacheOverview << md5((const char*)cut) << ";" << (const char*) cut << endl;
377	>	CacheOverview.close();
378	>	}
379	>
380	>	vector<TTree*> SkimmedTrees;
381	>	TTree *SkimmedTree[(int)allsamples.collection.size()];
382	>	for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
383		if(!allsamples.collection[isample].is_active) continue;
384		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
385		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
386		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
387		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
388	<	if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
389	<	if(!hashit) {
390	<	hashit=true;
391	<	completetree = SkimTree(isample)->CloneTree();
388	>	if(EdgeFitter::MarcoDebug) dout << "Considering : " << allsamples.collection[isample].samplename << endl;
389	>
390	>	string SkimName = GetSkimName(isample);
391	>	SkimmedTree[isample] = (TTree*)fEdgeCache->Get(SkimName.c_str());
392	>	if(!SkimmedTree[isample]) {
393	>	dout << "Need to generate tree for " << allsamples.collection[isample].filename << endl;
394	>	SkimmedTree[isample] = SkimTree(isample);
395	>	fEdgeCache->cd();
396	>	SkimmedTree[isample]->Write();
397		} else {
398	<	completetree->CopyEntries(SkimTree(isample));
398	>	dout << "Loaded tree for " << allsamples.collection[isample].filename << " from cache file " << endl;
399		}
400	<	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
400	>
401	>	SkimmedTrees.push_back(SkimmedTree[isample]);
402		}
403
404	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
404	>	TTree *completetree  = MergeTrees(SkimmedTrees);
405	>
406	>	if(EdgeFitter::MarcoDebug) dout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
407	>
408	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
409		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
410		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
411	<	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
412	<	RooRealVar weight("weight","weight",0,1000,"");
311	<	RooArgSet observables(mll,jzb,id1,id2,weight);
411	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
412	>	RooArgSet observables(mll,id1,id2,edgeWeight);
413
414		string title="CMS Data";
415		if(is_data!=1) title="CMS SIMULATION";
416	<	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
417	<	completetree->Write();
418	<	//  delete completetree;
419	<
420	<	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2");
421	<	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2");
321	<	EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
416	>	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
417	>	fEdgeCache->Close();
418	>	dout << "Edge cache closed." << endl;
419	>
420	>	EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
421	>	EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
422		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
423
424	<	eeSample->SetName("eeSample");
425	<	mmSample->SetName("mmSample");
326	<	emSample->SetName("emSample");
424	>	SFSample->SetName("SFSample");
425	>	OFSample->SetName("OFSample");
426		AllData->SetName("AllData");
427
428		if(EdgeFitter::MarcoDebug) {
429	<	cout << "Number of events in data sample = " << AllData->numEntries() << endl;
430	<	cout << "Number of events in ee sample = " << eeSample->numEntries() << endl;
431	<	cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
333	<	cout << "Number of events in em sample = " << emSample->numEntries() << endl;
429	>	dout << "Number of (weighted) events in full sample = " << AllData->sumEntries() << " (raw events : " << AllData->numEntries() << ")" <<  endl;
430	>	dout << "Number of (weighted) events in eemm sample = " << SFSample->sumEntries() << " (raw events : " << SFSample->numEntries() << ")" <<  endl;
431	>	dout << "Number of (weighted) events in em sample = " << OFSample->sumEntries() << " (raw events : " << OFSample->numEntries() << ")" <<  endl;
432		}
433	+
434	+	}
435	+
436	+	void EdgeFitter::DrawDatasetContent(int is_data) {
437	+	TCanvas* cSFdata = new TCanvas("cSFdata","cSFdata") ;
438	+
439	+	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
440	+	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
441	+	frame1->GetXaxis()->CenterTitle(1);
442	+	frame1->GetYaxis()->CenterTitle(1);
443	+	SFSample->plotOn(frame1,RooFit::Name("SFdata")) ;
444	+
445	+	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
446	+	frame2->GetXaxis()->CenterTitle(1);
447	+	frame2->GetYaxis()->CenterTitle(1);
448	+	OFSample->plotOn(frame2,RooFit::Name("OFdata")) ;
449	+
450	+	//  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);
451	+	//  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);
452	+
453	+	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);
454	+	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);
455	+
456	+	TH1F *sf_fit = (TH1F*)SFSample->createHistogram("sf_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
457	+	TH1F *of_fit = (TH1F*)OFSample->createHistogram("of_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
458	+
459	+	TLegend *leg = allsamples.allbglegend();
460	+	leg->AddEntry(sf_fit,"RooDataSet content","p");
461	+	leg->SetX1(0.58);
462	+
463	+	cSFdata->cd() ;
464	+	gPad->SetLeftMargin(0.15);
465	+	frame1->GetYaxis()->SetTitleOffset(1.4);
466	+	frame1->Draw();
467	+	sf_ref->Draw("histo,same");
468	+	sf_fit->Draw("e1,same");
469	+	leg->Draw("same");
470	+	if(is_data==data) DrawPrelim();
471	+	else DrawPrelim(PlottingSetup::luminosity,true);
472	+	CompleteSave(cSFdata,"Edge/SF_NoFit");
473	+
474	+	TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
475	+	cOFdata->cd() ;
476	+	gPad->SetLeftMargin(0.15);
477	+	frame2->SetMaximum(frame1->GetMaximum());
478	+	frame2->GetYaxis()->SetTitleOffset(1.4);
479	+	frame2->Draw();
480	+	of_ref->Draw("histo,same");
481	+	of_fit->Draw("e1,same");
482	+	leg->Draw("same");
483	+	if(is_data==data) DrawPrelim();
484	+	else DrawPrelim(PlottingSetup::luminosity,true);
485	+	CompleteSave(cOFdata,"Edge/OF_NoFit");
486	+
487	+	delete sf_fit;
488	+	delete of_fit;
489	+
490	+	delete sf_ref;
491	+	delete of_ref;
492	+
493	+	}
494	+
495	+	string WriteWithError(float central, float error, int digits) {
496	+	float ref=central;
497	+	if(ref<0) ref=-central;
498	+	int HighestSigDigit = 0;
499	+	if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
500	+	else HighestSigDigit = int(log(ref)/(log(10)));
501	+
502	+	float divider=pow(10.0,(double(HighestSigDigit-digits)));
503	+
504	+	stringstream result;
505	+	result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
506	+	return result.str();
507		}
508
509	+
510		string EdgeFitter::RandomStorageFile() {
511		TRandom3 *r = new TRandom3(0);
512		int rho = (int)r->Uniform(1,10000000);
#	Line 356 | Line 529 | Yield EdgeFitter::Get_T_estimate(float j
529		}
530
531		void EdgeFitter::DoFit(int is_data, float jzb_cut) {
532	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
532	>
533	>	stringstream prefix;
534	>	if(is_data==data) prefix << "data_";
535	>	if(is_data==mc) prefix << "mc_";
536	>	if(is_data==mcwithsignal) prefix << "mcwithS_";
537	>
538	>	prefix << EdgeFitter::Mode << "_" << jzb_cut;
539	>
540	>	if(EdgeFitter::AllowTriangle) prefix << "__H1";
541	>	else prefix << "__H0";
542	>
543	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
544	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
545		RooCategory sample("sample","sample") ;
546	<	sample.defineType("ee");
546	>	sample.defineType("SF");
547		//sample.defineType("mm");
548	<	sample.defineType("em");
549	<	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("ee",*eeSample),Import("mm",*mmSample),Import("em",*emSample));
550	<	RooDataSet combData("combData","combined data",mll,RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("em",*emSample));
548	>	sample.defineType("OF");
549	>
550	>	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
551	>	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
552	>
553
367	–
368	–
554		//First we make a fit to opposite flavor
555	<	RooRealVar fttbarem("fttbarem", "fttbarem", 100, 0, 10000);
556	<	RooRealVar par1ttbarem("par1ttbarem", "par1ttbarem", 1.6, 0.01, 4.0);
557	<	RooRealVar par2ttbarem("par2ttbarem", "par2ttbarem", 1.0);
558	<	RooRealVar par3ttbarem("par3ttbarem", "par3ttbarem", 0.028, 0.001, 1.0);
559	<	RooRealVar par4ttbarem("par4ttbarem", "par4ttbarem", 2.0);
560	<	RooSUSYBkgPdf ttbarem("ttbarem","ttbarem", mll , par1ttbarem, par2ttbarem, par3ttbarem, par4ttbarem);
561	<	RooAddPdf model_em("model_em","model_em", ttbarem, fttbarem);
555	>	RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 1.5*combData.sumEntries());
556	>	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
557	>	RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
558	>	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
559	>	RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
560	>	RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
561	>	RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
562		RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
563	<	simPdfOF.addPdf(model_em,"em");
564	<	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save());
563	>	simPdfOF.addPdf(model_OF,"OF");
564	>	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save(),RooFit::Extended(),RooFit::Minos(true));
565	>	dout << "============================= < OF results > =============================" << endl;
566		resultOF->Print();
567	<
568	<	RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarem");
569	<	float resultOFpar1 = resultOFpar1_->getVal();
570	<	//RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarem");
571	<	//float resultOFpar2 = resultOFpar2_->getVal();
572	<	//cout << "caca2.txt" << endl;
573	<
574	<	RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarem");
575	<	float resultOFpar3 = resultOFpar3_->getVal();
576	<
577	<	//RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarem");
578	<	//float resultOFpar4 = resultOFpar4_->getVal();
579	<	//cout << "caca4.txt" << endl;
567	>	dout << "============================= < /OF results > =============================" << endl;
568	>
569	>
570	>	RooPlot* frameO = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample"));
571	>	frameO->GetXaxis()->CenterTitle(1);
572	>	frameO->GetYaxis()->CenterTitle(1);
573	>	combData.plotOn(frameO,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
574	>	simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
575	>	simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
576	>
577	>	TCanvas* pof = new TCanvas("pof","pof") ;
578	>	pof->cd() ;
579	>	gPad->SetLeftMargin(0.15);
580	>	frameO->GetYaxis()->SetTitleOffset(1.4);
581	>	frameO->Draw();
582	>	if(is_data==data) DrawPrelim();
583	>	else DrawPrelim(PlottingSetup::luminosity,true);
584	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str()+"__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
585	>	else CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
586	>	delete pof;
587	>
588	>	if(resultOF->covQual()!=3) {
589	>	write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
590	>	dout << "covQual is " << resultOF->covQual() << endl;
591	>	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
592	>	else EdgeFitter::FixedMEdgeChi2_H0=-1;
593	>	if(EdgeFitter::RejectPointIfNoConvergence) return;
594	>	} else {
595	>	write_info(__FUNCTION__,"OF fit converged");
596	>	}
597	>
598	>	float StartingMedge=70;
599	>	if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
600
601
602	+	RooDataSet *ZDataSet = (RooDataSet*)EdgeFitter::AllData->reduce("id1==id2 && abs(mll-91.2)<20");
603	+
604	+	float maxZ = ZDataSet->sumEntries();
605	+	dout << "maxZ was set to " << maxZ << endl;
606	+	delete ZDataSet;
607	+
608	+
609		// Now same flavor
610	<	RooRealVar fzee("fzee", "fzee", 5, 0, 100000);
611	<	RooRealVar meanzee("meanzee", "meanzee", 91.1876, 89, 95);
612	<	//RooRealVar sigmazee("sigmazee", "sigmazee", 0.5);
613	<	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0, 100);
614	<	RooRealVar widthzee("widthzee", "widthzee", 2.94);
615	<
616	<	RooRealVar fttbaree("fttbaree", "fttbaree", 100, 0, 100000);
617	<	RooRealVar par1ttbaree("par1ttbaree", "par1ttbaree", resultOFpar1, 0, 100);
618	<	RooRealVar par2ttbaree("par2ttbaree", "par2ttbaree", 1.0);
619	<	RooRealVar par3ttbaree("par3ttbaree", "par3ttbaree", resultOFpar3, 0, 100);
620	<	RooRealVar par4ttbaree("par4ttbaree", "par4ttbaree", 2.0);
610	>	RooRealVar fzSF("fzSF", "fzSF", 39, 39, maxZ);
611	>	RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
612	>	//RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
613	>	RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
614	>	RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
615	>	widthzSF.setConstant(1);
616	>
617	>	/*RooRealVar fttbarSFx("fttbarSFx","fttbarSFx",0.2,1.8);
618	>	RooRealVar mRsfof("mRsfof","mRsfof",1.02);
619	>	RooRealVar wRsfof("wRsfof","wRsfof",0.07);
620	>
621	>	RooGaussian RSFOF("RSFOF","RSFOF",fttbarSFx,mRsfof,wRsfof);
622	>	RooProdPdf fttbarSF("fttbarSF","fttbarOF*RSFOF",RooArgList(RSFOF,fttbarOF));*/
623	>	RooRealVar fttbarSF("fttbarSF", "fttbarSF", fttbarOF.getVal(), 0.2*fttbarOF.getVal(), 1.5*fttbarOF.getVal());
624	>
625	>	RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", 1.02*par1ttbarOF.getVal(), (1.02-0.07)*par1ttbarOF.getVal(), (1.02+0.07)*par1ttbarOF.getVal());
626
627	<	RooRealVar fsignalee("fsignalee", "fsignalee", 10, 0, 400);
628	<	RooRealVar par1signalee("par1signalee", "par1signalee", 45, 20, 100);
629	<	RooRealVar par2signalee("par2signalee", "par2signalee", 2, 1, 10);
630	<	RooRealVar par3signalee("par3signalee", "par3signalee", 45, 0, 200);
627	>	RooRealVar fsignalSF("fsignalSF", "fsignalSF", 0, 0, 300);
628	>	RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
629	>	RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
630	>	RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
631
632	<	RooVoigtian zee("zee", "zee", mll, meanzee, widthzee, sigmazee);
632	>	RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
633
634	+	if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
635
636	<	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbaree, par2ttbaree, par3ttbaree, par4ttbaree);
637	<	//RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, par2signalee, par3signalee);
419	<	RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, sigmazee, par3signalee);
420	<
421	<	//RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
422	<	RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
423	<	RooAddPdf model_emu("model_emu","model_emu", RooArgList(ttbaree), RooArgList(fttbaree));
636	>	RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
637	>	RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
638
639	+	RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
640	+
641	+
642	+	if(!EdgeFitter::AllowTriangle) {
643	+	fsignalSF.setVal(0.0); // kill off the signal if we don't want the triangle
644	+	fsignalSF.setConstant(1);
645	+	par1signalSF.setConstant(1);
646	+	par2signalSF.setConstant(1);
647	+	par3signalSF.setConstant(1);
648	+	}
649
650		RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
651	<	simPdf.addPdf(model_ee,"ee") ;
652	<	simPdf.addPdf(model_emu,"em") ;
651	>	simPdf.addPdf(model_SF,"SF") ;
652	>	simPdf.addPdf(model_OF,"OF") ;
653
430	–	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save());
431	–	result->Print();
432	–
433	–	RooPlot* frame1 = mll.frame(RooFit::Bins(25),RooFit::Title("EE sample")) ;
434	–	combData.plotOn(frame1,RooFit::Cut("sample==sample::ee")) ;
435	–	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
436	–	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
437	–	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
438	–	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("signalee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
439	–
440	–	RooPlot* frame2 = mll.frame(RooFit::Bins(25),RooFit::Title("MM sample")) ;
441	–	combData.plotOn(frame2,RooFit::Cut("sample==sample::mm")) ;
442	–	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
443	–	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
444	–	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
445	–	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("signalmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
446	–
447	–
448	–	cout << "Result   : " << endl;
449	–	cout << "f signal : " << fsignalee.getVal() << " +/- " << fsignalee.getError() << endl;
450	–	cout << "f ttbar  : " << fttbaree.getVal() << " +/- " << fttbaree.getError() << endl;
451	–	cout << "f tt em  : " << fttbarem.getVal() << " +/- " << fttbarem.getError() << endl;
452	–	cout << "f z ee   : " << fzee.getVal() << " +/- " << fzee.getError() << endl;
654
655	<	// The same plot for the cointrol sample slice
455	<	RooPlot* frame3 = mll.frame(RooFit::Bins(25),RooFit::Title("EM sample")) ;
456	<	combData.plotOn(frame3,RooFit::Cut("sample==sample::em")) ;
457	<	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
458	<	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::Components("ttbarem"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
655	>	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
656
657	<	stringstream prefix;
658	<	if(is_data==data) prefix << "data_";
659	<	if(is_data==mc) prefix << "mc_";
660	<	if(is_data==mcwithsignal) prefix << "mcwithS_";
657	>	if(result->covQual()!=3) {
658	>	write_error(__FUNCTION__,"Full fit did not converge!!! Cannot continue!");
659	>	dout << "covQual is " << result->covQual() << endl;
660	>	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
661	>	else EdgeFitter::FixedMEdgeChi2_H0=-1;
662	>	if(EdgeFitter::RejectPointIfNoConvergence) return;
663	>	} else {
664	>	write_info(__FUNCTION__,"Full fit converged");
665	>	}
666	>
667	>	dout << "============================= < Full results > =============================" << endl;
668	>	result->Print();
669	>	dout << "============================= < /Full results > =============================" << endl;
670	>
671
672	<	prefix << "JZB_" << jzb_cut;
672	>	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
673	>	frame1->GetXaxis()->CenterTitle(1);
674	>	frame1->GetYaxis()->CenterTitle(1);
675	>	combData.plotOn(frame1,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::SF")) ;
676	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
677	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
678	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("DYSFonly"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
679	>	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));
680	>
681	>	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame1->chiSquare("FullFit", "SFdata", 3);
682	>	else EdgeFitter::FixedMEdgeChi2_H0=frame1->chiSquare("FullFit", "SFdata", 3);
683	>
684	>	dout << "Result   : " << endl;
685	>	if(EdgeFitter::AllowTriangle) dout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
686	>	//  dout << "f ttbar  : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
687	>	dout << "f ttbar  : " << fttbarSF.getVal() << " +/- NO ERROR CUZ ITS A PDF "<< endl;
688	>	dout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
689	>	dout << "f z SF   : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
690	>	if(EdgeFitter::AllowTriangle) dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
691	>	else dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
692
693	<	/*  cout << "fsignalee : " << fsignalee << endl;
694	<	cout << "fttbaree : " << fttbaree << endl;
695	<	cout << "fzee : " << fzee << endl;*/
693	>	// The same plot for the cointrol sample slice
694	>	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
695	>	frame3->GetXaxis()->CenterTitle(1);
696	>	frame3->GetYaxis()->CenterTitle(1);
697	>	frame3->SetMaximum(frame1->GetMaximum());
698	>	combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
699	>	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
700	>	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
701
702
703		TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
#	Line 476 | Line 707 | void EdgeFitter::DoFit(int is_data, floa
707		frame1->Draw();
708		if(is_data==data) DrawPrelim();
709		else DrawPrelim(PlottingSetup::luminosity,true);
710	<	CompleteSave(c,"Edge/"+prefix.str()+"eemm");
710	>	stringstream infotext;
711	>	infotext << "#splitline{Fit results (" << EdgeFitter::Mode << ">" << jzb_cut << "): }{#splitline{";
712	>	infotext << "N(Data) = " << EdgeFitter::SFSample->sumEntries() << "}{#splitline{";
713	>	infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
714	>	//infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
715	>	infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),0,3) << "}{#splitline{";
716	>	write_warning(any2string(__LINE__),"Don't have the error yet, need to complete this");
717	>	if(EdgeFitter::AllowTriangle) {
718	>	infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
719	>	infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
720	>	} else infotext << "}{}}}}}";
721	>
722	>	TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
723	>	infobox->SetNDC();
724	>	infobox->SetTextSize(0.03);
725	>	infobox->Draw();
726	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
727	>	else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
728		delete c;
729
482	–	TCanvas* d = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
483	–	d->cd() ;
484	–	gPad->SetLeftMargin(0.15);
485	–	frame2->GetYaxis()->SetTitleOffset(1.4);
486	–	frame2->Draw();
487	–	if(is_data==data) DrawPrelim();
488	–	else DrawPrelim(PlottingSetup::luminosity,true);
489	–	CompleteSave(d,"Edge/"+prefix.str()+"mm");
490	–	delete d;
491	–	//c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
492	–
730		TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
731		e->cd();
732		gPad->SetLeftMargin(0.15);
#	Line 497 | Line 734 | void EdgeFitter::DoFit(int is_data, floa
734		frame3->Draw();
735		if(is_data==data) DrawPrelim();
736		else DrawPrelim(PlottingSetup::luminosity,true);
737	<	CompleteSave(e,"Edge/"+prefix.str()+"emu");
737	>	if(EdgeFitter::FixedMEdge>=0) CompleteSave(e,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
738	>	else CompleteSave(e,"Edge/"+prefix.str()+"_OF",false,false);
739		delete e;
740
741	+
742	+
743	+
744		/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
745		f->cd();
746		gPad->SetLeftMargin(0.15);
#	Line 507 | Line 748 | void EdgeFitter::DoFit(int is_data, floa
748		frame4->Draw();
749		if(is_data==data) DrawPrelim();
750		else DrawPrelim(PlottingSetup::luminosity,true);
751	<	CompleteSave(f,"Edge/"+prefix.str()+"eemm");
751	>	CompleteSave(f,"Edge/"+prefix.str()+"_SF");
752		delete f;*/
753	<
753	>
754	>
755	>	/*
756	>	float maxZ=200;
757	>	RooWorkspace* wspace = new RooWorkspace();
758	>	stringstream mllvar;
759	>	mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
760	>	wspace->factory(mllvar.str().c_str());
761	>	wspace->var("mll")->setBins(30);
762	>	wspace->factory("nSig[1.,0.,100.]");
763	>	wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
764	>	wspace->factory("rME[1.12,1.05,1.19]");
765	>	wspace->factory("effUncert[1.]");
766	>	EdgeFitter::prepareLimits(wspace, true);
767	>	*/
768	>
769	>	write_warning(__FUNCTION__," A lot missing here to calculate limits");
770	>
771		}
772
773		void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
#	Line 522 | Line 780 | void EdgeFitter::DoEdgeFit(string mcjzb,
780
781		EdgeFitter::PrepareDatasets(is_data);
782
783	+	EdgeFitter::DrawDatasetContent(is_data);
784	+
785		RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
786		RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
787	+
788	+
789	+
790	+	EdgeFitter::AllowTriangle=false;
791		EdgeFitter::DoFit(is_data, jzb_cut);
792	+
793	+	write_info(__FUNCTION__,"TAKING SHORTCUT");return;
794	+
795	+	EdgeFitter::AllowTriangle=true;
796	+
797	+	bool ScanMassRange=false;
798	+	float ScanSteps=5.0;//GeV
799	+
800	+
801	+	if(ScanMassRange) {
802	+	TFile *fscan = new TFile("fscan.root","UPDATE");
803	+	TGraph *gr = new TGraph();
804	+	TGraph *Rgr = new TGraph();
805	+	stringstream GrName;
806	+	GrName << "ScanGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
807	+	stringstream RGrName;
808	+	RGrName << "ScanRatioGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
809	+	gr->SetName(GrName.str().c_str());
810	+	Rgr->SetName(RGrName.str().c_str());
811	+
812	+	int i=0;
813	+	for(float tempMedge=10;tempMedge<=300;tempMedge+=ScanSteps) {
814	+	write_info(__FUNCTION__,"Now testing Medge="+any2string(tempMedge)+" for "+EdgeFitter::Mode+">"+any2string(jzb_cut));
815	+	EdgeFitter::FixedMEdge=tempMedge;
816	+	EdgeFitter::DoFit(is_data, jzb_cut);
817	+	if(EdgeFitter::FixedMEdgeChi2_H1<0) continue;
818	+	gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1);
819	+	Rgr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1/EdgeFitter::FixedMEdgeChi2_H0);
820	+	i++;
821	+	}
822	+
823	+	TCanvas *ScanCan = new TCanvas("ScanCan","ScanCan",500,500);
824	+	gr->GetXaxis()->SetTitle("m_{edge}");
825	+	gr->GetXaxis()->CenterTitle();
826	+	gr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
827	+	gr->GetYaxis()->CenterTitle();
828	+	gr->GetYaxis()->SetTitleOffset(0.95);
829	+	gr->GetXaxis()->SetTitleOffset(0.9);
830	+	gr->SetLineColor(kBlue);
831	+	gr->SetTitle("");
832	+	gr->Draw("AL");
833	+	stringstream ScanCanSave;
834	+	ScanCanSave << "Edge/MEdgeScan_"+EdgeFitter::Mode+"_" << jzb_cut;
835	+	if(is_data) DrawPrelim();
836	+	else DrawMCPrelim();
837	+	CompleteSave(ScanCan,ScanCanSave.str());
838	+
839	+	Rgr->GetXaxis()->SetTitle("m_{edge}");
840	+	Rgr->GetXaxis()->CenterTitle();
841	+	Rgr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
842	+	Rgr->GetYaxis()->CenterTitle();
843	+	Rgr->GetYaxis()->SetTitleOffset(0.95);
844	+	Rgr->GetXaxis()->SetTitleOffset(0.9);
845	+	Rgr->SetLineColor(kBlue);
846	+	Rgr->SetTitle("");
847	+	Rgr->Draw("AL");
848	+	ScanCanSave.str("");
849	+	ScanCanSave << "Edge/MEdgeScan_Ratio_"+EdgeFitter::Mode+"_" << jzb_cut;
850	+	if(is_data) DrawPrelim();
851	+	else DrawMCPrelim();
852	+	CompleteSave(ScanCan,ScanCanSave.str());
853	+	fscan->cd();
854	+	gr->Write();
855	+	delete ScanCan;
856	+	fscan->Close();
857	+	} else {
858	+	EdgeFitter::DoFit(is_data, jzb_cut);
859	+	dout << "Chi^2 (H0) = " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
860	+	dout << "Chi^2 (H1) = " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
861	+	}
862	+
863	+
864		RooMsgService::instance().setGlobalKillBelow(msglevel);
865
530	–
866		f->Close();
867
868		}
869
870		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
871	<	for(int icut=0;icut<jzb_cut.size();icut++) {
871	>
872	>	EdgeFitter::Mode="JZB";
873	>	if(mcjzb=="met[4]") EdgeFitter::Mode="MET";
874	>
875	>	for(int icut=0;icut<(int)jzb_cut.size();icut++) {
876		stringstream addcut;
877		if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
878		if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";

Diff Legend

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