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.2 by buchmann, Tue Jul 19 07:04:41 2011 UTC vs.
Revision 1.46 by buchmann, Wed Nov 9 10:55:25 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 36 | Line 38 | int iplot=0;
38		int verbose=0;
39		string geqleq;
40		string mcjzbexpression;
41	+	bool automatized=false;//if we're running this fully automatized we don't want each function to flood the screen
42
43		TString geq_or_leq() {
44		if(geqleq=="geq") return TString(">=");
#	Line 43 | 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 74 | 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, 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(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(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(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(locmcjzbexpression)+")>>"+TString(hname),kBaseCut&&kMassCut&&JZBNegCut&&cutOSOF,"goff");
106	+
107	+	float obs=0;
108	+	float pred=0;
109	+	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
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(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(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(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(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;
143	+	return obs-pred;
144	+	}
145	+
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		// Define new histogram
#	Line 86 | Line 155 | TH1F* plotEff(TTree* events, TCut kbase,
155		nBins,jzbMin,jzbMax);
156		Float_t step = (jzbMax-jzbMin)/static_cast<Float_t>(nBins);
157
158	<	events->Draw(mcjzbexpression.c_str(),"genJZBSel>-400"&&kbase,"goff");
158	>	if(flipped==0) events->Draw(mcjzbexpression.c_str(),"genJZB>-400"&&kbase,"goff");
159	>	else events->Draw(("-"+mcjzbexpression).c_str(),"genJZB>-14000"&&kbase,"goff");
160		Float_t maxEff = events->GetSelectedRows();
161	<	if(verbose>0) std::cout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
161	>	if(verbose>0) dout << hname << " (" << informalname <<") " << maxEff <<  std::endl;
162
163	<	if(verbose>0) std::cout <<  "JZB max = " << jzbMax << std::endl;
163	>	if(verbose>0) dout <<  "JZB max = " << jzbMax << std::endl;
164		// Loop over steps to get efficiency curve
165		char cut[256];
166		for ( Int_t iBin = 0; iBin<nBins; ++iBin ) {
167	<	sprintf(cut,"genJZBSel>%3f",jzbMin+iBin*step);
167	>	sprintf(cut,"genJZB>%3f",jzbMin+iBin*step);
168		events->Draw(mcjzbexpression.c_str(),TCut(cut)&&kbase,"goff");
169		Float_t eff = static_cast<Float_t>(events->GetSelectedRows())/maxEff;
170	<	//     std::cout << "COUCOU " << __LINE__ << std::endl;
170	>	//     dout << "COUCOU " << __LINE__ << std::endl;
171		hJzbEff->SetBinContent(iBin+1,eff);
172		hJzbEff->SetBinError(iBin+1,TMath::Sqrt(eff*(1-eff)/maxEff));
173		}
#	Line 108 | Line 178 | TH1F* plotEff(TTree* events, TCut kbase,
178
179
180		//________________________________________________________________________________________
181	+	// Master Formula
182	+	void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
183	+
184	+	float x0 = eff[0];
185	+	float deltaPos = 0, deltaNeg = 0;
186	+	for(int k = 0; k < (eff.size()-1)/2; k++) {
187	+	float xneg = eff[2*k+2];
188	+	float xpos = eff[2*k+1];
189	+	if(xpos-x0>0 || xneg-x0>0) {
190	+	if(xpos-x0 > xneg-x0) {
191	+	deltaPos += (xpos-x0)*(xpos-x0);
192	+	} else {
193	+	deltaPos += (xneg-x0)*(xneg-x0);
194	+	}
195	+	}
196	+	if(x0-xpos>0 || x0-xneg>0) {
197	+	if(x0-xpos > x0-xneg) {
198	+	deltaNeg += (xpos-x0)*(xpos-x0);
199	+	} else {
200	+	deltaNeg += (xneg-x0)*(xneg-x0);
201	+	}
202	+	}
203	+	}
204	+	errHi = sqrt(deltaPos);
205	+	errLo = sqrt(deltaNeg);
206	+
207	+	}
208	+
209	+
210	+	//________________________________________________________________________________________
211	+	// Get normalization factor for the PDFs
212	+	float get_norm_pdf_factor(TTree *events, int k, string addcut) {
213	+
214	+	TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
215	+	char nameVar[20];
216	+	sprintf(nameVar, "pdfW[%d]", k);
217	+	events->Project("haux", nameVar, addcut.c_str());
218	+	float thisW = haux->Integral();
219	+	events->Project("haux", "pdfW[0]");
220	+	float normW = haux->Integral();
221	+
222	+	float factor=thisW/normW;
223	+
224	+	delete haux;
225	+
226	+	return factor;
227	+
228	+	}
229	+
230	+
231	+
232	+	//________________________________________________________________________________________
233		// Pile-up efficiency
234	<	float pileup(TTree *events, string informalname, Float_t myJzbMax = 140. ) {
234	>	float pileup(TTree *events, bool requireZ, string informalname, int flipped, string addcut="",Float_t myJzbMax = 140. ) {
235		nBins = 16;
236		jzbMax = myJzbMax;
237
238		// Acceptance cuts
239	<	TCut kbase("abs(genMllSel-91.2)<20&&pfJetGoodNum>2&&genZPtSel>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
239	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
240	>	TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
241	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
242
243	<	TH1F* hLM4 = plotEff(events,kbase,informalname);
243	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
244	>	TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
245		hLM4->SetMinimum(0.);
246
247		// Nominal function
#	Line 128 | Line 253 | float pileup(TTree *events, string infor
253		// Pimped-up function
254		TF1* funcUp = (TF1*)func->Clone();
255		funcUp->SetParameter( 0., func->GetParameter(0)/1.1); // 10% systematic error (up in sigma => 0.1 in erfc)
256	<	std::cout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
256	>	if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
257		<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
258
259	<	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
259	>	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
260
261		}
262
263		//____________________________________________________________________________________
264	+	// Effect of peak shifting
265	+	void PeakError(TTree *events,float &result, string mcjzb, float peakerr,int flipped,string addcut="") {
266	+	//Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
267	+	TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
268	+	TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
269	+	TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
270	+	TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
271	+	TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
272	+	TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
273	+	nBins = 1;
274	+	string informalname="PeakErrorCalculation";
275	+	float resup,resdown,rescent;
276	+	for(int i=0;i<3;i++) {
277	+	string poscut,negcut;
278	+	if(i==0) {
279	+	poscut=peakcentral;
280	+	negcut=npeakcentral;
281	+	} else if(i==1) {
282	+	poscut=peakdown;
283	+	negcut=npeakdown;
284	+	} else if(i==2) {
285	+	poscut=peakup;
286	+	negcut=npeakup;
287	+	}
288	+	float res;
289	+	if(addcut=="") res=allcontributionsplot(events,cutnJets,cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
290	+	else res=allcontributionsplot(events,cutnJets&&addcut.c_str(),cutmass,sidebandcut,poscut.c_str(),negcut.c_str(),flipped);
291	+	if(i==0) rescent=res;
292	+	else if(i==1) resdown=res;
293	+	else if(i==2) resup=res;
294	+	}
295	+	if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
296	+	else result=(TMath::Abs(rescent-resdown)/(float)rescent);
297	+	}
298	+
299	+
300	+	void  MCPartialefficiency(TTree *events,float &result, float &resulterr,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0, int type = 0) {
301	+	if(!events) {
302	+	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
303	+	result=0;resulterr=0;
304	+	return;
305	+	}
306	+
307	+	char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
308	+	// All acceptance cuts at gen. level
309	+	//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
310	+	TCut kbase("");
311	+	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
312	+	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
313	+	// Corresponding reco. cuts
314	+
315	+	TCut acceptance("genPt2 != 0");
316	+	TCut massId(cutmass&&cutOSSF);
317	+	TCut njets(cutnJets);
318	+	TCut jzbp;
319	+	TCut jzbn;
320	+	if(flipped==0) {
321	+	jzbp=TCut((TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
322	+	jzbn=TCut((TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
323	+	} else {
324	+	jzbp=TCut(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
325	+	jzbn=TCut(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
326	+	}
327	+	float ntotal = events->Draw("pt1", kbase, "goff");
328	+	TCut theCut;
329	+	switch(type) {
330	+	case 1:
331	+	theCut = kbase+acceptance;
332	+	break;
333	+	case 2:
334	+	theCut = kbase+massId;
335	+	break;
336	+	case 3:
337	+	theCut = kbase+massId+njets;
338	+	break;
339	+	case 4:
340	+	theCut = kbase+massId+njets+jzbn;
341	+	break;
342	+	default:
343	+	theCut = kbase+massId+njets+jzbn;
344	+	break;
345	+	}
346	+
347	+	string stheCut(theCut);
348	+	char var[20];
349	+	sprintf(var, "pdfW[%d]", k);
350	+
351	+	string svar(var);
352	+	string newtheCut;
353	+	if(k>0) newtheCut = "(" + stheCut + ")*" + svar;
354	+	else newtheCut = "(" + stheCut + ")"; // for k==0 or even k==-1 we don't need to evaluate PDFs
355	+
356	+	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
357	+	if(k>=0) events->Draw((mcjzbexpression+">>effh").c_str(), newtheCut.c_str(),"goff");
358	+	else events->Draw((mcjzbexpression+">>effh").c_str(), theCut,"goff");
359	+	Float_t sel = effh->Integral();
360	+	Float_t nsel=0;
361	+	//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.
362	+	float normFactor = 1;
363	+	if(k>=0) get_norm_pdf_factor(events, k, addcut);
364	+	sel = sel/normFactor;
365	+
366	+	result=(sel)/ntotal;
367	+	resulterr=TMath::Sqrt(sel/ntotal*(1+sel/ntotal)/ntotal);
368	+
369	+	}
370	+
371	+	//____________________________________________________________________________________
372		// Total selection efficiency (MC)
373	<	void MCefficiency(TTree *events,float &res, float &reserr,string mcjzb) {
373	>	//returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
374	>	Value MCefficiency(TTree *events,float &result, float &resulterr, int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
375	>	if(!events) {
376	>	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
377	>	result=0;
378	>	resulterr=0;
379	>	return Value(0,0);
380	>	}
381
382		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
383		// All acceptance cuts at gen. level
384	<	TCut kbase("abs(genMllSel-91.2)<20&&pfJetGoodNum>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
385	<	// Corresponding reco. cuts
146	<	TCut ksel("abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
384	>	//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
385	>	TCut kbase("");
386
387	<	events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
388	<	Float_t sel = events->GetSelectedRows();
389	<	events->Draw(mcjzbexpression.c_str(),kbase,"goff");
390	<	Float_t tot = events->GetSelectedRows();
387	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
388	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
389	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
390	>	// Corresponding reco. cuts
391
392	<	res=sel/tot;
393	<	reserr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
394	<	std::cout << "  MC efficiency: " << res << "+-" << reserr << std::endl;
392	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
393	>	TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
394	>	TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
395	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
396	>	if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
397	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
398	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
399	>	if(flipped>0) {
400	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
401	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
402	>	}
403	>	} else {
404	>	//for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
405	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
406	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
407	>	if(flipped>0) {
408	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
409	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
410	>	}
411	>	}
412	>
413	>	TCut posSide = kbase&&ksel;
414	>	TCut negSide = kbase&&ksel2;
415	>	string sposSide(posSide);
416	>	string snegSide(negSide);
417	>	char var[20];
418	>	sprintf(var, "pdfW[%d]", k);
419	>	if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
420	>	string svar(var);
421	>	string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
422	>	string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
423	>	string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
424	>	string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
425	>
426	>	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
427	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
428	>	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.
429	>	Float_t sel = effh->Integral();
430	>	Float_t nsel=0;
431	>
432	>	///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
433	>
434	>	if(ConsiderSignalContaminationForLimits) {
435	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
436	>	if(PlottingSetup::RestrictToMassPeak) {
437	>	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
438	>	nsel += effh->Integral();
439	>	} else {
440	>	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
441	>	nsel += effh->Integral();
442	>	events->Draw((mcjzbexpression+">>effh").c_str(), emnewPosSide.c_str(),"goff");
443	>	nsel += effh->Integral();
444	>	events->Draw((mcjzbexpression+">>effh").c_str(), emnewNegSide.c_str(),"goff");
445	>	nsel -= effh->Integral();
446	>	}
447	>	}
448	>
449	>	//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.
450	>	float normFactor = 1;
451	>	if(k>=0) get_norm_pdf_factor(events, k, addcut);
452	>	sel = sel/normFactor;
453	>	nsel = nsel/normFactor;
454	>
455	>	//      events->Draw(mcjzbexpression.c_str(),kbase,"goff");
456	>	//      Float_t tot = events->GetSelectedRows();
457	>	Float_t tot = Neventsinfile;
458	>
459	>	Value result_wo_signalcont;
460	>
461	>	if(ConsiderSignalContaminationForLimits) {
462	>	result=(sel-nsel)/tot;
463	>	resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
464	>	result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
465	>	} else {//no signal contamination considered:
466	>	result=(sel)/tot;
467	>	resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
468	>	result_wo_signalcont=Value(result,resulterr);
469	>	}
470	>	if(!automatized && k>0 ) dout << "PDF assessment [" << k << "] : ";
471	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
472	>	delete effh;
473	>	return result_wo_signalcont;
474		}
475
476	<	float JZBefficiency(TTree *events, string informalname) {
477	<	TCut kbase("abs(genMllSel-91.2)<20&&pfJetGoodNum>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
478	<	TH1F* hLM4 = plotEff(events,kbase,informalname);
476	>
477	>
478	>	//____________________________________________________________________________________
479	>	// Selection efficiency for one process (MC)
480	>	// not in use anymore.
481	>	/*
482	>	vector<float> processMCefficiency(TTree *events,int flipped, string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
483	>	vector<float> process_efficiencies;
484	>	for(int iprocess=0;iprocess<=10;iprocess++) {
485	>	float this_process_efficiency,efferr;
486	>	stringstream addcutplus;
487	>	addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
488	>	MCefficiency(events,this_process_efficiency, efferr,flipped,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
489	>	process_efficiencies.push_back(this_process_efficiency);
490	>	}
491	>	return process_efficiencies;
492	>	}
493	>	*/
494	>
495	>	void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, int flipped, bool requireZ, string addcut="") {
496	>	TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
497	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
498	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
499	>	TH1F* hLM4 = plotEff(events,kbase,informalname,flipped);
500		Int_t bin = hLM4->FindBin(jzbSel); // To get the error
501	<	std::cout << "  Efficiency at JZB==" << jzbSel  << std::endl;
502	<	std::cout << "    " << Interpolate(jzbSel,hLM4) << "+-" << hLM4->GetBinError(bin)  << std::endl;
503	<	return -1;
501	>	jzbeff=Interpolate(jzbSel,hLM4);
502	>	jzbefferr=hLM4->GetBinError(bin);
503	>	if(!automatized) dout << "  Efficiency at JZB==" << jzbSel  << std::endl;
504	>	if(!automatized) dout << "    " << jzbeff << "+-" << jzbefferr  << std::endl;
505		}
506
507		//________________________________________________________________________
508		// Effect of energy scale on efficiency
509	<	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
510	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0");
511	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
509	>	void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname, int flipped, bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
510	>	TCut kbase(genMassCut&&"genZPt>0");
511	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
512	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
513	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
514	>
515	>	TCut ksel(cutmass&&cutOSSF);
516		TCut nJets("pfJetGoodNum>2");
517		stringstream down,up;
518		down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
#	Line 177 | Line 521 | void JZBjetScale(TTree *events, float &j
521		TCut nJetsP(up.str().c_str());
522		TCut nJetsM(down.str().c_str());
523
180	–	if ( jzbSelection>0 ) jzbSel = jzbSelection;
181	–
524		if ( !(plotName.Length()>1) ) plotName = informalname;
525
526		nBins = 1; jzbMin = jzbSel*0.95; jzbMax = jzbSel*1.05;
527	<	TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname);
527	>	TH1F* hist = plotEff(events,(kbase&&ksel&&nJets),informalname,flipped);
528
529	<	TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname);
529	>	TH1F* histp = plotEff(events,(kbase&&ksel&&nJetsP),informalname,flipped);
530
531	<	TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname);
531	>	TH1F* histm = plotEff(events,(kbase&&ksel&&nJetsM),informalname,flipped);
532
533		// Dump some information
534		Float_t eff  = Interpolate(jzbSel,hist);
535		Float_t effp = Interpolate(jzbSel,histp);
536		Float_t effm = Interpolate(jzbSel,histm);
537	<	std::cout << "  Efficiency at JZB==" << jzbSel  << std::endl;
538	<	std::cout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
539	<	std::cout << "    central:  " << eff << std::endl;
540	<	std::cout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
541	<	jesup=(effp-eff)/eff*100.;
542	<	jesdown=(effm-eff)/eff*100.;
537	>	if(!automatized) dout << "  Efficiency at JZB==" << jzbSel  << std::endl;
538	>	if(!automatized) dout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
539	>	if(!automatized) dout << "    central:  " << eff << std::endl;
540	>	if(!automatized) dout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
541	>	jesup=(effp-eff)/eff;
542	>	jesdown=(effm-eff)/eff;
543		}
544
545		//________________________________________________________________________
546		// Effect of energy scale on JZB efficiency
547	<	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname) {
547	>	void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, int flipped, bool requireZ, string addcut) {
548
549	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPt>0&&pfJetGoodNum>2");
550	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
549	>	TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
550	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
551	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
552	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
553	>	TCut ksel(cutmass&&cutOSSF);
554
555		nBins =    50;
556		jzbMin =   0.5*jzbSel;
557		jzbMax =   2.0*jzbSel;
558
559	<	TH1F* hist = plotEff(events,kbase&&ksel,informalname);
559	>	TH1F* hist = plotEff(events,kbase&&ksel,informalname,flipped);
560
561		// Dump some information
562		Float_t eff  = Interpolate(jzbSel,hist);
563		Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
564		Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
565	<	std::cout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
566	<	std::cout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
567	<	std::cout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
568	<	up=((effp-eff)/eff)*100;
569	<	down=((effm-eff)/eff)*100;
565	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
566	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
567	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
568	>	up=((effp-eff)/eff);
569	>	down=((effm-eff)/eff);
570		}
571
572		//________________________________________________________________________
573		// JZB response (true/reco. vs. true)
574	<	void JZBresponse(TTree *events, bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
574	>	void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, int flipped, string addcut="", bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
575
576		jzbMin = 20;
577	<	TCut kbase("abs(genMllSel-91.2)<20&&genZPtSel>0&&pfJetGoodNum>2");
578	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
579	<
580	<	TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true",
581	<	nPeriods, jzbMin, myJzbMax, "" );
582	<
583	<	if (!isMET) events->Project("hJzbResp","("+TString(mcjzbexpression)+")/genJZBSel:genJZBSel",kbase&&ksel);
577	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
578	>	TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
579	>	if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
580	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
581	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
582	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
583	>	TCut ksel(cutmass&&cutOSSF);
584	>
585	>	TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
586	>
587	>	string locmcjzbexpression=mcjzbexpression;
588	>	if(flipped>0) locmcjzbexpression="-"+locmcjzbexpression;
589	>	string possibleminus="";
590	>	if(flipped>0) possibleminus="-";
591	>	if (!isMET) events->Project("hJzbResp","("+TString(locmcjzbexpression)+")/("+possibleminus+"genJZB):("+possibleminus+"genJZB)",kbase&&ksel);
592		else events->Project("hJzbResp","met[4]/genMET:genMET",kbase&&ksel);
593
594		hJzbResp->SetMaximum(1.2);
595		hJzbResp->SetMinimum(0.2);
596		hJzbResp->Fit("pol0","Q");
597		TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
598	<	cout << "  Response: " << fittedfunction->GetParameter(0) << " +/- " << fittedfunction->GetParError(0) << endl;
598	>	if(!fittedfunction) {
599	>	// in case there are not enough points passing our selection
600	>	cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
601	>	resp=1;
602	>	resperr=1;
603	>	} else {
604	>	resp=fittedfunction->GetParameter(0);
605	>	resperr=fittedfunction->GetParError(0);
606	>	if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
607	>	}
608	>	delete hJzbResp;
609		}
610
611
612	<	void do_systematics_for_one_file(TTree *events,string informalname, vector<vector<float> > &uncertainties,string mcjzb,string datajzb) {
612	>	//________________________________________________________________________________________
613	>	// PDF uncertainty
614	>	float get_pdf_uncertainty(TTree *events, int flipped, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
615	>	std::vector<float> efficiency;
616	>	for(int k = 1; k < NPdfs; k++) {
617	>	float result, resulterr;
618	>	Value flipval;
619	>	MCefficiency(events, result, resulterr, flipped, mcjzb, requireZ, Neventsinfile, addcut, k);
620	>	efficiency.push_back(result);
621	>	}
622	>	float errHi, errLow,err;
623	>	master_formula(efficiency, errHi, errLow);
624	>	err=errLow;
625	>	if(errHi>errLow) err=errHi;
626	>	if(!automatized) dout << "  Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
627	>	return err;
628	>
629	>	}
630	>
631	>	int get_npdfs(TTree *events) {
632	>	int NPDFs;
633	>	events->SetBranchAddress("NPdfs",&NPDFs);
634	>	events->GetEntry(1);
635	>	return NPDFs;
636	>	}
637
638	+
639	+	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) {
640		float JetEnergyScaleUncert=0.1;
641		float JZBScaleUncert=0.1;
642		mcjzbexpression=mcjzb;
643	+	float triggereff=5.0/100;// in range [0,1]
644	+	dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
645	+	float leptonseleff=2.0/100;// in range [0,1]
646	+	leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
647	+	dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
648
649	<	float triggereff=4;//percent!
650	<	cout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
257	<	float leptonseleff=2;//percent!
258	<	cout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
649	>	int NPdfs=0;
650	>	if(ismSUGRA) NPdfs = get_npdfs(events);
651
652	<	float mceff,mcefferr;
653	<	cout << "MC efficiencies:" << endl;
654	<	MCefficiency(events,mceff,mcefferr,mcjzb);
655	<	JZBefficiency(events,informalname);
652	>	float mceff,mcefferr,jzbeff,jzbefferr;
653	>	if(!automatized) dout << "MC efficiencies:" << endl;
654	>	Value flipefficiency;
655	>	Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,flipped,mcjzb,requireZ,Neventsinfile,addcut,-1);
656	>	if(!automatized) cout << "   Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
657	>
658	>	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,flipped,requireZ,addcut);
659	>	if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
660
661	<	std::cout << "Jet energy scale: " << std::endl;
661	>	if(!automatized) dout << "Error from Peak position:" << endl;
662	>	float sysfrompeak=0;
663	>	PeakError(events,sysfrompeak,mcjzb,peakerror,flipped,addcut);
664	>
665	>	if(!automatized) dout << "Jet energy scale: " << std::endl;
666		float jesup,jesdown;
667	<	JZBjetScale(events,jesdown,jesup,informalname,JetEnergyScaleUncert);
667	>	JZBjetScale(events,jesdown,jesup,informalname,flipped,requireZ,addcut,JetEnergyScaleUncert);
668
669	<	std::cout << "JZB scale: " << std::endl;
669	>	if(!automatized) dout << "JZB scale: " << std::endl;
670		float scaleup,scaledown,scalesyst;
671	<	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname);
671	>	doJZBscale(events,scaledown,scaleup,scalesyst,JZBScaleUncert,informalname,flipped,requireZ,addcut);
672
673	<	std::cout << "JZB response: " << std::endl;
674	<	JZBresponse(events);
675	<
676	<	std::cout << "Pileup: " << std::endl;
677	<	float resolution=pileup(events,informalname);
673	>	if(!automatized) dout << "JZB response: " << std::endl;
674	>	float resp,resperr;
675	>	if(PlottingSetup::computeJZBresponse) {
676	>	if(!automatized) dout << "JZB response: " << std::endl;
677	>	JZBresponse(events,requireZ,resp,resperr,flipped,addcut);
678	>	}
679	>
680	>	if(!automatized) dout << "Pileup: " << std::endl;
681	>	float resolution;
682	>	resolution=pileup(events,requireZ,informalname,flipped,addcut);
683	>
684	>	float PDFuncert=0;
685	>	if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
686	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, flipped, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
687	>
688	>	dout << "_______________________________________________" << endl;
689	>	dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
690	>	if(addcut!="") dout << "With additional cut: " << addcut;
691	>	dout << endl;
692	>	dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
693	>	dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
694	>	dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
695	>	dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
696	>	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
697	>	dout << "Resolution : " << resolution << endl; // in range [0,1]
698	>	dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
699	>	if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
700	>	if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
701	>	if(PlottingSetup::computeJZBresponse)dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
702
703	<	cout << "_______________________________________________" << endl;
704	<	cout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << endl;
705	<	cout << "Trigger efficiency: " << triggereff << endl;
706	<	cout << "Lepton Sel Eff: " << leptonseleff << endl;
707	<	cout << "For JZB>" << jzbSel << endl;
708	<	cout << "Jet energy scale: " << jesup << " " << jesdown << " --> suggesting: " << Round(0.5*(fabs(jesup)+fabs(jesdown)),1) << endl;
709	<	cout << "JZB Scale Uncert: " << scaledown << " " << scaleup << " --> suggesting: " << Round(0.5*(fabs(scaledown)+fabs(scaleup)),1) << endl;
710	<	cout << "Resolution : " << resolution << endl;
703	>	float toterr=0;
704	>	toterr+=(triggereff)*(triggereff);
705	>	toterr+=(leptonseleff)*(leptonseleff);
706	>	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
707	>	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
708	>	toterr+=(resolution*resolution);
709	>	toterr+=(sysfrompeak*sysfrompeak);
710	>	if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
711	>	dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
712	>	float systerr=TMath::Sqrt(toterr)*mceff;
713	>	toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
714
715	<	vector<float> uncert;
716	<	uncert.push_back(jzbSel);
717	<	uncert.push_back(triggereff);
291	<	uncert.push_back(leptonseleff);
292	<	uncert.push_back(0.5*(fabs(jesup)+fabs(jesdown)));
293	<	uncert.push_back(0.5*(fabs(scaledown)+fabs(scaleup)));
294	<	uncert.push_back(resolution);
715	>	dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst)   %" << endl;
716	>	dout << "     we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
717	>	dout << "_______________________________________________" << endl;
718
719	<	uncertainties.push_back(uncert);
720	<	}
721	<
722	<	vector<vector<float> > compute_systematics(string mcjzb, string datajzb, samplecollection &signalsamples, vector<float> bins) {
723	<	vector< vector<float> > uncertainties;
719	>	//Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
720	>	vector<float> res;
721	>	res.push_back(jzbSel);
722	>	res.push_back(mceff);
723	>	res.push_back(mcefferr);
724	>	res.push_back(toterr);
725	>	res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
726	>	if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
727	>	if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
728	>	res.push_back(fabs(resolution));
729	>	res.push_back(mceff_nosigcont.getValue());
730	>	res.push_back(mceff_nosigcont.getError());
731	>	if(ismSUGRA) res.push_back(PDFuncert);
732	>	results.push_back(res);
733	>	}
734	>
735	>	vector<vector<float> > compute_systematics(string mcjzb, float mcpeakerror, int flipped, string datajzb, samplecollection &signalsamples, vector<float> bins, bool requireZ=false) {
736	>	automatized=true;
737	>	vector< vector<float> > systematics;
738		for (int isignal=0; isignal<signalsamples.collection.size();isignal++) {
739	<	cout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
739	>	dout << "Looking at signal " << (signalsamples.collection)[isignal].filename << endl;
740		for(int ibin=0;ibin<bins.size();ibin++) {
741		jzbSel=bins[ibin];
742		geqleq="geq";
743	<	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].samplename,uncertainties,mcjzb,datajzb);
743	>	do_systematics_for_one_file((signalsamples.collection)[isignal].events,(signalsamples.collection)[isignal].Nentries,(signalsamples.collection)[isignal].samplename,systematics,flipped,mcjzb,datajzb,mcpeakerror,requireZ);
744		}//end of bin loop
745		}//end of signal loop
746	<	return uncertainties;
746	>	return systematics;
747		}

Diff Legend

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