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.3 by buchmann, Mon Jun 18 07:51:38 2012 UTC vs.
Revision 1.24 by buchmann, Thu Jun 20 19:38:37 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>
10	<	#include <RooRealVar.h>
11	<	#include <RooArgSet.h>
12	<	#include <RooDataSet.h>
10	>	#include <RooMCStudy.h>
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>
19	>	#include <RooVoigtian.h>
20	>	#include <RooMsgService.h>
21	>
22		#include <RooStats/ModelConfig.h>
23		#include "RooStats/ProfileLikelihoodCalculator.h"
24		#include "RooStats/LikelihoodInterval.h"
#	Line 15 | Line 28
28		#include "RooStats/HybridCalculatorOriginal.h"
29		#include "RooStats/HypoTestInverterOriginal.h"
30
31	+	//#include "ParametrizedEdge.C"
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 40 | 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);
63		void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
64		void PrepareDatasets(int);
65	<	void DoFit();
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 SetEdgeLimit(int , RooWorkspace *ws);
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
75		float jzbmax;
76		float mllmin;
#	Line 56 | Line 80 | namespace EdgeFitter {
80		RooDataSet* AllData;
81		RooDataSet* eeSample;
82		RooDataSet* mmSample;
83	<	RooDataSet* emSample;
83	>	RooDataSet* OFSample;
84
85	<	bool MarcoDebug;
85	>	bool MarcoDebug=true;
86	>
87	>	float FixedMEdge=-1;
88	>	float FixedMEdgeChi2_H1=-1;
89	>	float FixedMEdgeChi2_H0=-1;
90	>
91	>	bool RejectPointIfNoConvergence=false;
92	>
93	>	string Mode="UndefinedMode";
94	>
95	>	bool AllowTriangle=true;
96	>	}
97	>
98	>	TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
99	>	float massLM[1];
100	>	massLM[0]=mass;
101	>	float accEffLM[1];
102	>	accEffLM[0]=nEv/PlottingSetup::luminosity;
103	>	TGraph *lm = new TGraph(1, massLM, accEffLM);
104	>	lm->GetXaxis()->SetNoExponent(true);
105	>	lm->GetXaxis()->SetTitle("m_{cut} [GeV]");
106	>	lm->GetYaxis()->SetTitle("#sigma #times A [pb] 95% CL UL");
107	>	lm->GetXaxis()->SetLimits(0.,300.);
108	>	lm->GetYaxis()->SetRangeUser(0.,0.08);
109	>	lm->SetMarkerStyle(34);
110	>	lm->SetMarkerColor(kRed);
111	>	return lm;
112	>	}
113	>
114	>	vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
115	>	ws->ls();
116	>	ws->var("nSig")->setVal(0.);
117	>	ws->var("nSig")->setConstant(true);
118	>	RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
119	>	vector<RooDataSet> theToys;
120	>
121	>	RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
122	>	mcEE.generate(nToys,44,true);
123	>	RooMCStudy mcMM(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"MM"));
124	>	mcMM.generate(nToys,44,true);
125	>	RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
126	>	mcOSOF.generate(nToys,44,true);
127	>
128	>	RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
129	>	RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
130	>	RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
131	>	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
132	>	RooRealVar weight("weight","weight",0,1000,"");
133	>	RooArgSet observables(mll,jzb,id1,id2,weight);
134	>
135	>	for(int i=0;i<nToys;i++) {
136	>	RooDataSet* toyEE    = (RooDataSet*)mcEE.genData(i);
137	>	RooDataSet* toyMM    = (RooDataSet*)mcMM.genData(i);
138	>	RooDataSet* toyOSOF  = (RooDataSet*)mcOSOF.genData(i);
139	>	stringstream toyname;
140	>	toyname << "theToy_" << i;
141	>	write_warning(__FUNCTION__,"Problem while adding toys");
142	>	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));
143	>	theToys.push_back(toyData);
144	>	}
145	>	ws->var("nSig")->setVal(17.0);
146	>	ws->var("nSig")->setConstant(false);
147	>	return theToys;
148		}
149
150	<	float EdgeFitter::SetEdgeLimit(int is_data, RooWorkspace *ws) {
150	>	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
151	>	TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
152	>	vector<float> theLimits;
153	>	for(int itoy=0;itoy<(int)theToys.size();itoy++) {
154	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
155	>	if(theLimit > 0 ) gauLimit->Fill(theLimit);
156	>	}
157	>	const Int_t nQ = 4;
158	>	Double_t yQ[nQ] = {0.,0.,0.,0.};
159	>	Double_t xQ[nQ] = {0.022750132,0.1586552539,0.84134474609999998,0.977249868};
160	>	gauLimit->GetQuantiles(nQ,yQ,xQ);
161	>	median = gauLimit->GetMean();
162	>	//  median = median1(gauLimit);
163	>	twoSigmaDown = abs(yQ[0]-median);
164	>	sigmaDown = abs(yQ[1]-median);
165	>	sigmaUp = abs(yQ[2]-median);
166	>	twoSigmaUp = abs(yQ[3]-median);
167	>	dout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
168	>	}
169	>
170	>	void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
171	>	if(ComputeBands) {
172	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
173	>	vector<float> medVals;
174	>	vector<float> medLimits;
175	>	vector<float> sigmaLimitsDown;
176	>	vector<float> sigmaLimitsUp;
177	>	vector<float> twoSigmaLimitsDown;
178	>	vector<float> twoSigmaLimitsUp;
179	>	for(int i=20;i<=320;i+=40) {
180	>	ws->var("nSig")->setVal(10.0);
181	>	medVals.push_back((float)i);
182	>	ws->var("m0")->setVal((float)i);
183	>	ws->var("m0")->setConstant(true);
184	>	float Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp;
185	>	EdgeFitter::getMedianLimit(ws,theToys,Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp);
186	>	medLimits.push_back(Smedian);
187	>	sigmaLimitsDown.push_back(SsigmaDown);
188	>	sigmaLimitsUp.push_back(SsigmaUp);
189	>	twoSigmaLimitsDown.push_back(StwoSigmaDown);
190	>	twoSigmaLimitsUp.push_back(StwoSigmaUp);
191	>	}
192	>	write_warning(__FUNCTION__,"Still need to store limits");
193	>	} else {
194	>	vector<float> theVals;
195	>	vector<float> theLimits;
196	>	for(int i=20;i<300;i+=5) {
197	>	ws->var("nSig")->setVal(0.0);
198	>	theVals.push_back((float)i);
199	>	ws->var("m0")->setVal((float)i);
200	>	ws->var("m0")->setConstant(true);
201	>	//      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
202	>	write_error(__FUNCTION__,"Error while casting roo data set");
203	>	}
204	>
205	>	for(int i=0;i<(int)theLimits.size();i++) {
206	>	if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
207	>	dout << i << " : " << theVals[i] << endl;
208	>	theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
209	>	write_warning(__FUNCTION__,"Need to store limits");
210	>	}
211	>	write_warning(__FUNCTION__,"Need to store limits");
212	>	}
213	>	}
214	>	}
215	>
216	>
217	>	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
218	>	int numberOfToys=50;
219		RooRealVar mu("mu","nSig",0,10000,"");
220		RooArgSet poi = RooArgSet(mu);
221		RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
#	Line 70 | Line 224 | float EdgeFitter::SetEdgeLimit(int is_da
224		model->SetWorkspace(*ws);
225		model->SetPdf("combModel");
226		model->SetParametersOfInterest(poi);
227	<	RooAbsData* data = ws->data("data_obs");
227	>	//  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
228
229	<	RooStats::ProfileLikelihoodCalculator plc(*data, *model);
229	>	RooStats::ProfileLikelihoodCalculator plc(data, *model);
230		plc.SetNullParameters(*nullParams);
231		plc.SetTestSize(0.05);
232	+
233		RooStats::LikelihoodInterval* interval = plc.GetInterval();
234		RooStats::HypoTestResult *htr = plc.GetHypoTest();
235		double theLimit = interval->UpperLimit( mu );
236	<	cout << "Significance " << htr->Significance() << endl;
236	>	//  double significance = htr->Significance();
237
238		ws->defineSet("obs","nB");
239		ws->defineSet("poi","nSig");
#	Line 102 | Line 257 | RooAbsData* data = ws->data("data_obs");
257		slrts.SetAltParameters(*sb_model.GetSnapshot());
258		RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
259
260	<	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
260	>	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
261		hc.SetTestStatistic(2);
262	<	hc.SetNumberOfToys(50);
108	<
262	>	hc.SetNumberOfToys(numberOfToys);
263
264		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
265		hcInv.SetTestSize(0.05);
#	Line 114 | Line 268 | RooAbsData* data = ws->data("data_obs");
268		RooStats::HypoTestInverterResult* hcInt = hcInv.GetInterval();
269		float ulError = hcInt->UpperLimitEstimatedError();
270		theLimit = hcInt->UpperLimit();
271	<	cout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
271	>	dout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
272
273		return theLimit/PlottingSetup::luminosity;
274
275		}
276
277	+	string GetSkimName(int isample) {
278	+	return removefunnystring(allsamples.collection[isample].filename);
279	+	}
280	+
281		TTree* SkimTree(int isample) {
282	<	TTree* newTree = allsamples.collection[isample].events->CloneTree(0);
282	>	string TreeName = GetSkimName(isample);
283	>	cout << "About to skim " << TreeName << " for sample id " << isample << endl;
284	>	TTree* newTree = new TTree(TreeName.c_str(),TreeName.c_str());
285	>
286	>	float mll,edgeWeight;
287	>	int id1,id2;
288	>
289	>	newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
290	>	newTree->Branch("mll",&mll,"mll/F");
291	>	newTree->Branch("id1",&id1,"id1/I");
292	>	newTree->Branch("id2",&id2,"id2/I");
293	>	newTree->Branch("isample",&isample,"isample/I");
294	>
295		float xsweight=1.0;
296		if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
297		if(EdgeFitter::MarcoDebug) {
298	<	cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
299	<	cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
298	>	dout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
299	>	dout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
300		}
301	+
302	+	float tmll;
303	+	int tid1,tid2;
304	+	allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
305	+	allsamples.collection[isample].events->SetBranchAddress("id1",&tid1);
306	+	allsamples.collection[isample].events->SetBranchAddress("id2",&tid2);
307	+
308		TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
309	+	TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
310	+
311		float wgt=1.0;
133	–	allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
312		for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
313		allsamples.collection[isample].events->LoadTree(entry);
314		if (select->EvalInstance()) {
315		allsamples.collection[isample].events->GetEntry(entry);
316	<	wgt=wgt*xsweight;
316	>	mll=tmll;
317	>	id1=tid1;
318	>	id2=tid2;
319	>	wgt=Weight->EvalInstance();
320	>	edgeWeight=wgt*xsweight;
321		newTree->Fill();
322		}
323		}
324
325	<	if(EdgeFitter::MarcoDebug) cout << "     Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
325	>	if(EdgeFitter::MarcoDebug) dout << "   Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
326		return newTree;
327		}
328
#	Line 150 | Line 332 | void EdgeFitter::InitializeVariables(flo
332		jzbmax=_jzbmax;
333		cut=_cut;
334		}
335	<
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");
341	>	trees[0]->GetListOfClones()->Remove(newtree);
342	>	trees[0]->ResetBranchAddresses();
343	>	newtree->ResetBranchAddresses();
344	>
345	>	for(int itree=0;itree<trees.size();itree++) {
346	>	newtree->CopyAddresses(trees[itree]);
347	>	Long64_t nentries = trees[itree]->GetEntries();
348	>	for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
349	>	trees[itree]->GetEntry(iEntry);
350	>	newtree->Fill();
351	>	}
352	>	trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
353	>	if (newtree->GetTreeIndex()) {
354	>	newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
355	>	}
356	>	if (newtree && newtree->GetTreeIndex()) {
357	>	newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
358	>	}
359	>	}
360	>	return newtree;
361	>	}
362	>
363	>
364	>
365		void EdgeFitter::PrepareDatasets(int is_data) {
366	<	TTree *completetree;
367	<	bool hashit=0;
368	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
366	>	write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
367	>
368	>	ensure_directory_exists(PlottingSetup::cbafbasedir+"/EdgeCache/");
369	>
370	>	stringstream FileName;
371	>	FileName << PlottingSetup::cbafbasedir << "/EdgeCache/" << md5((const char*)cut) << ".root";
372	>
373	>	TFile *fEdgeCache = new TFile(FileName.str().c_str(),"UPDATE");
374	>	if(fEdgeCache->GetNkeys()==0) {
375	>	ofstream CacheOverview;
376	>	CacheOverview.open((PlottingSetup::cbafbasedir+"/EdgeCache/CacheOverview").c_str(),ios::app);
377	>	CacheOverview << md5((const char*)cut) << ";" << (const char*) cut << endl;
378	>	CacheOverview.close();
379	>	}
380	>
381	>	vector<TTree*> SkimmedTrees;
382	>	TTree *SkimmedTree[(int)allsamples.collection.size()];
383	>	for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
384		if(!allsamples.collection[isample].is_active) continue;
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==false) continue;//kick all samples that aren't data if we're looking for data.
387		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
388		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
389	<	if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
390	<	if(!hashit) {
391	<	hashit=true;
392	<	completetree = SkimTree(isample)->CloneTree();
389	>	if(EdgeFitter::MarcoDebug) dout << "Considering : " << allsamples.collection[isample].samplename << endl;
390	>
391	>	string SkimName = GetSkimName(isample);
392	>	SkimmedTree[isample] = (TTree*)fEdgeCache->Get(SkimName.c_str());
393	>	if(!SkimmedTree[isample]) {
394	>	dout << "Need to generate tree for " << allsamples.collection[isample].filename << endl;
395	>	SkimmedTree[isample] = SkimTree(isample);
396	>	fEdgeCache->cd();
397	>	SkimmedTree[isample]->Write();
398		} else {
399	<	completetree->CopyEntries(SkimTree(isample));
399	>	dout << "Loaded tree for " << allsamples.collection[isample].filename << " from cache file " << endl;
400		}
401	<	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
401	>
402	>	SkimmedTrees.push_back(SkimmedTree[isample]);
403		}
404
405	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
405	>	TTree *completetree  = MergeTrees(SkimmedTrees);
406	>
407	>	if(EdgeFitter::MarcoDebug) dout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
408	>
409	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
410		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
411		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
412	<	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
413	<	RooRealVar weight("weight","weight",0,1000,"");
178	<	RooArgSet observables(mll,jzb,id1,id2,weight);
412	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
413	>	RooArgSet observables(mll,id1,id2,edgeWeight);
414
415		string title="CMS Data";
416		if(is_data!=1) title="CMS SIMULATION";
417	<	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
418	<	completetree->Write();
419	<	//  delete completetree;
420	<
417	>	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
418	>	fEdgeCache->Close();
419	>	dout << "Edge cache closed." << endl;
420	>
421		EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==0");
422		EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==1");
423	<	EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
423	>	EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
424		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
425
426		eeSample->SetName("eeSample");
427		mmSample->SetName("mmSample");
428	<	emSample->SetName("emSample");
428	>	OFSample->SetName("OFSample");
429		AllData->SetName("AllData");
430
431		if(EdgeFitter::MarcoDebug) {
432	<	cout << "Number of events in data sample = " << AllData->numEntries() << endl;
433	<	cout << "Number of events in ee sample = " << eeSample->numEntries() << endl;
434	<	cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
435	<	cout << "Number of events in em sample = " << emSample->numEntries() << endl;
432	>	dout << "Number of (weighted) events in full sample = " << AllData->sumEntries() << " (raw events : " << AllData->numEntries() << ")" <<  endl;
433	>	dout << "Number of (weighted) events in ee sample = " << eeSample->sumEntries() << " (raw events : " << eeSample->numEntries() << ")" <<  endl;
434	>	dout << "Number of (weighted) events in mm sample = " << mmSample->sumEntries() << " (raw events : " << mmSample->numEntries() << ")" <<  endl;
435	>	dout << "Number of (weighted) events in em sample = " << OFSample->sumEntries() << " (raw events : " << OFSample->numEntries() << ")" <<  endl;
436		}
437	+
438	+	}
439	+
440	+	void EdgeFitter::DrawDatasetContent(int is_data) {
441	+	TCanvas* ceedata = new TCanvas("ceedata","ceedata") ;
442	+
443	+	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
444	+	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
445	+	frame1->GetXaxis()->CenterTitle(1);
446	+	frame1->GetYaxis()->CenterTitle(1);
447	+	eeSample->plotOn(frame1,RooFit::Name("eedata")) ;
448	+
449	+	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
450	+	frame2->GetXaxis()->CenterTitle(1);
451	+	frame2->GetYaxis()->CenterTitle(1);
452	+	OFSample->plotOn(frame2,RooFit::Name("OFdata")) ;
453	+
454	+	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);
455	+	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);
456	+	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);
457	+
458	+	TH1F *ee_fit = (TH1F*)eeSample->createHistogram("ee_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
459	+	TH1F *mm_fit = (TH1F*)eeSample->createHistogram("mm_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
460	+	TH1F *of_fit = (TH1F*)OFSample->createHistogram("of_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
461	+
462	+	TLegend *leg = allsamples.allbglegend();
463	+	leg->AddEntry(ee_fit,"RooDataSet content","p");
464	+	leg->SetX1(0.58);
465	+
466	+	ceedata->cd() ;
467	+	gPad->SetLeftMargin(0.15);
468	+	frame1->GetYaxis()->SetTitleOffset(1.4);
469	+	frame1->Draw();
470	+	ee_ref->Draw("histo,same");
471	+	ee_fit->Draw("e1,same");
472	+	leg->Draw("same");
473	+	if(is_data==data) DrawPrelim();
474	+	else DrawPrelim(PlottingSetup::luminosity,true);
475	+	CompleteSave(ceedata,"Edge/ee_NoFit");
476	+
477	+	ceedata->cd() ;
478	+	gPad->SetLeftMargin(0.15);
479	+	frame1->GetYaxis()->SetTitleOffset(1.4);
480	+	frame1->Draw();
481	+	mm_ref->Draw("histo,same");
482	+	mm_fit->Draw("e1,same");
483	+	leg->Draw("same");
484	+	if(is_data==data) DrawPrelim();
485	+	else DrawPrelim(PlottingSetup::luminosity,true);
486	+	CompleteSave(ceedata,"Edge/mm_NoFit");
487	+
488	+	TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
489	+	cOFdata->cd() ;
490	+	gPad->SetLeftMargin(0.15);
491	+	frame2->SetMaximum(frame1->GetMaximum());
492	+	frame2->GetYaxis()->SetTitleOffset(1.4);
493	+	frame2->Draw();
494	+	of_ref->Draw("histo,same");
495	+	of_fit->Draw("e1,same");
496	+	leg->Draw("same");
497	+	if(is_data==data) DrawPrelim();
498	+	else DrawPrelim(PlottingSetup::luminosity,true);
499	+	CompleteSave(cOFdata,"Edge/OF_NoFit");
500	+
501	+	delete ee_fit;
502	+	delete mm_fit;
503	+	delete of_fit;
504	+
505	+	delete ee_ref;
506	+	delete mm_ref;
507	+	delete of_ref;
508	+
509	+	}
510	+
511	+	string WriteWithError(float central, float error, int digits) {
512	+	float ref=central;
513	+	if(ref<0) ref=-central;
514	+	int HighestSigDigit = 0;
515	+	if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
516	+	else HighestSigDigit = int(log(ref)/(log(10)));
517	+
518	+	float divider=pow(10.0,(double(HighestSigDigit-digits)));
519	+
520	+	stringstream result;
521	+	result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
522	+	return result.str();
523		}
524
525	+
526		string EdgeFitter::RandomStorageFile() {
527		TRandom3 *r = new TRandom3(0);
528		int rho = (int)r->Uniform(1,10000000);
#	Line 213 | Line 535 | string EdgeFitter::RandomStorageFile() {
535		Yield EdgeFitter::Get_Z_estimate(float jzb_cut, int icut) {
536		if(MarcoDebug) write_error(__FUNCTION__,"Returning random Z yield");
537		Yield a(123,45,67); return a;
216	–
538		return PlottingSetup::allresults.predictions[icut].Zbkg;
539		}
540
#	Line 223 | Line 544 | Yield EdgeFitter::Get_T_estimate(float j
544		return PlottingSetup::allresults.predictions[icut].Flavorsym;
545		}
546
547	+	void EdgeFitter::DoFit(int is_data, float jzb_cut) {
548	+
549	+	stringstream prefix;
550	+	if(is_data==data) prefix << "data_";
551	+	if(is_data==mc) prefix << "mc_";
552	+	if(is_data==mcwithsignal) prefix << "mcwithS_";
553	+
554	+	prefix << EdgeFitter::Mode << "_" << jzb_cut;
555	+
556	+	if(EdgeFitter::AllowTriangle) prefix << "__H1";
557	+	else prefix << "__H0";
558	+
559	+	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
560	+	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
561	+	RooCategory sample("sample","sample") ;
562	+	sample.defineType("ee");
563	+	sample.defineType("mm");
564	+	sample.defineType("OF");
565	+
566	+	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));
567	+
568	+
569	+	//First we make a fit to opposite flavor
570	+	RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 1.5*combData.sumEntries());
571	+	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
572	+	RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
573	+	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
574	+	RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
575	+	RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
576	+	RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
577	+	RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
578	+	simPdfOF.addPdf(model_OF,"OF");
579	+	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save(),RooFit::Extended(),RooFit::Minos(true));
580	+	dout << "============================= < OF results > =============================" << endl;
581	+	resultOF->Print();
582	+	dout << "============================= < /OF results > =============================" << endl;
583	+
584	+
585	+	RooPlot* frameO = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample"));
586	+	frameO->GetXaxis()->CenterTitle(1);
587	+	frameO->GetYaxis()->CenterTitle(1);
588	+	combData.plotOn(frameO,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
589	+	simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
590	+	simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
591	+
592	+	TCanvas* pof = new TCanvas("pof","pof") ;
593	+	pof->cd() ;
594	+	gPad->SetLeftMargin(0.15);
595	+	frameO->GetYaxis()->SetTitleOffset(1.4);
596	+	frameO->Draw();
597	+	if(is_data==data) DrawPrelim();
598	+	else DrawPrelim(PlottingSetup::luminosity,true);
599	+	if(EdgeFitter::FixedMEdge>=0) CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str()+"__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
600	+	else CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
601	+	delete pof;
602	+
603	+	if(resultOF->covQual()!=3) {
604	+	write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
605	+	dout << "covQual is " << resultOF->covQual() << endl;
606	+	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
607	+	else EdgeFitter::FixedMEdgeChi2_H0=-1;
608	+	if(EdgeFitter::RejectPointIfNoConvergence) return;
609	+	} else {
610	+	write_info(__FUNCTION__,"OF fit converged");
611	+	}
612	+
613	+	float StartingMedge=70;
614	+	if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
615	+
616	+
617	+	RooDataSet *ZDataSet = (RooDataSet*)EdgeFitter::AllData->reduce("id1==id2 && abs(mll-91.2)<20");
618	+
619	+	float maxZ = ZDataSet->sumEntries();
620	+	dout << "maxZ was set to " << maxZ << endl;
621	+	delete ZDataSet;
622	+
623	+
624	+	// Now same flavor
625	+	RooRealVar fzee("fzee", "fzee", 39, 39, maxZ);
626	+	RooRealVar meanzee("meanzee", "meanzee", 91.1876, 89, 95);
627	+	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0.5, 5);
628	+	RooRealVar widthzee("widthzee", "widthzee", 2.94);
629	+	widthzee.setConstant(1);
630	+
631	+	RooRealVar fzmm("fzmm", "fzmm", 39, 39, maxZ);
632	+	RooRealVar meanzmm("meanzmm", "meanzmm", 91.1876, 89, 95);
633	+	RooRealVar sigmazmm("sigmazmm", "sigmazmm", 5, 0.5, 5);
634	+	RooRealVar widthzmm("widthzmm", "widthzmm", 2.94);
635	+	widthzmm.setConstant(1);
636	+
637	+	//for ttbar only the relative fraction are different - all other parameters are the same.
638	+	RooRealVar fttbaree("fttbaree", "fttbaree", fttbarOF.getVal(), 0.2*fttbarOF.getVal(), 1.5*fttbarOF.getVal());
639	+	RooRealVar fttbarmm("fttbarmm", "fttbarmm", fttbarOF.getVal(), 0.2*fttbarOF.getVal(), 1.5*fttbarOF.getVal());
640	+
641	+	RooRealVar par1ttbaree("par1ttbaree", "par1ttbaree", 1.02*par1ttbarOF.getVal(), (1.02-0.07)*par1ttbarOF.getVal(), (1.02+0.07)*par1ttbarOF.getVal());
642	+	RooRealVar par1ttbarmm("par1ttbarmm", "par1ttbarmm", 1.02*par1ttbarOF.getVal(), (1.02-0.07)*par1ttbarOF.getVal(), (1.02+0.07)*par1ttbarOF.getVal());
643	+
644	+	//for signal only the resolution and absolute fraction are different for ee / mm
645	+	RooRealVar fsignalee("fsignalee", "fsignalee", 0, 0, 300);
646	+	RooRealVar fsignalmm("fsignalmm", "fsignalmm", 0, 0, 300);
647	+
648	+	RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
649	+	RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
650	+	RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
651	+
652	+	RooVoigtian zee("zee", "zee", mll, meanzee, widthzee, sigmazee);
653	+	RooVoigtian zmm("zmm", "zmm", mll, meanzmm, widthzmm, sigmazmm);
654	+
655	+	if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
656	+
657	+	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbaree, par2ttbarOF, par3ttbarOF, par4ttbarOF);
658	+	RooSUSYTPdf signalee("signalee","signalee", mll , par1signalSF, sigmazee, par3signalSF);
659	+
660	+	RooSUSYBkgPdf ttbarmm("ttbarmm","ttbarmm", mll , par1ttbarmm, par2ttbarOF, par3ttbarOF, par4ttbarOF);
661	+	RooSUSYTPdf signalmm("signalmm","signalmm", mll , par1signalSF, sigmazmm, par3signalSF);
662	+
663	+	RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
664	+	RooAddPdf model_mm("model_mm","model_mm", RooArgList(zmm, ttbarmm, signalmm), RooArgList(fzmm, fttbarmm, fsignalmm));
665	+
666	+
667	+	if(!EdgeFitter::AllowTriangle) {
668	+	fsignalee.setVal(0.0); // kill off the signal if we don't want the triangle
669	+	fsignalee.setConstant(1);
670	+	fsignalmm.setVal(0.0); // kill off the signal if we don't want the triangle
671	+	fsignalmm.setConstant(1);
672	+	par1signalSF.setConstant(1);
673	+	par2signalSF.setConstant(1);
674	+	par3signalSF.setConstant(1);
675	+	}
676	+
677	+	RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
678	+	simPdf.addPdf(model_ee,"ee") ;
679	+	simPdf.addPdf(model_mm,"mm") ;
680	+	simPdf.addPdf(model_OF,"OF") ;
681	+
682	+
683	+	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
684	+
685	+	if(result->covQual()!=3) {
686	+	write_error(__FUNCTION__,"Full fit did not converge!!! Cannot continue!");
687	+	dout << "covQual is " << result->covQual() << endl;
688	+	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
689	+	else EdgeFitter::FixedMEdgeChi2_H0=-1;
690	+	if(EdgeFitter::RejectPointIfNoConvergence) return;
691	+	} else {
692	+	write_info(__FUNCTION__,"Full fit converged");
693	+	}
694	+
695	+	dout << "============================= < Full results > =============================" << endl;
696	+	result->Print();
697	+	dout << "============================= < /Full results > =============================" << endl;
698	+
699	+
700	+	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
701	+	frame1->GetXaxis()->CenterTitle(1);
702	+	frame1->GetYaxis()->CenterTitle(1);
703	+	combData.plotOn(frame1,RooFit::Name("eedata"),RooFit::Cut("sample==sample::ee")) ;
704	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
705	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
706	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("DYSFonly"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
707	+	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));
708	+
709	+	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame1->chiSquare("FullFit", "SFdata", 3);
710	+	else EdgeFitter::FixedMEdgeChi2_H0=frame1->chiSquare("FullFit", "SFdata", 3);
711	+
712	+	dout << "Result   : " << endl;
713	+	if(EdgeFitter::AllowTriangle) dout << "f signal : " << fsignalee.getVal() << " +/- " << fsignalee.getError() << endl;
714	+	if(EdgeFitter::AllowTriangle) dout << "f signal : " << fsignalmm.getVal() << " +/- " << fsignalmm.getError() << endl;
715	+	//  dout << "f ttbar  : " << fttbaree.getVal() << " +/- " << fttbaree.getError() << endl;
716	+	//  dout << "f ttbar  : " << fttbarmm.getVal() << " +/- " << fttbarmm.getError() << endl;
717	+	dout << "f ttbaree: " << fttbaree.getVal() << " +/- NO ERROR CUZ ITS A PDF "<< endl;
718	+	dout << "f ttbarmm: " << fttbarmm.getVal() << " +/- NO ERROR CUZ ITS A PDF "<< endl;
719	+	dout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
720	+	dout << "f z ee   : " << fzee.getVal() << " +/- " << fzee.getError() << endl;
721	+	dout << "f z mm   : " << fzmm.getVal() << " +/- " << fzmm.getError() << endl;
722	+	if(EdgeFitter::AllowTriangle) dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
723	+	else dout << "#Chi^{2}/NDF     : " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
724	+
725	+	// The same plot for the cointrol sample slice
726	+	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
727	+	frame3->GetXaxis()->CenterTitle(1);
728	+	frame3->GetYaxis()->CenterTitle(1);
729	+	frame3->SetMaximum(frame1->GetMaximum());
730	+	combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
731	+	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
732	+	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
733	+
734	+
735	+	TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
736	+	c->cd() ;
737	+	gPad->SetLeftMargin(0.15);
738	+	frame1->GetYaxis()->SetTitleOffset(1.4);
739	+	frame1->Draw();
740	+	if(is_data==data) DrawPrelim();
741	+	else DrawPrelim(PlottingSetup::luminosity,true);
742	+	stringstream infotext;
743	+	infotext << "#splitline{Fit results (" << EdgeFitter::Mode << ">" << jzb_cut << "): }{#splitline{";
744	+	infotext << "N(Data) = " << EdgeFitter::eeSample->sumEntries()+EdgeFitter::mmSample->sumEntries() << "}{#splitline{";
745	+	infotext << "N(Z+Jets) = " << WriteWithError(fzee.getVal()+fzmm.getVal(),sqrt(pow(fzee.getError(),2)+pow(fzmm.getError(),2)),3) << "}{#splitline{";
746	+	//infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
747	+	infotext << "N(t#bar{t}) = " << WriteWithError(fttbaree.getVal()+fttbarmm.getVal(),0,3) << "}{#splitline{";
748	+	write_warning(any2string(__LINE__),"Don't have the error yet, need to complete this");
749	+	if(EdgeFitter::AllowTriangle) {
750	+	infotext << "N(signal) = " << WriteWithError(fsignalee.getVal()+fsignalmm.getVal(),sqrt(pow(fsignalee.getError(),2)+pow(fsignalmm.getError(),2)),3) << "}{";
751	+	infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
752	+	} else infotext << "}{}}}}}";
753	+
754	+	TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
755	+	infobox->SetNDC();
756	+	infobox->SetTextSize(0.03);
757	+	infobox->Draw();
758	+	if(EdgeFitter::FixedMEdge>=0) CompleteSave(c,"Edge/"+prefix.str()+"_SF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
759	+	else CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
760	+	delete c;
761	+
762	+	TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
763	+	e->cd();
764	+	gPad->SetLeftMargin(0.15);
765	+	frame3->GetYaxis()->SetTitleOffset(1.4);
766	+	frame3->Draw();
767	+	if(is_data==data) DrawPrelim();
768	+	else DrawPrelim(PlottingSetup::luminosity,true);
769	+	if(EdgeFitter::FixedMEdge>=0) CompleteSave(e,"Edge/"+prefix.str()+"_OF__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge),false,false);
770	+	else CompleteSave(e,"Edge/"+prefix.str()+"_OF",false,false);
771	+	delete e;
772	+
773	+
774	+
775	+
776	+	/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
777	+	f->cd();
778	+	gPad->SetLeftMargin(0.15);
779	+	frame4->GetYaxis()->SetTitleOffset(1.4);
780	+	frame4->Draw();
781	+	if(is_data==data) DrawPrelim();
782	+	else DrawPrelim(PlottingSetup::luminosity,true);
783	+	CompleteSave(f,"Edge/"+prefix.str()+"_SF");
784	+	delete f;*/
785	+
786	+
787	+	/*
788	+	float maxZ=200;
789	+	RooWorkspace* wspace = new RooWorkspace();
790	+	stringstream mllvar;
791	+	mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
792	+	wspace->factory(mllvar.str().c_str());
793	+	wspace->var("mll")->setBins(30);
794	+	wspace->factory("nSig[1.,0.,100.]");
795	+	wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
796	+	wspace->factory("rME[1.12,1.05,1.19]");
797	+	wspace->factory("effUncert[1.]");
798	+	EdgeFitter::prepareLimits(wspace, true);
799	+	*/
800	+
801	+	write_warning(__FUNCTION__," A lot missing here to calculate limits");
802	+
803	+	}
804	+
805		void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
806
807	<	string storagefile=EdgeFitter::RandomStorageFile();
229	<	TFile *f = new TFile(storagefile.c_str(),"RECREATE");
230	<	EdgeFitter::InitializeVariables(iMllLow,iMllHigh,jzbHigh,cut);
231	<
232	<	Yield Zestimate=EdgeFitter::Get_Z_estimate(jzb_cut,icut);
233	<	Yield Testimate=EdgeFitter::Get_T_estimate(jzb_cut,icut);
234	<	cout << "Cut at JZB>" << jzb_cut << "; using estimates: " << endl;
235	<	cout << "       Z : " << Zestimate  << endl;
236	<	cout << "       T : " << Testimate  << endl;
807	>	TCut _cut(cut&&PlottingSetup::basiccut&&PlottingSetup::passtrig);
808
809	<	EdgeFitter::PrepareDatasets(is_data);
809	>	TFile *f = new TFile("workingfile.root","RECREATE");
810
811	<	EdgeFitter::eeSample->Write();
812	<	EdgeFitter::emSample->Write();
813	<	EdgeFitter::mmSample->Write();
814	<	EdgeFitter::AllData->Write();
811	>	EdgeFitter::InitializeVariables(PlottingSetup::iMllLow,PlottingSetup::iMllHigh,PlottingSetup::jzbHigh,_cut);
812	>
813	>	EdgeFitter::PrepareDatasets(is_data);
814	>
815	>	EdgeFitter::DrawDatasetContent(is_data);
816	>
817	>	RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
818	>	RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
819	>
820	>
821	>	EdgeFitter::AllowTriangle=false;
822	>	EdgeFitter::DoFit(is_data, jzb_cut);
823	>
824	>	write_info(__FUNCTION__,"TAKING SHORTCUT");return;
825	>
826	>	EdgeFitter::AllowTriangle=true;
827	>
828	>	bool ScanMassRange=false;
829	>	float ScanSteps=5.0;//GeV
830	>
831	>
832	>	if(ScanMassRange) {
833	>	TFile *fscan = new TFile("fscan.root","UPDATE");
834	>	TGraph *gr = new TGraph();
835	>	TGraph *Rgr = new TGraph();
836	>	stringstream GrName;
837	>	GrName << "ScanGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
838	>	stringstream RGrName;
839	>	RGrName << "ScanRatioGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
840	>	gr->SetName(GrName.str().c_str());
841	>	Rgr->SetName(RGrName.str().c_str());
842	>
843	>	int i=0;
844	>	for(float tempMedge=10;tempMedge<=300;tempMedge+=ScanSteps) {
845	>	write_info(__FUNCTION__,"Now testing Medge="+any2string(tempMedge)+" for "+EdgeFitter::Mode+">"+any2string(jzb_cut));
846	>	EdgeFitter::FixedMEdge=tempMedge;
847	>	EdgeFitter::DoFit(is_data, jzb_cut);
848	>	if(EdgeFitter::FixedMEdgeChi2_H1<0) continue;
849	>	gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1);
850	>	Rgr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1/EdgeFitter::FixedMEdgeChi2_H0);
851	>	i++;
852	>	}
853	>
854	>	TCanvas *ScanCan = new TCanvas("ScanCan","ScanCan",500,500);
855	>	gr->GetXaxis()->SetTitle("m_{edge}");
856	>	gr->GetXaxis()->CenterTitle();
857	>	gr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
858	>	gr->GetYaxis()->CenterTitle();
859	>	gr->GetYaxis()->SetTitleOffset(0.95);
860	>	gr->GetXaxis()->SetTitleOffset(0.9);
861	>	gr->SetLineColor(kBlue);
862	>	gr->SetTitle("");
863	>	gr->Draw("AL");
864	>	stringstream ScanCanSave;
865	>	ScanCanSave << "Edge/MEdgeScan_"+EdgeFitter::Mode+"_" << jzb_cut;
866	>	if(is_data) DrawPrelim();
867	>	else DrawMCPrelim();
868	>	CompleteSave(ScanCan,ScanCanSave.str());
869	>
870	>	Rgr->GetXaxis()->SetTitle("m_{edge}");
871	>	Rgr->GetXaxis()->CenterTitle();
872	>	Rgr->GetYaxis()->SetTitle("#Chi^{2} / NDF");
873	>	Rgr->GetYaxis()->CenterTitle();
874	>	Rgr->GetYaxis()->SetTitleOffset(0.95);
875	>	Rgr->GetXaxis()->SetTitleOffset(0.9);
876	>	Rgr->SetLineColor(kBlue);
877	>	Rgr->SetTitle("");
878	>	Rgr->Draw("AL");
879	>	ScanCanSave.str("");
880	>	ScanCanSave << "Edge/MEdgeScan_Ratio_"+EdgeFitter::Mode+"_" << jzb_cut;
881	>	if(is_data) DrawPrelim();
882	>	else DrawMCPrelim();
883	>	CompleteSave(ScanCan,ScanCanSave.str());
884	>	fscan->cd();
885	>	gr->Write();
886	>	delete ScanCan;
887	>	fscan->Close();
888	>	} else {
889	>	EdgeFitter::DoFit(is_data, jzb_cut);
890	>	dout << "Chi^2 (H0) = " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
891	>	dout << "Chi^2 (H1) = " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
892	>	}
893	>
894	>
895	>	RooMsgService::instance().setGlobalKillBelow(msglevel);
896	>
897		f->Close();
898
246	–	write_warning(__FUNCTION__,"Work continues here");
247	–
248	–	if(EdgeFitter::MarcoDebug) write_warning(__FUNCTION__,"Need to uncomment the next line (remove the output file)");
249	–	//  gSystem->Exec(("rm "+storagefile).c_str());
899		}
900
901		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
902	<	for(int icut=0;icut<jzb_cut.size();icut++) {
902	>
903	>	EdgeFitter::Mode="JZB";
904	>	if(mcjzb=="met[4]") EdgeFitter::Mode="MET";
905	>
906	>	for(int icut=0;icut<(int)jzb_cut.size();icut++) {
907		stringstream addcut;
908		if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
909		if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";
910		TCut jcut(addcut.str().c_str());
911
912	+
913		EdgeFitter::DoEdgeFit(mcjzb, datajzb, DataPeakError, MCPeakError, jzb_cut[icut], icut, is_data, jcut&&cut, signalevents);
914
915		}

Diff Legend

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