Transition to80 x

Makefile

@@ -41,7 +41,6 @@ LIB=-L$(BOOST)/lib -L/usr/lib64 # -L/usr/lib
 
 #### Make the list of modules from the list of .cc files in the SRC directory
 MODULES=$(shell ls $(SRCDIR)/*.cc | sed "s|.cc|.o|;s|$(SRCDIR)|$(OBJ_DIR)|g")
-MODULESZFitter=../../Calibration/ZFitter/lib/EnergyScaleCorrection_class.o
 #### Make the list of dependencies for a particular module
 
 default: signalPdf.exe $(BUILDDIR)/analysis
@@ -57,8 +56,6 @@ lib/%.o: $(SRCDIR)/%.cc
 	@echo "--> Making $@" 
 	@$(COMPILE.cc) $(CXXFLAGS) $(INCLUDE) $(MAKEDEPEND) -o $@ $<
 
-scales: ../../Calibration/ZFitter/lib/EnergyScaleCorrection_class.o
-../../Calibration/ZFitter/lib/EnergyScaleCorrection_class.o: ../../Calibration/ZFitter/src/EnergyScaleCorrection_class.cc
 	@echo "--> Making $@" 
 	@$(COMPILE.cc) $(CXXFLAGS) $(INCLUDE) $(MAKEDEPEND) -o $@ $< 
 
@@ -74,7 +71,7 @@ scales: ../../Calibration/ZFitter/lib/EnergyScaleCorrection_class.o
 $(BUILDDIR)/analysis: $(BUILDDIR)/analysis.cpp  $(MODULES) $(MODULESZFitter)
 	@echo "---> Making analysis $(COMPILE)"
 	@g++ $(CXXFLAGS) $(INCLUDE) $(MAKEDEPEND) -o $@ $< $(MODULES) $(MODULESZFitter) $(LIB) $(ROOT_LIB) $(ROOFIT_LIB) $(ROOSTAT_LIB) $(ROOT_FLAGS) \
-	-lboost_program_options -lTreePlayer $(CMSSW_BASE)/tmp/slc6_amd64_gcc491/src/FWCore/ParameterSet/src/FWCoreParameterSet/libFWCoreParameterSet.so
+	-lboost_program_options -lTreePlayer
 
 clean:
 	rm -f $(OBJ_DIR)/*.o

Systematics_To_Be_Evaluated.txt

@@ -1,4 +1,8 @@
 Smear_Muon_Scale
 Smear_Muon_ID_Iso
+Smear_Muon_Resolution
 Smear_Electron_Scale
+Smear_Electron_Resolution
 Smear_Jet_Scale
+Smear_Jet_Resolution
+

bin/analysis.cpp

@@ -27,13 +27,12 @@
 #include "FitRooDataSet.h"
 #include "rooFitFxns.h"
 #include "ToyThrower.h"
+#include "JetResolution.h"
 #include "analysisTools.h"
 #include "configReader.h"
 
 #include <unordered_set>
 
-#include "Calibration/ZFitter/interface/EnergyScaleCorrection_class.hh"
-
 #include <TStopwatch.h>
 #define _ENDSTRING std::string::npos
 //#define DEBUG
@@ -46,6 +45,7 @@ DY TANDP POWHEG AMC MAD POWINCL
 W
 WZ
 ZZ
+WW
 data TANDP
 */
 
@@ -59,7 +59,7 @@ class chainNames
 public:
 	chainNames(): ///< default constructor
 		all_modes(  // list of all possible modes
-	{"TT", "W", "WZ", "ZZ", "data", "DYPOWHEG", "DYMADHT", "DYAMC", "DYMAD", "DYPOWINCL", "signal"
+	{"TT", "W", "WZ", "ZZ", "WW", "data", "DYPOWHEG", "DYMADHT", "DYAMC", "DYMAD", "DYPOWINCL", "signal"
 	}
 	)
 	{
@@ -81,7 +81,7 @@ class chainNames
 		}
 		if(mode == "TT") {
 			//TTchainNames.push_back("TTJets_DiLept_v1");
-			TTchainNames.push_back("TTJets_DiLept_v2");
+			TTchainNames.push_back("TTJets");
 		} else if(mode.find("DY") != _ENDSTRING) {
 			//if(mode.Contains("TANDP") ) tree_channel = "_dytagandprobe";
 			std::string tagName = "";
@@ -91,45 +91,41 @@ class chainNames
 				std::cout << "ERROR looking for DY in EMu channel" << std::endl;
 				return TTchainNames;
 			}
-			if(mode.find("POWHEG") != _ENDSTRING) {
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_50to120");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_120to200");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_200to400");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_400to800");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_800to1400");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_1400to2300");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_2300to3500");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_3500to4500");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_4500to6000");
-				TTchainNames.push_back("DYTo" + tagName + "_powheg_6000toInf");
-			} else if(mode.find("AMC") != _ENDSTRING) {
+			if(mode.find("AMC") != _ENDSTRING) {
 				//amc at nlo inclusive sample gen dilepton mass greater than 50 GeV
 				TTchainNames.push_back("DYJets_amctnlo");
-			} else if(mode.find("MAD") != _ENDSTRING) {
-				//madgraph inclusive sample gen dilepton mass greater than 50 GeV
-				TTchainNames.push_back("DYJets_madgraph");
+			// } else if(mode.find("MAD") != _ENDSTRING) {
+			// 	//madgraph inclusive sample gen dilepton mass greater than 50 GeV
+			// 	TTchainNames.push_back("DYJets_madgraph");
 			} else if(mode.find("POWINCL") != _ENDSTRING && channel == Selector::EE) {
 				TTchainNames.push_back("DYToEE_powheg");
 			} else if(mode.find("MADHT") != _ENDSTRING) {
-				TTchainNames.push_back("DYJets_madgraph_ht100to200");
-				TTchainNames.push_back("DYJets_madgraph_ht200to400");
-				TTchainNames.push_back("DYJets_madgraph_ht400to600");
-				TTchainNames.push_back("DYJets_madgraph_ht600toInf");
+				TTchainNames.push_back("DYJets_HT_100to200");
+				TTchainNames.push_back("DYJets_HT_200to400");
+				TTchainNames.push_back("DYJets_HT_400to600");
+				TTchainNames.push_back("DYJets_HT_600to800");
+				TTchainNames.push_back("DYJets_HT_800to1200");
+				TTchainNames.push_back("DYJets_HT_1200to2500");
+				TTchainNames.push_back("DYJets_HT_2500toInf");
 			}
 		} else if(mode == "W") {
 			TTchainNames.push_back("WJetsLNu");
 		} else if(mode == "WZ") {
 			TTchainNames.push_back("WZ");
 		} else if(mode == "ZZ") {
 			TTchainNames.push_back("ZZ");
+		} else if(mode == "WW") {
+			TTchainNames.push_back("WW");
 		} else if(mode == "data") {
 			std::string dataTag = "";
 			if(channel == Selector::EMu)  dataTag = "MuEG";
 			if(channel == Selector::EE)   dataTag = "DoubleEG";
 			if(channel == Selector::MuMu) dataTag = "SingleMu";
+			//TTchainNames.push_back(dataTag + "_RunB");
 			TTchainNames.push_back(dataTag + "_RunC");
-			TTchainNames.push_back(dataTag + "_RunD_v3");
-			TTchainNames.push_back(dataTag + "_RunD_v4");
+			//TTchainNames.push_back(dataTag + "_RunD");
+			//TTchainNames.push_back(dataTag + "_RunE");
+			//TTchainNames.push_back(dataTag + "_RunF");
 		}
 		if(mode.find("WRto") != _ENDSTRING) {
 			TTchainNames.push_back(mode);
@@ -247,7 +243,7 @@ int main(int ac, char* av[])
 	}
 
 
