ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C (file contents):
Revision 1.3 by buchmann, Mon Jun 18 07:51:38 2012 UTC vs.
Revision 1.11 by buchmann, Wed Jun 12 05:43:20 2013 UTC

#	Line 1 | Line 1
1		#include <iostream>
2
3	+	#include <TVirtualIndex.h>
4	+
5		#include <RooRealVar.h>
6		#include <RooArgSet.h>
7		#include <RooDataSet.h>
8	<	#include <RooRealVar.h>
9	<	#include <RooArgSet.h>
10	<	#include <RooDataSet.h>
8	>	#include <RooMCStudy.h>
9	>	#include <RooCategory.h>
10	>
11	>	#include <RooPlot.h>
12	>	#include <RooSimultaneous.h>
13	>	#include <RooAddPdf.h>
14	>	#include <RooFitResult.h>
15	>	#include <RooVoigtian.h>
16	>	#include <RooMsgService.h>
17	>
18		#include <RooStats/ModelConfig.h>
19		#include "RooStats/ProfileLikelihoodCalculator.h"
20		#include "RooStats/LikelihoodInterval.h"
#	Line 15 | Line 24
24		#include "RooStats/HybridCalculatorOriginal.h"
25		#include "RooStats/HypoTestInverterOriginal.h"
26
27	+	//#include "ParametrizedEdge.C"
28	+	#include "EdgeModules/RooSUSYTPdf.cxx"
29	+	#include "EdgeModules/RooSUSYBkgPdf.cxx"
30	+
31
32		using namespace std;
33		using namespace PlottingSetup;
#	Line 40 | Line 53 | namespace EdgeFitter {
53
54		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree*);
55		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree*);
56	+	void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
57		void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
58		void PrepareDatasets(int);
59	<	void DoFit();
59	>	void DoFit(int is_data, float jzb_cut);
60		string RandomStorageFile();
61		Yield Get_Z_estimate(float,int);
62		Yield Get_T_estimate(float,int);
63	<	float SetEdgeLimit(int , RooWorkspace *ws);
63	>	float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
64	>	vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
65	>	void prepareLimits(RooWorkspace *ws, bool ComputeBands);
66	>	TGraph* prepareLM(float mass, float nEv);
67
68		float jzbmax;
69		float mllmin;
#	Line 54 | Line 71 | namespace EdgeFitter {
71		TCut cut;
72
73		RooDataSet* AllData;
74	<	RooDataSet* eeSample;
75	<	RooDataSet* mmSample;
76	<	RooDataSet* emSample;
74	>	RooDataSet* SFSample;
75	>	RooDataSet* OFSample;
76	>
77	>	bool MarcoDebug=true;
78	>
79	>	float FixedMEdge=-1;
80	>	float FixedMEdgeChi2=-1;
81
61	–	bool MarcoDebug;
82		}
83
84	<	float EdgeFitter::SetEdgeLimit(int is_data, RooWorkspace *ws) {
84	>	TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
85	>	float massLM[1];
86	>	massLM[0]=mass;
87	>	float accEffLM[1];
88	>	accEffLM[0]=nEv/PlottingSetup::luminosity;
89	>	TGraph *lm = new TGraph(1, massLM, accEffLM);
90	>	lm->GetXaxis()->SetNoExponent(true);
91	>	lm->GetXaxis()->SetTitle("m_{cut} [GeV]");
92	>	lm->GetYaxis()->SetTitle("#sigma #times A [pb] 95% CL UL");
93	>	lm->GetXaxis()->SetLimits(0.,300.);
94	>	lm->GetYaxis()->SetRangeUser(0.,0.08);
95	>	lm->SetMarkerStyle(34);
96	>	lm->SetMarkerColor(kRed);
97	>	return lm;
98	>	}
99	>
100	>	vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
101	>	ws->ls();
102	>	ws->var("nSig")->setVal(0.);
103	>	ws->var("nSig")->setConstant(true);
104	>	RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
105	>	vector<RooDataSet> theToys;
106	>
107	>	RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
108	>	mcEE.generate(nToys,44,true);
109	>	RooMCStudy mcMM(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"MM"));
110	>	mcMM.generate(nToys,44,true);
111	>	RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
112	>	mcOSOF.generate(nToys,44,true);
113	>
114	>	RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
115	>	RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
116	>	RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
117	>	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
118	>	RooRealVar weight("weight","weight",0,1000,"");
119	>	RooArgSet observables(mll,jzb,id1,id2,weight);
120	>
121	>	for(int i=0;i<nToys;i++) {
122	>	RooDataSet* toyEE    = (RooDataSet*)mcEE.genData(i);
123	>	RooDataSet* toyMM    = (RooDataSet*)mcMM.genData(i);
124	>	RooDataSet* toyOSOF  = (RooDataSet*)mcOSOF.genData(i);
125	>	stringstream toyname;
126	>	toyname << "theToy_" << i;
127	>	write_warning(__FUNCTION__,"Problem while adding toys");
128	>	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));
129	>	theToys.push_back(toyData);
130	>	}
131	>	ws->var("nSig")->setVal(17.0);
132	>	ws->var("nSig")->setConstant(false);
133	>	return theToys;
134	>	}
135	>
136	>	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
137	>	TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
138	>	vector<float> theLimits;
139	>	for(int itoy=0;itoy<(int)theToys.size();itoy++) {
140	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
141	>	if(theLimit > 0 ) gauLimit->Fill(theLimit);
142	>	}
143	>	const Int_t nQ = 4;
144	>	Double_t yQ[nQ] = {0.,0.,0.,0.};
145	>	Double_t xQ[nQ] = {0.022750132,0.1586552539,0.84134474609999998,0.977249868};
146	>	gauLimit->GetQuantiles(nQ,yQ,xQ);
147	>	median = gauLimit->GetMean();
148	>	//  median = median1(gauLimit);
149	>	twoSigmaDown = abs(yQ[0]-median);
150	>	sigmaDown = abs(yQ[1]-median);
151	>	sigmaUp = abs(yQ[2]-median);
152	>	twoSigmaUp = abs(yQ[3]-median);
153	>	cout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
154	>	}
155	>
156	>	void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
157	>	if(ComputeBands) {
158	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
159	>	vector<float> medVals;
160	>	vector<float> medLimits;
161	>	vector<float> sigmaLimitsDown;
162	>	vector<float> sigmaLimitsUp;
163	>	vector<float> twoSigmaLimitsDown;
164	>	vector<float> twoSigmaLimitsUp;
165	>	for(int i=20;i<=320;i+=40) {
166	>	ws->var("nSig")->setVal(10.0);
167	>	medVals.push_back((float)i);
168	>	ws->var("m0")->setVal((float)i);
169	>	ws->var("m0")->setConstant(true);
170	>	float Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp;
171	>	EdgeFitter::getMedianLimit(ws,theToys,Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp);
172	>	medLimits.push_back(Smedian);
173	>	sigmaLimitsDown.push_back(SsigmaDown);
174	>	sigmaLimitsUp.push_back(SsigmaUp);
175	>	twoSigmaLimitsDown.push_back(StwoSigmaDown);
176	>	twoSigmaLimitsUp.push_back(StwoSigmaUp);
177	>	}
178	>	write_warning(__FUNCTION__,"Still need to store limits");
179	>	} else {
180	>	vector<float> theVals;
181	>	vector<float> theLimits;
182	>	for(int i=20;i<300;i+=5) {
183	>	ws->var("nSig")->setVal(0.0);
184	>	theVals.push_back((float)i);
185	>	ws->var("m0")->setVal((float)i);
186	>	ws->var("m0")->setConstant(true);
187	>	//      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
188	>	write_error(__FUNCTION__,"Error while casting roo data set");
189	>	}
190	>
191	>	for(int i=0;i<(int)theLimits.size();i++) {
192	>	if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
193	>	cout << i << " : " << theVals[i] << endl;
194	>	theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
195	>	write_warning(__FUNCTION__,"Need to store limits");
196	>	}
197	>	write_warning(__FUNCTION__,"Need to store limits");
198	>	}
199	>	}
200	>	}
201	>
202	>
203	>	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
204	>	int numberOfToys=50;
205		RooRealVar mu("mu","nSig",0,10000,"");
206		RooArgSet poi = RooArgSet(mu);
207		RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
#	Line 70 | Line 210 | float EdgeFitter::SetEdgeLimit(int is_da
210		model->SetWorkspace(*ws);
211		model->SetPdf("combModel");
212		model->SetParametersOfInterest(poi);
213	<	RooAbsData* data = ws->data("data_obs");
213	>	//  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
214
215	<	RooStats::ProfileLikelihoodCalculator plc(*data, *model);
215	>	RooStats::ProfileLikelihoodCalculator plc(data, *model);
216		plc.SetNullParameters(*nullParams);
217		plc.SetTestSize(0.05);
218	+
219		RooStats::LikelihoodInterval* interval = plc.GetInterval();
220		RooStats::HypoTestResult *htr = plc.GetHypoTest();
221		double theLimit = interval->UpperLimit( mu );
222	<	cout << "Significance " << htr->Significance() << endl;
222	>	//  double significance = htr->Significance();
223
224		ws->defineSet("obs","nB");
225		ws->defineSet("poi","nSig");
#	Line 102 | Line 243 | RooAbsData* data = ws->data("data_obs");
243		slrts.SetAltParameters(*sb_model.GetSnapshot());
244		RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
245
246	<	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
246	>	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
247		hc.SetTestStatistic(2);
248	<	hc.SetNumberOfToys(50);
108	<
248	>	hc.SetNumberOfToys(numberOfToys);
249
250		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
251		hcInv.SetTestSize(0.05);
#	Line 128 | Line 268 | TTree* SkimTree(int isample) {
268		cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
269		cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
270		}
271	+	float edgeWeight;
272	+	newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
273	+	float tmll;
274	+	allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
275	+	//  int id1,id2;
276	+
277		TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
278	+	TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
279		float wgt=1.0;
280	<	allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
280	>	//  allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
281		for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
282		allsamples.collection[isample].events->LoadTree(entry);
283		if (select->EvalInstance()) {
284		allsamples.collection[isample].events->GetEntry(entry);
285	<	wgt=wgt*xsweight;
285	>	wgt=Weight->EvalInstance();
286	>	edgeWeight=wgt*xsweight;
287		newTree->Fill();
288		}
289		}
#	Line 150 | Line 298 | void EdgeFitter::InitializeVariables(flo
298		jzbmax=_jzbmax;
299		cut=_cut;
300		}
301	<
301	>
302	>	TTree* MergeTrees(vector<TTree*> trees) {
303	>	TTree * newtree = (TTree*)trees[0]->CloneTree();
304	>	trees[0]->GetListOfClones()->Remove(newtree);
305	>	trees[0]->ResetBranchAddresses();
306	>	newtree->ResetBranchAddresses();
307	>
308	>	for(int itree=1;itree<trees.size();itree++) {
309	>	newtree->CopyAddresses(trees[itree]);
310	>	Long64_t nentries = trees[itree]->GetEntries();
311	>	for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
312	>	trees[itree]->GetEntry(iEntry);
313	>	newtree->Fill();
314	>	}
315	>	trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
316	>	if (newtree->GetTreeIndex()) {
317	>	newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
318	>	}
319	>	if (newtree && newtree->GetTreeIndex()) {
320	>	newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
321	>	}
322	>	}
323	>	return newtree;
324	>	}
325	>
326	>
327	>
328		void EdgeFitter::PrepareDatasets(int is_data) {
329	<	TTree *completetree;
330	<	bool hashit=0;
331	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
329	>	write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
330	>	//  TFile *tempout = new TFile("tempout.root","RECREATE");
331	>	vector<TTree*> SkimmedTrees;
332	>	TTree *SkimmedTree[(int)allsamples.collection.size()];
333	>	for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
334		if(!allsamples.collection[isample].is_active) continue;
335		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
336		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
337		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
338		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
339		if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
340	<	if(!hashit) {
341	<	hashit=true;
342	<	completetree = SkimTree(isample)->CloneTree();
343	<	} else {
344	<	completetree->CopyEntries(SkimTree(isample));
345	<	}
346	<	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
340	>	SkimmedTrees.push_back(SkimTree(isample));
341	>	//      SkimmedTree[isample] = SkimTree(isample);
342	>	//    tempout->cd();
343	>	//    SkimmedTree[isample]->Write();
344	>	//    treelist->Add(SkimmedTree[isample]);
345	>	//treelist->Add(SkimTree(isample));
346	>	//    allsamples.collection[isample].tfile->Close();
347		}
348
349	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
349	>	TTree *completetree  = MergeTrees(SkimmedTrees);
350	>
351	>	//  for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
352	>	//    if(SkimmedTree[isample]) SkimmedTree[isample]->Delete();
353	>	//  }
354	>
355	>	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
356	>
357	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
358		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
359		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
360	<	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
361	<	RooRealVar weight("weight","weight",0,1000,"");
362	<	RooArgSet observables(mll,jzb,id1,id2,weight);
360	>	//RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
361	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
362	>	RooArgSet observables(mll,id1,id2,edgeWeight);
363
364		string title="CMS Data";
365		if(is_data!=1) title="CMS SIMULATION";
366	<	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
366	>	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
367		completetree->Write();
368	<	//  delete completetree;
368	>	delete completetree;
369	>	//  tempout->Close();
370
371	<	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==0");
372	<	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==1");
188	<	EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
371	>	EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
372	>	EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
373		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
374
375	<	eeSample->SetName("eeSample");
376	<	mmSample->SetName("mmSample");
193	<	emSample->SetName("emSample");
375	>	SFSample->SetName("SFSample");
376	>	OFSample->SetName("OFSample");
377		AllData->SetName("AllData");
378
379		if(EdgeFitter::MarcoDebug) {
380		cout << "Number of events in data sample = " << AllData->numEntries() << endl;
381	<	cout << "Number of events in ee sample = " << eeSample->numEntries() << endl;
382	<	cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
200	<	cout << "Number of events in em sample = " << emSample->numEntries() << endl;
381	>	cout << "Number of events in eemm sample = " << SFSample->numEntries() << endl;
382	>	cout << "Number of events in em sample = " << OFSample->numEntries() << endl;
383		}
384	+
385	+	}
386	+
387	+	string WriteWithError(float central, float error, int digits) {
388	+	float ref=central;
389	+	if(ref<0) ref=-central;
390	+	int HighestSigDigit = 0;
391	+	if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
392	+	else HighestSigDigit = int(log(ref)/(log(10)));
393	+
394	+	float divider=pow(10.0,(double(HighestSigDigit-digits)));
395	+
396	+	stringstream result;
397	+	result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
398	+	return result.str();
399		}
400
401	+
402		string EdgeFitter::RandomStorageFile() {
403		TRandom3 *r = new TRandom3(0);
404		int rho = (int)r->Uniform(1,10000000);
#	Line 213 | Line 411 | string EdgeFitter::RandomStorageFile() {
411		Yield EdgeFitter::Get_Z_estimate(float jzb_cut, int icut) {
412		if(MarcoDebug) write_error(__FUNCTION__,"Returning random Z yield");
413		Yield a(123,45,67); return a;
216	–
414		return PlottingSetup::allresults.predictions[icut].Zbkg;
415		}
416
#	Line 223 | Line 420 | Yield EdgeFitter::Get_T_estimate(float j
420		return PlottingSetup::allresults.predictions[icut].Flavorsym;
421		}
422
423	+	void EdgeFitter::DoFit(int is_data, float jzb_cut) {
424	+	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
425	+	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
426	+	RooCategory sample("sample","sample") ;
427	+	sample.defineType("SF");
428	+	//sample.defineType("mm");
429	+	sample.defineType("OF");
430	+
431	+	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
432	+	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
433	+
434	+
435	+	//First we make a fit to opposite flavor
436	+	RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 10000);
437	+	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
438	+	RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
439	+	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
440	+	RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
441	+	RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
442	+	RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
443	+	RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
444	+	simPdfOF.addPdf(model_OF,"OF");
445	+	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save(),RooFit::Extended());
446	+	resultOF->Print();
447	+
448	+	RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarOF");
449	+	float resultOFpar1 = resultOFpar1_->getVal();
450	+	//RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarOF");
451	+	//float resultOFpar2 = resultOFpar2_->getVal();
452	+	//cout << "caca2.txt" << endl;
453	+
454	+	RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarOF");
455	+	float resultOFpar3 = resultOFpar3_->getVal();
456	+
457	+	//RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarOF");
458	+	//float resultOFpar4 = resultOFpar4_->getVal();
459	+	//cout << "caca4.txt" << endl;
460	+
461	+	float StartingMedge=70;
462	+	if(EdgeFitter::FixedMEdge>0) StartingMedge=EdgeFitter::FixedMEdge;
463	+
464	+
465	+	// Now same flavor
466	+	RooRealVar fzSF("fzSF", "fzSF", 5, 0, 100000);
467	+	RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
468	+	//RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
469	+	RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
470	+	RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
471	+
472	+	RooRealVar fttbarSF("fttbarSF", "fttbarSF", 100, 0, 100000);
473	+	RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", resultOFpar1, 0, 100);
474	+	RooRealVar par2ttbarSF("par2ttbarSF", "par2ttbarSF", 1.0);
475	+	RooRealVar par3ttbarSF("par3ttbarSF", "par3ttbarSF", resultOFpar3, 0, 100);
476	+	RooRealVar par4ttbarSF("par4ttbarSF", "par4ttbarSF", 2.0);
477	+
478	+	RooRealVar fsignalSF("fsignalSF", "fsignalSF", 10, 0, 300);
479	+	RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
480	+	RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
481	+	RooRealVar par3signalSF("par3signalSF", "par3signalSF", StartingMedge, 0, 300);
482	+
483	+	RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
484	+
485	+	if(EdgeFitter::FixedMEdge>0) par3signalSF.setConstant();
486	+
487	+	RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarSF, par3ttbarSF, par4ttbarSF);
488	+	//RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, par2signalSF, par3signalSF);
489	+	RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
490	+
491	+	/*  par1ttbarSF.setConstant(true);
492	+	par2ttbarSF.setConstant(true);
493	+	par3ttbarSF.setConstant(true);
494	+	par4ttbarSF.setConstant(true);*/
495	+
496	+
497	+	//RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
498	+	RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
499	+	RooAddPdf model_em("model_em","model_em", RooArgList(ttbarSF), RooArgList(fttbarSF));
500	+
501	+
502	+	RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
503	+	simPdf.addPdf(model_SF,"SF") ;
504	+	simPdf.addPdf(model_em,"em") ;
505	+
506	+	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save(), RooFit::Extended());
507	+	result->Print();
508	+
509	+	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
510	+	frame1->GetXaxis()->CenterTitle(1);
511	+	combData.plotOn(frame1,RooFit::Name("SFdata"),RooFit::Cut("sample==sample::SF")) ;
512	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("FullFit"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
513	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("TTbarSFonly"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
514	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("DYSFonly"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
515	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Name("SignalSFonly"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
516	+
517	+	Double_t chi2 = frame1->chiSquare("FullFit", "SFdata", 3);
518	+
519	+
520	+	cout << "Result   : " << endl;
521	+	cout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
522	+	cout << "f ttbar  : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
523	+	cout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
524	+	cout << "f z SF   : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
525	+	cout << "Chi2     : " << chi2 << endl;
526	+
527	+	// The same plot for the cointrol sample slice
528	+	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
529	+	frame3->GetXaxis()->CenterTitle(1);
530	+	frame3->SetMaximum(frame1->GetMaximum());
531	+	combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
532	+	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
533	+	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
534	+
535	+
536	+	stringstream prefix;
537	+	if(is_data==data) prefix << "data_";
538	+	if(is_data==mc) prefix << "mc_";
539	+	if(is_data==mcwithsignal) prefix << "mcwithS_";
540	+
541	+	prefix << "JZB_" << jzb_cut;
542	+
543	+
544	+
545	+	TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
546	+	c->cd() ;
547	+	gPad->SetLeftMargin(0.15);
548	+	frame1->GetYaxis()->SetTitleOffset(1.4);
549	+	frame1->Draw();
550	+	if(is_data==data) DrawPrelim();
551	+	else DrawPrelim(PlottingSetup::luminosity,true);
552	+	stringstream infotext;
553	+	infotext << "#splitline{Fit results (JZB>" << jzb_cut << "): }{#splitline{";
554	+	infotext << "N(Data) = " << combData.sumEntries() << "}{#splitline{";
555	+	infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
556	+	infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
557	+	infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
558	+	infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
559	+
560	+	TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
561	+	infobox->SetNDC();
562	+	infobox->SetTextSize(0.03);
563	+	infobox->Draw();
564	+	CompleteSave(c,"Edge/"+prefix.str()+"_SF");
565	+	delete c;
566	+
567	+	TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
568	+	e->cd();
569	+	gPad->SetLeftMargin(0.15);
570	+	frame3->GetYaxis()->SetTitleOffset(1.4);
571	+	frame3->Draw();
572	+	if(is_data==data) DrawPrelim();
573	+	else DrawPrelim(PlottingSetup::luminosity,true);
574	+	CompleteSave(e,"Edge/"+prefix.str()+"_OF");
575	+	delete e;
576	+
577	+
578	+
579	+
580	+	/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
581	+	f->cd();
582	+	gPad->SetLeftMargin(0.15);
583	+	frame4->GetYaxis()->SetTitleOffset(1.4);
584	+	frame4->Draw();
585	+	if(is_data==data) DrawPrelim();
586	+	else DrawPrelim(PlottingSetup::luminosity,true);
587	+	CompleteSave(f,"Edge/"+prefix.str()+"_SF");
588	+	delete f;*/
589	+
590	+
591	+	/*
592	+	float maxZ=200;
593	+	RooWorkspace* wspace = new RooWorkspace();
594	+	stringstream mllvar;
595	+	mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
596	+	wspace->factory(mllvar.str().c_str());
597	+	wspace->var("mll")->setBins(30);
598	+	wspace->factory("nSig[1.,0.,100.]");
599	+	wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
600	+	wspace->factory("rME[1.12,1.05,1.19]");
601	+	wspace->factory("effUncert[1.]");
602	+	EdgeFitter::prepareLimits(wspace, true);
603	+	*/
604	+
605	+	write_warning(__FUNCTION__," A lot missing here to calculate limits");
606	+
607	+	}
608	+
609		void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
610
611	<	string storagefile=EdgeFitter::RandomStorageFile();
612	<	TFile *f = new TFile(storagefile.c_str(),"RECREATE");
613	<	EdgeFitter::InitializeVariables(iMllLow,iMllHigh,jzbHigh,cut);
614	<
615	<	Yield Zestimate=EdgeFitter::Get_Z_estimate(jzb_cut,icut);
233	<	Yield Testimate=EdgeFitter::Get_T_estimate(jzb_cut,icut);
234	<	cout << "Cut at JZB>" << jzb_cut << "; using estimates: " << endl;
235	<	cout << "       Z : " << Zestimate  << endl;
236	<	cout << "       T : " << Testimate  << endl;
611	>	TCut _cut(cut&&PlottingSetup::basiccut&&PlottingSetup::passtrig);
612	>
613	>	TFile *f = new TFile("workingfile.root","RECREATE");
614	>
615	>	EdgeFitter::InitializeVariables(PlottingSetup::iMllLow,PlottingSetup::iMllHigh,PlottingSetup::jzbHigh,_cut);
616
617		EdgeFitter::PrepareDatasets(is_data);
618	+
619	+	RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
620	+	RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
621	+	write_warning(__FUNCTION__,"Deactivated actual fitting procedure ATM");
622	+
623	+
624	+	bool ScanMassRange=false;
625	+
626	+	TGraph *gr = new TGraph();
627	+
628	+	if(ScanMassRange) {
629	+	int i=0;
630	+	for(float tempMedge=10;tempMedge<300;tempMedge+=0.1) {
631	+	EdgeFitter::FixedMEdge=tempMedge;
632	+	EdgeFitter::DoFit(is_data, jzb_cut);
633	+	gr->SetPoint(i,tempMedge,EdgeFitter::FixedMEdgeChi2);
634	+	}
635	+
636	+	TCanvas *ScanCan = new TCanvas("ScanCan","ScanCan",500,500);
637	+	gr->GetXaxis()->SetTitle("m_{edge}");
638	+	gr->GetXaxis()->CenterTitle();
639	+	gr->GetYaxis()->SetTitle("KS probability");
640	+	gr->GetYaxis()->CenterTitle();
641	+	gr->Draw("AP*");
642	+	stringstream ScanCanSave;
643	+	ScanCanSave << "Edge/MEdgeScan_JZB_" << jzb_cut;
644	+	ScanCan->SaveAs(ScanCanSave.str().c_str());
645	+	delete ScanCan;
646	+	} else {
647	+	EdgeFitter::DoFit(is_data, jzb_cut);
648	+	}
649	+
650	+
651	+
652	+
653	+
654	+	EdgeFitter::DoFit(is_data, jzb_cut);
655	+	RooMsgService::instance().setGlobalKillBelow(msglevel);
656	+
657
240	–	EdgeFitter::eeSample->Write();
241	–	EdgeFitter::emSample->Write();
242	–	EdgeFitter::mmSample->Write();
243	–	EdgeFitter::AllData->Write();
658		f->Close();
659
246	–	write_warning(__FUNCTION__,"Work continues here");
247	–
248	–	if(EdgeFitter::MarcoDebug) write_warning(__FUNCTION__,"Need to uncomment the next line (remove the output file)");
249	–	//  gSystem->Exec(("rm "+storagefile).c_str());
660		}
661
662		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
663	<	for(int icut=0;icut<jzb_cut.size();icut++) {
663	>	for(int icut=0;icut<(int)jzb_cut.size();icut++) {
664		stringstream addcut;
665		if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
666		if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";
667		TCut jcut(addcut.str().c_str());
668
669	+
670		EdgeFitter::DoEdgeFit(mcjzb, datajzb, DataPeakError, MCPeakError, jzb_cut[icut], icut, is_data, jcut&&cut, signalevents);
671
672		}

Diff Legend

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