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.42 by buchmann, Mon Oct 24 15:05:37 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 100 | Line 102 | float allcontributionsplot(TTree* events
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");
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
106	<	float obs = hossfp->Integral();
107	<	float pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
106	>	float obs=0;
107	>	float pred=0;
108	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
109	>	if(PlottingSetup::RestrictToMassPeak) {
110	>	hname=GetNumericHistoName();
111	>	TH1F* sbhossfp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
112	>	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSSF,"goff");
113	>	hname=GetNumericHistoName();
114	>	TH1F* sbhossfn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
115	>	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSSF,"goff");
116	>
117	>	hname=GetNumericHistoName();
118	>	TH1F* sbhosofp = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
119	>	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBPosCut&&cutOSOF,"goff");
120	>	hname=GetNumericHistoName();
121	>	TH1F* sbhosofn = new TH1F(hname.c_str(),hname.c_str(),1,-14000,14000);
122	>	events->Draw(TString(mcjzbexpression)+">>"+TString(hname),kBaseCut&&kSidebandCut&&JZBNegCut&&cutOSOF,"goff");
123	>
124	>	obs = hossfp->Integral();
125	>	pred= hossfn->Integral() + (1.0/3)*( hosofp->Integral() - hosofn->Integral() + sbhossfp->Integral() - sbhossfn->Integral() + sbhosofp->Integral() - sbhosofn->Integral());
126	>	delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
127	>	} else {
128	>	obs = hossfp->Integral();
129	>	pred= hossfn->Integral() + (hosofp->Integral() - hosofn->Integral());
130	>	}
131
132		delete hossfp,hossfn,hosofp,hosofn;
122	–	delete sbhossfp,sbhossfn,sbhosofp,sbhosofn;
133		return obs-pred;
134		}
135
#	Line 157 | Line 167 | TH1F* plotEff(TTree* events, TCut kbase,
167
168
169		//________________________________________________________________________________________
170	+	// Master Formula
171	+	void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
172	+
173	+	float x0 = eff[0];
174	+	float deltaPos = 0, deltaNeg = 0;
175	+	for(int k = 0; k < (eff.size()-1)/2; k++) {
176	+	float xneg = eff[2*k+2];
177	+	float xpos = eff[2*k+1];
178	+	if(xpos-x0>0 || xneg-x0>0) {
179	+	if(xpos-x0 > xneg-x0) {
180	+	deltaPos += (xpos-x0)*(xpos-x0);
181	+	} else {
182	+	deltaPos += (xneg-x0)*(xneg-x0);
183	+	}
184	+	}
185	+	if(x0-xpos>0 || x0-xneg>0) {
186	+	if(x0-xpos > x0-xneg) {
187	+	deltaNeg += (xpos-x0)*(xpos-x0);
188	+	} else {
189	+	deltaNeg += (xneg-x0)*(xneg-x0);
190	+	}
191	+	}
192	+	}
193	+	errHi = sqrt(deltaPos);
194	+	errLo = sqrt(deltaNeg);
195	+
196	+	}
197	+
198	+
199	+	//________________________________________________________________________________________
200	+	// Get normalization factor for the PDFs
201	+	float get_norm_pdf_factor(TTree *events, int k, string addcut) {
202	+
203	+	TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
204	+	char nameVar[20];
205	+	sprintf(nameVar, "pdfW[%d]", k);
206	+	events->Project("haux", nameVar, addcut.c_str());
207	+	float thisW = haux->Integral();
208	+	events->Project("haux", "pdfW[0]");
209	+	float normW = haux->Integral();
210	+
211	+	float factor=thisW/normW;
212	+
213	+	delete haux;
214	+
215	+	return factor;
216	+
217	+	}
218	+
219	+
220	+
221	+	//________________________________________________________________________________________
222		// Pile-up efficiency
223		float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
224		nBins = 16;
225		jzbMax = myJzbMax;
226
227		// Acceptance cuts
228	<	TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
228	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
229	>	TCut kbase(PlottingSetup::genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
230		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
231
232	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
232	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
233		TH1F* hLM4 = plotEff(events,kbase,informalname);
234		hLM4->SetMinimum(0.);
235
#	Line 182 | Line 245 | float pileup(TTree *events, bool require
245		if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
246		<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
247
248	<	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
248	>	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
249
250		}
251
252		//____________________________________________________________________________________
253		// Effect of peak shifting
254		void PeakError(TTree *events,float &result, string mcjzb, float peakerr,string addcut="") {
255	+	//Note: the cut used here is something like (JZBEXPRESSION+(peakerr)>50) without all the other cuts, to increase statistics (particularly for scans)
256		TString peakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
257		TString peakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+geq_or_leq()+TString(any2string(jzbSel)));
258		TString peakcentral("("+TString(mcjzb)+")"+geq_or_leq()+TString(any2string(jzbSel)));
259		TString npeakup("("+TString(mcjzb)+"+"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
260		TString npeakdown("("+TString(mcjzb)+"-"+TString(any2string(TMath::Abs(peakerr)))+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
261		TString npeakcentral("("+TString(mcjzb)+")"+ngeq_or_leq()+"-"+TString(any2string(jzbSel)));
198	–
262		nBins = 1;
263		string informalname="PeakErrorCalculation";
264		float resup,resdown,rescent;
#	Line 218 | Line 281 | void PeakError(TTree *events,float &resu
281		else if(i==1) resdown=res;
282		else if(i==2) resup=res;
283		}
284	<	if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent)*100;
285	<	else result=(TMath::Abs(rescent-resdown)/rescent)*100;
284	>	if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/(float)rescent);
285	>	else result=(TMath::Abs(rescent-resdown)/(float)rescent);
286		}
287
288		//____________________________________________________________________________________
289		// Total selection efficiency (MC)
290	<	void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="") {
290	>	//returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
291	>	Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
292	>	write_warning(__FUNCTION__,"Setting automatized to off!"); automatized=false;
293	>	if(!events) {
294	>	write_error(__FUNCTION__,"Tree passed for efficiency calculation is invalid!");
295	>	result=0;
296	>	resulterr=0;
297	>	return Value(0,0);
298	>	}
299
300		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
301		// All acceptance cuts at gen. level
302		//TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&genJZB"+geq_or_leq()+TString(jzbSelStr)+"&&genId1==-genId2");
303		TCut kbase("");
304	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
304	>
305	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
306	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
307		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
308		// Corresponding reco. cuts
309	<	TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
310	<	TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
311	<	events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
312	<	Float_t sel = events->GetSelectedRows();
313	<	events->Draw(mcjzbexpression.c_str(),kbase&&ksel2,"goff");
314	<	Float_t nsel = events->GetSelectedRows();
309	>
310	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
311	>	TCut ksel;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
312	>	TCut ksel2;//("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
313	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
314	>	if(PlottingSetup::RestrictToMassPeak||!ConsiderSignalContaminationForLimits) {
315	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
316	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&"id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
317	>	} else {
318	>	//for off peak analysis we don't use the OSSF condition here yet so we can recycle these two cuts for the em condition!
319	>	ksel=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+geq_or_leq()+TString(jzbSelStr)));
320	>	ksel2=TCut("pfJetGoodNum>2"&&cutmass&&(TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr)));
321	>	}
322	>
323	>	TCut posSide = kbase&&ksel;
324	>	TCut negSide = kbase&&ksel2;
325	>	string sposSide(posSide);
326	>	string snegSide(negSide);
327	>	char var[20];
328	>	sprintf(var, "pdfW[%d]", k);
329	>	if(k==-1) sprintf(var,"1.0");//case in which we don't want to evaluate PDFs
330	>	string svar(var);
331	>	string newPosSide = "((id1==id2)&&(" + sposSide + "))*" + svar;
332	>	string newNegSide = "((id1==id2)&&(" + snegSide + "))*" + svar;
333	>	string emnewPosSide = "((id1!=id2)&&(" + sposSide + "))*" + svar; // only used for off peak analysis
334	>	string emnewNegSide = "((id1!=id2)&&(" + snegSide + "))*" + svar; // only used for off peak analysis
335	>
336	>	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
337	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
338	>	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.
339	>	Float_t sel = effh->Integral();
340	>	Float_t nsel=0;
341	>
342	>	///----------------------------------------------- THIS PART REQUIRES STUDYING! -------------------------
343	>
344	>	if(ConsiderSignalContaminationForLimits) {
345	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
346	>	if(PlottingSetup::RestrictToMassPeak) {
347	>	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
348	>	nsel += effh->Integral();
349	>	} else {
350	>	events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
351	>	nsel += effh->Integral();
352	>	events->Draw((mcjzbexpression+">>effh").c_str(), emnewPosSide.c_str(),"goff");
353	>	nsel += effh->Integral();
354	>	events->Draw((mcjzbexpression+">>effh").c_str(), emnewNegSide.c_str(),"goff");
355	>	nsel -= effh->Integral();
356	>	}
357	>	}
358	>
359	>	//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.
360	>	float normFactor = 1;
361	>	if(k>=0) get_norm_pdf_factor(events, k, addcut);
362	>	sel = sel/normFactor;
363	>	nsel = nsel/normFactor;
364	>
365		//      events->Draw(mcjzbexpression.c_str(),kbase,"goff");
366		//      Float_t tot = events->GetSelectedRows();
367		Float_t tot = Neventsinfile;
368
369	<	result=(sel-nsel)/tot;
370	<	resulterr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
371	<	dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ")" << std::endl;
369	>	Value result_wo_signalcont;
370	>
371	>	if(ConsiderSignalContaminationForLimits) {
372	>	result=(sel-nsel)/tot;
373	>	resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
374	>	result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
375	>	} else {//no signal contamination considered:
376	>	result=(sel)/tot;
377	>	resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
378	>	result_wo_signalcont=Value(result,resulterr);
379	>	}
380	>	if(!automatized && k>0 ) dout << "PDF assessment: ";
381	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , signal contamination : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
382	>	delete effh;
383	>	return result_wo_signalcont;
384	>	}
385	>
386	>
387	>
388	>	//____________________________________________________________________________________
389	>	// Selection efficiency for one process (MC)
390	>	vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
391	>	vector<float> process_efficiencies;
392	>	for(int iprocess=0;iprocess<=10;iprocess++) {
393	>	float this_process_efficiency,efferr;
394	>	stringstream addcutplus;
395	>	addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
396	>	MCefficiency(events,this_process_efficiency, efferr,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
397	>	process_efficiencies.push_back(this_process_efficiency);
398	>	}
399	>	return process_efficiencies;
400		}
401	+
402
403		void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
404	<	TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
253	<	cout << "Getting started with JZB efficiency" << endl;
404	>	TCut kbase(genMassCut&&"genNjets>2&&genZPt>0"&&cutmass&&cutOSSF);
405		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
406	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
406	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
407		TH1F* hLM4 = plotEff(events,kbase,informalname);
408		Int_t bin = hLM4->FindBin(jzbSel); // To get the error
409		jzbeff=Interpolate(jzbSel,hLM4);
#	Line 264 | Line 415 | void JZBefficiency(TTree *events, string
415		//________________________________________________________________________
416		// Effect of energy scale on efficiency
417		void JZBjetScale(TTree *events, float &jesdown, float &jesup, string informalname,bool requireZ,string addcut="",float syst=0.1, Float_t jzbSelection=-1, TString plotName = "" ) {
418	<	TCut kbase("abs(genMll-91.2)<20&&genZPt>0");
418	>	TCut kbase(genMassCut&&"genZPt>0");
419		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
420	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
420	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
421	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
422
423	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
423	>	TCut ksel(cutmass&&cutOSSF);
424		TCut nJets("pfJetGoodNum>2");
425		stringstream down,up;
426		down << "pfJetGoodNum"<<30*(1-syst)<<">=3";
#	Line 294 | Line 446 | void JZBjetScale(TTree *events, float &j
446		if(!automatized) dout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
447		if(!automatized) dout << "    central:  " << eff << std::endl;
448		if(!automatized) dout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
449	<	jesup=(effp-eff)/eff*100.;
450	<	jesdown=(effm-eff)/eff*100.;
449	>	jesup=(effp-eff)/eff;
450	>	jesdown=(effm-eff)/eff;
451		}
452
453		//________________________________________________________________________
454		// Effect of energy scale on JZB efficiency
455		void doJZBscale(TTree *events, float &down, float &up, float &syst, float systematic, string informalname, bool requireZ, string addcut) {
456
457	<	TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
457	>	TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
458		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
459	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
460	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
459	>	flag_this_change(__FUNCTION__,__LINE__,true);//PlottingSetup::RestrictToMassPeak
460	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
461	>	TCut ksel(cutmass&&cutOSSF);
462
463		nBins =    50;
464		jzbMin =   0.5*jzbSel;
#	Line 317 | Line 470 | void doJZBscale(TTree *events, float &do
470		Float_t eff  = Interpolate(jzbSel,hist);
471		Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
472		Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
473	<	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
473	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
474		if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
475	<	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
476	<	up=((effp-eff)/eff)*100;
477	<	down=((effm-eff)/eff)*100;
475	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
476	>	up=((effp-eff)/eff);
477	>	down=((effm-eff)/eff);
478		}
479
480		//________________________________________________________________________
#	Line 329 | Line 482 | void doJZBscale(TTree *events, float &do
482		void JZBresponse(TTree *events, bool requireZ, float &resp, float &resperr, string addcut="",bool isMET = kFALSE, Float_t myJzbMax = 200., Int_t nPeriods = 9 ) {
483
484		jzbMin = 20;
485	<	TCut kbase("abs(genMll-91.2)<20&&genZPt>0&&genNjets>2");
485	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
486	>	TCut kbase(genMassCut&&"genZPt>0&&genNjets>2");
487		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
488	<	if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
489	<	TCut ksel("abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
488	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
489	>	if(requireZ&&PlottingSetup::RestrictToMassPeak) kbase=kbase&&"TMath::Abs(genMID)==23";
490	>	flag_this_change(__FUNCTION__,__LINE__,false);//PlottingSetup::RestrictToMassPeak
491	>	TCut ksel(cutmass&&cutOSSF);
492
493		TProfile* hJzbResp = new TProfile("hJzbResp","JZB response  ; JZB true (GeV/c); JZB reco. / JZB true", nPeriods, jzbMin, myJzbMax, "" );
494
#	Line 343 | Line 499 | void JZBresponse(TTree *events, bool req
499		hJzbResp->SetMinimum(0.2);
500		hJzbResp->Fit("pol0","Q");
501		TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
502	<	resp=fittedfunction->GetParameter(0);
503	<	resperr=fittedfunction->GetParError(0);
504	<	if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
502	>	if(!fittedfunction) {
503	>	// in case there are not enough points passing our selection
504	>	cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
505	>	resp=1;
506	>	resperr=1;
507	>	} else {
508	>	resp=fittedfunction->GetParameter(0);
509	>	resperr=fittedfunction->GetParError(0);
510	>	if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
511	>	}
512		delete hJzbResp;
513		}
514
515
516	<	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="") {
516	>	//________________________________________________________________________________________
517	>	// PDF uncertainty
518	>	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
519	>	std::vector<float> efficiency;
520	>	for(int k = 1; k < NPdfs; k++) {
521	>	float result, resulterr;
522	>	MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
523	>	efficiency.push_back(result);
524	>	}
525	>	float errHi, errLow,err;
526	>	master_formula(efficiency, errHi, errLow);
527	>	err=errLow;
528	>	if(errHi>errLow) err=errHi;
529	>	if(!automatized) dout << "  Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
530	>	return err;
531	>
532	>	}
533	>
534	>	int get_npdfs(TTree *events) {
535	>	int NPDFs;
536	>	events->SetBranchAddress("NPdfs",&NPDFs);
537	>	events->GetEntry(1);
538	>	return NPDFs;
539	>	}
540
541	+
542	+	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) {
543		float JetEnergyScaleUncert=0.1;
544		float JZBScaleUncert=0.1;
545		mcjzbexpression=mcjzb;
546	<	float triggereff=4;//percent!
546	>	float triggereff=5.0/100;// in range [0,1]
547		dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
548	<	float leptonseleff=2;//percent!
548	>	float leptonseleff=2.0/100;// in range [0,1]
549	>	leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
550		dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
551
552	+	int NPdfs=0;
553	+	if(ismSUGRA) NPdfs = get_npdfs(events);
554	+
555		float mceff,mcefferr,jzbeff,jzbefferr;
556		if(!automatized) dout << "MC efficiencies:" << endl;
557	<	MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut);
558	<	JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
557	>	Value mceff_nosigcont = MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
558	>	if(!automatized) cout << "   Without signal contamination, we find an efficiency of " << mceff_nosigcont << endl;
559	>
560	>	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
561		if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
562
563		if(!automatized) dout << "Error from Peak position:" << endl;
#	Line 380 | Line 574 | void do_systematics_for_one_file(TTree *
574
575		if(!automatized) dout << "JZB response: " << std::endl;
576		float resp,resperr;
577	<	JZBresponse(events,requireZ,resp,resperr,addcut);
577	>	if(PlottingSetup::computeJZBresponse) {
578	>	if(!automatized) dout << "JZB response: " << std::endl;
579	>	JZBresponse(events,requireZ,resp,resperr,addcut);
580	>	}
581
582		if(!automatized) dout << "Pileup: " << std::endl;
583	<	float resolution=pileup(events,requireZ,informalname,addcut);
584	<
583	>	float resolution;
584	>	resolution=pileup(events,requireZ,informalname,addcut);
585	>
586	>	float PDFuncert=0;
587	>	if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
588	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
589	>
590		dout << "_______________________________________________" << endl;
591		dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
592		if(addcut!="") dout << "With additional cut: " << addcut;
593		dout << endl;
594	<	dout << "MC efficiency: " << 100*mceff << "+/-" << 100*mcefferr << endl;
595	<	dout << "Trigger efficiency: " << triggereff << endl;
596	<	dout << "Lepton Sel Eff: " << leptonseleff << endl;
597	<	dout << "Jet energy scale: " << jesup << " " << jesdown << endl;
598	<	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl;
599	<	dout << "Resolution : " << resolution << endl;
600	<	dout << "From peak : " << sysfrompeak << endl;
601	<	dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl;
602	<	dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl;
594	>	dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
595	>	dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
596	>	dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
597	>	dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
598	>	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
599	>	dout << "Resolution : " << resolution << endl; // in range [0,1]
600	>	dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
601	>	if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
602	>	if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
603	>	if(PlottingSetup::computeJZBresponse)dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
604
605		float toterr=0;
606	<	toterr+=(triggereff/100)*(triggereff/100);
607	<	toterr+=(leptonseleff/100)*(leptonseleff/100);
608	<	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup/100)*(jesup/100); else toterr+=(jesdown/100)*(jesdown/100);
609	<	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup/100)*(scaleup/100); else toterr+=(scaledown/100)*(scaledown/100);
610	<	toterr+=(resolution/100)*(resolution/100);
611	<	toterr+=(sysfrompeak/100)*(sysfrompeak/100);
612	<	toterr=TMath::Sqrt(toterr);
613	<	dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*toterr << " (syst)   %" << endl;
614	<	dout << "     we thus use the sqrt of the sum of the squares which is : " << 100*TMath::Sqrt(mcefferr*mcefferr+(toterr*toterr)) << endl;
606	>	toterr+=(triggereff)*(triggereff);
607	>	toterr+=(leptonseleff)*(leptonseleff);
608	>	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
609	>	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
610	>	toterr+=(resolution*resolution);
611	>	toterr+=(sysfrompeak*sysfrompeak);
612	>	if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
613	>	dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
614	>	float systerr=TMath::Sqrt(toterr)*mceff;
615	>	toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
616	>
617	>	dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst)   %" << endl;
618	>	dout << "     we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
619	>	dout << "_______________________________________________" << endl;
620	>
621	>	//Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
622		vector<float> res;
623		res.push_back(jzbSel);
624		res.push_back(mceff);
625		res.push_back(mcefferr);
626		res.push_back(toterr);
627		res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
628	<	if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup/100)); else res.push_back(fabs(jesdown)/100);
629	<	if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)/100); else res.push_back(fabs(scaledown)/100);
630	<	res.push_back(fabs(resolution)/100);
628	>	if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
629	>	if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
630	>	res.push_back(fabs(resolution));
631	>	res.push_back(mceff_nosigcont.getValue());
632	>	res.push_back(mceff_nosigcont.getError());
633	>	if(ismSUGRA) res.push_back(PDFuncert);
634		results.push_back(res);
635		}
636

Diff Legend

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