-	EnergyScaleCorrection_class eSmearer("Calibration/ZFitter/data/scales_smearings/74X_Prompt_2015");
+	//EnergyScaleCorrection_class eSmearer("ExoAnalysis/Calibration/ZFitter/data/scales_smearings/74X_Prompt_2015");
 
 	//------------------------------ check if modes given in the command line are allowed
 	for(auto s : modes ) {
@@ -260,12 +256,12 @@ int main(int ac, char* av[])
 
 	Selector::tag_t cut_channel;
 	if(channel == Selector::EMu) {
-		cut_channel = Selector::getTag(channel_cut_str);
+	  //cut_channel = Selector::getTag(channel_cut_str);
 		//outFileTag += channel_cut_str;
 	} else
 		cut_channel = channel;
 
-	configReader myReader("configs/2015-v1.conf");
+	configReader myReader("configs/2016-v1.conf");
 	if(debug) std::cout << myReader << std::endl;
 
 	std::cout << "[INFO] Selected modes: \n";
@@ -296,6 +292,8 @@ int main(int ac, char* av[])
 
 	std::map< std::pair<Selector::tag_t,  int>, std::pair<int, int> > mass_cut = getMassCutMap();
 	std::vector<int> mass_vec = getMassVec();
+	TString dataPUfn = "MyDataPileupHistogramSingleMuonC.root";
+	std::map<float, double> pu_weights = PUreweight(dataPUfn);
 
 	std::string treeName = "miniTree" + chainNames_.getTreeName(channel, isTagAndProbe, isLowDiLepton);
 	unsigned long long zMass60to120EvtCount = 0;	///<count the number of evts from each dataset with 60 < dilepton_mass < 120 which pass loose selector cuts
@@ -307,8 +305,6 @@ int main(int ac, char* av[])
 
 	for(auto mode : modes) {
 		bool isData = chainNames_.isData(mode);
-		if(isData) eSmearer.doScale = true;
-		else eSmearer.doSmearings = true;
 
 		TChain *c = myReader.getMiniTreeChain(chainNames_.getChainNames(mode, channel, isTagAndProbe), treeName);
 #ifdef DEBUG
@@ -412,6 +408,10 @@ int main(int ac, char* av[])
 			std::cout << "the number of reco electrons in the event =\t" << nEle << std::endl;
 #endif
 
+			//apply JER
+			Rand.SetSeed(seed + 1);
+			JetResolution( &myEvent, Rand, isData);
+
 			if(nEle > 0) {
 				///if there are electrons in the event, then write the electron SF and SF errors into the miniTreeEvent object named myEvent
 				///before calling the Selector constructor
@@ -455,20 +455,6 @@ int main(int ac, char* av[])
 #ifdef DEBUGG
 				std::cout << "found an event passing preselection" << std::endl;
 #endif
-				for(unsigned int iEle = 0; iEle < nEle; ++iEle) {
-#ifdef DEBUGG
-					std::cout << "looping over reco electrons in the event passing preselection" << std::endl;
-#endif
-					TLorentzVector& p4 = (*myEvent.electrons_p4)[iEle];
-					if(isData && !(channel_str == "MuMu") ) { //only scales are corrected  HEEP and Reco ID SFs set to 1.0, errors set to 0., skip this step if the channel is not EE or EMu
-						(*myEvent.electron_scale)[iEle] = eSmearer.ScaleCorrection(myEvent.run, fabs(p4.Eta()) < 1.479, (*myEvent.electron_r9)[iEle], p4.Eta(), p4.Et());
-						//(*myEvent.electron_scale)[iEle] = 1.0;
-						(*myEvent.electron_smearing)[iEle] = 0.;
-					} else if(!isData && !(channel_str == "MuMu") ) { // only the smearings are corrected, skip this step if the channel != EE or EMu (this could be entered when running over DYToLL samples)
-						(*myEvent.electron_scale)[iEle] = 1.;
-						(*myEvent.electron_smearing)[iEle] = eSmearer.getSmearingSigma(myEvent.run, fabs(p4.Eta()) < 1.479, (*myEvent.electron_r9)[iEle], p4.Eta(), p4.Et(), EnergyScaleCorrection_class::kRho, 0);
-					}
-				}
 				myEventVector.push_back(myEvent);
 			}
 			myEvent.clear();
@@ -482,23 +468,21 @@ int main(int ac, char* av[])
 		bool loop_one = true;
 		int seed_i = seed + 1;
 
-		// electron systematics!
-		bool Flag_Smear_Electron_Scale = false;
-		for(unsigned int iii = 0; iii < List_Systematics.size(); iii++) {
-			if(List_Systematics[iii] == "Smear_Electron_Scale") Flag_Smear_Electron_Scale = true;
-		}
-
-
 		for(int i = 0; i < nToys + 1; ++i, ++seed_i) {
 
-			Rand.SetSeed(seed_i);
+			Rand.SetSeed(seed_i);			
+
 			//for central values, we take the central value of Mu ID/ISO efficiencies and dont smear for JES systematics
-			// Roch and Electron scales are smeared with a pre-defined seed(1), to give conistent results.
+			// Roch and Electron scales are smeared with a pre-defined seed(1), to give consistent results.
 			if(loop_one) {
 				Random_Numbers_for_Systematics_Up_Down[0] = 0.;//Mu Eff ID
 				Random_Numbers_for_Systematics_Up_Down[1] = 0.;//Mu Eff ISO
-				Random_Numbers_for_Systematics_Up_Down[2] = 0.; //Rand.Gaus(0.0, 1.);// Electron Scale(Data)
-				Random_Numbers_for_Systematics_Up_Down[3] = 0.;//JES
+				Random_Numbers_for_Systematics_Up_Down[2] = 0.;//Mu Res
+				Random_Numbers_for_Systematics_Up_Down[3] = 0.; //Electron Scale(Data)
+				Random_Numbers_for_Systematics_Up_Down[4] = 0.; //Electron Smear(MC)
+				Random_Numbers_for_Systematics_Up_Down[5] = 0.;//JES
+				Random_Numbers_for_Systematics_Up_Down[6] = 0.;//JER
+
 			} else {
 				for(int Rand_Up_Down_Iter = 0; Rand_Up_Down_Iter < Total_Number_of_Systematics_Up_Down; Rand_Up_Down_Iter++)
 					Random_Numbers_for_Systematics_Up_Down[Rand_Up_Down_Iter] = Rand.Gaus(0.0, 1.);
@@ -545,35 +529,6 @@ int main(int ac, char* av[])
 					Random_Numbers_for_Systematics_Smear[Rand_Smear_Iter] = Rand.Gaus(0.0, 1.);
 				ToyThrower( &myEventIt, Random_Numbers_for_Systematics_Smear, Random_Numbers_for_Systematics_Up_Down, seed_i, List_Systematics, isData);
 
-
-				unsigned int nEle = myEventIt.electrons_p4->size();
-				for(unsigned int iEle = 0; iEle < nEle; ++iEle) {
-					TLorentzVector& p4 = (*myEventIt.electrons_p4)[iEle];
-
-					if(isData) { //only scales are corrected
-#ifdef DEBUG
-						std::cout << p4.Et() << "-->";
-#endif
-						p4 *= (*myEventIt.electron_scale)[iEle];
-						// since the Random_Numbers_for_Systematics is the same for all the categories, then the uncertainties are considered fully correlated
-						if(loop_one == false && Flag_Smear_Electron_Scale) p4 *= 1 + Random_Numbers_for_Systematics_Up_Down[2] * eSmearer.ScaleCorrectionUncertainty(myEventIt.run, fabs(p4.Eta()) < 1.479, 0., p4.Eta(), p4.Et());
-#ifdef DEBUG
-						std::cout << p4.Et() << "\t" << p4.Eta() << "\t" << (*myEventIt.electron_r9)[iEle] << "\t" << (*myEventIt.electron_scale)[iEle] << "\t" << (eSmearer.ScaleCorrectionUncertainty(myEventIt.run, fabs(p4.Eta()) < 1.479, 0., p4.Eta(), p4.Et())) << std::endl;
-#endif
-					} else { // only the smearings are corrected
-						double r = (loop_one == false && Flag_Smear_Electron_Scale) ? Rand.Gaus(0., 1.) : 0;
-						(*myEventIt.electron_smearing)[iEle] = eSmearer.getSmearingSigma(myEventIt.run, fabs(p4.Eta()) < 1.479, 0., p4.Eta(), p4.Et(), EnergyScaleCorrection_class::kRho, r);
-#ifdef DEBUG
-						std::cout << p4.Et() << "-->";
-#endif
-						p4 *= Rand.Gaus(1., (*myEventIt.electron_smearing)[iEle]);
-#ifdef DEBUG
-						std::cout << p4.Et() << "\t" << p4.Eta() << "\t" << (*myEventIt.electron_r9)[iEle] << "\t" <<	(*myEventIt.electron_smearing)[iEle] << "\t" <<  eSmearer.getSmearingSigma(myEventIt.run, fabs(p4.Eta()) < 1.479, 0., p4.Eta(), p4.Et(), EnergyScaleCorrection_class::kRho, 0) << "\t" << eSmearer.getSmearingSigma(myEventIt.run, fabs(p4.Eta()) < 1.479, 0., p4.Eta(), p4.Et(), EnergyScaleCorrection_class::kRho, r) << std::endl;
-#endif
-
-					}
-				}
-
 				Selector tmp_selEvent(myEventIt);
 				selEvent = tmp_selEvent;
 				// Select events with one good WR candidate
@@ -582,22 +537,10 @@ int main(int ac, char* av[])
 				// 1 -- MuMuJJ Channel
 				// 2 -- EMuJJ Channel
 
-				// if(debug && selEvent.isPassing(channel)) {
-				// 	std::cout << std::endl << selEvent.dilepton_mass << std::endl;
-				// 	std::cout << "Mu " << ev  << std::endl;
-				// 	for(auto m : * (myEventIt.muons_p4))
-				// 		std::cout << m.Pt() << " " << m.Eta() << std::endl;
-				// 	std::cout << "Ele" << std::endl;
-				// 	for(auto m : * (myEventIt.electrons_p4))
-				// 		std::cout << m.Pt() << " " << m.Eta() << std::endl;
-				// 	std::cout << "Jet" << std::endl;
-				// 	for(auto m : * (myEventIt.jets_p4))
-				// 		std::cout << m.Pt() << " " << m.Eta() << std::endl;
-				// }
-
 				if(loop_one && selEvent.isPassingLooseCuts(channel)) {
 					if(isData == false) {
-						selEvent.weight *= myReader.getNorm1fb(selEvent.datasetName) * integratedLumi; // the weight is the event weight * single object weights
+					  selEvent.weight *= myReader.getNorm1fb(selEvent.datasetName) * integratedLumi * pu_weights[int(selEvent.nPU)]; // the weight is the event weight * single object weights
+
 #ifdef DEBUGG
 						std::cout << "PU weight=\t" << selEvent.pu_weight << std::endl;
 						std::cout << "num vertices=\t" << selEvent.nPV << std::endl;
@@ -638,7 +581,7 @@ int main(int ac, char* av[])
 
 
 					if(isData == false) {
-						selEvent.weight *= myReader.getNorm1fb(selEvent.datasetName) * integratedLumi; // the weight is the event weight * single object weights
+					  selEvent.weight *= myReader.getNorm1fb(selEvent.datasetName) * integratedLumi * pu_weights[int(selEvent.nPU)]; // the weight is the event weight * single object weights
 
 						//multiply by an additional weight when processing DY samples
 						if(mode.find("DY") != _ENDSTRING && !ignoreDyScaleFactors) {

configs/2016-v1.conf

@@ -0,0 +1,6 @@
+datasetFile=configs/datasets_80X.dat
+jsonFile=/afs/cern.ch/cms/CAF/CMSCOMM/COMM_DQM/certification/Collisions16/13TeV/Cert_271036-280385_13TeV_PromptReco_Collisions16_JSON.txt
+jsonName=246908-260627-Prompt_25ns-golden_silver-v1
+productionTAG=_WRv01
+skimProductionTAG=_WRv01
+unblind=false