ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C

Comparing UserCode/cbrown/AnalysisFramework/Plotting/Modules/Systematics.C (file contents):
Revision 1.42 by buchmann, Mon Oct 24 15:05:37 2011 UTC vs.
Revision 1.55 by buchmann, Thu Nov 24 08:15:39 2011 UTC

#	Line 85 | Line 85 | Double_t Interpolate(Double_t x, TH1 *hi
85
86		//____________________________________________________________________________________
87		// Plotting with all contributions, i.e. sidebands, peak, osof,ossf ... (for a systematic)
88	<	float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut) {
88	>	float allcontributionsplot(TTree* events, TCut kBaseCut, TCut kMassCut, TCut kSidebandCut, TCut JZBPosCut, TCut JZBNegCut, int flipped) {
89		iplot++;
90		int count=iplot;
91	+	string locmcjzbexpression=mcjzbexpression;
92		// Define new histogram
93		string hname=GetNumericHistoName();
94		TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
95	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
95	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
96		hname=GetNumericHistoName();
97		TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
98	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
98	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
99
100		hname=GetNumericHistoName();
101		TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
102	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
102	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
103		hname=GetNumericHistoName();
104		TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
105	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
105	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
106
107		float obs=0;
108		float pred=0;
#	Line 109 | Line 110 | float allcontributionsplot(TTree* events
110		if(PlottingSetup::RestrictToMassPeak) {
111		hname=GetNumericHistoName();
112		TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
113	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
113	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
114		hname=GetNumericHistoName();
115		TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
116	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
116	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
117
118		hname=GetNumericHistoName();
119		TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
120	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
120	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
121		hname=GetNumericHistoName();
122		TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
123	<	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
123	>	events->Draw("("+TString(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
124
125		obs = hossfp->Integral();
126		pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
127	+
128	+	if(flipped>0) {
129	+	obs = hossfn->Integral();
130	+	pred= hossfp->Integral() - (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
131	+	}
132		delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
133		} else {
134		obs = hossfp->Integral();
135		pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
136	+	if(flipped>0) {
137	+	obs = hossfn->Integral();
138	+	pred= hossfp->Integral() - (hosofp->Integral() - hosofn->Integral());;
139	+	}
140		}
141
142		delete hossfp,hossfn,hosofp,hosofn;
#	Line 136 | Line 146 | float allcontributionsplot(TTree* events
146
147		//____________________________________________________________________________________
148		// Efficiency plot
149	<	TH1F* plotEff(TTree* events, TCut kbase, TString informalname) {
149	>	TH1F* plotEff(TTree* events, TCut kbase, TString informalname, int flipped) {
150		iplot++;
151		int count=iplot;
152	+	iplot++;
153	+	int count2=iplot;
154		// Define new histogram
155		char hname[30]; sprintf(hname,"hJzbEff%d",count);
156	<	TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB (GeV/c); Efficiency",
157	<	nBins,jzbMin,jzbMax);
156	>	char hname2[30]; sprintf(hname2,"hJzbEff%d",count2);
157	>	TH1F* hJzbEff = new TH1F(hname,"JZB selection efficiency ; JZB [GeV]; Efficiency",nBins,jzbMin,jzbMax);
158	>	TH1F* hJzbEff2= new TH1F(hname2,"JZB selection efficiency ; JZB [GeV]; Efficiency",1,-14000,14000);
159		Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
160	<
161	<	events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
162	<	Float_t maxEff = events->GetSelectedRows();
160	>
161	>	if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
162	>	else events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),("genJZB>-400"&&kbase),"goff");
163	>	Float_t maxEff = hJzbEff2->Integral();
164		if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
165
166		if(verbose>0) dout <<  "JZB max = " << jzbMax << std::endl;
#	Line 154 | Line 168 | TH1F* plotEff(TTree* events, TCut kbase,
168		char cut[256];
169		for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
170		sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
171	<	events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
172	<	Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
173	<	//     dout << "COUCOU " << __LINE__ << std::endl;
171	>	if(flipped==0) events->Draw((mcjzbexpression+">>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
172	>	if(flipped>0) events->Draw(("(-"+mcjzbexpression+")>>"+(string)hname2).c_str(),(TCut(cut)&&kbase),"goff");
173	>	Float_t eff = static_cast<Float_t>(hJzbEff2->Integral())/maxEff;
174		hJzbEff->SetBinContent(iBin+1,eff);
175		hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
176		}
177	+	delete hJzbEff2;
178		return hJzbEff;
164	–
165	–
179		}
180
181
182	+
183		//________________________________________________________________________________________
184		// Master Formula
185		void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
#	Line 220 | Line 234 | float get_norm_pdf_factor(TTree *events,
234
235		//________________________________________________________________________________________
236		// Pile-up efficiency
237	<	float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
237	>	float pileup(TTree *events, bool requireZ, string informalname, int flipped, string addcut="",Float_t myJzbMax = 140. ) {
238		nBins = 16;
239		jzbMax = myJzbMax;
240
#	Line 230 | Line 244 | float pileup(TTree *events, bool require
244		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
245
246		if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
247	<	TH1F* hLM4 = plotEff(events,kbase,informalname);
247	>	TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
248		hLM4->SetMinimum(0.);
249
250		// Nominal function
#	Line 241 | Line 255 | float pileup(TTree *events, bool require
255
256		// Pimped-up function
257		TF1* funcUp = (TF1*)func->Clone();
258	<	funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
258	>	funcUp->SetParameter( 0, func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
259		if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
260		<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
261
#	Line 251 | Line 265 | float pileup(TTree *events, bool require
265
266		//____________________________________________________________________________________
267		// Effect of peak shifting
268	<	void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
268	>	void PeakError(TTree *events,float &result, string mcjzb, float peakerr,int flipped,string addcut="") {
269		//Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
270		TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
271		TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
#	Line 275 | Line 289 | void PeakError(TTree *events,float &resu
289		negcut=npeakup;
290		}
291		float res;
292	<	if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
293	<	else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
292	>	if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
293	>	else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
294		if(i==0) rescent=res;
295		else if(i==1) resdown=res;
296		else if(i==2) resup=res;
297		}
298		if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
299		else result=(TMath::Abs(rescent-resdown)/(float)rescent);
300	+	if(!automatized) cout << "    " << result << endl;
301	+	}
302	+
303	+
304	+	void  MCPartialefficiency(TTree *events,float &result, float &resulterr,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0, int type = 0) {
305	+	if(!events) {
306	+	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
307	+	result=0;resulterr=0;
308	+	return;
309	+	}
310	+
311	+	char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
312	+	char metSel[256]; sprintf(metSel, "met[4] > %f", jzbSel);
313	+	string metSelection(metSel);
314	+	// All acceptance cuts at gen. level
315	+	//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
316	+	TCut kbase("");
317	+	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
318	+	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
319	+	// Corresponding reco. cuts
320	+
321	+	TCut acceptance("genPt2 != 0");
322	+	TCut massId(cutmass&&cutOSSF);
323	+	TCut njets(cutnJets);
324	+	TCut jzbp;
325	+	TCut jzbn;
326	+	TCut met(("pfJetGoodNum > 1 && abs(mll-91.2) < 10.0 && id1 == id2 &&" + metSelection).c_str());
327	+	if(flipped==0) {
328	+	jzbp=TCut((TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
329	+	jzbn=TCut((TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
330	+	} else {
331	+	jzbp=TCut(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
332	+	jzbn=TCut(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
333	+	}
334	+	float ntotal = events->Draw("pt1", addcut.c_str(), "goff");
335	+	TCut theCut;
336	+	switch(type) {
337	+	case 1:
338	+	theCut = kbase+acceptance;
339	+	break;
340	+	case 2:
341	+	theCut = kbase+massId;
342	+	break;
343	+	case 3:
344	+	theCut = kbase+massId+njets;
345	+	break;
346	+	case 4:
347	+	theCut = kbase+massId+njets+jzbn;
348	+	break;
349	+	case 5:
350	+	theCut = kbase + met;
351	+	break;
352	+	default:
353	+	theCut = kbase+massId+njets+jzbn;
354	+	break;
355	+	}
356	+
357	+	string stheCut(theCut);
358	+	char var[20];
359	+	sprintf(var, "pdfW[%d]", k);
360	+
361	+	string svar(var);
362	+	string newtheCut;
363	+	if(k>0) newtheCut = "(" + stheCut + ")*" + svar;
364	+	else newtheCut = "(" + stheCut + ")"; // for k==0 or even k==-1 we don't need to evaluate PDFs
365	+
366	+	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
367	+	if(k>=0) events->Draw((mcjzbexpression+">>effh").c_str(), newtheCut.c_str(),"goff");
368	+	else events->Draw((mcjzbexpression+">>effh").c_str(), theCut,"goff");
369	+	Float_t sel = effh->Integral();
370	+	Float_t nsel=0;
371	+	//Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
372	+	float normFactor = 1;
373	+	if(k>=0) get_norm_pdf_factor(events, k, addcut);
374	+	sel = sel/normFactor;
375	+
376	+	result=(sel)/ntotal;
377	+	resulterr=TMath::Sqrt(sel/ntotal*(1+sel/ntotal)/ntotal);
378	+
379	+	delete effh;
380		}
381
382		//____________________________________________________________________________________
383		// Total selection efficiency (MC)
384		//returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
385	<	Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
292	<	write_warning(__FUNCTION__,"Setting automatized to off!"); automatized=false;
385	>	Value MCefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
386		if(!events) {
387		write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
388		result=0;
#	Line 314 | Line 407 | Value MCefficiency(TTree *events,float &
407		if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
408		ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
409		ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
410	+	if(flipped>0) {
411	+	ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
412	+	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
413	+	}
414		} else {
415		//for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
416		ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
417		ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
418	+	if(flipped>0) {
419	+	ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
420	+	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
421	+	}
422		}
423
424		TCut posSide = kbase&&ksel;
#	Line 334 | Line 435 | Value MCefficiency(TTree *events,float &
435		string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
436
437		TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
438	<	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
439	<	else events->Draw((mcjzbexpression+">>effh").c_str(), (sposSide+"&&(id1==id2)").c_str(),"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
438	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), TCut(newPosSide.c_str())*PlottingSetup::Weight,"goff");
439	>	else events->Draw((mcjzbexpression+">>effh").c_str(), TCut((sposSide+"&&(id1==id2)").c_str())*PlottingSetup::Weight,"goff");//the OSSF condition is added for the offpeak analysis, in onpeak case it's there already but doesn't change anything.
440	>
441		Float_t sel = effh->Integral();
442		Float_t nsel=0;
443
#	Line 344 | Line 446 | Value MCefficiency(TTree *events,float &
446		if(ConsiderSignalContaminationForLimits) {
447		flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
448		if(PlottingSetup::RestrictToMassPeak) {
449	<	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
449	>	events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*PlottingSetup::Weight,"goff");
450		nsel += effh->Integral();
451		} else {
452	<	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
452	>	events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(newNegSide.c_str()))*PlottingSetup::Weight,"goff");
453		nsel += effh->Integral();
454	<	events->Draw((mcjzbexpression+">>effh").c_str(), emnewPosSide.c_str(),"goff");
454	>	events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewPosSide.c_str()))*PlottingSetup::Weight,"goff");
455		nsel += effh->Integral();
456	<	events->Draw((mcjzbexpression+">>effh").c_str(), emnewNegSide.c_str(),"goff");
456	>	events->Draw((mcjzbexpression+">>effh").c_str(), (TCut(emnewNegSide.c_str()))*PlottingSetup::Weight,"goff");
457		nsel -= effh->Integral();
458		}
459		}
#	Line 362 | Line 464 | Value MCefficiency(TTree *events,float &
464		sel = sel/normFactor;
465		nsel = nsel/normFactor;
466
365	–	//      events->Draw(mcjzbexpression.c_str(),kbase,"goff");
366	–	//      Float_t tot = events->GetSelectedRows();
467		Float_t tot = Neventsinfile;
468
469		Value result_wo_signalcont;
#	Line 377 | Line 477 | Value MCefficiency(TTree *events,float &
477		resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
478		result_wo_signalcont=Value(result,resulterr);
479		}
480	<	if(!automatized && k>0 ) dout << "PDF assessment: ";
480	>	if(!automatized && k>0 ) dout << "PDF assessment [" << k << "] : ";
481		if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
482		delete effh;
483		return result_wo_signalcont;
#	Line 387 | Line 487 | Value MCefficiency(TTree *events,float &
487
488		//____________________________________________________________________________________
489		// Selection efficiency for one process (MC)
490	<	vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
490	>	// not in use anymore.
491	>	/*
492	>	vector<float> processMCefficiency(TTree *events,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
493		vector<float> process_efficiencies;
494		for(int iprocess=0;iprocess<=10;iprocess++) {
495		float this_process_efficiency,efferr;
496		stringstream addcutplus;
497		addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
498	<	MCefficiency(events,this_process_efficiency, efferr,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
498	>	MCefficiency(events,this_process_efficiency, efferr,flipped,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
499		process_efficiencies.push_back(this_process_efficiency);
500		}
501		return process_efficiencies;
502		}
503	<
503	>	*/
504
505	<	void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
505	>	void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, int flipped, bool requireZ, string addcut="") {
506		TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
507		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
508		if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
509	<	TH1F* hLM4 = plotEff(events,kbase,informalname);
509	>	TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
510		Int_t bin = hLM4->FindBin(jzbSel); // To get the error
511		jzbeff=Interpolate(jzbSel,hLM4);
512		jzbefferr=hLM4->GetBinError(bin);
#	Line 414 | Line 516 | void JZBefficiency(TTree *events, string
516
517		//________________________________________________________________________
518		// Effect of energy scale on efficiency
519	<	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
519	>	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname, int flipped, bool requireZ,string addcut="",Float_t jzbSelection=-1, TString plotName = "" ) {
520		TCut kbase(genMassCut&&"genZPt>0");
521		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
522		flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
#	Line 423 | Line 525 | void JZBjetScale(TTree *events, float &j
525		TCut ksel(cutmass&&cutOSSF);
526		TCut nJets("pfJetGoodNum>2");
527		stringstream down,up;
528	<	down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
529	<	up << "pfJetGoodNum"<<30*(1+syst)<<">=3";
528	>	down << "pfJetGoodNumn1sigma>=3";
529	>	up << "pfJetGoodNump1sigma>=3";
530
531		TCut nJetsP(up.str().c_str());
532		TCut nJetsM(down.str().c_str());
#	Line 432 | Line 534 | void JZBjetScale(TTree *events, float &j
534		if ( !(plotName.Length()>1) ) plotName = informalname;
535
536		nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
537	<	TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
538	<
539	<	TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
438	<
439	<	TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
537	>	TH1F* hist  = plotEff(events,(kbase&&ksel&&nJets),informalname,flipped);
538	>	TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname,flipped);
539	>	TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname,flipped);
540
541		// Dump some information
542		Float_t eff  = Interpolate(jzbSel,hist);
#	Line 452 | Line 552 | void JZBjetScale(TTree *events, float &j
552
553		//________________________________________________________________________
554		// Effect of energy scale on JZB efficiency
555	<	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
555	>	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, int flipped, bool requireZ, string addcut) {
556
557		TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
558		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
#	Line 464 | Line 564 | void doJZBscale(TTree *events, float &do
564		jzbMin =   0.5*jzbSel;
565		jzbMax =   2.0*jzbSel;
566
567	<	TH1F* hist = plotEff(events,kbase&&ksel,informalname);
567	>	TH1F* hist = plotEff(events,kbase&&ksel,informalname,flipped);
568
569		// Dump some information
570		Float_t eff  = Interpolate(jzbSel,hist);
#	Line 479 | Line 579 | void doJZBscale(TTree *events, float &do
579
580		//________________________________________________________________________
581		// JZB response (true/reco. vs. true)
582	<	void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
582	>	void JZBresponse(TTree *events, string name, bool requireZ, float &resp, float &resperr, int flipped, string addcut="", bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
583
584		jzbMin = 20;
585		flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
#	Line 492 | Line 592 | void JZBresponse(TTree *events, bool req
592
593		TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
594
595	<	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZB:genJZB",kbase&&ksel);
595	>	string locmcjzbexpression=mcjzbexpression;
596	>	if(flipped>0) locmcjzbexpression="-"+locmcjzbexpression;
597	>	string possibleminus="";
598	>	if(flipped>0) possibleminus="-";
599	>	if (!isMET) events->Project("hJzbResp","("+TString(locmcjzbexpression)+")/("+possibleminus+"genJZB):("+possibleminus+"genJZB)",kbase&&ksel);
600		else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
601
602		hJzbResp->SetMaximum(1.2);
#	Line 515 | Line 619 | void JZBresponse(TTree *events, bool req
619
620		//________________________________________________________________________________________
621		// PDF uncertainty
622	<	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
622	>	float get_pdf_uncertainty(TTree *events, int flipped, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
623		std::vector<float> efficiency;
624		for(int k = 1; k < NPdfs; k++) {
625		float result, resulterr;
626	<	MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
626	>	Value flipval;
627	>	MCefficiency(events, result, resulterr, flipped, mcjzb, requireZ, Neventsinfile, addcut, k);
628		efficiency.push_back(result);
629		}
630		float errHi, errLow,err;
#	Line 539 | Line 644 | int get_npdfs(TTree *events) {
644		}
645
646
647	<	void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
648	<	float JetEnergyScaleUncert=0.1;
544	<	float JZBScaleUncert=0.1;
647	>	void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,int flipped, string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
648	>	float JZBScaleUncert=0.05;
649		mcjzbexpression=mcjzb;
650	<	float triggereff=5.0/100;// in range [0,1]
650	>	float triggereff=2.0/100;// in range [0,1]
651		dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
652		float leptonseleff=2.0/100;// in range [0,1]
653		leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
#	Line 554 | Line 658 | void do_systematics_for_one_file(TTree *
658
659		float mceff,mcefferr,jzbeff,jzbefferr;
660		if(!automatized) dout << "MC efficiencies:" << endl;
661	<	Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
661	>	Value flipefficiency;
662	>	Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,addcut,-1);
663		if(!automatized) cout << "   Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
664
665	<	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
665	>	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,flipped,requireZ,addcut);
666		if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
667
668	<	if(!automatized) dout << "Error from Peak position:" << endl;
668	>	if(!automatized) dout << "Error from Peak position:";
669		float sysfrompeak=0;
670	<	PeakError(events,sysfrompeak,mcjzb,peakerror,addcut);
670	>	PeakError(events,sysfrompeak,mcjzb,peakerror,flipped,addcut);
671
672	<	if(!automatized) dout << "Jet energy scale: " << std::endl;
672	>	if(!automatized) dout << "Jet energy scale (JES): " << std::endl;
673		float jesup,jesdown;
674	<	JZBjetScale(events,jesdown,jesup,informalname,requireZ,addcut,JetEnergyScaleUncert);
675	<
674	>	JZBjetScale(events,jesdown,jesup,informalname,flipped,requireZ,addcut);
675	>
676		if(!automatized) dout << "JZB scale: " << std::endl;
677		float scaleup,scaledown,scalesyst;
678	<	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
678	>	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
679
680		if(!automatized) dout << "JZB response: " << std::endl;
681		float resp,resperr;
682		if(PlottingSetup::computeJZBresponse) {
683		if(!automatized) dout << "JZB response: " << std::endl;
684	<	JZBresponse(events,requireZ,resp,resperr,addcut);
684	>	if(!ismSUGRA) JZBresponse(events,informalname,requireZ,resp,resperr,flipped,addcut);
685		}
686
687		if(!automatized) dout << "Pileup: " << std::endl;
688	<	float resolution;
689	<	resolution=pileup(events,requireZ,informalname,addcut);
688	>	//  float resolution;
689	>	//resolution=pileup(events,requireZ,informalname,flipped,addcut);
690
691		float PDFuncert=0;
692		if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
693	<	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
693	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
694
695		dout << "_______________________________________________" << endl;
696		dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
#	Line 596 | Line 701 | void do_systematics_for_one_file(TTree *
701		dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
702		dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
703		dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
704	<	dout << "Resolution : " << resolution << endl; // in range [0,1]
704	>	//  dout << "Resolution : " << resolution << endl; // in range [0,1]
705		dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
706		if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
707		if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
#	Line 607 | Line 712 | void do_systematics_for_one_file(TTree *
712		toterr+=(leptonseleff)*(leptonseleff);
713		if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
714		if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
715	<	toterr+=(resolution*resolution);
715	>	//  toterr+=(resolution*resolution);
716		toterr+=(sysfrompeak*sysfrompeak);
717		if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
718		dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
#	Line 627 | Line 732 | void do_systematics_for_one_file(TTree *
732		res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
733		if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
734		if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
735	<	res.push_back(fabs(resolution));
735	>	//  res.push_back(fabs(resolution));
736	>	res.push_back(0.0);
737		res.push_back(mceff_nosigcont.getValue());
738		res.push_back(mceff_nosigcont.getError());
739		if(ismSUGRA) res.push_back(PDFuncert);
740		results.push_back(res);
741		}
742
743	<	vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
743	>	vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, int flipped, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
744		automatized=true;
745		vector< vector<float> > systematics;
746		for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
747	<	dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
748	<	for(int ibin=0;ibin<bins.size();ibin++) {
749	<	jzbSel=bins[ibin];
750	<	geqleq="geq";
751	<	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
752	<	}//end of bin loop
747	>	dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
748	>	for(int ibin=0;ibin<bins.size();ibin++) {
749	>	jzbSel=bins[ibin];
750	>	geqleq="geq";
751	>	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,flipped,mcjzb,datajzb,mcpeakerror,requireZ);
752	>	}//end of bin loop
753		}//end of signal loop
754		return systematics;
755		}

Diff Legend

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