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.5 by buchmann, Wed Jul 20 08:51:33 2011 UTC vs.
Revision 1.33 by buchmann, Mon Sep 26 15:09:41 2011 UTC

#	Line 1 | Line 1
1		#include <iostream>
2		#include <vector>
3		#include <sys/stat.h>
4	+	#include <algorithm>
5	+	#include <cmath>
6
7		#include <TMath.h>
8		#include <TColor.h>
#	Line 44 | Line 46 | TString geq_or_leq() {
46		return TString("GEQ_OR_LEQ_ERROR");
47		}
48
49	+	TString ngeq_or_leq() {
50	+	if(geqleq=="geq") return TString("<=");
51	+	if(geqleq=="leq") return TString(">=");
52	+	return TString("NGEQ_OR_LEQ_ERROR");
53	+	}
54	+
55		//______________________________________________________________________________
56		Double_t Interpolate(Double_t x, TH1 *histo)
57		{
#	Line 75 | Line 83 | Double_t Interpolate(Double_t x, TH1 *hi
83		}
84		}
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) {
89	+	iplot++;
90	+	int count=iplot;
91	+	// Define new histogram
92	+	string hname=GetNumericHistoName();
93	+	TH1F* hossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
94	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSSF,"goff");
95	+	hname=GetNumericHistoName();
96	+	TH1F* hossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
97	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSSF,"goff");
98	+
99	+	hname=GetNumericHistoName();
100	+	TH1F* hosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
101	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBPosCut&&cutOSOF,"goff");
102	+	hname=GetNumericHistoName();
103	+	TH1F* hosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
104	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
105	+
106	+	hname=GetNumericHistoName();
107	+	TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
108	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
109	+	hname=GetNumericHistoName();
110	+	TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
111	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
112	+
113	+	hname=GetNumericHistoName();
114	+	TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
115	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
116	+	hname=GetNumericHistoName();
117	+	TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
118	+	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
119	+
120	+	float obs = hossfp->Integral();
121	+	float pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
122	+
123	+	delete hossfp,hossfn,hosofp,hosofn;
124	+	delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
125	+	return obs-pred;
126	+	}
127	+
128
129		//____________________________________________________________________________________
130		// Efficiency plot
#	Line 87 | Line 137 | TH1F* plotEff(TTree* events, TCut kbase,
137		nBins,jzbMin,jzbMax);
138		Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
139
140	<	events->Draw(mcjzbexpression.c_str(),"genJZBSel>-400"&&kbase,"goff");
140	>	events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
141		Float_t maxEff = events->GetSelectedRows();
142		if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
143
#	Line 95 | Line 145 | TH1F* plotEff(TTree* events, TCut kbase,
145		// Loop over steps to get efficiency curve
146		char cut[256];
147		for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
148	<	sprintf(cut,"genJZBSel>%3f",jzbMin+iBin*step);
148	>	sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
149		events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
150		Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
151		//     dout << "COUCOU " << __LINE__ << std::endl;
#	Line 109 | Line 159 | TH1F* plotEff(TTree* events, TCut kbase,
159
160
161		//________________________________________________________________________________________
162	+	// Master Formula
163	+	void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
164	+
165	+	float x0 = eff[0];
166	+	float deltaPos = 0, deltaNeg = 0;
167	+	for(int k = 0; k < (eff.size()-1)/2; k++) {
168	+	float xneg = eff[2*k+2];
169	+	float xpos = eff[2*k+1];
170	+	if(xpos-x0>0 || xneg-x0>0) {
171	+	if(xpos-x0 > xneg-x0) {
172	+	deltaPos += (xpos-x0)*(xpos-x0);
173	+	} else {
174	+	deltaPos += (xneg-x0)*(xneg-x0);
175	+	}
176	+	}
177	+	if(x0-xpos>0 || x0-xneg>0) {
178	+	if(x0-xpos > x0-xneg) {
179	+	deltaNeg += (xpos-x0)*(xpos-x0);
180	+	} else {
181	+	deltaNeg += (xneg-x0)*(xneg-x0);
182	+	}
183	+	}
184	+	}
185	+	errHi = sqrt(deltaPos);
186	+	errLo = sqrt(deltaNeg);
187	+
188	+	}
189	+
190	+
191	+	//________________________________________________________________________________________
192	+	// Get normalization factor for the PDFs
193	+	float get_norm_pdf_factor(TTree *events, int k) {
194	+
195	+	TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
196	+	char nameVar[20];
197	+	sprintf(nameVar, "pdfW[%d]", k);
198	+	events->Project("haux", nameVar);
199	+	float thisW = haux->Integral();
200	+	events->Project("haux", "pdfW[0]");
201	+	float normW = haux->Integral();
202	+
203	+	float factor=thisW/normW;
204	+
205	+	delete haux;
206	+
207	+	return factor;
208	+
209	+	}
210	+
211	+
212	+
213	+	//________________________________________________________________________________________
214		// Pile-up efficiency
215	<	float pileup(TTree *events, string informalname, Float_t myJzbMax = 140. ) {
215	>	float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
216		nBins = 16;
217		jzbMax = myJzbMax;
218
219		// Acceptance cuts
220	<	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
220	>	TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
221	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
222
223	+	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
224		TH1F* hLM4 = plotEff(events,kbase,informalname);
225		hLM4->SetMinimum(0.);
226
#	Line 132 | Line 236 | float pileup(TTree *events, string infor
236		if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
237		<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
238
239	<	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
239	>	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
240
241		}
242
243		//____________________________________________________________________________________
244	+	// Effect of peak shifting
245	+	void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
246	+	TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
247	+	TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
248	+	TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
249	+	TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
250	+	TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
251	+	TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
252	+
253	+	nBins = 1;
254	+	string informalname="PeakErrorCalculation";
255	+	float resup,resdown,rescent;
256	+	for(int i=0;i<3;i++) {
257	+	string poscut,negcut;
258	+	if(i==0) {
259	+	poscut=peakcentral;
260	+	negcut=npeakcentral;
261	+	} else if(i==1) {
262	+	poscut=peakdown;
263	+	negcut=npeakdown;
264	+	} else if(i==2) {
265	+	poscut=peakup;
266	+	negcut=npeakup;
267	+	}
268	+	float res;
269	+	if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
270	+	else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str());
271	+	if(i==0) rescent=res;
272	+	else if(i==1) resdown=res;
273	+	else if(i==2) resup=res;
274	+	}
275	+	if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent);
276	+	else result=(TMath::Abs(rescent-resdown)/rescent);
277	+	}
278	+
279	+	//____________________________________________________________________________________
280		// Total selection efficiency (MC)
281	<	void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb) {
281	>	//returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
282	>	Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
283
284		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
285		// All acceptance cuts at gen. level
286	<	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
286	>	//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
287	>	TCut kbase("");
288	>	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
289	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
290		// Corresponding reco. cuts
291	<	TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
292	<
293	<	events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
294	<	Float_t sel = events->GetSelectedRows();
295	<	events->Draw(mcjzbexpression.c_str(),kbase,"goff");
296	<	Float_t tot = events->GetSelectedRows();
297	<
298	<	result=sel/tot;
299	<	resulterr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
300	<	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << std::endl;
291	>	TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
292	>	TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
293	>	TCut posSide = kbase&&ksel;
294	>	TCut negSide = kbase&&ksel2;
295	>	string sposSide(posSide);
296	>	string snegSide(negSide);
297	>	char var[20];
298	>	sprintf(var, "pdfW[%d]", k);
299	>	string svar(var);
300	>	string newPosSide = "(" + sposSide + ")*" + svar;
301	>	string newNegSide = "(" + snegSide + ")*" + svar;
302	>
303	>	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
304	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
305	>	else events->Draw((mcjzbexpression+">>effh").c_str(), sposSide.c_str(),"goff");
306	>	Float_t sel = effh->Integral();
307	>	Float_t nsel=0;
308	>	if(ConsiderSignalContaminationForLimits) {
309	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
310	>	else events->Draw((mcjzbexpression+">>effh").c_str(), snegSide.c_str(),"goff");
311	>	nsel = effh->Integral();
312	>	}
313	>	//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.
314	>	float normFactor = 1;
315	>	if(k>=0) get_norm_pdf_factor(events, k);
316	>	sel = sel/normFactor;
317	>	nsel = nsel/normFactor;
318	>
319	>	//      events->Draw(mcjzbexpression.c_str(),kbase,"goff");
320	>	//      Float_t tot = events->GetSelectedRows();
321	>	Float_t tot = Neventsinfile;
322	>
323	>	Value result_wo_signalcont;
324	>
325	>	if(ConsiderSignalContaminationForLimits) {
326	>	result=(sel-nsel)/tot;
327	>	resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
328	>	result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
329	>	} else {//no signal contamination considered:
330	>	result=(sel)/tot;
331	>	resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
332	>	result_wo_signalcont=Value(result,resulterr);
333	>	}
334	>	if(!automatized && k>0 ) dout << "PDF assessment: ";
335	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
336	>	delete effh;
337	>	}
338	>
339	>
340	>	//____________________________________________________________________________________
341	>	// Selection efficiency for one process (MC)
342	>	vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
343	>	vector<float> process_efficiencies;
344	>	for(int iprocess=0;iprocess<=10;iprocess++) {
345	>	float this_process_efficiency,efferr;
346	>	stringstream addcutplus;
347	>	addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
348	>	MCefficiency(events,this_process_efficiency, efferr,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
349	>	process_efficiencies.push_back(this_process_efficiency);
350	>	}
351	>	return process_efficiencies;
352		}
353	+
354
355	<	float JZBefficiency(TTree *events, string informalname) {
356	<	TCut kbase("abs(genMllSel-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
355	>	void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
356	>	TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
357	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
358	>	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
359		TH1F* hLM4 = plotEff(events,kbase,informalname);
360		Int_t bin = hLM4->FindBin(jzbSel); // To get the error
361	+	jzbeff=Interpolate(jzbSel,hLM4);
362	+	jzbefferr=hLM4->GetBinError(bin);
363		if(!automatized) dout << "  Efficiency at JZB==" << jzbSel  << std::endl;
364	<	if(!automatized) dout << "    " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin)  << std::endl;
165	<	return -1;
364	>	if(!automatized) dout << "    " << jzbeff << "+-" << jzbefferr  << std::endl;
365		}
366
367		//________________________________________________________________________
368		// Effect of energy scale on efficiency
369	<	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
370	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0");
369	>	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
370	>	TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
371	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
372	>	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
373	>
374		TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
375		TCut nJets("pfJetGoodNum>2");
376		stringstream down,up;
#	Line 178 | Line 380 | void JZBjetScale(TTree *events, float &j
380		TCut nJetsP(up.str().c_str());
381		TCut nJetsM(down.str().c_str());
382
181	–	if ( jzbSelection>0 ) jzbSel = jzbSelection;
182	–
383		if ( !(plotName.Length()>1) ) plotName = informalname;
384
385		nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
#	Line 197 | Line 397 | void JZBjetScale(TTree *events, float &j
397		if(!automatized) dout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
398		if(!automatized) dout << "    central:  " << eff << std::endl;
399		if(!automatized) dout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
400	<	jesup=(effp-eff)/eff*100.;
401	<	jesdown=(effm-eff)/eff*100.;
400	>	jesup=(effp-eff)/eff;
401	>	jesdown=(effm-eff)/eff;
402		}
403
404		//________________________________________________________________________
405		// Effect of energy scale on JZB efficiency
406	<	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname) {
406	>	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
407
408	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&genNjets>2");
408	>	TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
409	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
410	>	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
411		TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
412
413		nBins =    50;
#	Line 218 | Line 420 | void doJZBscale(TTree *events, float &do
420		Float_t eff  = Interpolate(jzbSel,hist);
421		Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
422		Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
423	<	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
423	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
424		if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
425	<	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
426	<	up=((effp-eff)/eff)*100;
427	<	down=((effm-eff)/eff)*100;
425	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
426	>	up=((effp-eff)/eff);
427	>	down=((effm-eff)/eff);
428		}
429
430		//________________________________________________________________________
431		// JZB response (true/reco. vs. true)
432	<	void JZBresponse(TTree *events, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
432	>	void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
433
434		jzbMin = 20;
435	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&genNjets>2");
435	>	TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
436	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
437	>	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
438		TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
439
440	<	TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true",
237	<	nPeriods, jzbMin, myJzbMax, "" );
440	>	TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
441
442	<	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZBSel:genJZBSel",kbase&&ksel);
442	>	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZB:genJZB",kbase&&ksel);
443		else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
444
445		hJzbResp->SetMaximum(1.2);
446		hJzbResp->SetMinimum(0.2);
447		hJzbResp->Fit("pol0","Q");
448		TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
449	<	if(!automatized) dout << "  Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
449	>	if(!fittedfunction) {
450	>	// in case there are not enough points passing our selection
451	>	cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
452	>	resp=1;
453	>	resperr=1;
454	>	} else {
455	>	resp=fittedfunction->GetParameter(0);
456	>	resperr=fittedfunction->GetParError(0);
457	>	if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
458	>	}
459		delete hJzbResp;
460		}
461
462
463	<	void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb) {
463	>	//________________________________________________________________________________________
464	>	// PDF uncertainty
465	>	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
466	>	std::vector<float> efficiency;
467	>	for(int k = 1; k < NPdfs; k++) {
468	>	float result, resulterr;
469	>	MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
470	>	efficiency.push_back(result);
471	>	}
472	>	float errHi, errLow,err;
473	>	master_formula(efficiency, errHi, errLow);
474	>	err=errLow;
475	>	if(errHi>errLow) err=errHi;
476	>	if(!automatized) dout << "  Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
477	>	return err;
478	>
479	>	}
480	>
481	>	int get_npdfs(TTree *events) {
482	>	int NPDFs;
483	>	events->SetBranchAddress("NPdfs",&NPDFs);
484	>	events->GetEntry(1);
485	>	return NPDFs;
486	>	}
487
488	+
489	+	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) {
490		float JetEnergyScaleUncert=0.1;
491		float JZBScaleUncert=0.1;
492		mcjzbexpression=mcjzb;
493	<
257	<	float triggereff=4;//percent!
493	>	float triggereff=5.0/100;// in range [0,1]
494		dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
495	<	float leptonseleff=2;//percent!
495	>	float leptonseleff=2.0/100;// in range [0,1]
496	>	leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
497		dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
498
499	<	float mceff,mcefferr;
499	>	int NPdfs=0;
500	>	if(ismSUGRA) NPdfs = get_npdfs(events);
501	>
502	>	float mceff,mcefferr,jzbeff,jzbefferr;
503		if(!automatized) dout << "MC efficiencies:" << endl;
504	<	MCefficiency(events,mceff,mcefferr,mcjzb);
505	<	JZBefficiency(events,informalname);
504	>	MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
505	>	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
506	>	if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
507
508	+	if(!automatized) dout << "Error from Peak position:" << endl;
509	+	float sysfrompeak=0;
510	+	PeakError(events,sysfrompeak,mcjzb,peakerror,addcut);
511	+
512		if(!automatized) dout << "Jet energy scale: " << std::endl;
513		float jesup,jesdown;
514	<	JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert);
514	>	JZBjetScale(events,jesdown,jesup,informalname,requireZ,addcut,JetEnergyScaleUncert);
515
516		if(!automatized) dout << "JZB scale: " << std::endl;
517		float scaleup,scaledown,scalesyst;
518	<	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname);
518	>	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,requireZ,addcut);
519
520		if(!automatized) dout << "JZB response: " << std::endl;
521	<	JZBresponse(events);
521	>	float resp,resperr;
522	>	if(PlottingSetup::computeJZBresponse) {
523	>	if(!automatized) dout << "JZB response: " << std::endl;
524	>	JZBresponse(events,requireZ,resp,resperr,addcut);
525	>	}
526
527		if(!automatized) dout << "Pileup: " << std::endl;
528	<	float resolution=pileup(events,informalname);
529	<
528	>	float resolution;
529	>	resolution=pileup(events,requireZ,informalname,addcut);
530	>
531	>	float PDFuncert=0;
532	>	if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
533	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
534	>
535		dout << "_______________________________________________" << endl;
536	<	dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
537	<	dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl;
538	<	dout << "Trigger efficiency: " << triggereff << endl;
539	<	dout << "Lepton Sel Eff: " << leptonseleff << endl;
540	<	dout << "For JZB>" << jzbSel << endl;
541	<	dout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
542	<	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
543	<	dout << "Resolution : " << resolution << endl;
544	<
536	>	dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
537	>	if(addcut!="") dout << "With additional cut: " << addcut;
538	>	dout << endl;
539	>	dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
540	>	dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
541	>	dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
542	>	dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
543	>	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
544	>	dout << "Resolution : " << resolution << endl; // in range [0,1]
545	>	dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
546	>	if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
547	>	if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
548	>	if(PlottingSetup::computeJZBresponse)dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
549
550		float toterr=0;
551	<	toterr+=(triggereff/100)*(triggereff/100);
552	<	toterr+=(leptonseleff/100)*(leptonseleff/100);
553	<	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup/100)*(jesup/100); else toterr+=(jesdown/100)*(jesdown/100);
554	<	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup/100)*(scaleup/100); else toterr+=(scaledown/100)*(scaledown/100);
555	<	toterr+=(resolution/100)*(resolution/100);
556	<	toterr=TMath::Sqrt(toterr);
557	<	dout << "FINAL RESULT : " << mceff << " +/- "<< mcefferr << " (stat) +/- " << 100*toterr << " (syst)" << endl;
558	<	dout << "     we thus use the sqrt of the sum of the squares which is : " << 100*TMath::Sqrt(mcefferr*mcefferr+(toterr*toterr)) << endl;
551	>	toterr+=(triggereff)*(triggereff);
552	>	toterr+=(leptonseleff)*(leptonseleff);
553	>	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
554	>	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
555	>	toterr+=(resolution*resolution);
556	>	toterr+=(sysfrompeak*sysfrompeak);
557	>	if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
558	>	dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
559	>	float systerr=TMath::Sqrt(toterr)*mceff;
560	>	toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
561	>
562	>	dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst)   %" << endl;
563	>	dout << "     we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
564	>
565	>	//Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
566		vector<float> res;
567		res.push_back(jzbSel);
568		res.push_back(mceff);
569		res.push_back(mcefferr);
570		res.push_back(toterr);
571		res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
572	<
572	>	if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
573	>	if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
574	>	res.push_back(fabs(resolution));
575	>	if(ismSUGRA) res.push_back(PDFuncert);
576		results.push_back(res);
577		}
578
579	<	vector<vector<float> > compute_systematics(string mcjzb, string datajzb, samplecollection &signalsamples, vector<float> bins) {
579	>	vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
580		automatized=true;
581		vector< vector<float> > systematics;
582		for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
#	Line 316 | Line 584 | vector<vector<float> > compute_systemati
584		for(int ibin=0;ibin<bins.size();ibin++) {
585		jzbSel=bins[ibin];
586		geqleq="geq";
587	<	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb);
587	>	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,mcjzb,datajzb,mcpeakerror,requireZ);
588		}//end of bin loop
589		}//end of signal loop
590		return systematics;

Diff Legend

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