ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/cbrown/Development/Plotting/Modules/EdgeLimit.C (file contents):
Revision 1.5 by buchmann, Wed Jun 27 08:53:23 2012 UTC vs.
Revision 1.10 by buchmann, Tue Jun 11 19:27:47 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>
#	Line 23 | Line 25
25		#include "RooStats/HypoTestInverterOriginal.h"
26
27		//#include "ParametrizedEdge.C"
28	<	#include "/shome/pablom/RooFit/Pdfs/RooSUSYTPdf.cxx"
29	<	#include "/shome/pablom/RooFit/Pdfs/RooSUSYBkgPdf.cxx"
28	>	#include "EdgeModules/RooSUSYTPdf.cxx"
29	>	#include "EdgeModules/RooSUSYBkgPdf.cxx"
30
31
32		using namespace std;
#	Line 51 | 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);
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(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 calcExclusion(RooWorkspace *ws, RooDataSet *data = NULL);
64	<	void prepareLimits(RooWorkspace *ws);
63	<	vector<RooDataSet*> generateToys(RooWorkspace *ws, int nToys);
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
#	Line 70 | Line 71 | namespace EdgeFitter {
71		TCut cut;
72
73		RooDataSet* AllData;
74	<	RooDataSet* eeSample;
75	<	RooDataSet* mmSample;
75	<	RooDataSet* emSample;
74	>	RooDataSet* SFSample;
75	>	RooDataSet* OFSample;
76
77	<	bool MarcoDebug;
77	>	bool MarcoDebug=true;
78		}
79
80		TGraph* EdgeFitter::prepareLM(float mass, float nEv) {
#	Line 93 | Line 93 | TGraph* EdgeFitter::prepareLM(float mass
93		return lm;
94		}
95
96	<	vector<RooDataSet*> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
96	>	vector<RooDataSet> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
97	>	ws->ls();
98		ws->var("nSig")->setVal(0.);
99		ws->var("nSig")->setConstant(true);
100		RooFitResult* fit = ws->pdf("combModel")->fitTo(*ws->data("data_obs"),RooFit::Save());
101	<	vector<RooDataSet*> theToys;
101	>	vector<RooDataSet> theToys;
102
103		RooMCStudy mcEE(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"EE"));
104		mcEE.generate(nToys,44,true);
#	Line 106 | Line 107 | vector<RooDataSet*> EdgeFitter::generate
107		RooMCStudy mcOSOF(*ws->pdf("combModel"),RooArgSet(*ws->var("inv")),RooFit::Slice(*ws->cat("cat"),"OSOF"));
108		mcOSOF.generate(nToys,44,true);
109
110	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
110	>	RooRealVar mll("m_{ll}","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
111		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
112		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
113		RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
#	Line 120 | Line 121 | vector<RooDataSet*> EdgeFitter::generate
121		stringstream toyname;
122		toyname << "theToy_" << i;
123		write_warning(__FUNCTION__,"Problem while adding toys");
124	<	//    RooDataSet *toyData = RooDataSet(toyname.str(),toyname.str(),observables,RooFit::Index(ws->cat("cat")),RooFit::Import("OSOF",*toyOSOF),RooFit::Import("EE",*toyEE),RooFit::Import("MM",*toyMM));
125	<	//    theToys.push_back(toyData);
124	>	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));
125	>	theToys.push_back(toyData);
126		}
127		ws->var("nSig")->setVal(17.0);
128		ws->var("nSig")->setConstant(false);
129		return theToys;
130		}
131
132	<	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
132	>	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
133		TH1F *gauLimit = new TH1F("gausLimit","gausLimit",60,0.,80./PlottingSetup::luminosity);
134		vector<float> theLimits;
135	<	for(int itoy=0;itoy<theToys.size();itoy++) {
136	<	float theLimit = calcExclusion(ws,theToys[itoy]);
135	>	for(int itoy=0;itoy<(int)theToys.size();itoy++) {
136	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
137		if(theLimit > 0 ) gauLimit->Fill(theLimit);
138		}
139		const Int_t nQ = 4;
#	Line 150 | Line 151 | void EdgeFitter::getMedianLimit(RooWorks
151
152		void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
153		if(ComputeBands) {
154	<	vector<RooDataSet*> theToys = EdgeFitter::generateToys(ws,50);
154	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
155		vector<float> medVals;
156		vector<float> medLimits;
157		vector<float> sigmaLimitsDown;
#	Line 179 | Line 180 | void EdgeFitter::prepareLimits(RooWorksp
180		theVals.push_back((float)i);
181		ws->var("m0")->setVal((float)i);
182		ws->var("m0")->setConstant(true);
183	<	theLimits.push_back(calcExclusion(ws));
183	>	//      theLimits.push_back(calcExclusion(ws,(RooDataSet)*ws->data("data_obs"),true));
184	>	write_error(__FUNCTION__,"Error while casting roo data set");
185		}
186
187	<	for(int i=0;i<theLimits.size();i++) {
187	>	for(int i=0;i<(int)theLimits.size();i++) {
188		if((theLimits[i]<2.0/PlottingSetup::luminosity)||(theLimits[i]>40.0/PlottingSetup::luminosity)) {
189		cout << i << " : " << theVals[i] << endl;
190		theLimits[i] = (theLimits[i+2]+theLimits[i-2])/2.0;
#	Line 194 | Line 196 | void EdgeFitter::prepareLimits(RooWorksp
196		}
197
198
199	<	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet *data) {
199	>	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet data, bool LoadDataObs) {
200	>	int numberOfToys=50;
201		RooRealVar mu("mu","nSig",0,10000,"");
202		RooArgSet poi = RooArgSet(mu);
203		RooArgSet *nullParams = (RooArgSet*)poi.snapshot();
#	Line 203 | Line 206 | float EdgeFitter::calcExclusion(RooWorks
206		model->SetWorkspace(*ws);
207		model->SetPdf("combModel");
208		model->SetParametersOfInterest(poi);
209	<	if(!data) data = (RooDataSet*)ws->data("data_obs");
209	>	//  if(LoadDataObs) data = (RooDataSet)*ws->data("data_obs");
210
211	<	RooStats::ProfileLikelihoodCalculator plc(*data, *model);
211	>	RooStats::ProfileLikelihoodCalculator plc(data, *model);
212		plc.SetNullParameters(*nullParams);
213		plc.SetTestSize(0.05);
214	+
215		RooStats::LikelihoodInterval* interval = plc.GetInterval();
216		RooStats::HypoTestResult *htr = plc.GetHypoTest();
217		double theLimit = interval->UpperLimit( mu );
218	<	cout << "Significance " << htr->Significance() << endl;
218	>	//  double significance = htr->Significance();
219
220		ws->defineSet("obs","nB");
221		ws->defineSet("poi","nSig");
#	Line 235 | Line 239 | float EdgeFitter::calcExclusion(RooWorks
239		slrts.SetAltParameters(*sb_model.GetSnapshot());
240		RooStats::ProfileLikelihoodTestStat profll = RooStats::ProfileLikelihoodTestStat(*b_model.GetPdf());
241
242	<	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(*data, sb_model, b_model,0,0);
242	>	RooStats::HybridCalculatorOriginal hc = RooStats::HybridCalculatorOriginal(data, sb_model, b_model,0,0);
243		hc.SetTestStatistic(2);
244	<	hc.SetNumberOfToys(50);
244	>	hc.SetNumberOfToys(numberOfToys);
245
246		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
247		hcInv.SetTestSize(0.05);
#	Line 260 | Line 264 | TTree* SkimTree(int isample) {
264		cout << "   Original tree contains " << allsamples.collection[isample].events->GetEntries() << endl;
265		cout << "   Going to reduce it with cut " << EdgeFitter::cut << endl;
266		}
267	+	float edgeWeight;
268	+	newTree->Branch("edgeWeight",&edgeWeight,"edgeWeight/F");
269	+	float tmll;
270	+	allsamples.collection[isample].events->SetBranchAddress("mll",&tmll);
271	+	//  int id1,id2;
272	+
273		TTreeFormula *select = new TTreeFormula("select", EdgeFitter::cut, allsamples.collection[isample].events);
274	+	TTreeFormula *Weight = new TTreeFormula("Weight", cutWeight, allsamples.collection[isample].events);
275		float wgt=1.0;
276	<	allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
276	>	//  allsamples.collection[isample].events->SetBranchAddress(cutWeight,&wgt);
277		for (Int_t entry = 0 ; entry < allsamples.collection[isample].events->GetEntries() ; entry++) {
278		allsamples.collection[isample].events->LoadTree(entry);
279		if (select->EvalInstance()) {
280		allsamples.collection[isample].events->GetEntry(entry);
281	<	wgt=wgt*xsweight;
281	>	wgt=Weight->EvalInstance();
282	>	edgeWeight=wgt*xsweight;
283		newTree->Fill();
284		}
285		}
#	Line 282 | Line 294 | void EdgeFitter::InitializeVariables(flo
294		jzbmax=_jzbmax;
295		cut=_cut;
296		}
297	<
297	>
298	>	TTree* MergeTrees(vector<TTree*> trees) {
299	>	TTree * newtree = (TTree*)trees[0]->CloneTree();
300	>	trees[0]->GetListOfClones()->Remove(newtree);
301	>	trees[0]->ResetBranchAddresses();
302	>	newtree->ResetBranchAddresses();
303	>
304	>	for(int itree=1;itree<trees.size();itree++) {
305	>	newtree->CopyAddresses(trees[itree]);
306	>	Long64_t nentries = trees[itree]->GetEntries();
307	>	for (Long64_t iEntry=0;iEntry<nentries;iEntry++) {
308	>	trees[itree]->GetEntry(iEntry);
309	>	newtree->Fill();
310	>	}
311	>	trees[itree]->ResetBranchAddresses(); // Disconnect from new tree
312	>	if (newtree->GetTreeIndex()) {
313	>	newtree->GetTreeIndex()->Append(trees[itree]->GetTreeIndex(),kTRUE);
314	>	}
315	>	if (newtree && newtree->GetTreeIndex()) {
316	>	newtree->GetTreeIndex()->Append(0,kFALSE); // Force the sorting
317	>	}
318	>	}
319	>	return newtree;
320	>	}
321	>
322	>
323	>
324		void EdgeFitter::PrepareDatasets(int is_data) {
287	–	TTree *completetree;
325		write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
326	<	bool hashit=0;
327	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
326	>	//  TFile *tempout = new TFile("tempout.root","RECREATE");
327	>	vector<TTree*> SkimmedTrees;
328	>	TTree *SkimmedTree[(int)allsamples.collection.size()];
329	>	for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
330		if(!allsamples.collection[isample].is_active) continue;
331		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
332		if(is_data==1&&allsamples.collection[isample].is_data==false) continue;//kick all samples that aren't data if we're looking for data.
333		if(is_data!=1&&allsamples.collection[isample].is_data==true) continue;//kick all data samples when looking for MC
334		if(is_data!=2&&allsamples.collection[isample].is_signal==true) continue;//remove signal sample if we don't want it.
335		if(EdgeFitter::MarcoDebug) cout << "Considering : " << allsamples.collection[isample].samplename << endl;
336	<	if(!hashit) {
337	<	hashit=true;
338	<	completetree = SkimTree(isample)->CloneTree();
339	<	} else {
340	<	completetree->CopyEntries(SkimTree(isample));
341	<	}
342	<	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
336	>	SkimmedTrees.push_back(SkimTree(isample));
337	>	//      SkimmedTree[isample] = SkimTree(isample);
338	>	//    tempout->cd();
339	>	//    SkimmedTree[isample]->Write();
340	>	//    treelist->Add(SkimmedTree[isample]);
341	>	//treelist->Add(SkimTree(isample));
342	>	//    allsamples.collection[isample].tfile->Close();
343		}
344
345	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
345	>	TTree *completetree  = MergeTrees(SkimmedTrees);
346	>
347	>	//  for(int isample=0;isample<(int)allsamples.collection.size();isample++) {
348	>	//    if(SkimmedTree[isample]) SkimmedTree[isample]->Delete();
349	>	//  }
350	>
351	>	if(EdgeFitter::MarcoDebug) cout << "Complete tree now contains " << completetree->GetEntries() << " entries " << endl;
352	>
353	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
354		RooRealVar id1("id1","id1",0,1,"GeV/c^{2}");
355		RooRealVar id2("id2","id2",0,1,"GeV/c^{2}");
356	<	RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
357	<	RooRealVar weight("weight","weight",0,1000,"");
358	<	RooArgSet observables(mll,jzb,id1,id2,weight);
356	>	//RooRealVar jzb("jzb","jzb",-jzbmax,jzbmax,"GeV/c");
357	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
358	>	RooArgSet observables(mll,id1,id2,edgeWeight);
359
360		string title="CMS Data";
361		if(is_data!=1) title="CMS SIMULATION";
362	<	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","weight");
362	>	RooDataSet LAllData("LAllData",title.c_str(),completetree,observables,"","edgeWeight");
363		completetree->Write();
364	<	//  delete completetree;
364	>	delete completetree;
365	>	//  tempout->Close();
366
367	<	EdgeFitter::eeSample = (RooDataSet*)LAllData.reduce("id1==id2");
368	<	EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2");
321	<	EdgeFitter::emSample = (RooDataSet*)LAllData.reduce("id1!=id2");
367	>	EdgeFitter::SFSample = (RooDataSet*)LAllData.reduce("id1==id2");
368	>	EdgeFitter::OFSample = (RooDataSet*)LAllData.reduce("id1!=id2");
369		EdgeFitter::AllData  = (RooDataSet*)LAllData.reduce("id1!=id2||id1==id2");
370
371	<	eeSample->SetName("eeSample");
372	<	mmSample->SetName("mmSample");
326	<	emSample->SetName("emSample");
371	>	SFSample->SetName("SFSample");
372	>	OFSample->SetName("OFSample");
373		AllData->SetName("AllData");
374
375		if(EdgeFitter::MarcoDebug) {
376		cout << "Number of events in data sample = " << AllData->numEntries() << endl;
377	<	cout << "Number of events in ee sample = " << eeSample->numEntries() << endl;
378	<	cout << "Number of events in mm sample = " << mmSample->numEntries() << endl;
333	<	cout << "Number of events in em sample = " << emSample->numEntries() << endl;
377	>	cout << "Number of events in eemm sample = " << SFSample->numEntries() << endl;
378	>	cout << "Number of events in em sample = " << OFSample->numEntries() << endl;
379		}
380	+
381	+	}
382	+
383	+	string WriteWithError(float central, float error, int digits) {
384	+	float ref=central;
385	+	if(ref<0) ref=-central;
386	+	int HighestSigDigit = 0;
387	+	if(ref>1) HighestSigDigit = int(log(ref)/log(10))+1;
388	+	else HighestSigDigit = int(log(ref)/(log(10)));
389	+
390	+	float divider=pow(10.0,(double(HighestSigDigit-digits)));
391	+
392	+	stringstream result;
393	+	result << divider*int(central/divider+0.5) << " #pm " << divider*int(error/divider+0.5);
394	+	return result.str();
395		}
396
397	+
398		string EdgeFitter::RandomStorageFile() {
399		TRandom3 *r = new TRandom3(0);
400		int rho = (int)r->Uniform(1,10000000);
#	Line 356 | Line 417 | Yield EdgeFitter::Get_T_estimate(float j
417		}
418
419		void EdgeFitter::DoFit(int is_data, float jzb_cut) {
420	<	RooRealVar mll("mll","mll",mllmin,mllmax,"GeV/c^{2}");
420	>	RooRealVar mll("mll","m_{ll}",mllmin,mllmax,"GeV/c^{2}");
421	>	RooRealVar edgeWeight("edgeWeight","edgeWeight",0,1000,"");
422		RooCategory sample("sample","sample") ;
423	<	sample.defineType("ee");
423	>	sample.defineType("SF");
424		//sample.defineType("mm");
425	<	sample.defineType("em");
426	<	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("ee",*eeSample),Import("mm",*mmSample),Import("em",*emSample));
427	<	RooDataSet combData("combData","combined data",mll,RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::Import("em",*emSample));
425	>	sample.defineType("OF");
426	>
427	>	//RooDataSet combData("combData","combined data",mll,Index(sample),Import("SF",*SFSample),Import("mm",*mmSample),Import("OF",*OFSample));
428	>	RooDataSet combData("combData","combined data",RooArgSet(mll,edgeWeight),RooFit::Index(sample),RooFit::Import("SF",*SFSample),RooFit::Import("OF",*OFSample),RooFit::WeightVar(edgeWeight));
429	>
430
367	–
368	–
431		//First we make a fit to opposite flavor
432	<	RooRealVar fttbarem("fttbarem", "fttbarem", 100, 0, 10000);
433	<	RooRealVar par1ttbarem("par1ttbarem", "par1ttbarem", 1.6, 0.01, 4.0);
434	<	RooRealVar par2ttbarem("par2ttbarem", "par2ttbarem", 1.0);
435	<	RooRealVar par3ttbarem("par3ttbarem", "par3ttbarem", 0.028, 0.001, 1.0);
436	<	RooRealVar par4ttbarem("par4ttbarem", "par4ttbarem", 2.0);
437	<	RooSUSYBkgPdf ttbarem("ttbarem","ttbarem", mll , par1ttbarem, par2ttbarem, par3ttbarem, par4ttbarem);
438	<	RooAddPdf model_em("model_em","model_em", ttbarem, fttbarem);
432	>	RooRealVar fttbarOF("fttbarOF", "fttbarOF", 100, 0, 10000);
433	>	RooRealVar par1ttbarOF("par1ttbarOF", "par1ttbarOF", 1.6, 0.01, 4.0);
434	>	RooRealVar par2ttbarOF("par2ttbarOF", "par2ttbarOF", 1.0);
435	>	RooRealVar par3ttbarOF("par3ttbarOF", "par3ttbarOF", 0.028, 0.001, 1.0);
436	>	RooRealVar par4ttbarOF("par4ttbarOF", "par4ttbarOF", 2.0);
437	>	RooSUSYBkgPdf ttbarOF("ttbarOF","ttbarOF", mll , par1ttbarOF, par2ttbarOF, par3ttbarOF, par4ttbarOF);
438	>	RooAddPdf model_OF("model_OF","model_OF", ttbarOF, fttbarOF);
439		RooSimultaneous simPdfOF("simPdfOF","simultaneous pdf", sample) ;
440	<	simPdfOF.addPdf(model_em,"em");
440	>	simPdfOF.addPdf(model_OF,"OF");
441		RooFitResult *resultOF = simPdfOF.fitTo(combData, RooFit::Save());
442		resultOF->Print();
443
444	<	RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarem");
444	>	RooRealVar* resultOFpar1_ = (RooRealVar*) resultOF->floatParsFinal().find("par1ttbarOF");
445		float resultOFpar1 = resultOFpar1_->getVal();
446	<	//RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarem");
446	>	//RooRealVar* resultOFpar2_ = (RooRealVar*) resultOF->floatParsFinal().find("par2ttbarOF");
447		//float resultOFpar2 = resultOFpar2_->getVal();
448		//cout << "caca2.txt" << endl;
449
450	<	RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarem");
450	>	RooRealVar* resultOFpar3_ = (RooRealVar*) resultOF->floatParsFinal().find("par3ttbarOF");
451		float resultOFpar3 = resultOFpar3_->getVal();
452
453	<	//RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarem");
453	>	//RooRealVar* resultOFpar4_ = (RooRealVar*) resultOF->floatParsFinal().find("par4ttbarOF");
454		//float resultOFpar4 = resultOFpar4_->getVal();
455		//cout << "caca4.txt" << endl;
456
457
458		// Now same flavor
459	<	RooRealVar fzee("fzee", "fzee", 5, 0, 100000);
460	<	RooRealVar meanzee("meanzee", "meanzee", 91.1876, 89, 95);
461	<	//RooRealVar sigmazee("sigmazee", "sigmazee", 0.5);
462	<	RooRealVar sigmazee("sigmazee", "sigmazee", 5, 0, 100);
463	<	RooRealVar widthzee("widthzee", "widthzee", 2.94);
464	<
465	<	RooRealVar fttbaree("fttbaree", "fttbaree", 100, 0, 100000);
466	<	RooRealVar par1ttbaree("par1ttbaree", "par1ttbaree", resultOFpar1, 0, 100);
467	<	RooRealVar par2ttbaree("par2ttbaree", "par2ttbaree", 1.0);
468	<	RooRealVar par3ttbaree("par3ttbaree", "par3ttbaree", resultOFpar3, 0, 100);
469	<	RooRealVar par4ttbaree("par4ttbaree", "par4ttbaree", 2.0);
459	>	RooRealVar fzSF("fzSF", "fzSF", 5, 0, 100000);
460	>	RooRealVar meanzSF("meanzSF", "meanzSF", 91.1876, 89, 95);
461	>	//RooRealVar sigmazSF("sigmazSF", "sigmazSF", 0.5);
462	>	RooRealVar sigmazSF("sigmazSF", "sigmazSF", 5, 0.5, 5);
463	>	RooRealVar widthzSF("widthzSF", "widthzSF", 2.94);
464	>
465	>	RooRealVar fttbarSF("fttbarSF", "fttbarSF", 100, 0, 100000);
466	>	RooRealVar par1ttbarSF("par1ttbarSF", "par1ttbarSF", resultOFpar1, 0, 100);
467	>	RooRealVar par2ttbarSF("par2ttbarSF", "par2ttbarSF", 1.0);
468	>	RooRealVar par3ttbarSF("par3ttbarSF", "par3ttbarSF", resultOFpar3, 0, 100);
469	>	RooRealVar par4ttbarSF("par4ttbarSF", "par4ttbarSF", 2.0);
470
471	<	RooRealVar fsignalee("fsignalee", "fsignalee", 10, 0, 400);
472	<	RooRealVar par1signalee("par1signalee", "par1signalee", 45, 20, 100);
473	<	RooRealVar par2signalee("par2signalee", "par2signalee", 2, 1, 10);
474	<	RooRealVar par3signalee("par3signalee", "par3signalee", 45, 0, 200);
471	>	RooRealVar fsignalSF("fsignalSF", "fsignalSF", 10, 0, 300);
472	>	RooRealVar par1signalSF("par1signalSF", "par1signalSF", 45, 20, 100);
473	>	RooRealVar par2signalSF("par2signalSF", "par2signalSF", 2, 1, 10);
474	>	RooRealVar par3signalSF("par3signalSF", "par3signalSF", 70, 0, 300);
475
476	<	RooVoigtian zee("zee", "zee", mll, meanzee, widthzee, sigmazee);
476	>	RooVoigtian zSF("zSF", "zSF", mll, meanzSF, widthzSF, sigmazSF);
477
478
479	<	RooSUSYBkgPdf ttbaree("ttbaree","ttbaree", mll , par1ttbaree, par2ttbaree, par3ttbaree, par4ttbaree);
480	<	//RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, par2signalee, par3signalee);
481	<	RooSUSYTPdf signalee("signalee","signalee", mll , par1signalee, sigmazee, par3signalee);
479	>	RooSUSYBkgPdf ttbarSF("ttbarSF","ttbarSF", mll , par1ttbarSF, par2ttbarSF, par3ttbarSF, par4ttbarSF);
480	>	//RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, par2signalSF, par3signalSF);
481	>	RooSUSYTPdf signalSF("signalSF","signalSF", mll , par1signalSF, sigmazSF, par3signalSF);
482	>
483	>	/*  par1ttbarSF.setConstant(true);
484	>	par2ttbarSF.setConstant(true);
485	>	par3ttbarSF.setConstant(true);
486	>	par4ttbarSF.setConstant(true);*/
487	>
488
489	<	//RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
490	<	RooAddPdf model_ee("model_ee","model_ee", RooArgList(zee, ttbaree, signalee), RooArgList(fzee, fttbaree, fsignalee));
491	<	RooAddPdf model_emu("model_emu","model_emu", RooArgList(ttbaree), RooArgList(fttbaree));
489	>	//RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
490	>	RooAddPdf model_SF("model_SF","model_SF", RooArgList(zSF, ttbarSF, signalSF), RooArgList(fzSF, fttbarSF, fsignalSF));
491	>	RooAddPdf model_em("model_em","model_em", RooArgList(ttbarSF), RooArgList(fttbarSF));
492
493
494		RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;
495	<	simPdf.addPdf(model_ee,"ee") ;
496	<	simPdf.addPdf(model_emu,"em") ;
495	>	simPdf.addPdf(model_SF,"SF") ;
496	>	simPdf.addPdf(model_em,"em") ;
497
498		RooFitResult *result = simPdf.fitTo(combData, RooFit::Save());
499		result->Print();
500
501	<	RooPlot* frame1 = mll.frame(RooFit::Bins(25),RooFit::Title("EE sample")) ;
502	<	combData.plotOn(frame1,RooFit::Cut("sample==sample::ee")) ;
503	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
504	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("ttbaree"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
505	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("zee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
506	<	simPdf.plotOn(frame1,RooFit::Slice(sample,"ee"),RooFit::Components("signalee"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
507	<
508	<	RooPlot* frame2 = mll.frame(RooFit::Bins(25),RooFit::Title("MM sample")) ;
441	<	combData.plotOn(frame2,RooFit::Cut("sample==sample::mm")) ;
442	<	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
443	<	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("ttbarmm"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
444	<	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("zmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
445	<	simPdf.plotOn(frame2,RooFit::Slice(sample,"mm"),RooFit::Components("signalmm"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
501	>	RooPlot* frame1 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("EE sample")) ;
502	>	frame1->GetXaxis()->CenterTitle(1);
503	>	combData.plotOn(frame1,RooFit::Cut("sample==sample::SF")) ;
504	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
505	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Components("ttbarSF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
506	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Components("zSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kRed));
507	>	simPdf.plotOn(frame1,RooFit::Slice(sample,"SF"),RooFit::Components("signalSF"), RooFit::ProjWData(sample, combData), RooFit::LineStyle(kDashed), RooFit::LineColor(kGreen));
508	>
509
510
511		cout << "Result   : " << endl;
512	<	cout << "f signal : " << fsignalee.getVal() << " +/- " << fsignalee.getError() << endl;
513	<	cout << "f ttbar  : " << fttbaree.getVal() << " +/- " << fttbaree.getError() << endl;
514	<	cout << "f tt em  : " << fttbarem.getVal() << " +/- " << fttbarem.getError() << endl;
515	<	cout << "f z ee   : " << fzee.getVal() << " +/- " << fzee.getError() << endl;
512	>	cout << "f signal : " << fsignalSF.getVal() << " +/- " << fsignalSF.getError() << endl;
513	>	cout << "f ttbar  : " << fttbarSF.getVal() << " +/- " << fttbarSF.getError() << endl;
514	>	cout << "f tt OF  : " << fttbarOF.getVal() << " +/- " << fttbarOF.getError() << endl;
515	>	cout << "f z SF   : " << fzSF.getVal() << " +/- " << fzSF.getError() << endl;
516
517		// The same plot for the cointrol sample slice
518	<	RooPlot* frame3 = mll.frame(RooFit::Bins(25),RooFit::Title("EM sample")) ;
519	<	combData.plotOn(frame3,RooFit::Cut("sample==sample::em")) ;
520	<	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
521	<	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"em"),RooFit::Components("ttbarem"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
518	>	RooPlot* frame3 = mll.frame(RooFit::Bins(int((mllmax-mllmin)/5.0)),RooFit::Title("OF sample")) ;
519	>	frame3->GetXaxis()->CenterTitle(1);
520	>	frame3->SetMaximum(frame1->GetMaximum());
521	>	combData.plotOn(frame3,RooFit::Cut("sample==sample::OF")) ;
522	>	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::ProjWData(sample,combData), RooFit::LineColor(kBlack)) ;
523	>	simPdfOF.plotOn(frame3,RooFit::Slice(sample,"OF"),RooFit::Components("ttbarOF"),RooFit::ProjWData(sample,combData),RooFit::LineStyle(kDashed)) ;
524	>
525
526		stringstream prefix;
527		if(is_data==data) prefix << "data_";
#	Line 464 | Line 530 | void EdgeFitter::DoFit(int is_data, floa
530
531		prefix << "JZB_" << jzb_cut;
532
467	–	/*  cout << "fsignalee : " << fsignalee << endl;
468	–	cout << "fttbaree : " << fttbaree << endl;
469	–	cout << "fzee : " << fzee << endl;*/
533
534
535		TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
#	Line 476 | Line 539 | void EdgeFitter::DoFit(int is_data, floa
539		frame1->Draw();
540		if(is_data==data) DrawPrelim();
541		else DrawPrelim(PlottingSetup::luminosity,true);
542	<	CompleteSave(c,"Edge/"+prefix.str()+"eemm");
542	>	stringstream infotext;
543	>	infotext << "#splitline{Fit results (JZB>" << jzb_cut << "): }{#splitline{";
544	>	infotext << "N(Data) = " << combData.sumEntries() << "}{#splitline{";
545	>	infotext << "N(Z+Jets) = " << WriteWithError(fzSF.getVal(),fzSF.getError(),3) << "}{#splitline{";
546	>	infotext << "N(t#bar{t}) = " << WriteWithError(fttbarSF.getVal(),fttbarSF.getError(),3) << "}{#splitline{";
547	>	infotext << "N(signal) = " << WriteWithError(fsignalSF.getVal(),fsignalSF.getError(),3) << "}{";
548	>	infotext << "m_{edge} = " << WriteWithError(par3signalSF.getVal(),par3signalSF.getError(),3) << "}}}}}";
549	>
550	>	TLatex *infobox = new TLatex(0.57,0.75,infotext.str().c_str());
551	>	infobox->SetNDC();
552	>	infobox->SetTextSize(0.03);
553	>	infobox->Draw();
554	>	CompleteSave(c,"Edge/"+prefix.str()+"_SF");
555		delete c;
556
482	–	TCanvas* d = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
483	–	d->cd() ;
484	–	gPad->SetLeftMargin(0.15);
485	–	frame2->GetYaxis()->SetTitleOffset(1.4);
486	–	frame2->Draw();
487	–	if(is_data==data) DrawPrelim();
488	–	else DrawPrelim(PlottingSetup::luminosity,true);
489	–	CompleteSave(d,"Edge/"+prefix.str()+"mm");
490	–	delete d;
491	–	//c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw();
492	–
557		TCanvas* e = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
558		e->cd();
559		gPad->SetLeftMargin(0.15);
#	Line 497 | Line 561 | void EdgeFitter::DoFit(int is_data, floa
561		frame3->Draw();
562		if(is_data==data) DrawPrelim();
563		else DrawPrelim(PlottingSetup::luminosity,true);
564	<	CompleteSave(e,"Edge/"+prefix.str()+"emu");
564	>	CompleteSave(e,"Edge/"+prefix.str()+"_OF");
565		delete e;
566
567	+
568	+
569	+
570		/*  TCanvas* f = new TCanvas("rf501_simultaneouspdfem","rf403_simultaneouspdfem") ;
571		f->cd();
572		gPad->SetLeftMargin(0.15);
#	Line 507 | Line 574 | void EdgeFitter::DoFit(int is_data, floa
574		frame4->Draw();
575		if(is_data==data) DrawPrelim();
576		else DrawPrelim(PlottingSetup::luminosity,true);
577	<	CompleteSave(f,"Edge/"+prefix.str()+"eemm");
577	>	CompleteSave(f,"Edge/"+prefix.str()+"_SF");
578		delete f;*/
579	<
579	>
580	>
581	>	/*
582	>	float maxZ=200;
583	>	RooWorkspace* wspace = new RooWorkspace();
584	>	stringstream mllvar;
585	>	mllvar << "mll[" << (mllmax-mllmin)/2 << "," << mllmin << "," << mllmax << "]";
586	>	wspace->factory(mllvar.str().c_str());
587	>	wspace->var("mll")->setBins(30);
588	>	wspace->factory("nSig[1.,0.,100.]");
589	>	wspace->factory(("nZ[0.04.,0.,"+any2string(maxZ)+"]").c_str());
590	>	wspace->factory("rME[1.12,1.05,1.19]");
591	>	wspace->factory("effUncert[1.]");
592	>	EdgeFitter::prepareLimits(wspace, true);
593	>	*/
594	>
595	>	write_warning(__FUNCTION__," A lot missing here to calculate limits");
596	>
597		}
598
599		void EdgeFitter::DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, float jzb_cut, int icut, int is_data, TCut cut, TTree *signalevents=0) {
#	Line 524 | Line 608 | void EdgeFitter::DoEdgeFit(string mcjzb,
608
609		RooFit::MsgLevel msglevel = RooMsgService::instance().globalKillBelow();
610		RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
611	+	write_warning(__FUNCTION__,"Deactivated actual fitting procedure ATM");
612		EdgeFitter::DoFit(is_data, jzb_cut);
613		RooMsgService::instance().setGlobalKillBelow(msglevel);
614
#	Line 533 | Line 618 | void EdgeFitter::DoEdgeFit(string mcjzb,
618		}
619
620		void DoEdgeFit(string mcjzb, string datajzb, float DataPeakError, float MCPeakError, vector<float> jzb_cut, int is_data, TCut cut, TTree *signalevents=0) {
621	<	for(int icut=0;icut<jzb_cut.size();icut++) {
621	>	for(int icut=0;icut<(int)jzb_cut.size();icut++) {
622		stringstream addcut;
623		if(is_data==1) addcut << "(" << datajzb << ">" << jzb_cut[icut] << ")";
624		if(is_data!=1) addcut << "(" << mcjzb << ">" << jzb_cut[icut] << ")";

Diff Legend

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