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.30 by buchmann, Thu Jun 27 06:02:20 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 "ExtendedMath.C"
33	>	#include "EdgeModules/RooSUSYTPdf.cxx"
34	>	#include "EdgeModules/RooSUSYBkgPdf.cxx"
35	>	#include "EdgeModules/RooSUSYCompleteBkgPdf.cxx"
36
37
38	+	#include "md5/md5.h"
39	+	#include "md5/md5.C"
40	+
41		using namespace std;
42		using namespace PlottingSetup;
43
#	Line 51 | Line 62 | namespace EdgeFitter {
62
63		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree*);
64		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree*);
65	<	void getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
65	>	void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
66		void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
67		void PrepareDatasets(int);
68	+	void DrawDatasetContent(int);
69		void DoFit(int is_data, float jzb_cut);
70		string RandomStorageFile();
71		Yield Get_Z_estimate(float,int);
72		Yield Get_T_estimate(float,int);
73	<	float calcExclusion(RooWorkspace *ws, RooDataSet *data = NULL);
74	<	void prepareLimits(RooWorkspace *ws);
63	<	vector<RooDataSet*> generateToys(RooWorkspace *ws, int nToys);
73	>	float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
74	>	vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
75		void prepareLimits(RooWorkspace *ws, bool ComputeBands);
76		TGraph* prepareLM(float mass, float nEv);
77
#	Line 72 | Line 83 | namespace EdgeFitter {
83		RooDataSet* AllData;
84		RooDataSet* eeSample;
85		RooDataSet* mmSample;
86	<	RooDataSet* emSample;
86	>	RooDataSet* OFSample;
87	>	RooDataSet* SFSample;
88	>
89	>	bool MarcoDebug=true;
90	>
91	>	float FixedMEdge=-1;
92	>	float FixedMEdgeChi2_H1=-1;
93	>	float FixedMEdgeChi2_H0=-1;
94
95	<	bool MarcoDebug;
95	>	bool RejectPointIfNoConvergence=false;
96	>
97	>	string Mode="UndefinedMode";
98	>
99	>	bool AllowTriangle=true;
100		}
101
102		TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
#	Line 93 | Line 115 | TGraph* EdgeFitter::prepareLM(float mass
115		return lm;
116		}
117
118	<	vector<RooDataSet*> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
118	>	vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
119	>	ws->ls();
120		ws->var("nSig")->setVal(0.);
121		ws->var("nSig")->setConstant(true);
122		RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
123	<	vector<RooDataSet*> theToys;
123	>	vector<RooDataSet> theToys;
124
125		RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
126		mcEE.generate(nToys,44,true);
#	Line 106 | Line 129 | vector<RooDataSet*> EdgeFitter::generate
129		RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
130		mcOSOF.generate(nToys,44,true);
131
132	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
132	>	RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
133		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
134		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
135		RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
#	Line 120 | Line 143 | vector<RooDataSet*> EdgeFitter::generate
143		stringstream toyname;
144		toyname << "theToy_" << i;
145		write_warning(__FUNCTION__,"Problem while adding toys");
146	<	//    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));
147	<	//    theToys.push_back(toyData);
146	>	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));
147	>	theToys.push_back(toyData);
148		}
149		ws->var("nSig")->setVal(17.0);
150		ws->var("nSig")->setConstant(false);
151		return theToys;
152		}
153
154	<	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
154	>	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
155		TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
156		vector<float> theLimits;
157	<	for(int itoy=0;itoy<theToys.size();itoy++) {
158	<	float theLimit = calcExclusion(ws,theToys[itoy]);
157	>	for(int itoy=0;itoy<(int)theToys.size();itoy++) {
158	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
159		if(theLimit > 0 ) gauLimit->Fill(theLimit);
160		}
161		const Int_t nQ = 4;
#	Line 145 | Line 168 | void EdgeFitter::getMedianLimit(RooWorks
168		sigmaDown = abs(yQ[1]-median);
169		sigmaUp = abs(yQ[2]-median);
170		twoSigmaUp = abs(yQ[3]-median);
171	<	cout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
171	>	dout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
172		}
173
174		void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
175		if(ComputeBands) {
176	<	vector<RooDataSet*> theToys = EdgeFitter::generateToys(ws,50);
176	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
177		vector<float> medVals;
178		vector<float> medLimits;
179		vector<float> sigmaLimitsDown;
#	Line 179 | Line 202 | void EdgeFitter::prepareLimits(RooWorksp
202		theVals.push_back((float)i);
203		ws->var("m0")->setVal((float)i);
204		ws->var("m0")->setConstant(true);
205	<	theLimits.push_back(calcExclusion(ws));
205	>	//      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
206	>	write_error(__FUNCTION__,"Error while casting roo data set");
207		}
208
209	<	for(int i=0;i<theLimits.size();i++) {
209	>	for(int i=0;i<(int)theLimits.size();i++) {
210		if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
211	<	cout << i << " : " << theVals[i] << endl;
211	>	dout << i << " : " << theVals[i] << endl;
212		theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
213		write_warning(__FUNCTION__,"Need to store limits");
214		}
#	Line 194 | Line 218 | void EdgeFitter::prepareLimits(RooWorksp
218		}
219
220
221	<	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet *data) {
221	>	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
222	>	int numberOfToys=50;
223		RooRealVar mu("mu","nSig",0,10000,"");
224		RooArgSet poi = RooArgSet(mu);
225		RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
#	Line 203 | Line 228 | float EdgeFitter::calcExclusion(RooWorks
228		model->SetWorkspace(*ws);
229		model->SetPdf("combModel");
230		model->SetParametersOfInterest(poi);
231	<	if(!data) data = (RooDataSet*)ws->data("data_obs");
231	>	//  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
232
233	<	RooStats::ProfileLikelihoodCalculator plc(*data, *model);
233	>	RooStats::ProfileLikelihoodCalculator plc(data, *model);
234		plc.SetNullParameters(*nullParams);
235		plc.SetTestSize(0.05);
236	+
237		RooStats::LikelihoodInterval* interval = plc.GetInterval();
238		RooStats::HypoTestResult *htr = plc.GetHypoTest();
239		double theLimit = interval->UpperLimit( mu );
240	<	cout << "Significance " << htr->Significance() << endl;
240	>	//  double significance = htr->Significance();
241
242		ws->defineSet("obs","nB");
243		ws->defineSet("poi","nSig");
#	Line 235 | Line 261 | float EdgeFitter::calcExclusion(RooWorks
261		slrts.SetAltParameters(*sb_model.GetSnapshot());
262		RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
263
264	<	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
264	>	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
265		hc.SetTestStatistic(2);
266	<	hc.SetNumberOfToys(50);
266	>	hc.SetNumberOfToys(numberOfToys);
267
268		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
269		hcInv.SetTestSize(0.05);
#	Line 246 | Line 272 | float EdgeFitter::calcExclusion(RooWorks
272		RooStats::HypoTestInverterResult* hcInt = hcInv.GetInterval();
273		float ulError = hcInt->UpperLimitEstimatedError();
274		theLimit = hcInt->UpperLimit();
275	<	cout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
275	>	dout << "Found upper limit : " << theLimit << "+/-" << ulError << endl;
276
277		return theLimit/PlottingSetup::luminosity;
278
279		}
280
281	+	string GetSkimName(int isample) {
282	+	return removefunnystring(allsamples.collection[isample].filename);
283	+	}
284	+
285		TTree* SkimTree(int isample) {
286	<	TTree* newTree = allsamples.collection[isample].events->CloneTree(0);
286	>	string TreeName = GetSkimName(isample);
287	>	cout << "About to skim " << TreeName << " for sample id " << isample << endl;
288	>	TTree* newTree = new TTree(TreeName.c_str(),TreeName.c_str());
289	>
290	>	float mll,edgeWeight;
291	>	int id1,id2;
292	>
293	>	newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
294	>	newTree->Branch("mll",&mll,"mll/F");
295	>	newTree->Branch("id1",&id1,"id1/I");
296	>	newTree->Branch("id2",&id2,"id2/I");
297	>	newTree->Branch("isample",&isample,"isample/I");
298	>
299		float xsweight=1.0;
300		if(allsamples.collection[isample].is_data==false) xsweight=luminosity*allsamples.collection[isample].weight;
301		if(EdgeFitter::MarcoDebug) {
302	<	cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
303	<	cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
302	>	dout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
303	>	dout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
304		}
305	+
306	+	float tmll;
307	+	int tid1,tid2;
308	+	allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
309	+	allsamples.collection[isample].events->SetBranchAddress("id1",&tid1);
310	+	allsamples.collection[isample].events->SetBranchAddress("id2",&tid2);
311	+
312		TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
313	+	TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
314	+
315		float wgt=1.0;
265	–	allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
316		for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
317		allsamples.collection[isample].events->LoadTree(entry);
318		if (select->EvalInstance()) {
319		allsamples.collection[isample].events->GetEntry(entry);
320	<	wgt=wgt*xsweight;
320	>	mll=tmll;
321	>	id1=tid1;
322	>	id2=tid2;
323	>	wgt=Weight->EvalInstance();
324	>	edgeWeight=wgt*xsweight;
325		newTree->Fill();
326		}
327		}
328
329	<	if(EdgeFitter::MarcoDebug) cout << "     Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
329	>	if(EdgeFitter::MarcoDebug) dout << "   Reduced tree contains " << newTree->GetEntries() << " entries " << endl;
330		return newTree;
331		}
332
#	Line 282 | Line 336 | void EdgeFitter::InitializeVariables(flo
336		jzbmax=_jzbmax;
337		cut=_cut;
338		}
339	<
339	>
340	>	TTree* MergeTrees(vector<TTree*> trees) {
341	>	assert(trees.size()>0);
342	>	TTree * newtree = (TTree*)trees[0]->CloneTree(0);
343	>	newtree->SetTitle("FullTree");
344	>	newtree->SetName("FullTree");
345	>	trees[0]->GetListOfClones()->Remove(newtree);
346	>	trees[0]->ResetBranchAddresses();
347	>	newtree->ResetBranchAddresses();
348	>
349	>	for(int itree=0;itree<trees.size();itree++) {
350	>	newtree->CopyAddresses(trees[itree]);
351	>	Long64_t nentries = trees[itree]->GetEntries();
352	>	for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
353	>	trees[itree]->GetEntry(iEntry);
354	>	newtree->Fill();
355	>	}
356	>	trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
357	>	if (newtree->GetTreeIndex()) {
358	>	newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
359	>	}
360	>	if (newtree && newtree->GetTreeIndex()) {
361	>	newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
362	>	}
363	>	}
364	>	return newtree;
365	>	}
366	>
367	>
368	>
369		void EdgeFitter::PrepareDatasets(int is_data) {
287	–	TTree *completetree;
370		write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
371	<	bool hashit=0;
372	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
371	>
372	>	ensure_directory_exists(PlottingSetup::cbafbasedir+"/EdgeCache/");
373	>
374	>	stringstream FileName;
375	>	FileName << PlottingSetup::cbafbasedir << "/EdgeCache/" << md5((const char*)cut) << ".root";
376	>
377	>
378	>	TFile *fEdgeCache;
379	>	if(PlottingSetup::do_CleanCache) fEdgeCache = new TFile(FileName.str().c_str(),"RECREATE");
380	>	else fEdgeCache = new TFile(FileName.str().c_str(),"UPDATE");
381	>
382	>	if(fEdgeCache->GetNkeys()==0) {
383	>	ofstream CacheOverview;
384	>	CacheOverview.open((PlottingSetup::cbafbasedir+"/EdgeCache/CacheOverview").c_str(),ios::app);
385	>	CacheOverview << md5((const char*)cut) << ";" << (const char*) cut << endl;
386	>	CacheOverview.close();
387	>	}
388	>
389	>	vector<TTree*> SkimmedTrees;
390	>	TTree *SkimmedTree[(int)allsamples.collection.size()];
391	>	for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
392		if(!allsamples.collection[isample].is_active) continue;
393		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
394		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
395		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
396		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
397	<	if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
398	<	if(!hashit) {
399	<	hashit=true;
400	<	completetree = SkimTree(isample)->CloneTree();
397	>	if(EdgeFitter::MarcoDebug) dout << "Considering : " << allsamples.collection[isample].samplename << endl;
398	>
399	>	string SkimName = GetSkimName(isample);
400	>	SkimmedTree[isample] = (TTree*)fEdgeCache->Get(SkimName.c_str());
401	>	if(!SkimmedTree[isample]) {
402	>	dout << "Need to generate tree for " << allsamples.collection[isample].filename << endl;
403	>	SkimmedTree[isample] = SkimTree(isample);
404	>	fEdgeCache->cd();
405	>	SkimmedTree[isample]->Write();
406		} else {
407	<	completetree->CopyEntries(SkimTree(isample));
407	>	dout << "Loaded tree for " << allsamples.collection[isample].filename << " from cache file " << endl;
408		}
409	<	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
409	>
410	>	SkimmedTrees.push_back(SkimmedTree[isample]);
411		}
412
413	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
413	>	TTree *completetree  = MergeTrees(SkimmedTrees);
414	>
415	>	if(EdgeFitter::MarcoDebug) dout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
416	>
417	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
418		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
419		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
420	<	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
421	<	RooRealVar weight("weight","weight",0,1000,"");
311	<	RooArgSet observables(mll,jzb,id1,id2,weight);
420	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
421	>	RooArgSet observables(mll,id1,id2,edgeWeight);
422
423		string title="CMS Data";
424		if(is_data!=1) title="CMS SIMULATION";
425	<	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
426	<	completetree->Write();
427	<	//  delete completetree;
428	<
429	<	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2");
430	<	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2");
431	<	EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
425	>	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
426	>	fEdgeCache->Close();
427	>	dout << "Edge cache closed." << endl;
428	>
429	>	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==0");
430	>	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==1");
431	>	EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
432	>	EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
433		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
434
435	+
436		eeSample->SetName("eeSample");
437		mmSample->SetName("mmSample");
438	<	emSample->SetName("emSample");
438	>	OFSample->SetName("OFSample");
439	>	SFSample->SetName("SFSample");
440		AllData->SetName("AllData");
441
442	+	eeSample->SetTitle("eeSample");
443	+	mmSample->SetTitle("mmSample");
444	+	OFSample->SetTitle("OFSample");
445	+	SFSample->SetTitle("SFSample");
446	+	AllData->SetTitle("AllData");
447	+
448		if(EdgeFitter::MarcoDebug) {
449	<	cout << "Number of events in data sample = " << AllData->numEntries() << endl;
450	<	cout << "Number of events in ee sample = " << eeSample->numEntries() << endl;
451	<	cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
452	<	cout << "Number of events in em sample = " << emSample->numEntries() << endl;
449	>	dout << "Number of (weighted) events in full sample = " << AllData->sumEntries() << " (raw events : " << AllData->numEntries() << ")" <<  endl;
450	>	dout << "Number of (weighted) events in ee sample = " << eeSample->sumEntries() << " (raw events : " << eeSample->numEntries() << ")" <<  endl;
451	>	dout << "Number of (weighted) events in mm sample = " << mmSample->sumEntries() << " (raw events : " << mmSample->numEntries() << ")" <<  endl;
452	>	dout << "Number of (weighted) events in SF sample = " << SFSample->sumEntries() << " (raw events : " << SFSample->numEntries() << ")" <<  endl;
453	>	dout << "Number of (weighted) events in em sample = " << OFSample->sumEntries() << " (raw events : " << OFSample->numEntries() << ")" <<  endl;
454		}
455	+
456		}
457
458	+	void EdgeFitter::DrawDatasetContent(int is_data) {
459	+	TCanvas* ceedata = new TCanvas("ceedata","ceedata") ;
460	+
461	+	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
462	+	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("ee sample")) ;
463	+	frame1->GetXaxis()->CenterTitle(1);
464	+	frame1->GetYaxis()->CenterTitle(1);
465	+	eeSample->plotOn(frame1,RooFit::Name("eedata"));
466	+
467	+	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("mm sample")) ;
468	+	frame2->GetXaxis()->CenterTitle(1);
469	+	frame2->GetYaxis()->CenterTitle(1);
470	+	mmSample->plotOn(frame2,RooFit::Name("mmdata"));
471	+
472	+	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
473	+	frame3->GetXaxis()->CenterTitle(1);
474	+	frame3->GetYaxis()->CenterTitle(1);
475	+	OFSample->plotOn(frame3,RooFit::Name("OFdata"));
476	+
477	+	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);
478	+	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);
479	+	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);
480	+
481	+	ee_ref->SetFillColor(TColor::GetColor("#3ADF00"));
482	+	mm_ref->SetFillColor(TColor::GetColor("#3ADF00"));
483	+	of_ref->SetFillColor(TColor::GetColor("#3ADF00"));
484	+
485	+	TH1F *ee_fit = (TH1F*)eeSample->createHistogram("ee_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
486	+	TH1F *mm_fit = (TH1F*)mmSample->createHistogram("mm_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
487	+	TH1F *of_fit = (TH1F*)OFSample->createHistogram("of_fit",mll,RooFit::Binning(int((mllmax-mllmin)/5)),RooFit::SumW2Error(true));
488	+
489	+	TLegend *leg = make_legend();
490	+	leg->AddEntry(ee_fit,"RooDataSet content","p");
491	+	leg->AddEntry(ee_ref,"CBAF cross-check","f");
492	+	leg->SetX1(0.5);
493	+	leg->SetY1(0.75);
494	+
495	+	ceedata->cd() ;
496	+	gPad->SetLeftMargin(0.15);
497	+	frame1->SetMaximum(1.2*frame1->GetMaximum());
498	+	frame1->GetYaxis()->SetTitleOffset(1.4);
499	+	frame1->Draw();
500	+	ee_ref->Draw("histo,same");
501	+	ee_fit->Draw("e1,same");
502	+	leg->Draw("same");
503	+	if(is_data==data) DrawPrelim();
504	+	else DrawPrelim(PlottingSetup::luminosity,true);
505	+	CompleteSave(ceedata,"Edge/PreFit_ee");
506	+
507	+	ceedata->cd() ;
508	+	gPad->SetLeftMargin(0.15);
509	+	frame2->SetMaximum(1.2*frame2->GetMaximum());
510	+	frame2->GetYaxis()->SetTitleOffset(1.4);
511	+	frame2->Draw();
512	+	mm_ref->Draw("histo,same");
513	+	mm_fit->Draw("e1,same");
514	+	leg->Draw("same");
515	+	if(is_data==data) DrawPrelim();
516	+	else DrawPrelim(PlottingSetup::luminosity,true);
517	+	CompleteSave(ceedata,"Edge/PreFit_mm");
518	+
519	+	TCanvas* cOFdata = new TCanvas("cOFdata","cOFdata") ;
520	+	cOFdata->cd() ;
521	+	gPad->SetLeftMargin(0.15);
522	+	frame3->SetMaximum(1.2*frame3->GetMaximum());
523	+	frame3->GetYaxis()->SetTitleOffset(1.4);
524	+	frame3->Draw();
525	+	of_ref->Draw("histo,same");
526	+	of_fit->Draw("e1,same");
527	+	leg->Draw("same");
528	+	if(is_data==data) DrawPrelim();
529	+	else DrawPrelim(PlottingSetup::luminosity,true);
530	+	CompleteSave(cOFdata,"Edge/PreFit_OF");
531	+
532	+	delete ee_fit;
533	+	delete mm_fit;
534	+	delete of_fit;
535	+
536	+	delete ee_ref;
537	+	delete mm_ref;
538	+	delete of_ref;
539	+	}
540	+
541	+	string WriteWithError(float central, float error, int digits) {
542	+	float ref=central;
543	+	if(ref<0) ref=-central;
544	+	int HighestSigDigit = 0;
545	+	if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
546	+	else HighestSigDigit = int(log(ref)/(log(10)));
547	+
548	+	float divider=pow(10.0,(double(HighestSigDigit-digits)));
549	+
550	+	stringstream result;
551	+	result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
552	+	return result.str();
553	+	}
554	+
555	+
556		string EdgeFitter::RandomStorageFile() {
557		TRandom3 *r = new TRandom3(0);
558		int rho = (int)r->Uniform(1,10000000);
#	Line 356 | Line 575 | Yield EdgeFitter::Get_T_estimate(float j
575		}
576
577		void EdgeFitter::DoFit(int is_data, float jzb_cut) {
578	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
578	>
579	>	stringstream prefix;
580	>	if(is_data==data) prefix << "data_";
581	>	if(is_data==mc) prefix << "mc_";
582	>	if(is_data==mcwithsignal) prefix << "mcwithS_";
583	>
584	>	prefix << EdgeFitter::Mode << "_" << jzb_cut;
585	>
586	>	if(!EdgeFitter::AllowTriangle && EdgeFitter::FixedMEdge>=0) prefix << "__MEdgeFix_"+any2string(EdgeFitter::FixedMEdge);
587	>
588	>	if(EdgeFitter::AllowTriangle) prefix << "__H1";
589	>	else prefix << "__H0";
590	>
591	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
592	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
593		RooCategory sample("sample","sample") ;
594		sample.defineType("ee");
595	<	//sample.defineType("mm");
596	<	sample.defineType("em");
597	<	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("ee",*eeSample),Import("mm",*mmSample),Import("em",*emSample));
598	<	RooDataSet combData("combData","combined data",mll,RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("em",*emSample));
595	>	sample.defineType("mm");
596	>	sample.defineType("OF");
597	>	sample.defineType("SF");
598	>
599	>	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("mm",*mmSample),RooFit::Import("OF",*OFSample),RooFit::Import("SF",*SFSample),RooFit::WeightVar(edgeWeight));
600	>
601	>	//------------------------------------------------------------------------------------------------------------------------------------
602	>	//    _____ _               __
603	>	//   / ____| |             /_ |
604	>	//  | (___ | |_ ___ _ __    | |
605	>	//   \___ \| __/ _ \ '_ \   | |
606	>	//   ____) | ||  __/ |_) |  | |
607	>	//  |_____/ \__\___| .__/   |_|
608	>	//                 | |
609	>	//                 |_|
610	>	//
611	>	// Fit OF to get good initial values for flavor-symmetric background
612	>
613	>
614
367	–
368	–
615		//First we make a fit to opposite flavor
616	<	RooRealVar fttbarem("fttbarem", "fttbarem", 100, 0, 10000);
617	<	RooRealVar par1ttbarem("par1ttbarem", "par1ttbarem", 1.6, 0.01, 4.0);
618	<	RooRealVar par2ttbarem("par2ttbarem", "par2ttbarem", 1.0);
619	<	RooRealVar par3ttbarem("par3ttbarem", "par3ttbarem", 0.028, 0.001, 1.0);
620	<	RooRealVar par4ttbarem("par4ttbarem", "par4ttbarem", 2.0);
621	<	RooSUSYBkgPdf ttbarem("ttbarem","ttbarem", mll , par1ttbarem, par2ttbarem, par3ttbarem, par4ttbarem);
622	<	RooAddPdf model_em("model_em","model_em", ttbarem, fttbarem);
616	>	RooRealVar fttbarOF("fttbarOF", "fttbarOF", 0.8*OFSample->sumEntries(), 0, 1.2*OFSample->sumEntries());
617	>	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.77, 0.01, 4.0);
618	>	RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
619	>	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.0272, 0.001, 1.0);
620	>	RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
621	>	RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
622	>	RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
623		RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
624	<	simPdfOF.addPdf(model_em,"em");
625	<	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save());
624	>	simPdfOF.addPdf(model_OF,"OF");
625	>	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save(),RooFit::Extended(),RooFit::Minos(true));
626	>	dout << "============================= < OF results > =============================" << endl;
627		resultOF->Print();
628	+	dout << "============================= < /OF results > =============================" << endl;
629	+
630	+
631	+	RooPlot* frameO = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample"));
632	+	frameO->GetXaxis()->CenterTitle(1);
633	+	frameO->GetYaxis()->CenterTitle(1);
634	+	combData.plotOn(frameO,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
635	+	simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
636	+	simPdfOF.plotOn(frameO,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
637	+
638	+	TCanvas* pof = new TCanvas("pof","pof") ;
639	+	pof->cd() ;
640	+	gPad->SetLeftMargin(0.15);
641	+	frameO->GetYaxis()->SetTitleOffset(1.4);
642	+	frameO->Draw();
643	+	if(is_data==data) DrawPrelim();
644	+	else DrawPrelim(PlottingSetup::luminosity,true);
645	+	CompleteSave(pof,"Edge/OF__OFFitonly_"+prefix.str(),false,false);
646	+	delete pof;
647	+
648	+	TCanvas* pof2 = new TCanvas("pof2","pof2") ;
649	+	pof2->cd();
650	+	gPad->SetLeftMargin(0.15);
651	+	frameO->GetYaxis()->SetTitleOffset(1.4);
652	+	frameO->Draw();
653	+	if(is_data==data) DrawPrelim();
654	+	else DrawPrelim(PlottingSetup::luminosity,true);
655	+	stringstream OFFitParameterSummary;
656	+	OFFitParameterSummary << "#splitline{ftbbarOF = " <<  WriteWithError(fttbarOF.getVal(),fttbarOF.getError(),3)  << "}{";
657	+	OFFitParameterSummary << "#splitline{par1ttbarOF = " <<  WriteWithError(par1ttbarOF.getVal(),par1ttbarOF.getError(),3)  << "}{";
658	+	OFFitParameterSummary << "#splitline{par2ttbarOF = " <<  WriteWithError(par2ttbarOF.getVal(),par2ttbarOF.getError(),3)  << "}{";
659	+	OFFitParameterSummary << "#splitline{par3ttbarOF = " <<  WriteWithError(par3ttbarOF.getVal(),par3ttbarOF.getError(),3)  << "}{";
660	+	OFFitParameterSummary << "par4ttbarOF = " <<  WriteWithError(par4ttbarOF.getVal(),par4ttbarOF.getError(),3);
661	+	OFFitParameterSummary << "}}}}";
662	+	TLatex *OFFitinfobox = new TLatex(0.57,0.75,OFFitParameterSummary.str().c_str());
663	+	OFFitinfobox->SetNDC();
664	+	OFFitinfobox->SetTextSize(0.03);
665	+	OFFitinfobox->Draw();
666
667	<	RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarem");
668	<	float resultOFpar1 = resultOFpar1_->getVal();
669	<	//RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarem");
670	<	//float resultOFpar2 = resultOFpar2_->getVal();
671	<	//cout << "caca2.txt" << endl;
672	<
673	<	RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarem");
674	<	float resultOFpar3 = resultOFpar3_->getVal();
675	<
676	<	//RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarem");
677	<	//float resultOFpar4 = resultOFpar4_->getVal();
678	<	//cout << "caca4.txt" << endl;
679	<
680	<
681	<	// Now same flavor
682	<	RooRealVar fzee("fzee", "fzee", 5, 0, 100000);
683	<	RooRealVar meanzee("meanzee", "meanzee", 91.1876, 89, 95);
684	<	//RooRealVar sigmazee("sigmazee", "sigmazee", 0.5);
685	<	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0, 100);
686	<	RooRealVar widthzee("widthzee", "widthzee", 2.94);
687	<
688	<	RooRealVar fttbaree("fttbaree", "fttbaree", 100, 0, 100000);
689	<	RooRealVar par1ttbaree("par1ttbaree", "par1ttbaree", resultOFpar1, 0, 100);
690	<	RooRealVar par2ttbaree("par2ttbaree", "par2ttbaree", 1.0);
691	<	RooRealVar par3ttbaree("par3ttbaree", "par3ttbaree", resultOFpar3, 0, 100);
692	<	RooRealVar par4ttbaree("par4ttbaree", "par4ttbaree", 2.0);
667	>
668	>	CompleteSave(pof2,"Edge/OF__OFFitonly_"+prefix.str()+"__INFO",false,false);
669	>
670	>	delete pof2;
671	>
672	>	if(resultOF->covQual()!=3) {
673	>	write_error(__FUNCTION__,"OF fit did not converge!!! Cannot continue!");
674	>	dout << "covQual is " << resultOF->covQual() << endl;
675	>	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
676	>	else EdgeFitter::FixedMEdgeChi2_H0=-1;
677	>	if(EdgeFitter::RejectPointIfNoConvergence) return;
678	>	} else {
679	>	write_info(__FUNCTION__,"OF fit converged");
680	>	}
681	>
682	>
683	>	//    _____ _               ___
684	>	//   / ____| |             |__ \
685	>	//  | (___ | |_ ___ _ __      ) |
686	>	//   \___ \| __/ _ \ '_ \    / /
687	>	//   ____) | ||  __/ |_) |  / /_
688	>	//  |_____/ \__\___| .__/  |____|
689	>	//                 | |
690	>	//                 |_|
691	>	//
692	>	// Set up all the models
693	>
694	>
695	>	// Step 2a: set up edge position (if fixed), and maximum Z yield
696	>	float StartingMedge=70;
697	>	if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
698
699	<	RooRealVar fsignalee("fsignalee", "fsignalee", 10, 0, 400);
410	<	RooRealVar par1signalee("par1signalee", "par1signalee", 45, 20, 100);
411	<	RooRealVar par2signalee("par2signalee", "par2signalee", 2, 1, 10);
412	<	RooRealVar par3signalee("par3signalee", "par3signalee", 45, 0, 200);
413	<
414	<	RooVoigtian zee("zee", "zee", mll, meanzee, widthzee, sigmazee);
699	>	RooDataSet *ZDataSet = (RooDataSet*)EdgeFitter::AllData->reduce("id1==id2 && abs(mll-91.2)<20");
700
701	+	float maxZ = ZDataSet->sumEntries();
702	+	dout << "maxZ was set to " << maxZ << endl;
703	+	delete ZDataSet;
704
705	<	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbaree, par2ttbaree, par3ttbaree, par4ttbaree);
418	<	//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));
705	>	RooRealVar JustOne("JustOne","Just One",1) ;
706
707
708	<	RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
709	<	simPdf.addPdf(model_ee,"ee") ;
710	<	simPdf.addPdf(model_emu,"em") ;
429	<
430	<	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save());
431	<	result->Print();
708	>	// Step 2b: set up the models!
709	>	RooRealVar eefrac("eefrac","eefrac",eeSample->sumEntries()/(eeSample->sumEntries()+mmSample->sumEntries()));
710	>	RooFormulaVar mmfrac("mmfrac","1-eefrac",RooArgSet(eefrac));
711
712	<	RooPlot* frame1 = mll.frame(RooFit::Bins(25),RooFit::Title("EE sample")) ;
713	<	combData.plotOn(frame1,RooFit::Cut("sample==sample::ee")) ;
714	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
715	<	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;
453	<
454	<	// 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)) ;
712	>	RooRealVar rOneOverEMFrac("rOneOverEMFrac","rOneOverEMFrac",1.02,1.02-0.07,1.02+0.07);
713	>	RooRealVar rmean("rmean","rmean",1.02);
714	>	RooRealVar rwidth("rwidth","rwidth",0.07);
715	>	RooGaussian OneOverEMFrac("OneOverEMFrac","OneOverEMFrac",rOneOverEMFrac,rmean,rwidth);
716
460	–	stringstream prefix;
461	–	if(is_data==data) prefix << "data_";
462	–	if(is_data==mc) prefix << "mc_";
463	–	if(is_data==mcwithsignal) prefix << "mcwithS_";
717
718	<	prefix << "JZB_" << jzb_cut;
718	>	// FIRST: Drell-Yan: different parameters for resolution of DY (and of course yield), but same mean and so on
719	>	RooRealVar widthz("widthz", "widthz", 2.94);
720	>	widthz.setConstant(1);
721	>	RooRealVar meanz("meanz", "meanz", 91.1876, 89, 93);
722	>	RooRealVar fzSF("fzSF","fzSF",39,0,maxZ);
723	>	RooFormulaVar fzee("fzee","eefrac*fzSF",RooArgSet(eefrac,fzSF));
724	>	RooFormulaVar fzmm("fzmm","mmfrac*fzSF",RooArgSet(mmfrac,fzSF));
725	>	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0.5, 5);
726	>	RooRealVar sigmazmm("sigmazmm", "sigmazmm", 5, 0.5, 5);
727	>
728	>	RooVoigtian zee("zee", "zee", mll, meanz, widthz, sigmazee);
729	>	RooVoigtian zmm("zmm", "zmm", mll, meanz, widthz, sigmazmm);
730	>
731	>
732	>
733	>	//SECOND: Flavor-Symmetry contribution. Only relative fraction is different, all other parameters are the same;
734	>	RooFormulaVar fttbaree("fttbaree","eefrac*fttbarOF*OneOverEMFrac",RooArgSet(eefrac,fttbarOF,OneOverEMFrac));
735	>	RooFormulaVar fttbarmm("fttbarmm","mmfrac*fttbarOF*OneOverEMFrac",RooArgSet(mmfrac,fttbarOF,OneOverEMFrac));
736	>
737	>	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
738	>	RooSUSYBkgPdf ttbarmm("ttbarmm","ttbarmm", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
739	>	RooAddPdf ttbarSF("ttbarSF","ttbarSF", RooArgList(ttbaree,ttbarmm), RooArgList(fttbaree,fttbarmm));
740	>
741	>
742	>	//for signal only the resolution and absolute fraction are different for ee / mm
743	>	RooRealVar fsignal("fsignal", "fsignal", 0, 0, (eeSample->sumEntries()+mmSample->sumEntries()-OFSample->sumEntries()));
744	>	RooFormulaVar fsignalee("fsignalee","eefrac*fsignal",RooArgSet(eefrac,fsignal));
745	>	RooFormulaVar fsignalmm("fsignalmm","mmfrac*fsignal",RooArgSet(mmfrac,fsignal));
746	>
747	>	RooRealVar par1signal("par1signal", "par1signal", 45, 20, 100);
748	>	par1signal.setConstant(1); // doesn't have any effect.
749	>	RooRealVar par3signal("par3signal", "par3signal", StartingMedge, 0, 300); // this is the edge position
750	>	if(EdgeFitter::FixedMEdge>0) par3signal.setConstant();
751	>
752	>	RooSUSYTPdf signalee("signalee","signalee", mll , par1signal, sigmazee, par3signal);
753	>	RooSUSYTPdf signalmm("signalmm","signalmm", mll , par1signal, sigmazmm, par3signal);
754	>	RooAddPdf signalSF("signalSF","signalSF",RooArgList(signalee,signalmm),RooArgList(fsignalee,fsignalmm));
755	>
756	>	if(!EdgeFitter::AllowTriangle) {
757	>	fsignal.setVal(0.0); // kill off the signal if we don't want the triangle
758	>	fsignal.setConstant(1);
759	>	par1signal.setConstant(1);
760	>	par3signal.setConstant(1);
761	>	}
762	>
763
467	–	/*  cout << "fsignalee : " << fsignalee << endl;
468	–	cout << "fttbaree : " << fttbaree << endl;
469	–	cout << "fzee : " << fzee << endl;*/
764
765	+	//COMPLETE MODEL:
766	+	RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
767	+	RooAddPdf model_mm("model_mm","model_mm", RooArgList(zmm, ttbarmm, signalmm), RooArgList(fzmm, fttbarmm, fsignalmm));
768	+	RooAddPdf model_SF("model_SF","model_SF",RooArgList(model_ee,model_mm),RooArgList(JustOne,JustOne));
769	+
770	+
771	+	//    _____ _               ____
772	+	//   / ____| |             |___ \
773	+	//  | (___ | |_ ___ _ __     __) |
774	+	//   \___ \| __/ _ \ '_ \   |__ <
775	+	//   ____) | ||  __/ |_) |  ___) |
776	+	//  |_____/ \__\___| .__/  |____/
777	+	//                 | |
778	+	//                 |_|
779	+	//
780	+	//  Fitting
781
782	+	RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
783	+	simPdf.addPdf(model_ee,"ee") ;
784	+	simPdf.addPdf(model_mm,"mm") ;
785	+	simPdf.addPdf(model_OF,"OF") ;
786	+
787	+	//RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true),RooFit::ExternalConstraints(OneOverEMFrac));
788	+	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended(),RooFit::Minos(true));
789	+
790	+	if(result->covQual()!=3) {
791	+	write_error(__FUNCTION__,"Full fit did not converge!!! Cannot continue!");
792	+	dout << "covQual is " << result->covQual() << endl;
793	+	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=-1;
794	+	else EdgeFitter::FixedMEdgeChi2_H0=-1;
795	+	if(EdgeFitter::RejectPointIfNoConvergence) return;
796	+	} else {
797	+	write_info(__FUNCTION__,"Full fit converged");
798	+	}
799	+
800	+	dout << "============================= < Full results > =============================" << endl;
801	+	result->Print();
802	+	dout << "============================= < /Full results > =============================" << endl;
803	+
804	+	//    _____ _               _  _
805	+	//   / ____| |             | || |
806	+	//  | (___ | |_ ___ _ __   | || |_
807	+	//   \___ \| __/ _ \ '_ \  |__   _|
808	+	//   ____) | ||  __/ |_) |    | |
809	+	//  |_____/ \__\___| .__/     |_|
810	+	//                 | |
811	+	//                 |_|
812	+	//
813	+	//   Present results
814	+
815	+	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
816	+	frame1->GetXaxis()->CenterTitle(1);
817	+	frame1->GetYaxis()->CenterTitle(1);
818	+	combData.plotOn(frame1,RooFit::Name("eedata"),RooFit::Cut("sample==sample::ee")) ;
819	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("FullFitEE"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
820	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("TTbarEEonly"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
821	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Name("DYEEonly"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
822	+	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));
823	+
824	+
825	+	RooPlot* frame2 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("mm sample")) ;
826	+	frame2->GetXaxis()->CenterTitle(1);
827	+	frame2->GetYaxis()->CenterTitle(1);
828	+	combData.plotOn(frame2,RooFit::Name("mmdata"),RooFit::Cut("sample==sample::mm")) ;
829	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("FullFitMM"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
830	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("TTbarMMonly"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
831	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Name("DYMMonly"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
832	+	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));
833	+
834	+
835	+	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
836	+	frame3->GetXaxis()->CenterTitle(1);
837	+	frame3->GetYaxis()->CenterTitle(1);
838	+	combData.plotOn(frame3,RooFit::Name("OFdata"),RooFit::Cut("sample==sample::OF")) ;
839	+	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("FullFitOF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
840	+	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("TTbarOFonly"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
841	+	simPdf.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Name("DYOFonly"),RooFit::Components("zOF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
842	+
843	+	RooAbsPdf *ee_result = simPdf.getPdf("ee");
844	+	RooAbsPdf *mm_result = simPdf.getPdf("mm");
845	+	RooAddPdf SF_result("SF_result","SF_result",RooArgList(*ee_result,*mm_result));
846	+
847	+
848	+	TCanvas* g = new TCanvas("g","g") ;
849	+
850	+	RooPlot* frame4 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("SF sample")) ;
851	+	frame4->GetXaxis()->CenterTitle(1);
852	+	frame4->GetYaxis()->CenterTitle(1);
853	+	EdgeFitter::SFSample->plotOn(frame4,RooFit::Name("SFdata")) ;
854	+	SF_result.plotOn(frame4,RooFit::Name("FullFitSF"),RooFit::ProjWData(*EdgeFitter::SFSample),RooFit::LineColor(kBlack));
855	+	SF_result.plotOn(frame4,RooFit::Name("TTbarSFonly"),RooFit::ProjWData(*EdgeFitter::SFSample),RooFit::Components("ttbaree,ttbarmm"),RooFit::LineColor(kBlue),RooFit::LineStyle(kDashed));
856	+	SF_result.plotOn(frame4,RooFit::Name("DYSFonly"),RooFit::ProjWData(*EdgeFitter::SFSample),RooFit::Components("zmm,zee"), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
857	+	if(EdgeFitter::AllowTriangle) SF_result.plotOn(frame4,RooFit::Name("SignalSFonly"),RooFit::Components("signalee,signalmm"), RooFit::ProjWData(*EdgeFitter::SFSample), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
858	+
859	+	/// ********************************************************************************************************************************
860		TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
861		c->cd() ;
862		gPad->SetLeftMargin(0.15);
#	Line 476 | Line 864 | void EdgeFitter::DoFit(int is_data, floa
864		frame1->Draw();
865		if(is_data==data) DrawPrelim();
866		else DrawPrelim(PlottingSetup::luminosity,true);
867	<	CompleteSave(c,"Edge/"+prefix.str()+"eemm");
868	<	delete c;
867	>	CompleteSave(c,"Edge/"+prefix.str()+"_ee",false,false);
868	>	stringstream eeFitParameterSummary;
869	>	eeFitParameterSummary << "#splitline{t#bar{t} : " <<  WriteWithError((OneOverEMFrac.getVal()*fttbarOF.getVal()*eefrac.getVal()),(OneOverEMFrac.getVal()*fttbarOF.getVal()*eefrac.getVal())*sqrt(pow(fttbarOF.getError()/fttbarOF.getVal(),2)+pow(eefrac.getError()/eefrac.getVal(),2)),3)  << "}{";//no taking uncertainty on R(SF/OF) into consideration!
870	>	eeFitParameterSummary << "#splitline{DY    : " <<  WriteWithError((fzSF.getVal()*eefrac.getVal()),(fzSF.getVal()*eefrac.getVal())*sqrt(pow(fzSF.getError()/fzSF.getVal(),2)+pow(eefrac.getError()/eefrac.getVal(),2)),3)  << "}{";
871	>	eeFitParameterSummary << "#splitline{signal    : " <<  WriteWithError((fsignal.getVal()*eefrac.getVal()),(fsignal.getVal()*eefrac.getVal())*sqrt(pow(fsignal.getError()/fsignal.getVal(),2)+pow(eefrac.getError()/eefrac.getVal(),2)),3)  << "}{";
872	>	eeFitParameterSummary << "}}}";
873	>	TLatex *eeFitinfobox = new TLatex(0.57,0.75,eeFitParameterSummary.str().c_str());
874	>	eeFitinfobox->SetNDC();
875	>	eeFitinfobox->SetTextSize(0.03);
876	>	eeFitinfobox->Draw();
877	>	CompleteSave(c,"Edge/"+prefix.str()+"_ee__INFO",false,false);
878
879	<	TCanvas* d = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
483	<	d->cd() ;
879	>	c->cd() ;
880		gPad->SetLeftMargin(0.15);
881		frame2->GetYaxis()->SetTitleOffset(1.4);
882		frame2->Draw();
883		if(is_data==data) DrawPrelim();
884		else DrawPrelim(PlottingSetup::luminosity,true);
885	<	CompleteSave(d,"Edge/"+prefix.str()+"mm");
886	<	delete d;
887	<	//c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
885	>	CompleteSave(c,"Edge/"+prefix.str()+"_mm",false,false);
886	>	stringstream mmFitParameterSummary;
887	>	mmFitParameterSummary << "#splitline{t#bar{t} : " <<  WriteWithError((OneOverEMFrac.getVal()*fttbarOF.getVal()*(1-eefrac.getVal())),(OneOverEMFrac.getVal()*fttbarOF.getVal()*(1-eefrac.getVal()))*sqrt(pow(fttbarOF.getError()/fttbarOF.getVal(),2)+pow(eefrac.getError()/(1-eefrac.getVal()),2)),3)  << "}{";//no taking uncertainty on R(SF/OF) into consideration!
888	>	mmFitParameterSummary << "#splitline{DY    : " <<  WriteWithError((fzSF.getVal()*(1-eefrac.getVal())),(fzSF.getVal()*(1-eefrac.getVal()))*sqrt(pow(fzSF.getError()/fzSF.getVal(),2)+pow(eefrac.getError()/(1-eefrac.getVal()),2)),3)  << "}{";
889	>	mmFitParameterSummary << "#splitline{signal    : " <<  WriteWithError((fsignal.getVal()*(1-eefrac.getVal())),(fsignal.getVal()*(1-eefrac.getVal()))*sqrt(pow(fsignal.getError()/fsignal.getVal(),2)+pow(eefrac.getError()/(1-eefrac.getVal()),2)),3)  << "}{";
890	>	mmFitParameterSummary << "}}}";
891	>	TLatex *mmFitinfobox = new TLatex(0.57,0.75,mmFitParameterSummary.str().c_str());
892	>	mmFitinfobox->SetNDC();
893	>	mmFitinfobox->SetTextSize(0.03);
894	>	mmFitinfobox->Draw();
895	>	CompleteSave(c,"Edge/"+prefix.str()+"_mm__INFO",false,false);
896
897	<	TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
494	<	e->cd();
897	>	c->cd() ;
898		gPad->SetLeftMargin(0.15);
899		frame3->GetYaxis()->SetTitleOffset(1.4);
900		frame3->Draw();
901		if(is_data==data) DrawPrelim();
902		else DrawPrelim(PlottingSetup::luminosity,true);
903	<	CompleteSave(e,"Edge/"+prefix.str()+"emu");
904	<	delete e;
905	<
503	<	/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
504	<	f->cd();
903	>	CompleteSave(c,"Edge/"+prefix.str()+"_OF",false,false);
904	>
905	>	c->cd() ;
906		gPad->SetLeftMargin(0.15);
907		frame4->GetYaxis()->SetTitleOffset(1.4);
908		frame4->Draw();
909		if(is_data==data) DrawPrelim();
910		else DrawPrelim(PlottingSetup::luminosity,true);
911	<	CompleteSave(f,"Edge/"+prefix.str()+"eemm");
912	<	delete f;*/
913	<
911	>	CompleteSave(c,"Edge/"+prefix.str()+"_SF",false,false);
912	>	stringstream SFFitParameterSummary;
913	>	SFFitParameterSummary << "#splitline{t#bar{t} : " <<  WriteWithError((OneOverEMFrac.getVal()*fttbarOF.getVal()),(OneOverEMFrac.getVal()*fttbarOF.getVal())*sqrt(pow(fttbarOF.getError()/fttbarOF.getVal(),2)),3)  << "}{";//no taking uncertainty on R(SF/OF) into consideration!
914	>	SFFitParameterSummary << "#splitline{DY    : " <<  WriteWithError((fzSF.getVal()),(fzSF.getVal())*sqrt(pow(fzSF.getError()/fzSF.getVal(),2)),3)  << "}{";
915	>	SFFitParameterSummary << "#splitline{signal    : " <<  WriteWithError(fsignal.getVal(),(fsignal.getVal())*sqrt(pow(fsignal.getError()/fsignal.getVal(),2)),3)  << "}{";
916	>	SFFitParameterSummary << "}}}";
917	>	TLatex *SFFitinfobox = new TLatex(0.57,0.75,SFFitParameterSummary.str().c_str());
918	>	SFFitinfobox->SetNDC();
919	>	SFFitinfobox->SetTextSize(0.03);
920	>	SFFitinfobox->Draw();
921	>	CompleteSave(c,"Edge/"+prefix.str()+"_SF__INFO",false,false);
922	>
923	>	delete c;
924	>	/// ********************************************************************************************************************************
925	>
926	>	if(result->covQual()==3) {
927	>	if(EdgeFitter::AllowTriangle) EdgeFitter::FixedMEdgeChi2_H1=frame4->chiSquare("FullFitSF", "SFdata", 3);
928	>	else EdgeFitter::FixedMEdgeChi2_H0=frame4->chiSquare("FullFitSF", "SFdata", 3);
929	>	}
930	>
931		}
932
933		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 940 | void EdgeFitter::DoEdgeFit(string mcjzb,
940
941		EdgeFitter::PrepareDatasets(is_data);
942
943	+	EdgeFitter::DrawDatasetContent(is_data);
944	+
945		RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
946		RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
947	+
948	+
949	+	EdgeFitter::AllowTriangle=false;
950		EdgeFitter::DoFit(is_data, jzb_cut);
951	+
952	+	EdgeFitter::AllowTriangle=true;
953	+
954	+	bool ScanMassRange=true;
955	+	float ScanSteps=5.0;//GeV
956	+
957	+
958	+	if(ScanMassRange) {
959	+	TFile *fscan = new TFile("fscan.root","UPDATE");
960	+	TGraph *gr = new TGraph();
961	+	TGraph *Rgr = new TGraph();
962	+	stringstream GrName;
963	+	GrName << "ScanGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
964	+	stringstream RGrName;
965	+	RGrName << "ScanRatioGraphFor_" << EdgeFitter::Mode << "_" << jzb_cut;
966	+	gr->SetName(GrName.str().c_str());
967	+	Rgr->SetName(RGrName.str().c_str());
968	+
969	+	int i=0;
970	+	for(float tempMedge=10;tempMedge<=300;tempMedge+=ScanSteps) {
971	+	write_info(__FUNCTION__,"Now testing Medge="+any2string(tempMedge)+" for "+EdgeFitter::Mode+">"+any2string(jzb_cut));
972	+	EdgeFitter::FixedMEdge=tempMedge;
973	+	EdgeFitter::DoFit(is_data, jzb_cut);
974	+	if(EdgeFitter::FixedMEdgeChi2_H1<0) continue;
975	+	gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1);
976	+	Rgr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2_H1/EdgeFitter::FixedMEdgeChi2_H0);
977	+	i++;
978	+	}
979	+
980	+	TCanvas *ScanCan = new TCanvas("ScanCan","ScanCan",500,500);
981	+	gr->GetXaxis()->SetTitle("m_{edge}");
982	+	gr->GetXaxis()->CenterTitle();
983	+	gr->GetYaxis()->SetTitle("#chi^{2}");
984	+	gr->GetYaxis()->CenterTitle();
985	+	gr->GetYaxis()->SetTitleOffset(0.95);
986	+	gr->GetXaxis()->SetTitleOffset(0.9);
987	+	gr->SetLineColor(kBlue);
988	+	gr->SetTitle("");
989	+	gr->Draw("AL");
990	+	stringstream ScanCanSave;
991	+	ScanCanSave << "Edge/MEdgeScan_"+EdgeFitter::Mode+"_" << jzb_cut;
992	+	if(is_data) DrawPrelim();
993	+	else DrawMCPrelim();
994	+	CompleteSave(ScanCan,ScanCanSave.str());
995	+
996	+	Rgr->GetXaxis()->SetTitle("m_{edge}");
997	+	Rgr->GetXaxis()->CenterTitle();
998	+	Rgr->GetYaxis()->SetTitle("#frac{#chi^{2}(H_{1})}{#chi^{2}(H_{0})}");
999	+	Rgr->GetYaxis()->CenterTitle();
1000	+	Rgr->GetYaxis()->SetTitleOffset(0.95);
1001	+	Rgr->GetXaxis()->SetTitleOffset(0.9);
1002	+	Rgr->SetLineColor(kBlue);
1003	+	Rgr->SetTitle("");
1004	+	Rgr->Draw("AL");
1005	+	ScanCanSave.str("");
1006	+	ScanCanSave << "Edge/MEdgeScan_Ratio_"+EdgeFitter::Mode+"_" << jzb_cut;
1007	+	if(is_data) DrawPrelim();
1008	+	else DrawMCPrelim();
1009	+	CompleteSave(ScanCan,ScanCanSave.str());
1010	+	fscan->cd();
1011	+	gr->Write();
1012	+	delete ScanCan;
1013	+	fscan->Close();
1014	+	} else {
1015	+	EdgeFitter::DoFit(is_data, jzb_cut);
1016	+	dout << "Chi^2 (H0) = " << EdgeFitter::FixedMEdgeChi2_H0 << endl;
1017	+	dout << "Chi^2 (H1) = " << EdgeFitter::FixedMEdgeChi2_H1 << endl;
1018	+	}
1019	+
1020	+
1021		RooMsgService::instance().setGlobalKillBelow(msglevel);
1022
530	–
1023		f->Close();
1024
1025		}
1026
1027		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
1028	<	for(int icut=0;icut<jzb_cut.size();icut++) {
1028	>
1029	>	EdgeFitter::Mode="JZB";
1030	>	if(mcjzb=="met[4]") EdgeFitter::Mode="MET";
1031	>
1032	>	for(int icut=0;icut<(int)jzb_cut.size();icut++) {
1033		stringstream addcut;
1034		if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
1035		if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";

Diff Legend

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