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.23 by pablom, Wed Aug 31 10:21:17 2011 UTC vs.
Revision 1.33 by buchmann, Mon Sep 26 15:09:41 2011 UTC

#	Line 1 | Line 1
1		#include <iostream>
2		#include <vector>
3		#include <sys/stat.h>
4	+	#include <algorithm>
5	+	#include <cmath>
6
7		#include <TMath.h>
8		#include <TColor.h>
#	Line 190 | Line 192 | void master_formula(std::vector<float> e
192		// Get normalization factor for the PDFs
193		float get_norm_pdf_factor(TTree *events, int k) {
194
195	<	TH1F *haux = new TH1F("haux", "", 1000, 0, 15);
195	>	TH1F *haux = new TH1F("haux", "", 10000, 0, 5);
196		char nameVar[20];
197		sprintf(nameVar, "pdfW[%d]", k);
198		events->Project("haux", nameVar);
199	<	float factor = haux->Integral();
199	>	float thisW = haux->Integral();
200	>	events->Project("haux", "pdfW[0]");
201	>	float normW = haux->Integral();
202	>
203	>	float factor=thisW/normW;
204
205		delete haux;
206
#	Line 272 | Line 278 | void PeakError(TTree *events,float &resu
278
279		//____________________________________________________________________________________
280		// Total selection efficiency (MC)
281	<	void MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
281	>	//returns the efficiency WITHOUT signal contamination, and the result and resulterr contain the result and the corresponding error
282	>	Value MCefficiency(TTree *events,float &result, float &resulterr,string mcjzb,bool requireZ,int Neventsinfile, string addcut="", int k = 0) {
283
284		char jzbSelStr[256]; sprintf(jzbSelStr,"%f",jzbSel);
285		// All acceptance cuts at gen. level
#	Line 293 | Line 300 | void MCefficiency(TTree *events,float &r
300		string newPosSide = "(" + sposSide + ")*" + svar;
301		string newNegSide = "(" + snegSide + ")*" + svar;
302
303	<	events->Draw(mcjzbexpression.c_str(), newPosSide.c_str(),"goff");
304	<	Float_t sel = events->GetSelectedRows();
303	>	TH1F *effh= new TH1F("effh","effh",1,-14000,14000);
304	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newPosSide.c_str(),"goff");
305	>	else events->Draw((mcjzbexpression+">>effh").c_str(), sposSide.c_str(),"goff");
306	>	Float_t sel = effh->Integral();
307		Float_t nsel=0;
308		if(ConsiderSignalContaminationForLimits) {
309	<	events->Draw(mcjzbexpression.c_str(), newNegSide.c_str(),"goff");
310	<	nsel = events->GetSelectedRows();
309	>	if(k>=0)events->Draw((mcjzbexpression+">>effh").c_str(), newNegSide.c_str(),"goff");
310	>	else events->Draw((mcjzbexpression+">>effh").c_str(), snegSide.c_str(),"goff");
311	>	nsel = effh->Integral();
312		}
313		//Corrections due to normalization in the PDF. This has to be applied as well to the number of events in a file if the definition changes at some point.
314	<	float normFactor = get_norm_pdf_factor(events, k);
314	>	float normFactor = 1;
315	>	if(k>=0) get_norm_pdf_factor(events, k);
316		sel = sel/normFactor;
317		nsel = nsel/normFactor;
318
#	Line 309 | Line 320 | void MCefficiency(TTree *events,float &r
320		//      Float_t tot = events->GetSelectedRows();
321		Float_t tot = Neventsinfile;
322
323	+	Value result_wo_signalcont;
324	+
325		if(ConsiderSignalContaminationForLimits) {
326		result=(sel-nsel)/tot;
327		resulterr=(1.0/tot)*TMath::Sqrt(sel+nsel+(sel-nsel)*(sel-nsel)/tot);
328	+	result_wo_signalcont=Value(sel/tot,TMath::Sqrt(sel/tot*(1+sel/tot)/tot));
329		} else {//no signal contamination considered:
330		result=(sel)/tot;
331		resulterr=TMath::Sqrt(sel/tot*(1+sel/tot)/tot);
332	+	result_wo_signalcont=Value(result,resulterr);
333		}
334	<	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ")" << std::endl;
334	>	if(!automatized && k>0 ) dout << "PDF assessment: ";
335	>	if(!automatized) dout << "  MC efficiency: " << result << "+-" << resulterr << "  ( JZB>" << jzbSel << " : " << sel << " , JZB<-" << jzbSel << " : " << nsel << " and nevents=" << tot << ") with normFact=" << normFactor << std::endl;
336	>	delete effh;
337	>	}
338	>
339	>
340	>	//____________________________________________________________________________________
341	>	// Selection efficiency for one process (MC)
342	>	vector<float> processMCefficiency(TTree *events,string mcjzb,bool requireZ,int Neventsinfile, string addcut) {
343	>	vector<float> process_efficiencies;
344	>	for(int iprocess=0;iprocess<=10;iprocess++) {
345	>	float this_process_efficiency,efferr;
346	>	stringstream addcutplus;
347	>	addcutplus<<addcut<<"&&(process=="<<iprocess<<")";
348	>	MCefficiency(events,this_process_efficiency, efferr,mcjzb,requireZ,Neventsinfile, addcutplus.str(),-1);
349	>	process_efficiencies.push_back(this_process_efficiency);
350	>	}
351	>	return process_efficiencies;
352		}
353	+
354
355		void JZBefficiency(TTree *events, string informalname, float &jzbeff, float &jzbefferr, bool requireZ, string addcut="") {
356		TCut kbase("abs(genMll-91.2)<20&&genNjets>2&&genZPt>0&&abs(mll-91.2)<20&&((id1+1)*(id2+1)*ch1*ch2)!=-2");
#	Line 387 | Line 420 | void doJZBscale(TTree *events, float &do
420		Float_t eff  = Interpolate(jzbSel,hist);
421		Float_t effp = Interpolate(jzbSel*(1.+systematic),hist);
422		Float_t effm = Interpolate(jzbSel*(1.-systematic),hist);
423	<	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<syst*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
423	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.+systematic)  << "(-"<<systematic*100<<"%) : " << effp << " (" << ((effp-eff)/eff)*100. << "%)"  << std::endl;
424		if(!automatized) dout << "  efficiency at JZB==" << jzbSel  << ": " << eff << std::endl;
425	<	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<syst*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
425	>	if(!automatized) dout << "  efficiency at JZB==" << jzbSel*(1.-systematic)  << "(-"<<systematic*100<<"%) : " << effm << " (" << ((effm-eff)/eff)*100. << "%)"  << std::endl;
426		up=((effp-eff)/eff);
427		down=((effm-eff)/eff);
428		}
#	Line 413 | Line 446 | void JZBresponse(TTree *events, bool req
446		hJzbResp->SetMinimum(0.2);
447		hJzbResp->Fit("pol0","Q");
448		TF1 *fittedfunction = hJzbResp->GetFunction("pol0");
449	<	resp=fittedfunction->GetParameter(0);
450	<	resperr=fittedfunction->GetParError(0);
451	<	if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
449	>	if(!fittedfunction) {
450	>	// in case there are not enough points passing our selection
451	>	cout << "OOPS response function invalid, assuming 100% error !!!!" << endl;
452	>	resp=1;
453	>	resperr=1;
454	>	} else {
455	>	resp=fittedfunction->GetParameter(0);
456	>	resperr=fittedfunction->GetParError(0);
457	>	if(!automatized) dout << "  Response: " << resp << " +/- " << resperr << endl;
458	>	}
459		delete hJzbResp;
460		}
461
462
463		//________________________________________________________________________________________
464		// PDF uncertainty
465	<	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, string addcut="") {
426	<
427	<	int numberOfPDFs = 44;
465	>	float get_pdf_uncertainty(TTree *events, string mcjzb, bool requireZ, int Neventsinfile, int NPdfs, string addcut="") {
466		std::vector<float> efficiency;
467	<	for(int k = 1; k < numberOfPDFs; k++) {
467	>	for(int k = 1; k < NPdfs; k++) {
468		float result, resulterr;
469		MCefficiency(events, result, resulterr, mcjzb, requireZ, Neventsinfile, addcut, k);
470		efficiency.push_back(result);
471		}
472	<	float errHi, errLow;
472	>	float errHi, errLow,err;
473		master_formula(efficiency, errHi, errLow);
474	<	if(errHi>errLow) return errHi;
475	<	return errLow;
474	>	err=errLow;
475	>	if(errHi>errLow) err=errHi;
476	>	if(!automatized) dout << "  Uncertainty from PDF: " << errLow << " (low) and " << errHi << "(high) ---> Picked " << err << endl;
477	>	return err;
478
479		}
480
481	<
481	>	int get_npdfs(TTree *events) {
482	>	int NPDFs;
483	>	events->SetBranchAddress("NPdfs",&NPDFs);
484	>	events->GetEntry(1);
485	>	return NPDFs;
486	>	}
487	>
488
489		void do_systematics_for_one_file(TTree *events,int Neventsinfile,string informalname, vector<vector<float> > &results,string mcjzb,string datajzb,float peakerror,bool requireZ=false, string addcut="", bool ismSUGRA=false) {
444	–
490		float JetEnergyScaleUncert=0.1;
491		float JZBScaleUncert=0.1;
492		mcjzbexpression=mcjzb;
493	<	float triggereff=4.0/100;// in range [0,1]
493	>	float triggereff=5.0/100;// in range [0,1]
494		dout << "Trigger efficiency not implemented in this script  yet, still using external one" << endl;
495		float leptonseleff=2.0/100;// in range [0,1]
496	+	leptonseleff=TMath::Sqrt(leptonseleff*leptonseleff+leptonseleff*leptonseleff); // because the 2% is per lepton
497		dout << "Lepton selection efficiency not implemented in this script  yet, still using external one" << endl;
498
499	+	int NPdfs=0;
500	+	if(ismSUGRA) NPdfs = get_npdfs(events);
501	+
502		float mceff,mcefferr,jzbeff,jzbefferr;
503		if(!automatized) dout << "MC efficiencies:" << endl;
504	<	MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut);
505	<	JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
504	>	MCefficiency(events,mceff,mcefferr,mcjzb,requireZ,Neventsinfile,addcut,-1);
505	>	if(PlottingSetup::computeJZBefficiency) JZBefficiency(events,informalname,jzbeff,jzbefferr,requireZ,addcut);
506		if(!automatized) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << endl;
507
508		if(!automatized) dout << "Error from Peak position:" << endl;
#	Line 470 | Line 519 | void do_systematics_for_one_file(TTree *
519
520		if(!automatized) dout << "JZB response: " << std::endl;
521		float resp,resperr;
522	<	JZBresponse(events,requireZ,resp,resperr,addcut);
522	>	if(PlottingSetup::computeJZBresponse) {
523	>	if(!automatized) dout << "JZB response: " << std::endl;
524	>	JZBresponse(events,requireZ,resp,resperr,addcut);
525	>	}
526
527		if(!automatized) dout << "Pileup: " << std::endl;
528	<	float resolution=pileup(events,requireZ,informalname,addcut);
528	>	float resolution;
529	>	resolution=pileup(events,requireZ,informalname,addcut);
530
531		float PDFuncert=0;
532	<	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, addcut);
532	>	if(!automatized) dout << "Assessing PDF uncertainty: " << std::endl;
533	>	if(ismSUGRA) PDFuncert = get_pdf_uncertainty(events, mcjzb, requireZ, Neventsinfile, NPdfs, addcut);
534
535		dout << "_______________________________________________" << endl;
536		dout << "                 SUMMARY FOR " << informalname << " with JZB>" << jzbSel << "  (all in %) ";
#	Line 489 | Line 543 | void do_systematics_for_one_file(TTree *
543		dout << "JZB Scale Uncert: " << scaledown << " " << scaleup << endl; // in range [0,1]
544		dout << "Resolution : " << resolution << endl; // in range [0,1]
545		dout << "From peak : " << sysfrompeak << endl; // in range [0,1]
546	<	dout << "PDF uncertainty  : " << PDFuncert << " (not yet included below) " << endl; // in range [0,1]
547	<	dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
548	<	dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
546	>	if(ismSUGRA) dout << "PDF uncertainty  : " << PDFuncert << endl; // in range [0,1]
547	>	if(PlottingSetup::computeJZBefficiency) dout << "JZB efficiency: " << jzbeff << "+/-" << jzbefferr << " (not yet included below) " << endl; // in range [0,1]
548	>	if(PlottingSetup::computeJZBresponse)dout << "JZB response  : " << resp << " +/-" << resperr << " (not yet included below) " << endl; // in range [0,1]
549
550		float toterr=0;
551		toterr+=(triggereff)*(triggereff);

Diff Legend

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