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

Diff Legend

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