diff --git a/README.md b/README.md
index 2038a7f..0ca5b69 100644
--- a/README.md
+++ b/README.md
@@ -37,7 +37,7 @@ you can use the below commands for the setup.
 This also sets up the environment for the batch system at UCSB
 
 ~~~~bash
-git clone --recurse-submodules git@github.com:richstu/draw_pico.git
+git clone --recursive git@github.com:richstu/draw_pico.git
 cd draw_pico
 source set_env.sh
 ./compile.py or scons
diff --git a/inc/zgamma/KinZfitter.h b/inc/zgamma/KinZfitter.h
new file mode 100644
index 0000000..361113a
--- /dev/null
+++ b/inc/zgamma/KinZfitter.h
@@ -0,0 +1,121 @@
+/*************************************************************************
+*  Authors:   Tongguang CHeng(IHEP, Beijing) Hualin Mei(UF)
+*************************************************************************/
+#ifndef KinZfitter_h
+#define KinZfitter_h
+
+// C++ includes
+#include <iostream>
+#include <complex>
+#include <cstdlib>
+#include <vector>
+#include <fstream>
+#include <sstream>
+#include <cmath>
+#include <map>
+
+// ROOT includes
+#include "TString.h"
+#include "TLorentzVector.h"
+
+// ROOFIT
+#include "RooRealVar.h"
+#include "RooArgSet.h"
+#include "RooGaussian.h"
+#include "RooBreitWigner.h"
+#include "RooProdPdf.h"
+#include "RooDataSet.h"
+#include "RooGlobalFunc.h"
+#include "RooDataHist.h"
+#include "RooHistPdf.h"
+#include "RooCBShape.h"
+#include "RooMinuit.h"
+#include "RooFormulaVar.h"
+#include "RooAddPdf.h"
+#include "RooGenericPdf.h"
+#include "RooFFTConvPdf.h"
+#include "RooFitResult.h"
+
+// fit result covariance matrix
+#include <TMatrixDSym.h>
+
+using namespace std;
+
+class KinZfitter {
+ public:
+
+  KinZfitter(bool isData);
+
+  /// Kinematic fit of lepton momenta
+  void Setup(std::map<unsigned int, TLorentzVector> selectedLeptons, std::map<unsigned int, TLorentzVector> selectedFsrPhotons, std::map<unsigned int, double> errorLeptons);
+
+  ///
+  void KinRefitZ1();
+
+  int  PerZ1Likelihood(double & l1, double & l2, double & lph1, double & lph2);
+  void SetZ1Result(double l1, double l2, double lph1, double lph2);
+  double pterr(TLorentzVector fsrPhoton);
+  double masserror(std::vector<TLorentzVector> p4s, std::vector<double> pTErrs);
+  double masserrorFullCov(std::vector<TLorentzVector> p4s, TMatrixDSym covMatrix);
+
+  // result wrappers
+  double GetRefitM4l();
+  double GetM4l();
+  double GetRefitMZ1();
+
+  double GetMZ1Err();
+  double GetRefitM4lErr();
+  double GetM4lErr();
+  double GetRefitM4lErrFullCov();
+
+  std::vector<TLorentzVector> GetRefitP4s();
+  std::vector<TLorentzVector> GetP4s();
+
+  ////////////////
+
+  // Need to check in deep whether this function is doable
+  // RooFormulaVar p1DOTp2(RooRealVar pT1, RooRealVar theta1, RooRealVar phi1, RooRealVar m1, TString index1, RooRealVar pT2, RooRealVar theta2, RooRealVar phi2, RooRealVar m2, TString index2);
+
+ private:
+
+  /// True mZ/mZ1 shape, final states
+  TString PDFName_, fs_;
+
+  /// debug flag
+  bool debug_;
+
+  /// whether use correction for pT error
+  bool isCorrPTerr_;
+  /// whether use data or mc correction
+  bool isData_;
+
+  void initZs(std::map<unsigned int, TLorentzVector> selectedLeptons, std::map<unsigned int, TLorentzVector> selectedFsrPhoton, std::map<unsigned int, double> errorLeptons);
+
+  /// lepton ids for Z1
+  std::vector<int> idsZ1_;
+  /// lepton ids that fsr photon associated to
+  std::vector<int> idsFsrZ1_;
+  /// (Four) TLorentzVectors that form the Higgs Candidate
+  std::vector<TLorentzVector> p4sZ1_, p4sZ1ph_;
+  std::vector<TLorentzVector> p4sZ1REFIT_, p4sZ1phREFIT_;
+
+  /// pTerr vector
+  std::vector<double> pTerrsZ1_, pTerrsZ1ph_;
+  std::vector<double> pTerrsZ1REFIT_, pTerrsZ1phREFIT_;
+
+  // covariance matrix
+  // what directly coming from Refit
+  TMatrixDSym covMatrixZ1_;
+
+  // refit energy scale with respect to reco pT
+  double lZ1_l1_, lZ1_l2_;
+  double lZ1_ph1_, lZ1_ph2_;
+
+  // True mZ1 shape parameters
+  // double sgVal_, aVal_, nVal_, fVal_, meanVal_, sigmaVal_, f1Val_;
+  double meanCB_, sigmaCB_, alphaCB_, nCB_, meanGauss1_, sigmaGauss1_, f1_, meanGauss2_, sigmaGauss2_, f2_, meanGauss3_, sigmaGauss3_, f3_;
+  // double aVal_, meanVal_, nVal_, sgVal_, bwsigVal_;
+
+};
+
+#endif
diff --git a/inc/zgamma/zg_utilities.hpp b/inc/zgamma/zg_utilities.hpp
index 951b891..1e11587 100644
--- a/inc/zgamma/zg_utilities.hpp
+++ b/inc/zgamma/zg_utilities.hpp
@@ -32,5 +32,9 @@ namespace ZgUtilities {
   double cos_theta(const Baby &b, bool gen = false);
   double cos_Theta(const Baby &b, bool gen = false);
   double Getphi(const Baby &b, bool gen = false);
+  double KinRefit(const Baby &b);
+  std::vector<TLorentzVector> RefitP4(const Baby &b);
+  double AssignL1Error(const Baby &b);
+  double AssignL2Error(const Baby &b);
 }
 #endif
