propagate systematics updates

MakeAllDCsyst.C

@@ -197,6 +197,8 @@ class KSystProcessor {
 			auto ssyst = MakeStat(sname, stat_yield);
 			hsyst.push_back(ssyst);
 			pctDiffMap.Add("MCStatErr",fabs(1-stat_yield/nominal_yield)*100);
+			//average rather than max stat err
+			pctDiffMap.Add("MCStatOverallErr",MakeStatOverall(nominal)*100);
 
 			//genMHT correction and unc for fastsim
 			if(genMHT){
@@ -258,6 +260,7 @@ class KSystProcessor {
 			);
 		}
 		//functions to be specialized
+		double MakeStatOverall(T* hist);
 		double Integral(T* hist);
 		void DivideBound(T* hist, const string& binname);
 		T* MakeStat(const string& sname, double& yield);
@@ -271,6 +274,18 @@ class KSystProcessor {
 		vector<T*> hsyst;
 };
 
+template <> double KSystProcessor<TH1F>::MakeStatOverall(TH1F* hist) {
+	double err = 0.;
+	double hint = hist->IntegralAndError(0,hist->GetNbinsX()+1,err);
+	return err/hint;
+}
+
+template <> double KSystProcessor<TH2F>::MakeStatOverall(TH2F* hist) {
+	double err = 0.;
+	double hint = hist->IntegralAndError(0,hist->GetNbinsX()+1,0,hist->GetNbinsY()+1,err);
+	return err/hint;
+}
+
 template<> double KSystProcessor<TH1F>::Integral(TH1F* hist){ return hist->Integral(0,hist->GetNbinsX()+1); }
 template<> double KSystProcessor<TH2F>::Integral(TH2F* hist){ return hist->Integral(0,hist->GetNbinsX()+1,0,hist->GetNbinsY()+1); }
 

README.md

@@ -260,16 +260,16 @@ root -b -l -q 'KPlotDriver.C+("root://cmseos.fnal.gov//store/user/lpcsusyhad/SVJ
 To create histograms for datacards and do post-processing:
 ```
 cd batch
-./DCsub_svj.sh -t SVJData,SVJBkg -y Run2 -v trig5/sigfull
+./DCsub_svj.sh -t SVJData,SVJBkg -y Run2 -v trig6/sigfull
 (wait for jobs to finish)
-../finalizeDatacardsSVJ.sh -d /store/user/pedrok/SVJ2017/Datacards/trig5/sigfull -k
-./DCsub_svj.sh -t SVJScan -y Run2 -v trig5/sigfull
+../finalizeDatacardsSVJ.sh -d /store/user/pedrok/SVJ2017/Datacards/trig6/sigfull -k
+./DCsub_svj.sh -t SVJScan -y Run2 -v trig6/sigfull
 ```
 
 Signal systematics studies:
 ```
 cd batch
-./HSsub.sh -c block,scan -d /store/user/pedrok/SVJ2017/Datacards/trig5/sigfull -S scan -k -n 1
+./HSsub.sh -c block,scan -d /store/user/pedrok/SVJ2017/Datacards/trig6/sigfull -S scan -k -n 1
 cd ..
-python summarizeSyst.py -d /store/user/pedrok/SVJ2017/Datacards/trig5/sigfull -c "mDark==20.&&rinv==0.3&&alpha==-2."
+python summarizeSyst.py -d /store/user/pedrok/SVJ2017/Datacards/trig6/sigfull -c "mDark==20.&&rinv==0.3&&alpha==-2."
 ```

batch/DCsub_svj.sh

@@ -79,7 +79,7 @@ for TYPE in ${TYPES[@]}; do
 	fi
 
 	if [[ -z "$NOSYST" && ("$TYPE" == SVJ  ||  "$TYPE" == "SVJScan"*) ]]; then
-		SYSTS=${NOMINAL},trigfnuncUp,trigfnuncDown,puuncUp,puuncDown,pdfalluncUp,pdfalluncDown,psisruncUp,psisruncDown,psfsruncUp,psfsruncDown
+		SYSTS=${NOMINAL},trigfnuncUp,trigfnuncDown,puuncUp,puuncDown,pdfuncUp,pdfuncDown,psisruncUp,psisruncDown,psfsruncUp,psfsruncDown
 		VARS=JECup,JECdown,JERup,JERdown,JESup,JESdown
 	else
 		SYSTS=${NOMINAL}

input/input_selection_svj_syst_2016.txt

@@ -49,6 +49,12 @@ pdfalluncUp
 pdfalluncDown
 	MCWeight	in:options[input/input_svj_MCWeight_options_2016.txt]	int:pdfallunc[-1]
 	common
+pdfuncUp
+	MCWeight	in:options[input/input_svj_MCWeight_options_2016.txt]	int:pdfunc[1]
+	common
+pdfuncDown
+	MCWeight	in:options[input/input_svj_MCWeight_options_2016.txt]	int:pdfunc[-1]
+	common
 psisruncUp
 	MCWeight	in:options[input/input_svj_MCWeight_options_2016.txt]	b:psnorm[0]	int:psisrunc[1]
 	common

input/input_selection_svj_syst_2017.txt

@@ -49,6 +49,12 @@ pdfalluncUp
 pdfalluncDown
 	MCWeight	in:options[input/input_svj_MCWeight_options_2017.txt]	int:pdfallunc[-1]
 	common
+pdfuncUp
+	MCWeight	in:options[input/input_svj_MCWeight_options_2017.txt]	int:pdfunc[1]
+	common
+pdfuncDown
+	MCWeight	in:options[input/input_svj_MCWeight_options_2017.txt]	int:pdfunc[-1]
+	common
 psisruncUp
 	MCWeight	in:options[input/input_svj_MCWeight_options_2017.txt]	b:psnorm[0]	int:psisrunc[1]
 	common

input/input_selection_svj_syst_2018POST.txt

@@ -50,6 +50,12 @@ pdfalluncUp
 pdfalluncDown
 	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	int:pdfallunc[-1]
 	common
+pdfuncUp
+	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	int:pdfunc[1]
+	common
+pdfuncDown
+	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	int:pdfunc[-1]
+	common
 psisruncUp
 	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	b:psnorm[0]	int:psisrunc[1]
 	common

input/input_selection_svj_syst_2018PRE.txt

@@ -49,6 +49,12 @@ pdfalluncUp
 pdfalluncDown
 	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	int:pdfallunc[-1]
 	common
+pdfuncUp
+	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	int:pdfunc[1]
+	common
+pdfuncDown
+	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	int:pdfunc[-1]
+	common
 psisruncUp
 	MCWeight	in:options[input/input_svj_MCWeight_options_2018.txt]	b:psnorm[0]	int:psisrunc[1]
 	common

input/input_svj_ext_sig_syst.txt

@@ -36,9 +36,9 @@ hist	mc	JERup	s:legname[JER up]	c:color[kBlue]	i:linestyle[7]
 hist	mc	JERdown	s:legname[JER down]	c:color[kRed]	i:linestyle[7]
 	base	ext	SVJ_3100_20_0.3_peak_MC2017	vs:exthisto_in[SVJ_mZprime3100_mDark20_rinv03_alphapeak_MC2017JERDown]	vs:exthisto_out[MTAK8]
 hist	mc	PDFup	s:legname[PDF up]	c:color[kBlue]	i:linestyle[7]
-	base	ext	SVJ_3100_20_0.3_peak_MC2017	vs:exthisto_in[SVJ_mZprime3100_mDark20_rinv03_alphapeak_MCRun2pdfalluncUp]	vs:exthisto_out[MTAK8]
+	base	ext	SVJ_3100_20_0.3_peak_MC2017	vs:exthisto_in[SVJ_mZprime3100_mDark20_rinv03_alphapeak_MCRun2pdfuncUp]	vs:exthisto_out[MTAK8]
 hist	mc	PDFdown	s:legname[PDF down]	c:color[kRed]	i:linestyle[7]
-	base	ext	SVJ_3100_20_0.3_peak_MC2017	vs:exthisto_in[SVJ_mZprime3100_mDark20_rinv03_alphapeak_MCRun2pdfalluncDown]	vs:exthisto_out[MTAK8]
+	base	ext	SVJ_3100_20_0.3_peak_MC2017	vs:exthisto_in[SVJ_mZprime3100_mDark20_rinv03_alphapeak_MCRun2pdfuncDown]	vs:exthisto_out[MTAK8]
 hist	mc	PUup	s:legname[PU up]	c:color[kBlue]	i:linestyle[7]
 	base	ext	SVJ_3100_20_0.3_peak_MC2017	vs:exthisto_in[SVJ_mZprime3100_mDark20_rinv03_alphapeak_MC2017puuncUp]	vs:exthisto_out[MTAK8]
 hist	mc	PUdown	s:legname[PU down]	c:color[kRed]	i:linestyle[7]

scripts/plotSyst2D.C

@@ -16,6 +16,8 @@
 #include <iostream>
 #include <iomanip>
 #include <sstream>
+#include <utility>
+#include <tuple>
 
 using namespace std;
 
@@ -24,6 +26,7 @@ map<string,string> legnames = {
 	{"JEC", "JEC uncertainty [%]"},
 	{"JER", "JER uncertainty [%]"},
 	{"MCStatErr", "MC statistical uncertainty [%]"},
+	{"MCStatOverallErr", "MC statistical uncertainty [%]"},
 	{"MHTSyst", "H_{T}^{miss} modeling uncertainty [%]"},
 	{"puaccunc", "pileup acceptance uncertainty [%]"},
 	{"puunc", "pileup reweighting uncertainty [%]"},
@@ -33,6 +36,7 @@ map<string,string> legnames = {
 	{"prefireunc","L1 prefiring weight uncertainty [%]"},
 	{"hemvetounc","HEM veto uncertainty [%]"},
 	{"contam","SL contamination [%]"},
+	{"pdfunc", "PDF uncertainty [%]"},
 	{"pdfallunc", "PDF uncertainty [%]"},
 	{"psisrunc", "PS ISR uncertainty [%]"},
 	{"psfsrunc", "PS FSR uncertainty [%]"},
@@ -49,36 +53,59 @@ map<string,string> legnames = {
 	{"alpha", "#alpha_{dark}"},
 };
 
-void getBinsQty(string rqty, TTree* tree, double& min, double& max, int& nbins, bool isMassPlane){
+vector<double> getVariableBins(const vector<double>& vals){
+	vector<double> bins; bins.reserve(vals.size()+1);
+	for(unsigned i = 1; i < vals.size(); ++i){
+		//midpoint between bin centers
+		double midpt = (vals[i]+vals[i-1])/2.;
+		//starting point for first bin
+		if(i==1) bins.push_back(2*vals[i-1]-midpt);
+		bins.push_back(midpt);
+		//ending point for last bin
+		if(i==vals.size()-1) bins.push_back(2*vals[i]-midpt);
+	}
+	return bins;
+}
+
+pair<bool,vector<double>> getBinsQty(string rqty, TTree* tree, bool isMassPlane){
+	//fixed-width case: min, max, nbins
+	//variable-width case: bin boundaries (length nbins+1)
+	vector<double> bin_info;
+	bool fixed = true;
 	if(rqty=="mZprime"){
-		min = 1400; max = 5200; nbins = 19;
+		bin_info = {1400,5200,19};
 	}
 	else if(rqty=="mDark"){
-		min = -4; max = 106; nbins = 11;
+		vector<double> vals{1,5,10,20,30,40,50,60,70,80,90,100};
+		bin_info = getVariableBins(vals);
+		fixed = false;
 	}
 	else if(rqty=="rinv"){
-		min = -0.05; max = 1.05; nbins = 11;
+		vector<double> vals{0,0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0};
+		bin_info = getVariableBins(vals);
+		fixed = false;
 	}
 	else if(rqty=="alpha"){
-		min = -3.5; max = -0.5; nbins = 3;
+		bin_info = {-3.5,-0.5,3};
 	}
 	else {
-		min = 0; max = tree->GetMaximum(rqty.c_str()); nbins = 40;
+		bin_info = {0, tree->GetMaximum(rqty.c_str()), 40};
 	}
 
 	if(isMassPlane){
-		min = tree->GetMinimum(rqty.c_str());
-		nbins = (max-min)/25;
+		bin_info[0] = tree->GetMinimum(rqty.c_str());
+		bin_info[2] = (bin_info[1]-bin_info[0])/25;
 	}
+
+	return make_pair(fixed,bin_info);
 }
 
 //make 2D plot
 void plotSyst(TTree* tree, string xqty, string yqty, string model, string outdir, vector<string> printformats){
 	//get ranges (use mMother for derived quantity deltaM)
 	string rxqty = (xqty=="deltaM") ? "mMother" : xqty;
 	string ryqty = (yqty=="deltaM") ? "mMother" : yqty;
-	int nbinsx = 0, nbinsy = 0;
-	double xmin = 0, xmax = 0, ymin = 0, ymax = 0;
+	vector<double> binsx, binsy;
 	bool isPlane = false;
 	bool isMassPlane = false;
 	string cutPlane = "";
@@ -106,14 +133,31 @@ void plotSyst(TTree* tree, string xqty, string yqty, string model, string outdir
 		cutPlane = "*(rinv==0.3&&mDark==20)";
 	}
 
-	getBinsQty(rxqty, tree, xmin, xmax, nbinsx, isMassPlane);
-	getBinsQty(ryqty, tree, ymin, ymax, nbinsy, isMassPlane);
+	bool fixedx,fixedy;
+	tie(fixedx,binsx) = getBinsQty(rxqty, tree, isMassPlane);
+	tie(fixedy,binsy) = getBinsQty(ryqty, tree, isMassPlane);
 
 	//make name
 	string hname = "plotSyst2D_"+yqty+"_vs_"+xqty+"_"+model;
 
 	//initialize histos
-	TH2F* hbase = new TH2F("hbase","",nbinsx,xmin,xmax,nbinsy,ymin,ymax);
+	TH2F* hbase = nullptr;
+	if(fixedx){
+		if(fixedy){
+			hbase = new TH2F("hbase","",binsx[2],binsx[0],binsx[1],binsy[2],binsy[0],binsy[1]);
+		}
+		else {
+			hbase = new TH2F("hbase","",binsx[2],binsx[0],binsx[1],binsy.size()-1,binsy.data());
+		}
+	}
+	else {
+		if(fixedy){
+			hbase = new TH2F("hbase","",binsx.size()-1,binsx.data(),binsy[2],binsy[0],binsy[1]);
+		}
+		else {
+			hbase = new TH2F("hbase","",binsx.size()-1,binsx.data(),binsy.size()-1,binsy.data());
+		}
+	}
 	if(isPlane){
 		hbase->GetXaxis()->SetTitle(legnames[rxqty].c_str());
 		hbase->GetYaxis()->SetTitle(legnames[ryqty].c_str());

scripts/processDatacardsSVJ.py

@@ -192,6 +192,8 @@ def write(self):
 
     # hardcoded list of correlated systematics
     correl_systs = set([
+        "pdfuncUp",
+        "pdfuncDown",
         "pdfalluncUp",
         "pdfalluncDown",
         "psisruncUp",

summarizeSyst.py

@@ -10,6 +10,7 @@
     "jetidunc": "Jet ID",
     "lumiunc": "Luminosity",
     "MCStatErr": "MC statistical",
+    "MCStatOverallErr": "MC statistical",
     "MHTSyst": "\\MHT modeling",
     "prefireunc": "Prefiring weight",
     "puunc": "Pileup reweighting",
@@ -18,6 +19,7 @@
     "trigsystunc": "Trigger systematic",
 # SVJ uncs below
     "pdfallunc": "PDF",
+    "pdfunc": "PDF",
     "psfsrunc": "FSR (parton shower)",
     "psisrunc": "ISR (parton shower)",
     "trigfnunc": "Trigger statistical",
@@ -124,7 +126,7 @@ def findMinMax(fname,dir,exclude,include,cut,flat,output):
     parser.add_argument("-d", "--dir", dest="dir", type=str, default="/store/user/pedrok/SUSY2015/Analysis/Datacards/Run2ProductionV12", help="directory (LFN) of systematics files")
     parser.add_argument("-k", "--skip", dest="skip", type=str, default=[], nargs='*', help="list of models to skip")
     parser.add_argument("-n", "--include", dest="include", type=str, default=[], nargs='*', help="list of systematics to include in calculations (empty = all)")
-    parser.add_argument("-x", "--exclude", dest="exclude", type=str, default=["contam"], help="list of systematics to exclude from calculations")
+    parser.add_argument("-x", "--exclude", dest="exclude", type=str, default=["contam"], nargs='*', help="list of systematics to exclude from calculations")
     parser.add_argument("-m", "--minimum", dest="minimum", type=float, default=0.01, help="minimum value to display, smaller values rounded to 0")
     parser.add_argument("-g", "--grep", dest="grep", type=str, default="", help="select only models matching this string")
     parser.add_argument("-c", "--cut", dest="cut", type=str, default="", help="apply cut when computing syst min/max")