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

#	Line 23 | Line 23
23		#include "RooStats/HypoTestInverterOriginal.h"
24
25		//#include "ParametrizedEdge.C"
26	<	#include "/shome/pablom/RooFit/Pdfs/RooSUSYTPdf.cxx"
27	<	#include "/shome/pablom/RooFit/Pdfs/RooSUSYBkgPdf.cxx"
26	>	#include "EdgeModules/RooSUSYTPdf.cxx"
27	>	#include "EdgeModules/RooSUSYBkgPdf.cxx"
28
29
30		using namespace std;
#	Line 51 | 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);
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(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 calcExclusion(RooWorkspace *ws, RooDataSet *data = NULL);
62	<	void prepareLimits(RooWorkspace *ws);
63	<	vector<RooDataSet*> generateToys(RooWorkspace *ws, int nToys);
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
#	Line 74 | 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) {
#	Line 93 | Line 92 | TGraph* EdgeFitter::prepareLM(float mass
92		return lm;
93		}
94
95	<	vector<RooDataSet*> EdgeFitter::generateToys(RooWorkspace *ws, int nToys) {
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;
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);
#	Line 120 | Line 119 | vector<RooDataSet*> EdgeFitter::generate
119		stringstream toyname;
120		toyname << "theToy_" << i;
121		write_warning(__FUNCTION__,"Problem while adding toys");
122	<	//    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));
123	<	//    theToys.push_back(toyData);
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	<	void EdgeFitter::getMedianLimit(RooWorkspace *ws,vector<RooDataSet*> theToys,float &median,float &sigmaDown, float &sigmaUp, float &twoSigmaDown, float &twoSigmaUp) {
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]);
134	>	float theLimit = calcExclusion(ws,theToys[itoy],false);
135		if(theLimit > 0 ) gauLimit->Fill(theLimit);
136		}
137		const Int_t nQ = 4;
#	Line 150 | Line 149 | void EdgeFitter::getMedianLimit(RooWorks
149
150		void EdgeFitter::prepareLimits(RooWorkspace *ws, bool ComputeBands) {
151		if(ComputeBands) {
152	<	vector<RooDataSet*> theToys = EdgeFitter::generateToys(ws,50);
152	>	vector<RooDataSet> theToys = EdgeFitter::generateToys(ws,50);
153		vector<float> medVals;
154		vector<float> medLimits;
155		vector<float> sigmaLimitsDown;
#	Line 179 | Line 178 | void EdgeFitter::prepareLimits(RooWorksp
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));
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<theLimits.size();i++) {
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;
#	Line 194 | Line 194 | void EdgeFitter::prepareLimits(RooWorksp
194		}
195
196
197	<	float EdgeFitter::calcExclusion(RooWorkspace *ws, RooDataSet *data) {
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 203 | Line 204 | float EdgeFitter::calcExclusion(RooWorks
204		model->SetWorkspace(*ws);
205		model->SetPdf("combModel");
206		model->SetParametersOfInterest(poi);
207	<	if(!data) data = (RooDataSet*)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 235 | Line 237 | float EdgeFitter::calcExclusion(RooWorks
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);
242	>	hc.SetNumberOfToys(numberOfToys);
243
244		RooStats::HypoTestInverterOriginal hcInv =  RooStats::HypoTestInverterOriginal(hc,*ws->var("nSig"));
245		hcInv.SetTestSize(0.05);
#	Line 260 | 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 267 | 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 286 | Line 294 | void EdgeFitter::InitializeVariables(flo
294		void EdgeFitter::PrepareDatasets(int is_data) {
295		TTree *completetree;
296		write_warning(__FUNCTION__,"Need to make this function ready for scans as well (use signal from scan samples)");
297	<	bool hashit=0;
298	<	for(int isample=0;isample<allsamples.collection.size();isample++) {
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) {
298	<	hashit=true;
299	<	completetree = SkimTree(isample)->CloneTree();
300	<	} else {
301	<	completetree->CopyEntries(SkimTree(isample));
302	<	}
303	<	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");
324		EdgeFitter::mmSample = (RooDataSet*)LAllData.reduce("id1==id2");
#	Line 332 | 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 357 | Line 362 | Yield EdgeFitter::Get_T_estimate(float j
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",mll,RooFit::Index(sample),RooFit::Import("ee",*eeSample),RooFit::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
367	–
368	–
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);
#	Line 464 | Line 470 | void EdgeFitter::DoFit(int is_data, floa
470
471		prefix << "JZB_" << jzb_cut;
472
467	–	/*  cout << "fsignalee : " << fsignalee << endl;
468	–	cout << "fttbaree : " << fttbaree << endl;
469	–	cout << "fzee : " << fzee << endl;*/
473
474
475		TCanvas* c = new TCanvas("rf501_simultaneouspdf","rf403_simultaneouspdf") ;
#	Line 500 | Line 503 | void EdgeFitter::DoFit(int is_data, floa
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);
#	Line 509 | Line 515 | void EdgeFitter::DoFit(int is_data, floa
515		else DrawPrelim(PlottingSetup::luminosity,true);
516		CompleteSave(f,"Edge/"+prefix.str()+"eemm");
517		delete f;*/
518	<
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) {
#	Line 524 | Line 538 | void EdgeFitter::DoEdgeFit(string mcjzb,
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
#	Line 533 | Line 548 | void EdgeFitter::DoEdgeFit(string mcjzb,
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] << ")";

Diff Legend

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