diff --git a/src/zgamma/KinZfitter.cpp b/src/zgamma/KinZfitter.cpp
new file mode 100644
index 0000000..763e482
--- /dev/null
+++ b/src/zgamma/KinZfitter.cpp
@@ -0,0 +1,633 @@
+/*************************************************************************
+ *  Authors:   Tongguang Cheng
+ *************************************************************************/
+#ifndef KinZfitter_cpp
+#define KinZfitter_cpp
+
+/// KinFitter header
+#include "zgamma/KinZfitter.h"
+
+KinZfitter::KinZfitter(bool isData) {
+
+  PDFName_ = "HZg_ggF_125GeV_ext1_M125_13TeV_powheg2_pythia8";
+  debug_ = true;
+  if(debug_) std::cout << "KinZfitter. The debug flag is ON with "<<PDFName_<< std::endl;
+  isData_ = isData;
+
+}
+
+
+double KinZfitter::masserrorFullCov(std::vector<TLorentzVector> p4s, TMatrixDSym covMatrix){
+
+  int ndim = 3*p4s.size();
+  if(debug_) cout<<""<<endl;
+  TMatrixD jacobian(1,ndim);
+  double e = 0; double mass = 0;
+  double px = 0; double py = 0; double pz = 0;
+  for (unsigned int ip = 0; ip < p4s.size(); ip++)
+    {
+      e = e + p4s[ip].E();
+      px = px + p4s[ip].Px();
+      py = py + p4s[ip].Py();
+      pz = pz + p4s[ip].Pz();
+    }
+  mass = TMath::Sqrt(e*e-px*px-py*py-pz*pz);
+  for (unsigned int i = 0, o = 0; i < p4s.size(); i++, o += 3)
+    {
+      double pxi = p4s[i].Px();
+      double pyi = p4s[i].Py();
+      double pzi = p4s[i].Pz();
+      double ei = p4s[i].E();
+      jacobian(0, o+0) = (e*(pxi/ei) - px)/mass;
+      jacobian(0, o+1) = (e*(pyi/ei) - py)/mass;
+      jacobian(0, o+2) = (e*(pzi/ei) - pz)/mass;
+    }
+  TMatrixDSym massCov = covMatrix.Similarity(jacobian);
+  double dm2 = massCov(0,0);
+  return (dm2 > 0 ? std::sqrt(dm2) : 0.0);
+
+}
+
+
+double KinZfitter::masserror( std::vector<TLorentzVector> Lep, std::vector<double> pterr){
+
+  TLorentzVector compositeParticle ;
+  for(unsigned int i=0; i<Lep.size(); i++)
+    {
+      compositeParticle+=Lep[i];
+    }
+  double mass  =  compositeParticle.M();
+  double masserr = 0;
+  for(unsigned int i=0; i<Lep.size(); i++)
+    {
+      TLorentzVector variedLep; // = Lep[i];
+      variedLep.SetPtEtaPhiM(Lep[i].Pt()+ pterr[i], Lep[i].Eta(), Lep[i].Phi(), Lep[i].M());
+      TLorentzVector compositeParticleVariation ;
+      for(unsigned int j=0; j<Lep.size(); j++)
+	{
+	  if(i!=j)compositeParticleVariation+=Lep[j];
+	  else compositeParticleVariation+=variedLep;
+	}
+      masserr += (compositeParticleVariation.M()-mass)*(compositeParticleVariation.M()-mass);
+    }
+  return sqrt(masserr);
+
+}
+
+
+double KinZfitter::pterr(TLorentzVector ph){
+
+  double C, S, N;
+  if (abs(ph.Eta()) < 1.48) {
+    C = 0.35 / 100;
+    S = 5.51 / 100;
+    N = 98. / 1000.;
+  } else {
+    C = 0;
+    S = 12.8 / 100;
+    N = 440. / 1000.;
+  }
+  double pherr = sqrt(C * C * ph.Energy() * ph.Energy() + S * S * ph.Energy() + N * N);
+  return pherr;
+
+}
+
+
+void KinZfitter::Setup(std::map<unsigned int, TLorentzVector> selectedLeptons, std::map<unsigned int, TLorentzVector> selectedFsrPhotons, std::map<unsigned int, double> errorLeptons) {
+
+  // reset everything for each event
+  idsZ1_.clear();
+
+  p4sZ1_.clear();
+  p4sZ1ph_.clear();
+  p4sZ1REFIT_.clear();
+  p4sZ1phREFIT_.clear();
+
+  pTerrsZ1_.clear();
+  pTerrsZ1ph_.clear();
+  pTerrsZ1REFIT_.clear();
+  pTerrsZ1phREFIT_.clear();
+
+  initZs(selectedLeptons, selectedFsrPhotons, errorLeptons);
+
+  if(debug_) cout<<"list ids"<<endl;
+  if(debug_) cout<<"IDs[0] "<<idsZ1_[0]<<" IDs[1] "<<idsZ1_[1]<<endl;
+
+  fs_="";
+  if(abs(idsZ1_[0])==11) fs_="2e";
+  if(abs(idsZ1_[0])==13) fs_="2mu";
+
+  if(debug_) cout<<"fs is "<<fs_<<endl;
+
+  /////////////
+  string paramZ1_dummy("/homes/psiddire/draw_pico/data/dummy.txt");
+
+  TString paramZ1 = TString( paramZ1_dummy.substr(0, paramZ1_dummy.length() - 9));
+
+  paramZ1+=PDFName_;
+  paramZ1+="_";
+  paramZ1+=+fs_;
+  paramZ1+=".txt";
+
+  if(debug_) cout<<"paramZ1 in "<<paramZ1<<endl;
+
+  std::ifstream input(paramZ1);
+  std::string line;
+  while (!input.eof() && std::getline(input,line))
+    {
+      std::istringstream iss(line);
+      string p; double val;
+      if(iss >> p >> val) {
+        if(p=="meanCB")      { meanCB_ = val;}
+        if(p=="sigmaCB")     { sigmaCB_ = val;}
+        if(p=="alphaCB")     { alphaCB_ = val;}
+        if(p=="nCB")         { nCB_ = val;}
+        if(p=="meanGauss1")  { meanGauss1_ = val;}
+        if(p=="sigmaGauss1") { sigmaGauss1_ = val;}
+        if(p=="f1")          { f1_ = val;}
+        if(p=="meanGauss2")  { meanGauss2_ = val; }
+        if(p=="sigmaGauss2") { sigmaGauss2_ = val;}
+        if(p=="f2")          { f2_ = val;}
+        if(p=="meanGauss3")  { meanGauss3_ = val;}
+        if(p=="sigmaGauss3") { sigmaGauss3_ = val;}
+        if(p=="f3")          { f3_ = val;}
+      }
+    }
+}
+
+///----------------------------------------------------------------------------------------------
+void KinZfitter::initZs(std::map<unsigned int, TLorentzVector> selectedLeptons, std::map<unsigned int, TLorentzVector> selectedFsrPhotons, std::map<unsigned int, double> errorLeptons) {
+
+  if(debug_) cout<<"init leptons"<<endl;
+
+  for(unsigned int il = 0; il < selectedLeptons.size(); il++)
+    {
+      double pTerr = 0;
+
+      TLorentzVector p4 = selectedLeptons[il];
+
+      int pdgId;
+      if (p4.M() > 0.05) {
+	pdgId = 13;
+	pTerr = errorLeptons[il];
+      }
+      else {
+	pdgId = 11;
+	pTerr = errorLeptons[il];
+	// pTerr = errorLeptons[il]*(p4.Pt()/p4.Et());
+      }
+      if(debug_) cout<<"pdg id "<<pdgId<<endl;
+      if(debug_) cout<<" pt err is "<<pTerr<<endl;
+
+      idsZ1_.push_back(pdgId);
+      pTerrsZ1_.push_back(pTerr);
+      p4sZ1_.push_back(p4);
+
+    }
+
+  if(debug_) cout<<"init fsr photons"<<endl;
+
+  for(unsigned int ifsr = 0; ifsr < selectedFsrPhotons.size(); ifsr++)
+    {
+
+      TLorentzVector p4 = selectedFsrPhotons[ifsr];
+      if(selectedFsrPhotons[ifsr].Pt()==0) continue;
+
+      if(debug_) cout<<"ifsr "<<ifsr<<endl;
+
+      double pTerr = 0;
+
+      pTerr = pterr(p4); //,isData_);
+
+      if(debug_) cout<<" pt err is "<<pTerr<<endl;
+
+      if(debug_) cout<<"for fsr Z1 photon"<<endl;
+
+      pTerrsZ1ph_.push_back(pTerr);
+      p4sZ1ph_.push_back(p4);
+
+    }
+
+  if(debug_) cout<<"p4sZ1ph_ "<<p4sZ1ph_.size()<<endl;
+
+}
+
+void KinZfitter::SetZ1Result(double l1, double l2, double lph1, double lph2) {
+
+  if(debug_) cout<<"start set Z1 result"<<endl;
+
+  // pT scale after refitting w.r.t. reco pT
+  lZ1_l1_ = l1; lZ1_l2_ = l2;
+  if(debug_) cout<<"l1 "<<l1<<" l2 "<<l2<<endl;
+  lZ1_ph1_ = lph1; lZ1_ph2_ = lph2;
+
+  TLorentzVector Z1_1 = p4sZ1_[0]; TLorentzVector Z1_2 = p4sZ1_[1];
+
+  TLorentzVector Z1_1_True(0,0,0,0);
+  Z1_1_True.SetPtEtaPhiM(lZ1_l1_*Z1_1.Pt(),Z1_1.Eta(),Z1_1.Phi(),Z1_1.M());
+  TLorentzVector Z1_2_True(0,0,0,0);
+  Z1_2_True.SetPtEtaPhiM(lZ1_l2_*Z1_2.Pt(),Z1_2.Eta(),Z1_2.Phi(),Z1_2.M());
+
+  p4sZ1REFIT_.push_back(Z1_1_True); p4sZ1REFIT_.push_back(Z1_2_True);
+
+  for(unsigned int ifsr1 = 0; ifsr1<p4sZ1ph_.size(); ifsr1++) {
+
+    TLorentzVector Z1ph = p4sZ1ph_[ifsr1];
+    TLorentzVector Z1phTrue(0,0,0,0);
+
+    double l = 1.0;
+    if(ifsr1==0) l = lZ1_ph1_; 
+    if(ifsr1==1) l = lZ1_ph2_;
+
+    Z1phTrue.SetPtEtaPhiM(l*Z1ph.Pt(),Z1ph.Eta(),Z1ph.Phi(),Z1ph.M());
+
+    p4sZ1phREFIT_.push_back(Z1phTrue);
+
+  }
+
+  if(debug_) cout<<"end set Z1 result"<<endl;
+
+}
+
+
+double KinZfitter::GetRefitMZ1()
+{
+
+  vector<TLorentzVector> p4s = GetRefitP4s();
+
+  TLorentzVector pZ1(0,0,0,0);
+
+  pZ1 = p4s[0] + p4s[1];
+
+  return pZ1.M();
+
+}
+
+
+double KinZfitter::GetMZ1Err()
+{
+
+  vector<TLorentzVector> p4s;
+  vector<double> pTErrs;
+
+  p4s.push_back(p4sZ1_[0]);
+  p4s.push_back(p4sZ1_[1]);
+
+  pTErrs.push_back(pTerrsZ1_[0]);
+  pTErrs.push_back(pTerrsZ1_[1]);
+
+  for(unsigned int ifsr1 = 0; ifsr1<p4sZ1ph_.size(); ifsr1++) {
+
+    p4s.push_back(p4sZ1ph_[ifsr1]);
+    pTErrs.push_back(pTerrsZ1ph_[ifsr1]);
+
+  }
+
+  return masserror(p4s,pTErrs);
+
+}
+
+
+vector<TLorentzVector> KinZfitter::GetRefitP4s()
+{
+
+  TLorentzVector Z1_1 = p4sZ1REFIT_[0];
+  TLorentzVector Z1_2 = p4sZ1REFIT_[1];
+
+  /// fsr photons
+  for(unsigned int ifsr1 = 0; ifsr1<p4sZ1phREFIT_.size(); ifsr1++) {
+
+    TLorentzVector Z1ph = p4sZ1phREFIT_[ifsr1];
+    if(ifsr1==0) Z1_1 = Z1_1 + Z1ph;
+    if(ifsr1==1) Z1_2 = Z1_2 + Z1ph;
+
+  }
+
+  vector<TLorentzVector> p4s;
+  p4s.push_back(Z1_1); p4s.push_back(Z1_2);
+
+  return p4s;
+
+}
+
+vector<TLorentzVector> KinZfitter::GetP4s()
+{
+
+  TLorentzVector Z1_1 = p4sZ1_[0];
+  TLorentzVector Z1_2 = p4sZ1_[1];
+
+  /// fsr photons
+  for(unsigned int ifsr1 = 0; ifsr1<p4sZ1ph_.size(); ifsr1++) {
+
+    TLorentzVector Z1ph = p4sZ1ph_[ifsr1];
+    if(ifsr1==0) Z1_1 = Z1_1 + Z1ph;
+    if(ifsr1==1) Z1_2 = Z1_2 + Z1ph;
+
+  }
+
+  vector<TLorentzVector> p4s;
+  p4s.push_back(Z1_1);
+  p4s.push_back(Z1_2);
+
+  return p4s;
+
+}
+
+void KinZfitter::KinRefitZ1()
+{
+
+  double l1, l2, lph1, lph2;
+  l1 = 1.0; l2 = 1.0; lph1 = 1.0; lph2 = 1.0;
+
+  PerZ1Likelihood(l1, l2, lph1, lph2);
+  if(debug_) cout<<"l1 "<<l1<<"; l2 "<<l2<<" lph1 "<<lph1<<" lph2 "<<lph2<<endl;
+  SetZ1Result(l1, l2, lph1, lph2);
+  if(debug_) cout<<"Z1 refit done"<<endl;
+
+}
+
+
+int KinZfitter::PerZ1Likelihood(double & l1, double & l2, double & lph1, double & lph2)
+{
+
+  l1= 1.0; l2 = 1.0;
+  lph1 = 1.0; lph2 = 1.0;
+
+  if(debug_) cout<<"start Z1 refit"<<endl;
+
+  TLorentzVector Z1_1 = p4sZ1_[0];
+  TLorentzVector Z1_2 = p4sZ1_[1];
+
+  double RECOpT1 = Z1_1.Pt();
+  double RECOpT2 = Z1_2.Pt();
+
+  double pTerrZ1_1 = pTerrsZ1_[0];
+  double pTerrZ1_2 = pTerrsZ1_[1];
+
+  if(debug_)cout<<"pT1 "<<RECOpT1<<" pTerrZ1_1 "<<pTerrZ1_1<<endl;
+  if(debug_)cout<<"pT2 "<<RECOpT2<<" pTerrZ1_2 "<<pTerrZ1_2<<endl;
+
+  //////////////
+
+  TLorentzVector Z1_ph1, Z1_ph2;
+  double pTerrZ1_ph1, pTerrZ1_ph2;
+  double RECOpTph1, RECOpTph2;
+
+  TLorentzVector nullFourVector(0, 0, 0, 0);
+  Z1_ph1=nullFourVector; Z1_ph2=nullFourVector;
+  RECOpTph1 = 0; RECOpTph2 = 0;
+  pTerrZ1_ph1 = 0; pTerrZ1_ph2 = 0;
+
+  if(p4sZ1ph_.size()>=1) {
+    Z1_ph1 = p4sZ1ph_[0]; pTerrZ1_ph1 = pTerrsZ1ph_[0];
+    RECOpTph1 = Z1_ph1.Pt();
+    if(debug_) cout<<"put in Z1 fsr photon 1 pT "<<RECOpTph1<<" pT err "<<pTerrZ1_ph1<<endl;
+  }
+  if(p4sZ1ph_.size()==2) {
+    //if(debug_) cout<<"put in Z1 fsr photon 2"<<endl;
+    Z1_ph2 = p4sZ1ph_[1]; pTerrZ1_ph2 = pTerrsZ1ph_[1];
+    RECOpTph2 = Z1_ph2.Pt();
+  }
+
+  RooRealVar* pT1RECO = new RooRealVar("pT1RECO", "pT1RECO", RECOpT1, 5, 500);
+  RooRealVar* pT2RECO = new RooRealVar("pT2RECO", "pT2RECO", RECOpT2, 5, 500);
+
+  double RECOpT1min = max(5.0, RECOpT1-2*pTerrZ1_1);
+  double RECOpT2min = max(5.0, RECOpT2-2*pTerrZ1_2);
+
+  RooRealVar* pTph1RECO = new RooRealVar("pTph1RECO", "pTph1RECO", RECOpTph1, 5, 500);
+  RooRealVar* pTph2RECO = new RooRealVar("pTph2RECO", "pTph2RECO", RECOpTph2, 5, 500);
+
+  double RECOpTph1min = max(0.5, RECOpTph1-2*pTerrZ1_ph1);
+  double RECOpTph2min = max(0.5, RECOpTph2-2*pTerrZ1_ph2);
+
+  // observables pT1,2,ph1,ph2
+  RooRealVar* pT1 = new RooRealVar("pT1", "pT1FIT", RECOpT1, RECOpT1min, RECOpT1+2*pTerrZ1_1 );
+  RooRealVar* pT2 = new RooRealVar("pT2", "pT2FIT", RECOpT2, RECOpT2min, RECOpT2+2*pTerrZ1_2 );
+
+  RooRealVar* m1 = new RooRealVar("m1", "m1", Z1_1.M());
+  RooRealVar* m2 = new RooRealVar("m2", "m2", Z1_2.M());
+
+  if(debug_) cout<<"m1 "<<m1->getVal()<<" m2 "<<m2->getVal()<<endl;
+
+  double Vtheta1, Vphi1, Vtheta2, Vphi2;
+  Vtheta1 = (Z1_1).Theta(); Vtheta2 = (Z1_2).Theta();
+  Vphi1 = (Z1_1).Phi(); Vphi2 = (Z1_2).Phi();
+
+  RooRealVar* theta1 = new RooRealVar("theta1", "theta1", Vtheta1);
+  RooRealVar* phi1   = new RooRealVar("phi1", "phi1", Vphi1);
+  RooRealVar* theta2 = new RooRealVar("theta2", "theta2", Vtheta2);
+  RooRealVar* phi2   = new RooRealVar("phi2", "phi2", Vphi2);
+
+  // dot product to calculate (p1+p2+ph1+ph2).M()
+  RooFormulaVar E1("E1", "TMath::Sqrt((@0*@0)/((TMath::Sin(@1))*(TMath::Sin(@1)))+@2*@2)",
+		   RooArgList(*pT1, *theta1, *m1));
+  RooFormulaVar E2("E2", "TMath::Sqrt((@0*@0)/((TMath::Sin(@1))*(TMath::Sin(@1)))+@2*@2)",
+		   RooArgList(*pT2, *theta2, *m2));
+  if(debug_) cout<<"E1 "<<E1.getVal()<<"; E2 "<<E2.getVal()<<endl;
+
+  /////
+
+  RooRealVar* pTph1 = new RooRealVar("pTph1", "pTph1FIT", RECOpTph1, RECOpTph1min, RECOpTph1+2*pTerrZ1_ph1 );
+  RooRealVar* pTph2 = new RooRealVar("pTph2", "pTph2FIT", RECOpTph2, RECOpTph2min, RECOpTph2+2*pTerrZ1_ph2 );
+
+  double Vthetaph1, Vphiph1, Vthetaph2, Vphiph2;
+  Vthetaph1 = (Z1_ph1).Theta(); Vthetaph2 = (Z1_ph2).Theta();
+  Vphiph1 = (Z1_ph1).Phi(); Vphiph2 = (Z1_ph2).Phi();
+
+  RooRealVar* thetaph1 = new RooRealVar("thetaph1", "thetaph1", Vthetaph1);
+  RooRealVar* phiph1   = new RooRealVar("phiph1", "phiph1", Vphiph1);
+  RooRealVar* thetaph2 = new RooRealVar("thetaph2", "thetaph2", Vthetaph2);
+  RooRealVar* phiph2   = new RooRealVar("phiph2", "phi2", Vphiph2);
+
+  RooFormulaVar Eph1("Eph1", "TMath::Sqrt((@0*@0)/((TMath::Sin(@1))*(TMath::Sin(@1))))",
+		     RooArgList(*pTph1, *thetaph1));
+  RooFormulaVar Eph2("Eph2", "TMath::Sqrt((@0*@0)/((TMath::Sin(@1))*(TMath::Sin(@1))))",
+		     RooArgList(*pTph2, *thetaph2));
+
+  //// dot products of 4-vectors
+
+  // 3-vector DOT
+  RooFormulaVar* p1v3D2 = new RooFormulaVar("p1v3D2",
+					    "@0*@1*( ((TMath::Cos(@2))*(TMath::Cos(@3)))/((TMath::Sin(@2))*(TMath::Sin(@3)))+(TMath::Cos(@4-@5)))",
+					    RooArgList(*pT1, *pT2, *theta1, *theta2, *phi1, *phi2));
+  if(debug_) cout<<"p1 DOT p2 is "<<p1v3D2->getVal()<<endl;
+  // 4-vector DOT metric 1 -1 -1 -1
+  RooFormulaVar p1D2("p1D2", "@0*@1-@2", RooArgList(E1, E2, *p1v3D2));
+
+  //lep DOT fsrPhoton1
+
+  // 3-vector DOT
+  RooFormulaVar* p1v3Dph1 = new RooFormulaVar("p1v3Dph1",
+					      "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@4-@5))",
+					      RooArgList(*pT1, *pTph1, *theta1, *thetaph1, *phi1, *phiph1));
+
+  // 4-vector DOT metric 1 -1 -1 -1
+  RooFormulaVar p1Dph1("p1Dph1", "@0*@1-@2", RooArgList(E1, Eph1, *p1v3Dph1));
+
+  // 3-vector DOT
+  RooFormulaVar* p2v3Dph1 = new RooFormulaVar("p2v3Dph1",
+					      "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@4-@5))",
+					      RooArgList(*pT2, *pTph1, *theta2, *thetaph1, *phi2, *phiph1));
+  // 4-vector DOT metric 1 -1 -1 -1
+  RooFormulaVar p2Dph1("p2Dph1", "@0*@1-@2", RooArgList(E2, Eph1, *p2v3Dph1));
+
+  // lep DOT fsrPhoton2
+
+  // 3-vector DOT
+  RooFormulaVar* p1v3Dph2 = new RooFormulaVar("p1v3Dph2",
+					      "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@4-@5))",
+					      RooArgList(*pT1, *pTph2, *theta1, *thetaph2, *phi1, *phiph2));
+
+  // 4-vector DOT metric 1 -1 -1 -1
+  RooFormulaVar p1Dph2("p1Dph2", "@0*@1-@2", RooArgList(E1, Eph2, *p1v3Dph2));
+
+  // 3-vector DOT
+  RooFormulaVar* p2v3Dph2 = new RooFormulaVar("p2v3Dph2",
+					      "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@4-@5))",
+					      RooArgList(*pT2, *pTph2, *theta2, *thetaph2, *phi2, *phiph2));
+  // 4-vector DOT metric 1 -1 -1 -1
+  RooFormulaVar p2Dph2("p2Dph2", "@0*@1-@2", RooArgList(E2, Eph2, *p2v3Dph2));
+
+  // fsrPhoton1 DOT fsrPhoton2
+
+  // 3-vector DOT
+  RooFormulaVar* ph1v3Dph2 = new RooFormulaVar("ph1v3Dph2",
+						 "@0*@1*( (TMath::Cos(@2)*TMath::Cos(@3))/(TMath::Sin(@2)*TMath::Sin(@3))+TMath::Cos(@4-@5))",
+					       RooArgList(*pTph1, *pTph2, *thetaph1, *thetaph2, *phiph1, *phiph2));
+  // 4-vector DOT metric 1 -1 -1 -1
+  RooFormulaVar ph1Dph2("ph1Dph2", "@0*@1-@2", RooArgList(Eph1, Eph2, *ph1v3Dph2));
+
+  // mZ1
+  RooFormulaVar* mZ1;
+  mZ1 = new RooFormulaVar("mZ1", "TMath::Sqrt(2*@0+@1*@1+@2*@2)", RooArgList(p1D2, *m1, *m2));
+  if(p4sZ1ph_.size()==1)
+    mZ1 = new RooFormulaVar("mZ1", "TMath::Sqrt(2*@0+2*@1+2*@2+@3*@3+@4*@4)",
+			    RooArgList(p1D2, p1Dph1, p2Dph1, *m1, *m2));
+  if(p4sZ1ph_.size()==2)
+    mZ1 = new RooFormulaVar("mZ1", "TMath::Sqrt(2*@0+2*@1+2*@2+2*@3+2*@4+2*@5+@6*@6+@7*@7)",
+			    RooArgList(p1D2, p1Dph1, p2Dph1, p1Dph2, p2Dph2, ph1Dph2, *m1, *m2));
+
+  if(debug_) cout<<"mZ1 is "<<mZ1->getVal()<<endl;
+
+  // pTerrs, 1, 2, ph1, ph2
+  RooRealVar sigmaZ1_1("sigmaZ1_1", "sigmaZ1_1", pTerrZ1_1);
+  RooRealVar sigmaZ1_2("sigmaZ1_2", "sigmaZ1_2", pTerrZ1_2);
+
+  RooRealVar sigmaZ1_ph1("sigmaZ1_ph1", "sigmaZ1_ph1", pTerrZ1_ph1);
+  RooRealVar sigmaZ1_ph2("sigmaZ1_ph2", "sigmaZ1_ph2", pTerrZ1_ph2);
+
+  // resolution for decay products
+  RooGaussian gauss1("gauss1", "gaussian PDF", *pT1RECO, *pT1, sigmaZ1_1);
+  RooGaussian gauss2("gauss2", "gaussian PDF", *pT2RECO, *pT2, sigmaZ1_2);
+
+  RooGaussian gaussph1("gaussph1", "gaussian PDF", *pTph1RECO, *pTph1, sigmaZ1_ph1);
+  RooGaussian gaussph2("gaussph2", "gaussian PDF", *pTph2RECO, *pTph2, sigmaZ1_ph2);
+
+  RooRealVar meanCB("meanCB", "", meanCB_);
+  RooRealVar sigmaCB("sigmaCB", "", sigmaCB_);
+  RooRealVar alphaCB("alphaCB", "", alphaCB_);
+  RooRealVar nCB("nCB", "", nCB_);
+  RooRealVar meanGauss1("meanGauss1", "", meanGauss1_);
+  RooRealVar sigmaGauss1("sigmaGauss1", "", sigmaGauss1_);
+  RooRealVar f1("f1", "", f1_);
+  RooRealVar meanGauss2("meanGauss2", "", meanGauss2_);
+  RooRealVar sigmaGauss2("sigmaGauss2", "", sigmaGauss2_);
+  RooRealVar f2("f2", "", f2_);
+  RooRealVar meanGauss3("meanGauss3", "", meanGauss3_);
+  RooRealVar sigmaGauss3("sigmaGauss3", "", sigmaGauss3_);
+  RooRealVar f3("f3", "", f3_);
+
+  RooCBShape* singleCB = new RooCBShape("singleCB", "", *mZ1, meanCB, sigmaCB, alphaCB, nCB);
+  RooGaussian* gaussShape1 = new RooGaussian("gaussShape1", "", *mZ1, meanGauss1, sigmaGauss1);
+  RooAddPdf* CBplusGauss = new RooAddPdf("CBplusGauss", "", *singleCB, *gaussShape1, f1);
+  RooGaussian* gaussShape2 = new RooGaussian("gaussShape2", "", *mZ1, meanGauss2, sigmaGauss2);
+  RooAddPdf* CBplusGaussplusGauss = new RooAddPdf("CBplusGaussplusGauss", "", *CBplusGauss, *gaussShape2, f2);
+  RooGaussian* gaussShape3 = new RooGaussian("gaussShape3", "", *mZ1, meanGauss3, sigmaGauss3);
+  RooAddPdf* CBplusGaussplusGaussplusGauss = new RooAddPdf("CBplusGaussplusGaussplusGauss", "", *CBplusGaussplusGauss, *gaussShape3, f3);
+
+  RooProdPdf *model;
+  model = new RooProdPdf("model", "model", RooArgList(gauss1, gauss2, *CBplusGaussplusGaussplusGauss) );
+  if(p4sZ1ph_.size()==1)
+    model = new RooProdPdf("model", "model", RooArgList(gauss1, gauss2, gaussph1, *CBplusGaussplusGaussplusGauss) );
+  if(p4sZ1ph_.size()==2)
+    model = new RooProdPdf("model", "model", RooArgList(gauss1, gauss2, gaussph1, gaussph2, *CBplusGaussplusGaussplusGauss) );
+
+  // observable set
+  RooArgSet *rastmp;
+  rastmp = new RooArgSet(*pT1RECO, *pT2RECO);
+  if(p4sZ1ph_.size()==1)
+    rastmp = new RooArgSet(*pT1RECO, *pT2RECO, *pTph1RECO);
+  if(p4sZ1ph_.size()>=2)
+    rastmp = new RooArgSet(*pT1RECO, *pT2RECO, *pTph1RECO, *pTph2RECO);
+
+  RooDataSet* pTs = new RooDataSet("pTs", "pTs", *rastmp);
+  pTs->add(*rastmp);
+
+  //RooAbsReal* nll;
+  //nll = model->createNLL(*pTs);
+  //RooMinuit(*nll).migrad();
+
+  RooFitResult* r = model->fitTo(*pTs, RooFit::Save(), RooFit::PrintLevel(-1));
+  const TMatrixDSym& covMatrix = r->covarianceMatrix();
+
+  const RooArgList& finalPars = r->floatParsFinal();
+  for (int i=0 ; i<finalPars.getSize(); i++) {
+    // TString name = TString(((RooRealVar*)finalPars.at(i))->GetName());
+    TString name = "Cast Error";
+
+    if(debug_) cout<<"name list of RooRealVar for covariance matrix "<<name<<endl;
+
+  }
+
+  int size = covMatrix.GetNcols();
+  //TMatrixDSym covMatrixTest_(size);
+  covMatrixZ1_.ResizeTo(size, size);
+  covMatrixZ1_ = covMatrix;
+
+  if(debug_) cout<<"save the covariance matrix"<<endl;
+
+  l1 = pT1->getVal()/RECOpT1;
+  l2 = pT2->getVal()/RECOpT2;
+  double pTerrZ1REFIT1 = pT1->getError();
+  double pTerrZ1REFIT2 = pT2->getError();
+
+  pTerrsZ1REFIT_.push_back(pTerrZ1REFIT1);
+  pTerrsZ1REFIT_.push_back(pTerrZ1REFIT2);
+
+  if(p4sZ1ph_.size()>=1) {
+
+    if(debug_) cout<<"set refit result for Z1 fsr photon 1"<<endl;
+
+    lph1 = pTph1->getVal()/RECOpTph1;
+    double pTerrZ1phREFIT1 = pTph1->getError();
+    if(debug_) cout<<"scale "<<lph1<<" pterr "<<pTerrZ1phREFIT1<<endl;
+
+    pTerrsZ1phREFIT_.push_back(pTerrZ1phREFIT1);
+
+  }
+  if(p4sZ1ph_.size()==2) {
+
+    lph2 = pTph2->getVal()/RECOpTph2;
+    double pTerrZ1phREFIT2 = pTph2->getError();
+    pTerrsZ1phREFIT_.push_back(pTerrZ1phREFIT2);
+
+  }
+
+  //delete nll;
+  delete r;
+  delete mZ1;
+  delete pT1; delete pT2; delete pTph1; delete pTph2;
+  delete pT1RECO; delete pT2RECO; delete pTph1RECO; delete pTph2RECO;
+  delete ph1v3Dph2; delete p1v3Dph1; delete p2v3Dph1; delete p1v3Dph2; delete p2v3Dph2;
+  delete model;
+  delete pTs;
+  delete rastmp;
+
+  if(debug_) cout<<"end Z1 refit"<<endl;
+
+  return 0;
+
+}
+
+#endif
diff --git a/src/zgamma/plot_KinFit.cxx b/src/zgamma/plot_KinFit.cxx
new file mode 100644
index 0000000..865e841
--- /dev/null
+++ b/src/zgamma/plot_KinFit.cxx
@@ -0,0 +1,113 @@
+#include <iostream>
+#include <string>
+#include <vector>
+#include <memory>
+#include <bitset>
+#include <unistd.h>
+#include <getopt.h>
+#include "TError.h"
+#include "TColor.h"
+#include "TVector3.h"
+#include "TLorentzVector.h"
+#include "core/baby.hpp"
+#include "core/process.hpp"
+#include "core/named_func.hpp"
+#include "core/plot_maker.hpp"
+#include "core/plot_opt.hpp"
+#include "core/palette.hpp"
+#include "core/table.hpp"
+#include "core/hist1d.hpp"
+#include "core/utilities.hpp"
+#include "zgamma/zg_utilities.hpp"
+using namespace std;
+using namespace PlotOptTypes;
+using namespace ZgUtilities;
+
+int main() {
+  gErrorIgnoreLevel = 6000;
+  Palette colors("txt/colors.txt","default");
+  Process::Type sig  = Process::Type::signal;
+  string bfolder("/net/cms17/cms17r0/pico/");
+  string ul_path(bfolder+"NanoAODv9/zgamma_signal/2016/signal/merged_zgmc_llg/");
+
+  NamedFunc kinFit("kinFit", [](const Baby &b) -> NamedFunc::ScalarType{
+      double mZ_fit = KinRefit(b);
+      return mZ_fit;
+    });
+  NamedFunc kinFitH("kinFitH", [](const Baby &b) -> NamedFunc::ScalarType{
+      std::vector<TLorentzVector> reFit = RefitP4(b);
+      TLorentzVector H = reFit[0] + reFit[1] + AssignGamma(b);
+      return H.M();
+    });
+  NamedFunc wgt("weight",[](const Baby &b) -> NamedFunc::ScalarType{
+      double weight = b.w_lumi();
+      return weight*100;
+    });
+
+  NamedFunc el_trigs("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_DZ");
+  NamedFunc mu_trigs("HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ || HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ");
+  NamedFunc trigs(el_trigs || mu_trigs);
+  NamedFunc mass_cuts("ll_m[llphoton_ill[0]] > 50 && "
+  		      "llphoton_m[0]+ll_m[llphoton_ill[0]]>=185 && "
+  		      "llphoton_m[0] > 100 && llphoton_m[0] < 180 && "
+  		      "photon_pt[llphoton_iph[0]]/llphoton_m[0] >= 15./110");
+  NamedFunc baseline("nllphoton > 0");
+  vector<NamedFunc> lep = {"ll_lepid[llphoton_ill[0]] == 11 && "
+  			   "el_pt[ll_i1[llphoton_ill[0]]] > 25 && "
+  			   "el_pt[ll_i2[llphoton_ill[0]]] > 15 && "
+  			   "el_sig[ll_i1[llphoton_ill[0]]] && "
+  			   "el_sig[ll_i2[llphoton_ill[0]]]",
+                           "ll_lepid[llphoton_ill[0]] == 13 && "
+  			   "mu_pt[ll_i1[llphoton_ill[0]]] > 20 && "
+  			   "mu_pt[ll_i2[llphoton_ill[0]]] > 10 && "
+  			   "mu_sig[ll_i1[llphoton_ill[0]]] && "
+  			   "mu_sig[ll_i2[llphoton_ill[0]]]"};
+  NamedFunc pho("photon_pt[llphoton_iph[0]] > 15 && "
+  		"photon_drmin[llphoton_iph[0]] > 0.4 && "
+  		"photon_sig[llphoton_iph[0]]");
+
+  PlotOpt log_lumi("txt/plot_styles.txt","CMSPaper");
+  log_lumi.Title(TitleType::info)
+          .YAxis(YAxisType::log)
+          .Stack(StackType::shapes)
+          .Overflow(OverflowType::none)
+          .YTitleOffset(1.75)
+          .AutoYAxis(false)
+          .UseCMYK(false)
+          .LeftMargin(0.17)
+          .LegendColumns(1)
+          .CanvasWidth(800)
+          .CanvasHeight(800)
+          .FileExtensions({"pdf"});
+
+  PlotOpt lin_lumi = log_lumi().YAxis(YAxisType::linear);
+  PlotOpt lin_stack = lin_lumi().Stack(StackType::signal_overlay);
+  PlotOpt log_stack = log_lumi().Stack(StackType::signal_overlay);
+  vector<PlotOpt> ops = {lin_stack};
+
+  auto proc_Red     = Process::MakeShared<Baby_pico>("HToZ#gamma (UL)", sig, 
+						     TColor::GetColor("#ff0000"), {ul_path+"*.root"}, trigs);
+
+  auto proc_Blue    = Process::MakeShared<Baby_pico>("HToZ#gamma (UL)", sig,
+						     TColor::GetColor("#0000ff"), {ul_path+"*.root"}, trigs);
+
+  proc_Red->SetLineWidth(3);
+  proc_Blue->SetLineWidth(3);
+
+  vector<shared_ptr<Process>> procs_Red = {proc_Red};
+  vector<shared_ptr<Process>> procs_Blue = {proc_Blue};
+
+  PlotMaker pm;
+  string lepName[2] = {"el/", "mu/"};
+  for(int i(0); i < 2; i++) {
+    NamedFunc selection = baseline && lep.at(i) && pho && mass_cuts;
+    pm.Push<Hist1D>(Axis(80, 100, 180, "llphoton_m[0]", "m_{ll#gamma} [GeV]", {}), selection, procs_Red, ops).Weight(wgt).Tag(lepName[i]+"llphoton_m");
+    pm.Push<Hist1D>(Axis(100, 50, 150, "ll_m[llphoton_ill[0]]", "m_{ll} [GeV]", {}), selection, procs_Red, ops).Weight(wgt).Tag(lepName[i]+"ll_m");
+    pm.Push<Hist1D>(Axis(80, 100, 180, kinFitH, "m_{ll#gamma} [GeV]", {}), selection, procs_Blue, ops).Weight(wgt).Tag(lepName[i]+"llphoton_m_Refit");
+    pm.Push<Hist1D>(Axis(100, 50, 150, kinFit, "m_{ll} [GeV]", {}), selection, procs_Blue, ops).Weight(wgt).Tag(lepName[i]+"ll_m_Refit");
+  }
+  pm.min_print_ = true;
+  pm.multithreaded_ = false;
+  pm.MakePlots(35.9);
+}
+
diff --git a/src/zgamma/zg_utilities.cpp b/src/zgamma/zg_utilities.cpp
index 4e5d767..b6044f9 100644
--- a/src/zgamma/zg_utilities.cpp
+++ b/src/zgamma/zg_utilities.cpp
@@ -5,6 +5,7 @@
 #include "core/utilities.hpp"
 #include "core/palette.hpp"
 #include "core/baby.hpp"
