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.8 by buchmann, Mon Jul 16 02:31:01 2012 UTC

#	Line 3 | Line 3
3		#include <RooRealVar.h>
4		#include <RooArgSet.h>
5		#include <RooDataSet.h>
6	<	#include <RooRealVar.h>
7	<	#include <RooArgSet.h>
8	<	#include <RooDataSet.h>
6	>	#include <RooMCStudy.h>
7	>	#include <RooCategory.h>
8	>
9	>	#include <RooPlot.h>
10	>	#include <RooSimultaneous.h>
11	>	#include <RooAddPdf.h>
12	>	#include <RooFitResult.h>
13	>	#include <RooVoigtian.h>
14	>	#include <RooMsgService.h>
15	>
16		#include <RooStats/ModelConfig.h>
17		#include "RooStats/ProfileLikelihoodCalculator.h"
18		#include "RooStats/LikelihoodInterval.h"
#	Line 15 | Line 22
22		#include "RooStats/HybridCalculatorOriginal.h"
23		#include "RooStats/HypoTestInverterOriginal.h"
24
25	+	//#include "ParametrizedEdge.C"
26	+	#include "EdgeModules/RooSUSYTPdf.cxx"
27	+	#include "EdgeModules/RooSUSYBkgPdf.cxx"
28	+
29
30		using namespace std;
31		using namespace PlottingSetup;
#	Line 40 | Line 51 | namespace EdgeFitter {
51
52		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree*);
53		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree*);
54	+	void getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp);
55		void InitializeVariables(float _mllmin, float _mllmax, float _jzbmax, TCut _cut);
56		void PrepareDatasets(int);
57	<	void DoFit();
57	>	void DoFit(int is_data, float jzb_cut);
58		string RandomStorageFile();
59		Yield Get_Z_estimate(float,int);
60		Yield Get_T_estimate(float,int);
61	<	float SetEdgeLimit(int , RooWorkspace *ws);
61	>	float calcExclusion(RooWorkspace *ws, RooDataSet data, bool calcExclusion);
62	>	vector<RooDataSet> generateToys(RooWorkspace *ws, int nToys);
63	>	void prepareLimits(RooWorkspace *ws, bool ComputeBands);
64	>	TGraph* prepareLM(float mass, float nEv);
65
66		float jzbmax;
67		float mllmin;
#	Line 58 | Line 73 | namespace EdgeFitter {
73		RooDataSet* mmSample;
74		RooDataSet* emSample;
75
76	<	bool MarcoDebug;
76	>	bool MarcoDebug=true;
77	>	}
78	>
79	>	TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
80	>	float massLM[1];
81	>	massLM[0]=mass;
82	>	float accEffLM[1];
83	>	accEffLM[0]=nEv/PlottingSetup::luminosity;
84	>	TGraph *lm = new TGraph(1, massLM, accEffLM);
85	>	lm->GetXaxis()->SetNoExponent(true);
86	>	lm->GetXaxis()->SetTitle("m_{cut} [GeV]");
87	>	lm->GetYaxis()->SetTitle("#sigma #times A [pb] 95% CL UL");
88	>	lm->GetXaxis()->SetLimits(0.,300.);
89	>	lm->GetYaxis()->SetRangeUser(0.,0.08);
90	>	lm->SetMarkerStyle(34);
91	>	lm->SetMarkerColor(kRed);
92	>	return lm;
93	>	}
94	>
95	>	vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
96	>	ws->var("nSig")->setVal(0.);
97	>	ws->var("nSig")->setConstant(true);
98	>	RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
99	>	vector<RooDataSet> theToys;
100	>
101	>	RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
102	>	mcEE.generate(nToys,44,true);
103	>	RooMCStudy mcMM(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"MM"));
104	>	mcMM.generate(nToys,44,true);
105	>	RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
106	>	mcOSOF.generate(nToys,44,true);
107	>
108	>	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
109	>	RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
110	>	RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
111	>	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
112	>	RooRealVar weight("weight","weight",0,1000,"");
113	>	RooArgSet observables(mll,jzb,id1,id2,weight);
114	>
115	>	for(int i=0;i<nToys;i++) {
116	>	RooDataSet* toyEE    = (RooDataSet*)mcEE.genData(i);
117	>	RooDataSet* toyMM    = (RooDataSet*)mcMM.genData(i);
118	>	RooDataSet* toyOSOF  = (RooDataSet*)mcOSOF.genData(i);
119	>	stringstream toyname;
120	>	toyname << "theToy_" << i;
121	>	write_warning(__FUNCTION__,"Problem while adding toys");
122	>	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));
123	>	theToys.push_back(toyData);
124	>	}
125	>	ws->var("nSig")->setVal(17.0);
126	>	ws->var("nSig")->setConstant(false);
127	>	return theToys;
128		}
129
130	<	float EdgeFitter::SetEdgeLimit(int is_data, RooWorkspace *ws) {
130	>	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
131	>	TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
132	>	vector<float> theLimits;
133	>	for(int itoy=0;itoy<theToys.size();itoy++) {
134	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
135	>	if(theLimit > 0 ) gauLimit->Fill(theLimit);
136	>	}
137	>	const Int_t nQ = 4;
138	>	Double_t yQ[nQ] = {0.,0.,0.,0.};
139	>	Double_t xQ[nQ] = {0.022750132,0.1586552539,0.84134474609999998,0.977249868};
140	>	gauLimit->GetQuantiles(nQ,yQ,xQ);
141	>	median = gauLimit->GetMean();
142	>	//  median = median1(gauLimit);
143	>	twoSigmaDown = abs(yQ[0]-median);
144	>	sigmaDown = abs(yQ[1]-median);
145	>	sigmaUp = abs(yQ[2]-median);
146	>	twoSigmaUp = abs(yQ[3]-median);
147	>	cout << median * PlottingSetup::luminosity << " " << sigmaUp * PlottingSetup::luminosity << endl;
148	>	}
149	>
150	>	void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
151	>	if(ComputeBands) {
152	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
153	>	vector<float> medVals;
154	>	vector<float> medLimits;
155	>	vector<float> sigmaLimitsDown;
156	>	vector<float> sigmaLimitsUp;
157	>	vector<float> twoSigmaLimitsDown;
158	>	vector<float> twoSigmaLimitsUp;
159	>	for(int i=20;i<=320;i+=40) {
160	>	ws->var("nSig")->setVal(10.0);
161	>	medVals.push_back((float)i);
162	>	ws->var("m0")->setVal((float)i);
163	>	ws->var("m0")->setConstant(true);
164	>	float Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp;
165	>	EdgeFitter::getMedianLimit(ws,theToys,Smedian,SsigmaDown,SsigmaUp,StwoSigmaDown,StwoSigmaUp);
166	>	medLimits.push_back(Smedian);
167	>	sigmaLimitsDown.push_back(SsigmaDown);
168	>	sigmaLimitsUp.push_back(SsigmaUp);
169	>	twoSigmaLimitsDown.push_back(StwoSigmaDown);
170	>	twoSigmaLimitsUp.push_back(StwoSigmaUp);
171	>	}
172	>	write_warning(__FUNCTION__,"Still need to store limits");
173	>	} else {
174	>	vector<float> theVals;
175	>	vector<float> theLimits;
176	>	for(int i=20;i<300;i+=5) {
177	>	ws->var("nSig")->setVal(0.0);
178	>	theVals.push_back((float)i);
179	>	ws->var("m0")->setVal((float)i);
180	>	ws->var("m0")->setConstant(true);
181	>	//      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
182	>	write_error(__FUNCTION__,"Error while casting roo data set");
183	>	}
184	>
185	>	for(int i=0;i<(int)theLimits.size();i++) {
186	>	if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
187	>	cout << i << " : " << theVals[i] << endl;
188	>	theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
189	>	write_warning(__FUNCTION__,"Need to store limits");
190	>	}
191	>	write_warning(__FUNCTION__,"Need to store limits");
192	>	}
193	>	}
194	>	}
195	>
196	>
197	>	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
198	>	int numberOfToys=50;
199		RooRealVar mu("mu","nSig",0,10000,"");
200		RooArgSet poi = RooArgSet(mu);
201		RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
#	Line 70 | Line 204 | float EdgeFitter::SetEdgeLimit(int is_da
204		model->SetWorkspace(*ws);
205		model->SetPdf("combModel");
206		model->SetParametersOfInterest(poi);
207	<	RooAbsData* data = ws->data("data_obs");
207	>	//  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
208
209	<	RooStats::ProfileLikelihoodCalculator plc(*data, *model);
209	>	RooStats::ProfileLikelihoodCalculator plc(data, *model);
210		plc.SetNullParameters(*nullParams);
211		plc.SetTestSize(0.05);
212	+
213		RooStats::LikelihoodInterval* interval = plc.GetInterval();
214		RooStats::HypoTestResult *htr = plc.GetHypoTest();
215		double theLimit = interval->UpperLimit( mu );
216	<	cout << "Significance " << htr->Significance() << endl;
216	>	//  double significance = htr->Significance();
217
218		ws->defineSet("obs","nB");
219		ws->defineSet("poi","nSig");
#	Line 102 | Line 237 | RooAbsData* data = ws->data("data_obs");
237		slrts.SetAltParameters(*sb_model.GetSnapshot());
238		RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
239
240	<	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
240	>	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
241		hc.SetTestStatistic(2);
242	<	hc.SetNumberOfToys(50);
108	<
242	>	hc.SetNumberOfToys(numberOfToys);
243
244		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
245		hcInv.SetTestSize(0.05);
#	Line 128 | Line 262 | TTree* SkimTree(int isample) {
262		cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
263		cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
264		}
265	+	float edgeWeight;
266	+	newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
267	+	float tmll;
268	+	allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
269	+	//  int id1,id2;
270	+
271		TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
272		float wgt=1.0;
273		allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
#	Line 135 | Line 275 | TTree* SkimTree(int isample) {
275		allsamples.collection[isample].events->LoadTree(entry);
276		if (select->EvalInstance()) {
277		allsamples.collection[isample].events->GetEntry(entry);
278	<	wgt=wgt*xsweight;
278	>	edgeWeight=wgt*xsweight;
279		newTree->Fill();
280		}
281		}
#	Line 153 | Line 293 | void EdgeFitter::InitializeVariables(flo
293
294		void EdgeFitter::PrepareDatasets(int is_data) {
295		TTree *completetree;
296	<	bool hashit=0;
297	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
296	>	write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
297	>	TList *treelist = new TList;
298	>	for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
299		if(!allsamples.collection[isample].is_active) continue;
300		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
301		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
302		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
303		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
304		if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
305	<	if(!hashit) {
165	<	hashit=true;
166	<	completetree = SkimTree(isample)->CloneTree();
167	<	} else {
168	<	completetree->CopyEntries(SkimTree(isample));
169	<	}
170	<	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
305	>	treelist->Add(SkimTree(isample));
306		}
307	+	completetree = TTree::MergeTrees(treelist);
308	+	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
309
310		RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
311		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
312		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
313	<	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
314	<	RooRealVar weight("weight","weight",0,1000,"");
315	<	RooArgSet observables(mll,jzb,id1,id2,weight);
313	>	//RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
314	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
315	>	RooArgSet observables(mll,id1,id2,edgeWeight);
316
317		string title="CMS Data";
318		if(is_data!=1) title="CMS SIMULATION";
319	<	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
319	>	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
320		completetree->Write();
321	<	//  delete completetree;
321	>	delete completetree;
322
323	<	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==0");
324	<	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2&&id1==1");
323	>	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2");
324	>	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2");
325		EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
326		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
327
#	Line 199 | Line 336 | void EdgeFitter::PrepareDatasets(int is_
336		cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
337		cout << "Number of events in em sample = " << emSample->numEntries() << endl;
338		}
339	+
340		}
341
342		string EdgeFitter::RandomStorageFile() {
#	Line 213 | Line 351 | string EdgeFitter::RandomStorageFile() {
351		Yield EdgeFitter::Get_Z_estimate(float jzb_cut, int icut) {
352		if(MarcoDebug) write_error(__FUNCTION__,"Returning random Z yield");
353		Yield a(123,45,67); return a;
216	–
354		return PlottingSetup::allresults.predictions[icut].Zbkg;
355		}
356
#	Line 223 | Line 360 | Yield EdgeFitter::Get_T_estimate(float j
360		return PlottingSetup::allresults.predictions[icut].Flavorsym;
361		}
362
363	+	void EdgeFitter::DoFit(int is_data, float jzb_cut) {
364	+	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
365	+	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
366	+	RooCategory sample("sample","sample") ;
367	+	sample.defineType("ee");
368	+	//sample.defineType("mm");
369	+	sample.defineType("em");
370	+
371	+	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("ee",*eeSample),Import("mm",*mmSample),Import("em",*emSample));
372	+	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("em",*emSample),RooFit::WeightVar(edgeWeight));
373	+
374	+
375	+	//First we make a fit to opposite flavor
376	+	RooRealVar fttbarem("fttbarem", "fttbarem", 100, 0, 10000);
377	+	RooRealVar par1ttbarem("par1ttbarem", "par1ttbarem", 1.6, 0.01, 4.0);
378	+	RooRealVar par2ttbarem("par2ttbarem", "par2ttbarem", 1.0);
379	+	RooRealVar par3ttbarem("par3ttbarem", "par3ttbarem", 0.028, 0.001, 1.0);
380	+	RooRealVar par4ttbarem("par4ttbarem", "par4ttbarem", 2.0);
381	+	RooSUSYBkgPdf ttbarem("ttbarem","ttbarem", mll , par1ttbarem, par2ttbarem, par3ttbarem, par4ttbarem);
382	+	RooAddPdf model_em("model_em","model_em", ttbarem, fttbarem);
383	+	RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
384	+	simPdfOF.addPdf(model_em,"em");
385	+	RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save());
386	+	resultOF->Print();
387	+
388	+	RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarem");
389	+	float resultOFpar1 = resultOFpar1_->getVal();
390	+	//RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarem");
391	+	//float resultOFpar2 = resultOFpar2_->getVal();
392	+	//cout << "caca2.txt" << endl;
393	+
394	+	RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarem");
395	+	float resultOFpar3 = resultOFpar3_->getVal();
396	+
397	+	//RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarem");
398	+	//float resultOFpar4 = resultOFpar4_->getVal();
399	+	//cout << "caca4.txt" << endl;
400	+
401	+
402	+	// Now same flavor
403	+	RooRealVar fzee("fzee", "fzee", 5, 0, 100000);
404	+	RooRealVar meanzee("meanzee", "meanzee", 91.1876, 89, 95);
405	+	//RooRealVar sigmazee("sigmazee", "sigmazee", 0.5);
406	+	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0, 100);
407	+	RooRealVar widthzee("widthzee", "widthzee", 2.94);
408	+
409	+	RooRealVar fttbaree("fttbaree", "fttbaree", 100, 0, 100000);
410	+	RooRealVar par1ttbaree("par1ttbaree", "par1ttbaree", resultOFpar1, 0, 100);
411	+	RooRealVar par2ttbaree("par2ttbaree", "par2ttbaree", 1.0);
412	+	RooRealVar par3ttbaree("par3ttbaree", "par3ttbaree", resultOFpar3, 0, 100);
413	+	RooRealVar par4ttbaree("par4ttbaree", "par4ttbaree", 2.0);
414	+
415	+	RooRealVar fsignalee("fsignalee", "fsignalee", 10, 0, 400);
416	+	RooRealVar par1signalee("par1signalee", "par1signalee", 45, 20, 100);
417	+	RooRealVar par2signalee("par2signalee", "par2signalee", 2, 1, 10);
418	+	RooRealVar par3signalee("par3signalee", "par3signalee", 45, 0, 200);
419	+
420	+	RooVoigtian zee("zee", "zee", mll, meanzee, widthzee, sigmazee);
421	+
422	+
423	+	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbaree, par2ttbaree, par3ttbaree, par4ttbaree);
424	+	//RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, par2signalee, par3signalee);
425	+	RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, sigmazee, par3signalee);
426	+
427	+	//RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
428	+	RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
429	+	RooAddPdf model_emu("model_emu","model_emu", RooArgList(ttbaree), RooArgList(fttbaree));
430	+
431	+
432	+	RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
433	+	simPdf.addPdf(model_ee,"ee") ;
434	+	simPdf.addPdf(model_emu,"em") ;
435	+
436	+	RooFitResult *result = simPdf.fitTo(combData, RooFit::Save());
437	+	result->Print();
438	+
439	+	RooPlot* frame1 = mll.frame(RooFit::Bins(25),RooFit::Title("EE sample")) ;
440	+	combData.plotOn(frame1,RooFit::Cut("sample==sample::ee")) ;
441	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
442	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
443	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
444	+	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("signalee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
445	+
446	+	RooPlot* frame2 = mll.frame(RooFit::Bins(25),RooFit::Title("MM sample")) ;
447	+	combData.plotOn(frame2,RooFit::Cut("sample==sample::mm")) ;
448	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
449	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
450	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
451	+	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("signalmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
452	+
453	+
454	+	cout << "Result   : " << endl;
455	+	cout << "f signal : " << fsignalee.getVal() << " +/- " << fsignalee.getError() << endl;
456	+	cout << "f ttbar  : " << fttbaree.getVal() << " +/- " << fttbaree.getError() << endl;
457	+	cout << "f tt em  : " << fttbarem.getVal() << " +/- " << fttbarem.getError() << endl;
458	+	cout << "f z ee   : " << fzee.getVal() << " +/- " << fzee.getError() << endl;
459	+
460	+	// The same plot for the cointrol sample slice
461	+	RooPlot* frame3 = mll.frame(RooFit::Bins(25),RooFit::Title("EM sample")) ;
462	+	combData.plotOn(frame3,RooFit::Cut("sample==sample::em")) ;
463	+	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
464	+	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::Components("ttbarem"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
465	+
466	+	stringstream prefix;
467	+	if(is_data==data) prefix << "data_";
468	+	if(is_data==mc) prefix << "mc_";
469	+	if(is_data==mcwithsignal) prefix << "mcwithS_";
470	+
471	+	prefix << "JZB_" << jzb_cut;
472	+
473	+
474	+
475	+	TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
476	+	c->cd() ;
477	+	gPad->SetLeftMargin(0.15);
478	+	frame1->GetYaxis()->SetTitleOffset(1.4);
479	+	frame1->Draw();
480	+	if(is_data==data) DrawPrelim();
481	+	else DrawPrelim(PlottingSetup::luminosity,true);
482	+	CompleteSave(c,"Edge/"+prefix.str()+"eemm");
483	+	delete c;
484	+
485	+	TCanvas* d = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
486	+	d->cd() ;
487	+	gPad->SetLeftMargin(0.15);
488	+	frame2->GetYaxis()->SetTitleOffset(1.4);
489	+	frame2->Draw();
490	+	if(is_data==data) DrawPrelim();
491	+	else DrawPrelim(PlottingSetup::luminosity,true);
492	+	CompleteSave(d,"Edge/"+prefix.str()+"mm");
493	+	delete d;
494	+	//c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
495	+
496	+	TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
497	+	e->cd();
498	+	gPad->SetLeftMargin(0.15);
499	+	frame3->GetYaxis()->SetTitleOffset(1.4);
500	+	frame3->Draw();
501	+	if(is_data==data) DrawPrelim();
502	+	else DrawPrelim(PlottingSetup::luminosity,true);
503	+	CompleteSave(e,"Edge/"+prefix.str()+"emu");
504	+	delete e;
505	+
506	+
507	+
508	+
509	+	/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
510	+	f->cd();
511	+	gPad->SetLeftMargin(0.15);
512	+	frame4->GetYaxis()->SetTitleOffset(1.4);
513	+	frame4->Draw();
514	+	if(is_data==data) DrawPrelim();
515	+	else DrawPrelim(PlottingSetup::luminosity,true);
516	+	CompleteSave(f,"Edge/"+prefix.str()+"eemm");
517	+	delete f;*/
518	+
519	+	//  RooWorkspace* wspace = new RooWorkspace();
520	+	//  stringstream mllvar;
521	+	//  mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax;
522	+	//  wspace->factory(mllvar.str().c_str());
523	+	//  wspace->var("mll")->setBins(30);
524	+	//  a lot of stuff missing here :-)
525	+	//  we need to store data_obs, weight, and so on in this space.
526	+	//  EdgeFitter::prepareLimits(wspace, true);
527	+	}
528	+
529		void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
530
531	<	string storagefile=EdgeFitter::RandomStorageFile();
532	<	TFile *f = new TFile(storagefile.c_str(),"RECREATE");
533	<	EdgeFitter::InitializeVariables(iMllLow,iMllHigh,jzbHigh,cut);
534	<
535	<	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;
531	>	TCut _cut(cut&&PlottingSetup::basiccut&&PlottingSetup::passtrig);
532	>
533	>	TFile *f = new TFile("workingfile.root","RECREATE");
534	>
535	>	EdgeFitter::InitializeVariables(PlottingSetup::iMllLow,PlottingSetup::iMllHigh,PlottingSetup::jzbHigh,_cut);
536
537		EdgeFitter::PrepareDatasets(is_data);
538	+
539	+	RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
540	+	RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
541	+	write_warning(__FUNCTION__,"Deactivated actual fitting procedure ATM");
542	+	EdgeFitter::DoFit(is_data, jzb_cut);
543	+	RooMsgService::instance().setGlobalKillBelow(msglevel);
544	+
545
240	–	EdgeFitter::eeSample->Write();
241	–	EdgeFitter::emSample->Write();
242	–	EdgeFitter::mmSample->Write();
243	–	EdgeFitter::AllData->Write();
546		f->Close();
547
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());
548		}
549
550		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
551	<	for(int icut=0;icut<jzb_cut.size();icut++) {
551	>	for(int icut=0;icut<(int)jzb_cut.size();icut++) {
552		stringstream addcut;
553		if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
554		if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";
555		TCut jcut(addcut.str().c_str());
556
557	+
558		EdgeFitter::DoEdgeFit(mcjzb, datajzb, DataPeakError, MCPeakError, jzb_cut[icut], icut, is_data, jcut&&cut, signalevents);
559
560		}

Diff Legend

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