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

Diff Legend

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