+#include "zgamma/KinZfitter.h"
 
 namespace ZgUtilities {
   using std::string;
@@ -33,7 +34,7 @@ namespace ZgUtilities {
         else                                        il = b.ll_i2()->at(0);
         l1.SetPtEtaPhiM(b.el_pt() ->at(il), 
                         b.el_eta()->at(il),
-                        b.el_phi()->at(il), 0.00511);
+                        b.el_phi()->at(il), 0.000511);
       }
       else if(b.ll_lepid()->at(0) == 13){
         if(b.mu_charge()->at(b.ll_i1()->at(0)) < 0) il = b.ll_i1()->at(0);
@@ -68,7 +69,7 @@ namespace ZgUtilities {
         else                                        il = b.ll_i1()->at(0);
         l2.SetPtEtaPhiM(b.el_pt() ->at(il), 
                         b.el_eta()->at(il),
-                        b.el_phi()->at(il), 0.00511);
+                        b.el_phi()->at(il), 0.000511);
       }
       else if(b.ll_lepid()->at(0) == 13){
         if(b.mu_charge()->at(b.ll_i1()->at(0)) < 0) il = b.ll_i2()->at(0);
@@ -81,6 +82,86 @@ namespace ZgUtilities {
     return l2;
   }
 
+  // Returns negative lepton 4-momentum error
+  double AssignL1Error(const Baby &b) {
+    double l1Err = 0.0;
+    double ecalEnergy = 0.0;
+    double perr = 0.0;
+    double err2 = 0.0;
+    int il(-1);
+    TLorentzVector l1;
+    if(b.ll_lepid()->at(0) == 11){
+      if(b.el_charge()->at(b.ll_i1()->at(0)) < 0) il = b.ll_i1()->at(0);
+      else                                        il = b.ll_i2()->at(0);
+      l1.SetPtEtaPhiM(b.el_pt() ->at(il),
+                      b.el_eta()->at(il),
+                      b.el_phi()->at(il), 0.000511);
+      if (b.el_ecal()->at(il)) {
+        ecalEnergy = (b.el_etPt()->at(il) + 1.0) * b.el_pt()->at(il);
+        if (abs(b.el_eta()->at(il)) < 1.5) {
+          err2 += (5.24e-02*5.24e-02)/ecalEnergy;
+          err2 += (2.01e-01*2.01e-01)/(ecalEnergy*ecalEnergy);
+          err2 += 1.00e-02*1.00e-02;
+        } else {
+          err2 += (1.46e-01*1.46e-01)/ecalEnergy;
+          err2 += (9.21e-01*9.21e-01)/(ecalEnergy*ecalEnergy);
+          err2 += 1.94e-03*1.94e-03;
+        }
+        perr = ecalEnergy * sqrt(err2);
+      }
+      else {
+        perr = b.el_energyErr()->at(il);
+      }
+      l1Err = perr * l1.Pt() / l1.P();
+    }
+    else if(b.ll_lepid()->at(0) == 13){
+      if(b.mu_charge()->at(b.ll_i1()->at(0)) < 0) il = b.ll_i1()->at(0);
+      else                                        il = b.ll_i2()->at(0);
+      l1Err = b.mu_ptErr()->at(il);
+    }
+    return l1Err;
+  }
+
+  // Returns positive lepton 4-momentum error
+  double AssignL2Error(const Baby &b) {
+    double l2Err = 0.0;
+    double ecalEnergy = 0.0;
+    double perr = 0.0;
+    double err2 = 0.0;
+    int il(-1);
+    TLorentzVector l2;
+    if(b.ll_lepid()->at(0) == 11){
+      if(b.el_charge()->at(b.ll_i1()->at(0)) < 0) il = b.ll_i2()->at(0);
+      else                                        il = b.ll_i1()->at(0);
+      l2.SetPtEtaPhiM(b.el_pt() ->at(il),
+                      b.el_eta()->at(il),
+                      b.el_phi()->at(il), 0.000511);
+      if (b.el_ecal()->at(il)) {
+        ecalEnergy = (b.el_etPt()->at(il) + 1.0) * b.el_pt()->at(il);
+        if (abs(b.el_eta()->at(il)) < 1.5) {
+          err2 += (5.24e-02*5.24e-02)/ecalEnergy;
+          err2 += (2.01e-01*2.01e-01)/(ecalEnergy*ecalEnergy);
+          err2 += 1.00e-02*1.00e-02;
+        } else {
+          err2 += (1.46e-01*1.46e-01)/ecalEnergy;
+          err2 += (9.21e-01*9.21e-01)/(ecalEnergy*ecalEnergy);
+          err2 += 1.94e-03*1.94e-03;
+        }
+        perr = ecalEnergy * sqrt(err2);
+      }
+      else {
+        perr = b.el_energyErr()->at(il);
+      }
+      l2Err = perr * l2.Pt() / l2.P();
+    }
+    else if(b.ll_lepid()->at(0) == 13){
+      if(b.mu_charge()->at(b.ll_i1()->at(0)) < 0) il = b.ll_i2()->at(0);
+      else                                        il = b.ll_i1()->at(0);
+      l2Err = b.mu_ptErr()->at(il);
+    }
+    return l2Err;
+  }
+
   // Returns Z 4-momentum
   TLorentzVector AssignZ(const Baby &b, bool gen) {
     TLorentzVector ll;
@@ -234,6 +315,53 @@ namespace ZgUtilities {
     else           phi = acos(cosphi);
     return phi;
   }
+
+  double KinRefit(const Baby &b) {
+    KinZfitter *kinZfitter;
+    bool isData = false;
+    kinZfitter = new KinZfitter(isData);
+
+    std::map<unsigned int, TLorentzVector> selectedLeptons;
+    selectedLeptons[0] = AssignL1(b);
+    selectedLeptons[1] = AssignL2(b);
+    std::map<unsigned int, double> errorLeptons;
+    errorLeptons[0] = AssignL1Error(b);
+    errorLeptons[1] = AssignL2Error(b);
+    std::map<unsigned int, TLorentzVector> selectedFsrMap;
+    // TLorentzVector nullFourVector(0, 0, 0, 0);  
+    // selectedFsrMap[0] = nullFourVector;
+    // selectedFsrMap[1] = nullFourVector;
+
+    // Need to updated selectedFsrMap to incorporate FSR contribution
+    kinZfitter->Setup(selectedLeptons, selectedFsrMap, errorLeptons);
+    kinZfitter->KinRefitZ1();
+    double massZ1REFIT = kinZfitter->GetRefitMZ1();
+    return massZ1REFIT;
+  }
+
+  std::vector<TLorentzVector> RefitP4(const Baby &b) {
+    KinZfitter *kinZfitter;
+    bool isData = false;
+    kinZfitter = new KinZfitter(isData);
+
+    std::map<unsigned int, TLorentzVector> selectedLeptons;
+    selectedLeptons[0] = AssignL1(b);
+    selectedLeptons[1] = AssignL2(b);
+    std::map<unsigned int, double> errorLeptons;
+    errorLeptons[0] = AssignL1Error(b);
+    errorLeptons[1] = AssignL2Error(b);
+    std::map<unsigned int, TLorentzVector> selectedFsrMap;
+    // TLorentzVector nullFourVector(0, 0, 0, 0);  
+    // selectedFsrMap[0] = nullFourVector;
+    // selectedFsrMap[1] = nullFourVector;
+
+    // Need to updated selectedFsrMap to incorporate FSR contribution
+    kinZfitter->Setup(selectedLeptons, selectedFsrMap, errorLeptons);
+    kinZfitter->KinRefitZ1();
+    std::vector<TLorentzVector> reFit = kinZfitter->GetRefitP4s();
+    return reFit;
+  }
+
 }
 
 
diff --git a/txt/variables/pico b/txt/variables/pico
index 0cd9371..302f0fc 100644
--- a/txt/variables/pico
+++ b/txt/variables/pico
@@ -6,6 +6,7 @@ int type
 bool stitch
 bool stitch_ht
 bool stitch_htmet
+bool stitch_photon
 bool stitch_dy
 
 ####################   Global   ###################
@@ -56,12 +57,21 @@ std::vector<float> llphoton_pt
 std::vector<float> llphoton_eta
 std::vector<float> llphoton_phi
 std::vector<float> llphoton_m
+std::vector<float> llphoton_l1_masserr
+std::vector<float> llphoton_l2_masserr
+std::vector<float> llphoton_ph_masserr
 std::vector<float> llphoton_dr
 std::vector<float> llphoton_dphi
 std::vector<float> llphoton_deta
 std::vector<float> llphoton_costhj
+std::vector<float> llphoton_costheta
+std::vector<float> llphoton_cosTheta
+std::vector<float> llphoton_psi
 std::vector<int>   llphoton_iph
 std::vector<int>   llphoton_ill
+std::vector<float> llphoton_dijet_dphi
+std::vector<float> llphoton_dijet_balance
+std::vector<float> llphoton_pTt2
 
 ###################   Jets   ##################
 int njet
@@ -86,6 +96,7 @@ std::vector<float> jet_puid_disc
 std::vector<float> jet_ne_emef
 std::vector<float> jet_qgl
 std::vector<bool>  jet_islep
+std::vector<bool>  jet_isvlep
 std::vector<bool>  jet_isphoton
 std::vector<bool>  jet_isgood
 std::vector<bool>  jet_id
@@ -117,6 +128,7 @@ std::vector<float> fjet_msoftdrop
 std::vector<float> fjet_deep_md_hbb_btv
 std::vector<float> fjet_mva_hbb_btv
 std::vector<float> fjet_deep_md_hbb_jme
+std::vector<float> fjet_deep_md_tvsqcd
 std::vector<float> fjet_deep_tvsqcd
 std::vector<int>   fjet_subjet_idx1
 std::vector<int>   fjet_subjet_idx2
@@ -160,6 +172,7 @@ int nvmu
 std::vector<float> mu_pt
 std::vector<float> mu_eta
 std::vector<float> mu_phi
+std::vector<float> mu_ptErr
 std::vector<float> mu_miniso
 std::vector<float> mu_reliso
 std::vector<float> mu_dz
@@ -167,15 +180,24 @@ std::vector<float> mu_dxy
 std::vector<float> mu_ip3d
 std::vector<float> mu_sip3d
 std::vector<bool>  mu_id
+std::vector<bool>  mu_mediumid
+std::vector<bool>  mu_tightid
+std::vector<char>  mu_highptid
 std::vector<bool>  mu_sig
 std::vector<int>   mu_charge
 std::vector<int>   mu_pflavor
+std::vector<int>   mu_fsrphotonid
 
 int nel
 int nvel
 std::vector<float> el_pt
 std::vector<float> el_eta
+std::vector<float> el_etasc
 std::vector<float> el_phi
+std::vector<float> el_etPt
+std::vector<float> el_eminusp
+std::vector<float> el_energyErr
+std::vector<int>   el_phidx
 std::vector<float> el_miniso
 std::vector<float> el_reliso
 std::vector<float> el_dz
@@ -184,10 +206,14 @@ std::vector<float> el_ip3d
 std::vector<float> el_sip3d
 std::vector<float> el_idmva
 std::vector<bool>  el_id
+std::vector<bool>  el_id80
+std::vector<bool>  el_id90
+std::vector<bool>  el_idLoose
 std::vector<bool>  el_sig
 std::vector<bool>  el_ispf
 std::vector<int>   el_charge
 std::vector<int>   el_pflavor
+std::vector<bool>  el_ecal
 
 ##################   Di-Leptons   ################
 int nll
@@ -195,12 +221,15 @@ std::vector<float> ll_pt
 std::vector<float> ll_eta
 std::vector<float> ll_phi
 std::vector<float> ll_m
+std::vector<int>   ll_l1_masserr
+std::vector<int>   ll_l2_masserr
 std::vector<float> ll_dr
 std::vector<float> ll_dphi
 std::vector<float> ll_deta
 std::vector<int>   ll_lepid
 std::vector<int>   ll_i1
 std::vector<int>   ll_i2
+std::vector<int>   ll_charge
 
 ##################   Photons   ###################
 int nphoton
@@ -217,10 +246,23 @@ std::vector<int>   photon_pflavor
 std::vector<bool>  photon_id
 std::vector<bool>  photon_elveto
 std::vector<bool>  photon_sig
+std::vector<bool>  photon_id80
+std::vector<float> photon_sieie
+std::vector<int>   photon_elidx
+std::vector<float> photon_jet_mindr
+std::vector<float> photon_zeppenfeld
+
+#################   FsrPhotons  ##################
+int nfsrphoton
+std::vector<float> fsrphoton_droveret2
+std::vector<float> fsrphoton_eta
+std::vector<float> fsrphoton_phi
+std::vector<float> fsrphoton_pt
+std::vector<float> fsrphoton_reliso
+std::vector<int>   fsrphoton_muonidx
 
 ###################   Tracks    ##################
 int ntk
-
 std::vector<float> tk_pt
 std::vector<float> tk_eta
 std::vector<float> tk_phi
@@ -230,7 +272,8 @@ std::vector<int>   tk_pdgid
 std::vector<float> tk_miniso_chg
 std::vector<float> tk_reliso_chg
 std::vector<float> tk_mt
-
+std::vector<bool>  tk_nano_electron
+std::vector<bool>  tk_nano_muon
 
 ###################   Quality   ##################
 bool pass_jets
@@ -245,9 +288,10 @@ bool pass_mubadtrk
 bool pass_badchhad
 bool pass_muon_jet
 # pass_ra2_badmu is for higgsino_angeles compatibility
-bool pass_ra2_badmu 
+bool pass_ra2_badmu
 bool pass_low_neutral_jet
 bool pass_htratio_dphi_tight
+bool pass_ecalnoisejet
 bool pass_badcalib
 bool pass_boosted
 bool pass
@@ -394,6 +438,9 @@ bool HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL
 bool HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL
 bool HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ
 bool HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ
+bool HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass3p8
+bool HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_Mass8
+bool HLT_Photon175
 
 ################  Systematic variations  - to be implemented as needed  ##################