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.27 by buchmann, Thu Sep 1 15:41:32 2011 UTC

#	Line 157 | Line 157 | TH1F* plotEff(TTree* events, TCut kbase,
157
158
159		//________________________________________________________________________________________
160	+	// Master Formula
161	+	void master_formula(std::vector<float> eff, float &errHi, float &errLo) {
162	+
163	+	float x0 = eff[0];
164	+	float deltaPos = 0, deltaNeg = 0;
165	+	for(int k = 0; k < (eff.size()-1)/2; k++) {
166	+	float xneg = eff[2*k+2];
167	+	float xpos = eff[2*k+1];
168	+	if(xpos-x0>0 || xneg-x0>0) {
169	+	if(xpos-x0 > xneg-x0) {
170	+	deltaPos += (xpos-x0)*(xpos-x0);
171	+	} else {
172	+	deltaPos += (xneg-x0)*(xneg-x0);
173	+	}
174	+	}
175	+	if(x0-xpos>0 || x0-xneg>0) {
176	+	if(x0-xpos > x0-xneg) {
177	+	deltaNeg += (xpos-x0)*(xpos-x0);
178	+	} else {
179	+	deltaNeg += (xneg-x0)*(xneg-x0);
180	+	}
181	+	}
182	+	}
183	+	errHi = sqrt(deltaPos);
184	+	errLo = sqrt(deltaNeg);
185	+
186	+	}
187	+
188	+
189	+	//________________________________________________________________________________________
190	+	// Get normalization factor for the PDFs
191	+	float get_norm_pdf_factor(TTree *events, int k) {
192	+
193	+	TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
194	+	char nameVar[20];
195	+	sprintf(nameVar, "pdfW[%d]", k);
196	+	events->Project("haux", nameVar);
197	+	float thisW = haux->Integral();
198	+	events->Project("haux", "pdfW[0]");
199	+	float normW = haux->Integral();
200	+
201	+	float factor=thisW/normW;
202	+
203	+	delete haux;
204	+
205	+	return factor;
206	+
207	+	}
208	+
209	+
210	+
211	+	//________________________________________________________________________________________
212		// Pile-up efficiency
213		float pileup(TTree *events, bool requireZ, string informalname, string addcut="",Float_t myJzbMax = 140. ) {
214		nBins = 16;
#	Line 182 | Line 234 | float pileup(TTree *events, bool require
234		if(!automatized) dout << "  PU: " << funcUp->Eval(jzbSel) << " " <<  func->Eval(jzbSel)
235		<< "(" << (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100. << "%)" << std::endl;
236
237	<	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel)*100.;
237	>	return (funcUp->Eval(jzbSel)-func->Eval(jzbSel))/func->Eval(jzbSel);
238
239		}
240
#	Line 218 | Line 270 | void PeakError(TTree *events,float &resu
270		else if(i==1) resdown=res;
271		else if(i==2) resup=res;
272		}
273	<	if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent)*100;
274	<	else result=(TMath::Abs(rescent-resdown)/rescent)*100;
273	>	if(TMath::Abs(rescent-resup)>TMath::Abs(rescent-resdown)) result=(TMath::Abs(rescent-resup)/rescent);
274	>	else result=(TMath::Abs(rescent-resdown)/rescent);
275		}
276
277		//____________________________________________________________________________________
278		// Total selection efficiency (MC)
279	<	void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="") {
279	>	void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
280
281		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
282		// All acceptance cuts at gen. level
#	Line 235 | Line 287 | void MCefficiency(TTree *events,float &r
287		// Corresponding reco. cuts
288		TCut ksel("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+geq_or_leq()+TString(jzbSelStr));
289		TCut ksel2("pfJetGoodNum>2&&abs(mll-91.2)<20&&id1==id2&&"+TString(mcjzb)+ngeq_or_leq()+TString("-")+TString(jzbSelStr));
290	<	events->Draw(mcjzbexpression.c_str(),kbase&&ksel,"goff");
291	<	Float_t sel = events->GetSelectedRows();
292	<	events->Draw(mcjzbexpression.c_str(),kbase&&ksel2,"goff");
293	<	Float_t nsel = events->GetSelectedRows();
290	>	TCut posSide = kbase&&ksel;
291	>	TCut negSide = kbase&&ksel2;
292	>	string sposSide(posSide);
293	>	string snegSide(negSide);
294	>	char var[20];
295	>	sprintf(var, "pdfW[%d]", k);
296	>	string svar(var);
297	>	string newPosSide = "(" + sposSide + ")*" + svar;
298	>	string newNegSide = "(" + snegSide + ")*" + svar;
299	>
300	>	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
301	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
302	>	else events->Draw((mcjzbexpression+">>effh").c_str(), sposSide.c_str(),"goff");
303	>	Float_t sel = effh->Integral();
304	>	Float_t nsel=0;
305	>	if(ConsiderSignalContaminationForLimits) {
306	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
307	>	else events->Draw((mcjzbexpression+">>effh").c_str(), snegSide.c_str(),"goff");
308	>	nsel = effh->Integral();
309	>	}
310	>	//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.
311	>	float normFactor = 1;
312	>	if(k>=0) get_norm_pdf_factor(events, k);
313	>	sel = sel/normFactor;
314	>	nsel = nsel/normFactor;
315	>
316		//      events->Draw(mcjzbexpression.c_str(),kbase,"goff");
317		//      Float_t tot = events->GetSelectedRows();
318		Float_t tot = Neventsinfile;
319
320	<	result=(sel-nsel)/tot;
321	<	resulterr=TMath::Sqrt(sel/tot*(1-sel/tot)/tot);
322	<	dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ")" << std::endl;
320	>	if(ConsiderSignalContaminationForLimits) {
321	>	result=(sel-nsel)/tot;
322	>	resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
323	>	} else {//no signal contamination considered:
324	>	result=(sel)/tot;
325	>	resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
326	>	}
327	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
328	>	delete effh;
329		}
330
331		void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
332		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;
333		if(addcut!="") kbase=kbase&&addcut.c_str();//this is mostly for SUSY scans (adding requirements on masses)
334		if(requireZ) kbase=kbase&&"TMath::Abs(genMID)==23";
335		TH1F* hLM4 = plotEff(events,kbase,informalname);
#	Line 294 | Line 373 | void JZBjetScale(TTree *events, float &j
373		if(!automatized) dout << "    JESup: " << effp << " (" << (effp-eff)/eff*100. << "%)" << std::endl;
374		if(!automatized) dout << "    central:  " << eff << std::endl;
375		if(!automatized) dout << "    JESdown: " << effm << " (" << (effm-eff)/eff*100. << "%)" << std::endl;
376	<	jesup=(effp-eff)/eff*100.;
377	<	jesdown=(effm-eff)/eff*100.;
376	>	jesup=(effp-eff)/eff;
377	>	jesdown=(effm-eff)/eff;
378		}
379
380		//________________________________________________________________________
#	Line 320 | Line 399 | void doJZBscale(TTree *events, float &do
399		if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
400		if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
401		if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
402	<	up=((effp-eff)/eff)*100;
403	<	down=((effm-eff)/eff)*100;
402	>	up=((effp-eff)/eff);
403	>	down=((effm-eff)/eff);
404		}
405
406		//________________________________________________________________________
#	Line 350 | Line 429 | void JZBresponse(TTree *events, bool req
429		}
430
431
432	<	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="") {
432	>	//________________________________________________________________________________________
433	>	// PDF uncertainty
434	>	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
435	>	std::vector<float> efficiency;
436	>	for(int k = 1; k < NPdfs; k++) {
437	>	float result, resulterr;
438	>	MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
439	>	efficiency.push_back(result);
440	>	}
441	>	float errHi, errLow,err;
442	>	master_formula(efficiency, errHi, errLow);
443	>	err=errLow;
444	>	if(errHi>errLow) err=errHi;
445	>	if(!automatized) dout << "  Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
446	>	return err;
447	>
448	>	}
449	>
450	>	int get_npdfs(TTree *events) {
451	>	int NPDFs;
452	>	events->SetBranchAddress("NPdfs",&NPDFs);
453	>	events->GetEntry(1);
454	>	return NPDFs;
455	>	}
456
457	+
458	+	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) {
459		float JetEnergyScaleUncert=0.1;
460		float JZBScaleUncert=0.1;
461		mcjzbexpression=mcjzb;
462	<	float triggereff=4;//percent!
462	>	float triggereff=4.0/100;// in range [0,1]
463		dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
464	<	float leptonseleff=2;//percent!
464	>	float leptonseleff=2.0/100;// in range [0,1]
465		dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
466
467	+	int NPdfs=0;
468	+	if(ismSUGRA) NPdfs = get_npdfs(events);
469	+
470		float mceff,mcefferr,jzbeff,jzbefferr;
471		if(!automatized) dout << "MC efficiencies:" << endl;
472	<	MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut);
472	>	MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
473		JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
474		if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
475
#	Line 384 | Line 491 | void do_systematics_for_one_file(TTree *
491
492		if(!automatized) dout << "Pileup: " << std::endl;
493		float resolution=pileup(events,requireZ,informalname,addcut);
494	<
494	>
495	>	float PDFuncert=0;
496	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
497	>
498		dout << "_______________________________________________" << endl;
499		dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
500		if(addcut!="") dout << "With additional cut: " << addcut;
501		dout << endl;
502	<	dout << "MC efficiency: " << 100*mceff << "+/-" << 100*mcefferr << endl;
503	<	dout << "Trigger efficiency: " << triggereff << endl;
504	<	dout << "Lepton Sel Eff: " << leptonseleff << endl;
505	<	dout << "Jet energy scale: " << jesup << " " << jesdown << endl;
506	<	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl;
507	<	dout << "Resolution : " << resolution << endl;
508	<	dout << "From peak : " << sysfrompeak << endl;
509	<	dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl;
510	<	dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl;
502	>	dout << "MC efficiency: " << mceff << "+/-" << mcefferr << endl; // in range [0,1]
503	>	dout << "Trigger efficiency: " << triggereff << endl; // in range [0,1]
504	>	dout << "Lepton Sel Eff: " << leptonseleff << endl; // in range [0,1]
505	>	dout << "Jet energy scale: " << jesup << " " << jesdown << endl; // in range [0,1]
506	>	dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
507	>	dout << "Resolution : " << resolution << endl; // in range [0,1]
508	>	dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
509	>	if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
510	>	dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
511	>	dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
512
513		float toterr=0;
514	<	toterr+=(triggereff/100)*(triggereff/100);
515	<	toterr+=(leptonseleff/100)*(leptonseleff/100);
516	<	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup/100)*(jesup/100); else toterr+=(jesdown/100)*(jesdown/100);
517	<	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup/100)*(scaleup/100); else toterr+=(scaledown/100)*(scaledown/100);
518	<	toterr+=(resolution/100)*(resolution/100);
519	<	toterr+=(sysfrompeak/100)*(sysfrompeak/100);
520	<	toterr=TMath::Sqrt(toterr);
521	<	dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*toterr << " (syst)   %" << endl;
522	<	dout << "     we thus use the sqrt of the sum of the squares which is : " << 100*TMath::Sqrt(mcefferr*mcefferr+(toterr*toterr)) << endl;
514	>	toterr+=(triggereff)*(triggereff);
515	>	toterr+=(leptonseleff)*(leptonseleff);
516	>	if(fabs(jesup)>fabs(jesdown)) toterr+=(jesup*jesup); else toterr+=(jesdown*jesdown);
517	>	if(fabs(scaleup)>fabs(scaledown)) toterr+=(scaleup*scaleup); else toterr+=(scaledown*scaledown);
518	>	toterr+=(resolution*resolution);
519	>	toterr+=(sysfrompeak*sysfrompeak);
520	>	if(ismSUGRA) toterr+=(PDFuncert*PDFuncert);
521	>	dout << "TOTAL SYSTEMATICS: " << TMath::Sqrt(toterr) << " --> " << TMath::Sqrt(toterr)*mceff << endl;
522	>	float systerr=TMath::Sqrt(toterr)*mceff;
523	>	toterr=TMath::Sqrt(toterr*mceff*mceff+mcefferr*mcefferr);//also includes stat err!
524	>
525	>	dout << "FINAL RESULT : " << 100*mceff << " +/- "<< 100*mcefferr << " (stat) +/- " << 100*systerr << " (syst)   %" << endl;
526	>	dout << "     we thus use the sqrt of the sum of the squares of the stat & syst err, which is : " << 100*toterr << endl;
527	>
528	>	//Do not modify the lines below or mess with the order; this order is expected by all limit calculating functions!
529		vector<float> res;
530		res.push_back(jzbSel);
531		res.push_back(mceff);
532		res.push_back(mcefferr);
533		res.push_back(toterr);
534		res.push_back(TMath::Sqrt((mcefferr)*(mcefferr)+(toterr*toterr)));
535	<	if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup/100)); else res.push_back(fabs(jesdown)/100);
536	<	if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)/100); else res.push_back(fabs(scaledown)/100);
537	<	res.push_back(fabs(resolution)/100);
535	>	if(fabs(jesup)>fabs(jesdown)) res.push_back(fabs(jesup)); else res.push_back(fabs(jesdown));
536	>	if(fabs(scaleup)>fabs(scaledown)) res.push_back(fabs(scaleup)); else res.push_back(fabs(scaledown));
537	>	res.push_back(fabs(resolution));
538	>	if(ismSUGRA) res.push_back(PDFuncert);
539		results.push_back(res);
540		}
541

Diff Legend

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