diff --git a/TopAnalysis/bin/analysisWrapper.cc b/TopAnalysis/bin/analysisWrapper.cc
index 6bf2fbf102c84d9b929c9f5bbca186e5fd0f0fc8..7290f734039de0f9259639bb783e4ef7c0cb0b31 100644
--- a/TopAnalysis/bin/analysisWrapper.cc
+++ b/TopAnalysis/bin/analysisWrapper.cc
@@ -22,7 +22,8 @@ void printHelp()
        << "\t --systVar  - specify single systematic variation" << endl
        << "\t --era      - era directory to use for corrections, uncertainties" << endl
        << "\t --normTag  - normalization tag" << endl
-       << "\t --method   - method to run" << endl;
+       << "\t --method   - method to run" << endl
+       << "\t --mvatree   - store selected events in a tree for mva" << endl;
 }
 
 //
@@ -32,7 +33,7 @@ int main(int argc, char* argv[])
   TString in(""),out(""),era(""),normTag(""),method("");
   std::string systVar("");
   bool runSysts(false);
-  bool debug(false);
+  bool debug(false), skimtree(false);
   int channel(0),charge(0),flag(0);
   for(int i=1;i<argc;i++){
     string arg(argv[i]);
@@ -45,6 +46,7 @@ int main(int argc, char* argv[])
     else if(arg.find("--in")!=string::npos && i+1<argc)       { in=argv[i+1]; i++;}
     else if(arg.find("--out")!=string::npos && i+1<argc)      { out=argv[i+1]; i++;}
     else if(arg.find("--debug")!=string::npos)                { debug=true; }
+    else if(arg.find("--mvatree")!=string::npos)                { skimtree=true; }
     else if(arg.find("--normTag")!=string::npos && i+1<argc)  { normTag=argv[i+1]; i++;}
     else if(arg.find("--era")!=string::npos && i+1<argc)      { era=argv[i+1]; i++;}
     else if(arg.find("--method")!=string::npos && i+1<argc)   { method=argv[i+1]; i++;}
@@ -85,7 +87,10 @@ int main(int argc, char* argv[])
 
   //check method to run
   if(method=="ExclusiveTop::RunExclusiveTop")          RunExclusiveTop(in,out,channel,charge,normH,puH,era,debug);
-  else if(method=="VBFVectorBoson::RunVBFVectorBoson") RunVBFVectorBoson(in,out,flag,normH,puH,era,debug);
+  else if(method=="VBFVectorBoson::RunVBFVectorBoson") {
+    VBFVectorBoson myVBF(in,out,flag,normH,puH,era,debug,skimtree);
+    myVBF.RunVBFVectorBoson();
+  }
   else {
       cout << "Check method=" << method <<endl;
       printHelp();
diff --git a/TopAnalysis/data/era2017/pileupWgts.root b/TopAnalysis/data/era2017/pileupWgts.root
index c854232f8adae7d85d77d0a5bf9eebb983f83737..d2c0f6a4f67a93c12b19b49e06e49e573a09348d 100644
Binary files a/TopAnalysis/data/era2017/pileupWgts.root and b/TopAnalysis/data/era2017/pileupWgts.root differ
diff --git a/TopAnalysis/interface/SelectionTools.h b/TopAnalysis/interface/SelectionTools.h
index 6ca12a59378edc3b146293b70baefb94abb03e00..00d32028b52e5cd24b0c8fe0aeee1bdc3c3ec32b 100644
--- a/TopAnalysis/interface/SelectionTools.h
+++ b/TopAnalysis/interface/SelectionTools.h
@@ -37,7 +37,7 @@ class SelectionTool {
   std::vector<Particle> &getSelLeptons()  { return leptons_; }
   std::vector<Particle> &getVetoLeptons() { return vetoLeptons_; }
   std::vector<Particle> flaggedPhotons(MiniEvent_t &ev);
-  std::vector<Particle> selPhotons(std::vector<Particle> &flaggedPhotons,int qualBit=LOOSE, double minPt = 30., double maxEta = 2.5, std::vector<Particle> veto = {});
+  std::vector<Particle> selPhotons(std::vector<Particle> &flaggedPhotons,int qualBit=LOOSE, std::vector<Particle> leptons = {}, double minPt = 30., double maxEta = 2.5, std::vector<Particle> veto = {});
   std::vector<Particle> &getSelPhotons()  { return photons_; }
   std::vector<Jet>      getGoodJets(MiniEvent_t &ev, double minPt = 30., double maxEta = 4.7, std::vector<Particle> leptons = {},std::vector<Particle> photons = {}); // changed for the moment to VBF
   std::vector<Jet>      &getJets()        { return jets_; }
diff --git a/TopAnalysis/interface/VBFVectorBoson.h b/TopAnalysis/interface/VBFVectorBoson.h
index 75dbad926f694534dbcf157840176bd9f2f0499b..53fecd556806e5378e6ccc07ba8de88a026f0447 100644
--- a/TopAnalysis/interface/VBFVectorBoson.h
+++ b/TopAnalysis/interface/VBFVectorBoson.h
@@ -4,15 +4,194 @@
 #include "TLorentzVector.h"
 #include "TopLJets2015/TopAnalysis/interface/ObjectTools.h"
 #include "TopLJets2015/TopAnalysis/interface/SelectionTools.h"
+#include <TFile.h>
+#include <TROOT.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TSystem.h>
+#include <TGraph.h>
+#include <TGraphAsymmErrors.h>
+
+#include "TopLJets2015/TopAnalysis/interface/MiniEvent.h"
+#include "TopLJets2015/TopAnalysis/interface/CommonTools.h"
+#include "TopLJets2015/TopAnalysis/interface/VBFVectorBoson.h"
+#include "TopLJets2015/TopAnalysis/interface/EfficiencyScaleFactorsWrapper.h"
+
+#include "PhysicsTools/CandUtils/interface/EventShapeVariables.h"
+
+#include <vector>
+#include <set>
+#include <iostream>
+#include <algorithm>
+#include "TRandom3.h"
+#include "TMath.h"
+using namespace std;
 
 //Vector boson will be either Z or photon at the moment
 
-void RunVBFVectorBoson(TString filename,
-                       TString outname,
-                       Int_t anFlag,
-                       TH1F *normH, 
-                       TH1F *genPU,
-                       TString era,
-                       Bool_t debug=false);
+struct Category{
+  float MM,A,VBF,HighPt,HighPtVBF,V1J;
+  Category(float * cat){
+    MM = cat[0];
+    A = cat[1];
+    VBF = cat[2];
+    HighPt = cat[3];
+    HighPtVBF = cat[4];
+    V1J = cat[5];
+  };
+  Category(){
+    MM = 0;
+    A = 0;
+    VBF = 0;
+    HighPt = 0;
+    HighPtVBF = 0; 
+    V1J = 0; 
+  };
+  void set(float * cat){
+    MM = cat[0];
+    A = cat[1];
+    VBF = cat[2];
+    HighPt = cat[3];
+    HighPtVBF = cat[4];
+    V1J = cat[5];
+  };
+};
+
+class VBFVectorBoson{
+public:
+	VBFVectorBoson(TString filename_,
+                       TString outname_,
+                       Int_t anFlag_,
+                       TH1F *normH_, 
+                       TH1F *genPU_,
+                       TString era_,
+                       Bool_t debug_=false, Bool_t skimtree_=false):
+	filename(filename_),outname(outname_),anFlag(anFlag_), era(era_), debug(debug_), skimtree(skimtree_)
+	{
+	  normH = (TH1F*)normH_->Clone("normH_c");
+	  genPU = (TH1F*)genPU_->Clone("genPu_c");
+	  fMVATree = NULL;
+	  newTree = NULL;
+	  init();
+	  setXsecs();
+          rnd.SetSeed(123456789);
+	};
+	~VBFVectorBoson(){}
+	void init(){
+		this->readTree();
+		cout << "...producing " << outname << " from " << nentries << " events" << endl;
+		this->prepareOutput();
+		this->bookHistograms();
+		this->loadCorrections();
+		if(skimtree){
+			this->addMVAvars();
+		}
+		selector = new SelectionTool(filename, debug, triggerList,SelectionTool::VBF);
+		std::cout << "init done" << std::endl;
+	}
+
+	void saveHistos();
+	void readTree();
+	void prepareOutput();
+	void bookHistograms();
+	void setGammaZPtWeights();
+	void loadCorrections();
+	void addMVAvars();
+	void fill(MiniEvent_t ev, TLorentzVector boson, std::vector<Jet> jets, std::vector<double> cplotwgts, TString c);
+	void RunVBFVectorBoson();
+
+
+
+private:
+	TString filename, outname, baseName;
+	Int_t anFlag, nentries;
+	TH1F * normH, * genPU;
+	TString era;
+	TH1* triggerList;
+	Bool_t debug, skimtree, isQCDEMEnriched;
+	TFile * f /*inFile*/, *fMVATree, *fOut;
+	TTree * t /*inTree*/, *newTree /*MVA*/;
+  	HistTool * ht;
+	MiniEvent_t ev;
+        float sihih,chiso,r9,hoe, ystar;
+	double mindrl;
+
+	TRandom3 rnd;
+	
+  	//LUMINOSITY+PILEUP
+	LumiTools * lumi;
+  
+  	//LEPTON EFFICIENCIES
+	EfficiencyScaleFactorsWrapper * gammaEffWR;
+  
+  	//JEC/JER
+  	JECTools * jec;
+	//Photon/Z pt weights
+  	std::map<TString,TGraph *> photonPtWgts;
+  	std::map<TString,std::pair<double,double> > photonPtWgtCtr;
+
+	//Variables to be added to the MVA Tree
+	float centraleta, forwardeta, jjetas, centjy, ncentjj, dphivj0, dphivj1, dphivj2, dphivj3;
+	float evtWeight, mjj, detajj , dphijj ,jjpt;
+	float isotropy, circularity,sphericity,	aplanarity, C, D;
+	float scalarht,balance, mht, training;
+
+	/////////////////////////////////////
+	// Categorie for VBF:              //
+	//   MM:A:VBF:HighPt:HighPtVBF:V1J // 
+	/////////////////////////////////////
+	Category category;
+
+	SelectionTool * selector;
+
+	/////////////////////////////////////////
+	// To select events for training       //
+	/////////////////////////////////////////
+
+	int useForTraining(){
+	  double myRnd = rnd.Rndm();
+	  if (myRnd >= 0.5) return 1;
+	  return 0;
+	}
+
+	/////////////////////////////////////////
+	// Quick and ugly cross section getter //
+	// To be updated such as to use the    //
+	// json file                           //
+	/////////////////////////////////////////
+	std::map<TString, float> xsecRefs;
 
+	void setXsecs(){
+	  xsecRefs["MC13TeV_TTJets"        ] = 832;
+	  xsecRefs["MC13TeV_ZZ"            ] = 0.5644; 
+	  xsecRefs["MC13TeV_WZ"            ] = 47.13;  
+	  xsecRefs["MC13TeV_WW"            ] = 12.178; 
+	  xsecRefs["MC13TeV_SingleTbar_tW" ] = 35.85;  
+	  xsecRefs["MC13TeV_SingleT_tW"    ] = 35.85;  
+	  xsecRefs["MC13TeV_DY50toInf"     ] = 5765.4; 
+	  xsecRefs["MC13TeV_QCDEM_15to20"  ] = 2302200;
+	  xsecRefs["MC13TeV_QCDEM_20to30"  ] = 5352960;
+	  xsecRefs["MC13TeV_QCDEM_30to50"  ] = 9928000;
+	  xsecRefs["MC13TeV_QCDEM_50to80"  ] = 2890800;
+	  xsecRefs["MC13TeV_QCDEM_80to120" ] = 350000; 
+	  xsecRefs["MC13TeV_QCDEM_120to170"] = 62964;  
+	  xsecRefs["MC13TeV_QCDEM_170to300"] = 18810;  
+	  xsecRefs["MC13TeV_QCDEM_300toInf"] = 1350;  
+	  xsecRefs["MC13TeV_GJets_HT40to100" ] = 20790 ;
+	  xsecRefs["MC13TeV_GJets_HT100to200"] = 9238;  
+	  xsecRefs["MC13TeV_GJets_HT200to400"] = 2305;  
+	  xsecRefs["MC13TeV_GJets_HT400to600"] = 274.4; 
+	  xsecRefs["MC13TeV_GJets_HT600toInf"] = 93.46; 
+	  xsecRefs["MC13TeV_EWKZJJ"        ] = 4.32;
+	}
+	
+	float getXsec(){
+	  for (auto const& x : xsecRefs){
+	      if (filename.Contains(x.first))
+		return x.second;
+	  }
+	  return 1;
+	}
+	
+};
 #endif
diff --git a/TopAnalysis/python/Plot.py b/TopAnalysis/python/Plot.py
index 02cb734f6d360d7d3b66185acb1bfbd9aa687288..0205a832f0ce37c0109552c983e7e3350c0161c0 100644
--- a/TopAnalysis/python/Plot.py
+++ b/TopAnalysis/python/Plot.py
@@ -149,7 +149,6 @@ class Plot(object):
                 pass
 
     def show(self, outDir,lumi,noStack=False,saveTeX=False,extraText=None,noRatio=False):
-
         if len(self.mc)<1 and self.dataH is None:
             print '%s has 0 or 1 MC!' % self.name
             return
@@ -184,7 +183,6 @@ class Plot(object):
                 p1.SetTopMargin(0.05)
             p1.SetBottomMargin(0.12)
         p1.Draw()
-
         p1.SetGridx(False)
         p1.SetGridy(False) #True)
         self._garbageList.append(p1)
@@ -213,6 +211,7 @@ class Plot(object):
             if self.data is None: self.finalize()
             leg.AddEntry( self.data, self.data.GetTitle(),'p')
             nlegCols += 1
+
         for h in self.mc:
 
             #compare
@@ -224,7 +223,6 @@ class Plot(object):
                     self.mc[h].SetTitle('#splitline{%s}{#chi^{2}=%3.1f (p-val: %3.3f)}'%(self.mc[h].GetTitle(),chi2,pval))
                 else:
                     refH.SetLineWidth(2)
-
             leg.AddEntry(self.mc[h], self.mc[h].GetTitle(), 'f')
             nlegCols += 1
         if nlegCols ==0 :
@@ -238,11 +236,11 @@ class Plot(object):
         totalMC = None
         stack = ROOT.THStack('mc','mc')
         for h in reversed(self.mc):
-
+            print "Integral is %f" %self.mc[h].Integral()
             if noStack:
                 self.mc[h].SetFillStyle(0)
                 self.mc[h].SetLineColor(self.mc[h].GetFillColor())
-
+                if lumi == 1 and self.mc[h].Integral() != 1 and self.mc[h].Integral() !=0 : self.mc[h].Scale(1./self.mc[h].Integral())
             stack.Add(self.mc[h],'hist')
 
             try:
@@ -365,7 +363,6 @@ class Plot(object):
             leg.AddEntry(self.spimpose[m],self.spimpose[m].GetTitle(),'l')
         if self.data is not None : self.data.Draw('p')
 
-
         if (totalMC is not None and totalMC.GetMaximumBin() > totalMC.GetNbinsX()/2.):
             inix = 0.15
             leg.SetX1(inix)
@@ -400,7 +397,6 @@ class Plot(object):
         except:
             pass
 
-
         #holds the ratio
         c.cd()
         if not noRatio and self.dataH and len(self.mc)>0 :
@@ -487,7 +483,6 @@ class Plot(object):
                     gr.Draw('p')
                 except:
                     pass
-
         #all done
         if p1: p1.RedrawAxis()
         c.cd()
@@ -495,10 +490,12 @@ class Plot(object):
         c.Update()
 
         #save
+        if lumi == 1: frame.GetYaxis().SetRangeUser(0,0.8)
         for ext in self.plotformats : c.SaveAs(os.path.join(outDir, self.name+'.'+ext))
         if self.savelog:
             p1.cd()
             frame.GetYaxis().SetRangeUser(1,maxY*50)
+            if lumi == 1: frame.GetYaxis().SetRangeUser(0.001,1)
             p1.SetLogy()
             c.cd()
             c.Modified()
@@ -508,6 +505,8 @@ class Plot(object):
         if saveTeX : self.convertToTeX(outDir=outDir)
 
 
+
+
     def convertToTeX(self, outDir):
         if len(self.mc)==0:
             print '%s is empty' % self.name
diff --git a/TopAnalysis/scripts/plotter.py b/TopAnalysis/scripts/plotter.py
index 8bf7ea5f1eadfe1ba0e51091210a66a74c4f39f9..db04a653bfdeb51f778ff63038bc068b7510ec7f 100644
--- a/TopAnalysis/scripts/plotter.py
+++ b/TopAnalysis/scripts/plotter.py
@@ -6,9 +6,11 @@ from collections import OrderedDict
 
 from TopLJets2015.TopAnalysis.Plot import *
 
+
 """
 steer the script
 """
+	
 def main():
 
     #configuration
@@ -28,7 +30,7 @@ def main():
     parser.add_option(      '--onlyData',    dest='onlyData' ,   help='only plots containing data',     default=False,             action='store_true')
     parser.add_option(      '--saveTeX',     dest='saveTeX' ,    help='save as tex file as well',       default=False,             action='store_true')
     parser.add_option(      '--rebin',       dest='rebin',       help='rebin factor',                   default=1,                 type=int)
-    parser.add_option('-l', '--lumi',        dest='lumi' ,       help='lumi to print out',              default=12900,              type=float)
+    parser.add_option('-l', '--lumi',        dest='lumi' ,       help='lumi to print out, if == 1 draw normalized',              default=12900,              type=float)
     parser.add_option(      '--lumiSpecs',   dest='lumiSpecs',   help='lumi specifications for some channels [tag:lumi,tag2:lumi2,...]', default=None,       type=str)
     parser.add_option(      '--only',        dest='only',        help='plot only these (csv)',          default='',                type='string')
     parser.add_option(      '--skip',        dest='skip',        help='skip these samples (csv)',          default='MC13TeV_TTJets_cflip',                type='string')
@@ -64,6 +66,8 @@ def main():
     
     skipList=opt.skip.split(',')
 
+
+
     #lumi specifications per tag
     lumiSpecs={}
     if opt.lumiSpecs:
@@ -129,7 +133,6 @@ def main():
 
                 for tkey in fIn.GetListOfKeys():
                     keyIsSyst=False
-
                     try:
                         key=tkey.GetName()
 
@@ -140,7 +143,6 @@ def main():
                                 keep=True
                                 break
                         if not keep: continue
-
                         histos = []
                         obj=fIn.Get(key)
                         if obj.InheritsFrom('TH2'):
@@ -170,7 +172,8 @@ def main():
 
                         for hist in histos:
                             if not isData and not '(data)' in sp[1]: 
-
+                                print hist.GetName()
+                                print hist.GetTitle()
                                 #check if a special scale factor needs to be applied
                                 sfVal=1.0                            
                                 for procToScale in procSF:
@@ -186,7 +189,6 @@ def main():
                                     if not tag+'_' in key: continue
                                     lumi=lumiSpecs[tag]
                                     break
-                                            
                                 hist.Scale(xsec*lumi*puNormSF*sfVal)                    
                             
                             #rebin if needed
@@ -209,7 +211,8 @@ def main():
     else:                outDir = opt.outDir
     os.system('mkdir -p %s' % outDir)
     os.system('rm %s/%s'%(outDir,opt.outName))
-    for p in plots : 
+    for p in plots:
+        print plots[p].name
         plots[p].mcUnc=opt.mcUnc
         if opt.saveLog    : plots[p].savelog=True
         skipPlot=False
@@ -220,7 +223,8 @@ def main():
             if not tag+'_' in p: continue
             lumi=lumiSpecs[tag]
             break
-        if not skipPlot : plots[p].show(outDir=outDir,lumi=lumi,noStack=opt.noStack,saveTeX=opt.saveTeX)
+        #continue
+        plots[p].show(outDir=outDir,lumi=lumi,noStack=opt.noStack,saveTeX=opt.saveTeX)
         plots[p].appendTo('%s/%s'%(outDir,opt.outName))
         plots[p].reset()
 
diff --git a/TopAnalysis/scripts/runLocalAnalysis.py b/TopAnalysis/scripts/runLocalAnalysis.py
index 058493644d5cae910c18c2e3e27ff00950b34a53..55aa7e949e4a5e78d0331623e69e2c2553f61d80 100755
--- a/TopAnalysis/scripts/runLocalAnalysis.py
+++ b/TopAnalysis/scripts/runLocalAnalysis.py
@@ -15,18 +15,20 @@ Wrapper to be used when run in parallel
 """
 def RunMethodPacked(args):
 
-    method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug=args
+    method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug,mvatree=args
     print 'Running ',method,' on ',inF
     print 'Output file',outF
     print 'Selection ch=',channel,' charge=',charge,' flag=',flag,' systs=',runSysts
     print 'Normalization applied from tag=',tag
     print 'Corrections will be retrieved for era=',era
+    print 'Make the mva tree? ', mvatree
 
     try:
         cmd='analysisWrapper --era %s --normTag %s --in %s --out %s --method %s --charge %d --channel %d --flag %d --systVar %s'\
             %(era, tag, inF, outF, method, charge, channel, flag, systVar)
         if runSysts : cmd += ' --runSysts'
         if debug : cmd += ' --debug'
+        if mvatree : cmd += ' --mvatree'
         print(cmd)
         os.system(cmd)
 
@@ -63,6 +65,7 @@ def main():
     parser.add_option(      '--exactonly',   dest='exactonly',   help='match only exact sample tags to process  [%default]',    default=False,      action='store_true')
     parser.add_option(      '--outputonly',        dest='outputonly',        help='filter job submission for a csv list of output files  [%default]',             default=None,       type='string')
     parser.add_option(      '--farmappendix',        dest='farmappendix',        help='Appendix to condor FARM directory [%default]',             default='',       type='string')
+    parser.add_option(      '--mvatree',       dest='mvatree',       help='make mva tree  [%default]',                            default=False,      action='store_true')
     (opt, args) = parser.parse_args()
 
     #parse selection lists
@@ -126,7 +129,7 @@ def main():
         for systVar in varList:
             outF=opt.output
             if systVar != 'nominal' and not systVar in opt.output: outF=opt.output[:-5]+'_'+systVar+'.root'
-            task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,opt.tag,opt.debug) )
+            task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,opt.tag,opt.debug, opt.mvatree) )
     else:
 
         inputTags=getEOSlslist(directory=opt.input,prepend='')
@@ -166,7 +169,7 @@ def main():
                         continue
                     if (len(outputOnlyList) > 1 and not outF in outputOnlyList):
                         continue
-                    task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,tag,opt.debug) )
+                    task_list.append( (opt.method,inF,outF,opt.channel,opt.charge,opt.flag,opt.runSysts,systVar,opt.era,tag,opt.debug, opt.mvatree) )
                 if (opt.skipexisting and nexisting): print '--skipexisting: %s - skipping %d of %d tasks as files already exist'%(systVar,nexisting,len(input_list))
 
     #run the analysis jobs
@@ -195,7 +198,7 @@ def main():
             condor.write('+JobFlavour = "{0}"\n'.format(opt.queue))
 
             jobNb=0
-            for method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug in task_list:
+            for method,inF,outF,channel,charge,flag,runSysts,systVar,era,tag,debug,mvatree in task_list:
 
                 jobNb+=1
                 cfgFile='%s'%(os.path.splitext(os.path.basename(outF))[0])
@@ -216,6 +219,7 @@ def main():
                         %(inF, localOutF, charge, channel, era, tag, flag, method, systVar)
                     if runSysts : runOpts += ' --runSysts'
                     if debug :    runOpts += ' --debug'
+                    if mvatree :    runOpts += ' --mvatree'
                     cfg.write('python %s/src/TopLJets2015/TopAnalysis/scripts/runLocalAnalysis.py %s\n'%(cmsswBase,runOpts))
                     if '/store' in outF:
                         cfg.write('xrdcp ${WORKDIR}/%s root://eoscms//eos/cms/%s\n'%(localOutF,outF))
diff --git a/TopAnalysis/src/MiniEvent.cc b/TopAnalysis/src/MiniEvent.cc
index a8fe1f92c7d3a3dab965cae7a805d889455ab1d4..436c1deccf5acd23faae8a37429d35eec118e080 100644
--- a/TopAnalysis/src/MiniEvent.cc
+++ b/TopAnalysis/src/MiniEvent.cc
@@ -292,6 +292,28 @@ void attachToMiniEventTree(TTree *t,MiniEvent_t &ev,bool full)
   t->SetBranchAddress("j_id",       ev.j_id);
   t->SetBranchAddress("j_csv",      ev.j_csv);
   t->SetBranchAddress("j_btag",     ev.j_btag);
+  t->SetBranchAddress("j_c1_00",    ev.j_c1_00);
+  t->SetBranchAddress("j_c1_02",    ev.j_c1_02);
+  t->SetBranchAddress("j_c1_05",    ev.j_c1_05);
+  t->SetBranchAddress("j_c1_10",    ev.j_c1_10);
+  t->SetBranchAddress("j_c1_20",    ev.j_c1_20);
+  t->SetBranchAddress("j_c2_00",    ev.j_c2_00);
+  t->SetBranchAddress("j_c2_02",    ev.j_c2_02);
+  t->SetBranchAddress("j_c2_05",    ev.j_c2_05);
+  t->SetBranchAddress("j_c2_10",    ev.j_c2_10);
+  t->SetBranchAddress("j_c2_20",    ev.j_c2_20);
+  t->SetBranchAddress("j_c3_00",    ev.j_c3_00);
+  t->SetBranchAddress("j_c3_02",    ev.j_c3_02);
+  t->SetBranchAddress("j_c3_05",    ev.j_c3_05);
+  t->SetBranchAddress("j_c3_10",    ev.j_c3_10);
+  t->SetBranchAddress("j_c3_20",    ev.j_c3_20);
+  t->SetBranchAddress("j_zg",       ev.j_zg);
+  t->SetBranchAddress("j_mult",     ev.j_mult);
+  t->SetBranchAddress("j_gaptd",    ev.j_gaptd);
+  t->SetBranchAddress("j_gawidth",  ev.j_gawidth);
+  t->SetBranchAddress("j_gathrust", ev.j_gathrust);
+  t->SetBranchAddress("j_tau32",    ev.j_tau32);
+  t->SetBranchAddress("j_tau21",    ev.j_tau21);
   t->SetBranchAddress("j_deepcsv",     ev.j_deepcsv);  
   t->SetBranchAddress("j_vtxpx",    ev.j_vtxpx);
   t->SetBranchAddress("j_vtxpy",    ev.j_vtxpy);
diff --git a/TopAnalysis/src/SelectionTools.cc b/TopAnalysis/src/SelectionTools.cc
index bbceac2a963dba821233665c28374cef31533bfb..af7554cf7078fee4faf784c98ce87534c5ff7243 100644
--- a/TopAnalysis/src/SelectionTools.cc
+++ b/TopAnalysis/src/SelectionTools.cc
@@ -41,7 +41,7 @@ TString SelectionTool::flagFinalState(MiniEvent_t &ev, std::vector<Particle> pre
 
   //decide the channel based on the lepton multiplicity and set lepton collections
   std::vector<Particle> tightLeptons( selLeptons(preselLeptons,TIGHT) );
-  std::vector<Particle> tightPhotons( selPhotons(preselPhotons,TIGHT) );
+  std::vector<Particle> tightPhotons( selPhotons(preselPhotons,TIGHT, tightLeptons) );
 
   TString chTag("");
   if(anType_==TOP)
@@ -70,12 +70,12 @@ TString SelectionTool::flagFinalState(MiniEvent_t &ev, std::vector<Particle> pre
           if( ch==13*13 && fabs(mll-91)<15 && (tightLeptons[0].pt()>30 || tightLeptons[1].pt()>30)) chTag="MM";          
           leptons_=tightLeptons;
         }
-      if(tightPhotons.size()>=1) {
+     else if(tightPhotons.size()>=1) {
         chTag="A";
         photons_=tightPhotons;
       }
     }
-
+  
   //select jets based on the leptons and photon candidates
   float maxJetEta(2.4);
   if(anType_==VBF) maxJetEta=4.7;
@@ -147,6 +147,10 @@ TString SelectionTool::flagFinalState(MiniEvent_t &ev, std::vector<Particle> pre
   if(chTag=="A")
     {
       if(!hasPhotonTrigger) chTag="";
+	  if((hasEETrigger || hasETrigger) && chTag == "A"){
+		cout<< "----------------- This is EE in fact!" <<endl;
+		chTag = "";
+	  }
       if(ev.isData && !isPhotonPD_) chTag="";
     }
       
@@ -293,7 +297,7 @@ std::vector<Particle> SelectionTool::flaggedPhotons(MiniEvent_t &ev)
 }
 
 //
-std::vector<Particle> SelectionTool::selPhotons(std::vector<Particle> &photons,int qualBit,double minPt, double maxEta,std::vector<Particle> veto){
+std::vector<Particle> SelectionTool::selPhotons(std::vector<Particle> &photons,int qualBit, std::vector<Particle> leptons,double minPt, double maxEta,std::vector<Particle> veto){
   std::vector<Particle> selPhotons;
   for(size_t i =0; i<photons.size(); i++)
     {
@@ -312,7 +316,14 @@ std::vector<Particle> SelectionTool::selPhotons(std::vector<Particle> &photons,i
       }
       if(skipThisPhoton) continue;
             
-      //lepton is selected
+      // cross-cleaning with leptos
+      bool overlapsWithPhysicsObject(false);
+      for (auto& lepton : leptons) {
+	if(photons[i].p4().DeltaR(lepton.p4())<0.4) overlapsWithPhysicsObject=true;
+      }
+      
+      if(overlapsWithPhysicsObject) continue;
+      //photon is selected
       selPhotons.push_back(photons[i]);
     }
 
diff --git a/TopAnalysis/src/VBFVectorBoson.cc b/TopAnalysis/src/VBFVectorBoson.cc
index d4f665dc4a2c0e7db7e914db3991db16e8d15abc..41c251950729b58c91db7be4fca2e5cc7ea5fa8a 100644
--- a/TopAnalysis/src/VBFVectorBoson.cc
+++ b/TopAnalysis/src/VBFVectorBoson.cc
@@ -18,7 +18,6 @@
 #include <set>
 #include <iostream>
 #include <algorithm>
-
 #include "TMath.h"
 
 using namespace std;
@@ -28,142 +27,37 @@ using namespace std;
 //
 
 //
-void RunVBFVectorBoson(TString filename,
-                       TString outname,
-                       Int_t anFlag,
-                       TH1F *normH,
-                       TH1F *genPU,
-                       TString era,
-                       Bool_t debug)
+void VBFVectorBoson::RunVBFVectorBoson()
 {
-  /////////////////////
-  // INITIALIZATION //
-  ///////////////////
-  
-  //PREPARE OUTPUT
-  TString baseName=gSystem->BaseName(outname); 
-  TString dirName=gSystem->DirName(outname);
-  TFile *fOut=TFile::Open(dirName+"/"+baseName,"RECREATE");
-  fOut->cd();
-
-  //read photon/Z pT weights if required
-  std::map<TString,TGraph *> photonPtWgts;
-  std::map<TString,std::pair<double,double> > photonPtWgtCtr;
-  if(anFlag>0) {
-    TString wgtUrl("${CMSSW_BASE}/src/TopLJets2015/TopAnalysis/test/analysis/VBFVectorBoson/raw/plots/ratio_plotter.root");
-    gSystem->ExpandPathName(wgtUrl);
-    TFile *wgtF=TFile::Open(wgtUrl);
-    if(wgtF) {
-      cout << "Reading photon/Z pT weights" << endl;
-      TString pfix(baseName.Contains("Data13TeV_") ? "" : "_mc_MC");
-      photonPtWgts["VBFA"]      = new TGraph((TH1* )wgtF->Get("VBFA_vectorbosonPt_ratio/VBFA_vectorbosonPt"+pfix));
-      photonPtWgts["HighPtA"]   = new TGraph((TH1* )wgtF->Get("HighPtA_vectorbosonPt_ratio/HighPtA_vectorbosonPt"+pfix));
-      photonPtWgtCtr["VBFA"]    = std::pair<double,double>(0.0,0.0);
-      photonPtWgtCtr["HighPtA"] = std::pair<double,double>(0.0,0.0);
-      wgtF->Close();
-    } else {
-      cout << "Requested to reweight photon spectrum but could not find " << wgtUrl << endl
-           << "Proceeding without" << endl;
-    }
-  }
-
-
-  //READ TREE FROM FILE
-  MiniEvent_t ev;
-  TFile *f = TFile::Open(filename);
-  TH1 *triggerList=(TH1 *)f->Get("analysis/triggerList");
-  TTree *t = (TTree*)f->Get("analysis/data");
-  attachToMiniEventTree(t,ev,true);
-  Int_t nentries(t->GetEntriesFast());
-  if (debug) nentries = 10000; //restrict number of entries for testing
-  t->GetEntry(0);
-
-  bool isQCDEMEnriched( filename.Contains("MC13TeV_QCDEM") );
-
-  cout << "...producing " << outname << " from " << nentries << " events" << endl;
-
-  //LUMINOSITY+PILEUP
-  LumiTools lumi(era,genPU);
-  
-  //LEPTON EFFICIENCIES
-  EfficiencyScaleFactorsWrapper gammaEffWR(filename.Contains("Data13TeV"),era);
-
-  //B-TAG CALIBRATION
-  BTagSFUtil btvSF(era,"DeepCSV",BTagEntry::OperatingPoint::OP_MEDIUM,"",0);
-  
-  //JEC/JER
-  JECTools jec(era);
-  
-  //BOOK HISTOGRAMS
-  HistTool ht;
-  ht.setNsyst(0);
-  ht.addHist("puwgtctr",      new TH1F("puwgtctr",         ";Weight sums;Events",2,0,2));  
-  ht.addHist("qscale",        new TH1F("qscale",           ";Q^{2} scale;Events",100,0,2000));  
-  ht.addHist("nvtx",          new TH1F("nvtx",             ";Vertex multiplicity;Events",100,-0.5,99.5));
-  ht.addHist("vpt", 	      new TH1F("vectorbosonPt",    ";Boson p_{T}[GeV];Events",25,0,500));
-  ht.addHist("vy", 	      new TH1F("vectorbosony",     ";Boson rapidity;Events",25,0,3));
-  ht.addHist("mindrl", 	      new TH1F("mindrl",           ";min #Delta R(boson,lepton);Events",25,0,6));
-  ht.addHist("sihih", 	      new TH1F("sihih",            ";#sigma(i#eta,i#eta);Events",50,0,0.1));
-  ht.addHist("hoe", 	      new TH1F("hoe",              ";h/e;Events",25,0,0.1));
-  ht.addHist("r9", 	      new TH1F("r9",               ";r9;Events",25,0,1.0));
-  ht.addHist("chiso", 	      new TH1F("chiso",            ";Charged isolation [GeV];Events",25,0,0.10));
-  ht.addHist("vystar",        new TH1F("vectorbosonystar", ";y-(1/2)(y_{j1}+y_{j2});Events",25,0,5));
-  ht.addHist("njets",         new TH1F("njets",            ";Jet multiplicity;Events",10,-0.5,9.5));
-  ht.addHist("mjj", 	      new TH1F("mjj",              ";Dijet invariant mass [GeV];Events",40,0,4000));
-  ht.addHist("leadpt",        new TH1F("leadpt",           ";Leading jet p_{T} [GeV];Events",25,0,500));
-  ht.addHist("subleadpt",     new TH1F("subleadpt"   ,     ";Sub-leading jet p_{T} [GeV];Events",25,0,500));
-  ht.addHist("drj1b",         new TH1F("drj1b",            ";#DeltaR(j_{1},boson);Events",25,0,8));
-  ht.addHist("drj2b",         new TH1F("drj2b"   ,         ";#DeltaR(j_{2},boson);Events",25,0,8));
-  ht.addHist("leadpumva",     new TH1F("leadpumva",        ";Pileup MVA;Events",25,-1,1));
-  ht.addHist("subleadpumva",  new TH1F("subleadpumva"   ,  ";Pileup MVA;Events",25,-1,1));
-  ht.addHist("centraleta",    new TH1F("centraleta",       ";Most central jet |#eta|;Events",25,0,5));
-  ht.addHist("forwardeta",    new TH1F("forwardeta",       ";Most forward jet |#eta|;Events",25,0,5));
-  ht.addHist("dijetpt",       new TH1F("dijetpt",          ";Dijet p_{T} [GeV];Events",20,0,1000));
-  ht.addHist("detajj",        new TH1F("detajj" ,          ";#Delta#eta(J,J);Events",20,0,8));
-  ht.addHist("dphijj",        new TH1F("dphijj" ,          ";#Delta#phi(J,J) [rad];Events",20,-3.15,3.15));
-  ht.addHist("ht",            new TH1F("ht",               ";H_{T} [GeV];Events",20,0,4000));
-  ht.addHist("mht",           new TH1F("mht",              ";Missing H_{T} [GeV];Events",20,0,500));
-  ht.addHist("balance",       new TH1F("balance",          ";System p_{T} balance [GeV];Events",20,0,300));
-  ht.addHist("sphericity",    new TH1F("sphericity",       ";Sphericity;Events",20,0,1.0));
-  ht.addHist("aplanarity",    new TH1F("aplanarity",       ";Aplanarity;Events",20,0,1.0));
-  ht.addHist("C",             new TH1F("C",                ";C;Events",20,0,1.0));
-  ht.addHist("D",             new TH1F("D",                ";D;Events",20,0,1.0));
-  ht.addHist("isotropy",      new TH1F("isotropy",         ";Isotropy;Events",20,0,1.0));
-  ht.addHist("circularity",   new TH1F("circularity",      "Circularity;;Events",20,0,1.0));
-
-  std::cout << "init done" << std::endl;
-
   ///////////////////////
   // LOOP OVER EVENTS //
   /////////////////////
-  
-  //EVENT SELECTION WRAPPER
-  SelectionTool selector(filename, debug, triggerList,SelectionTool::VBF);
-  
+  float xsec = getXsec();
   for (Int_t iev=0;iev<nentries;iev++)
     {
       t->GetEntry(iev);
       if(iev%10000==0) printf ("\r [%3.0f%%] done", 100.*(float)iev/(float)nentries);
-
       std::vector<double>plotwgts(1,1.0);
-      ht.fill("qscale",ev.g_qscale,plotwgts);
+      ht->fill("qscale",ev.g_qscale,plotwgts);
       
       //assign randomly a run period
-      TString period = lumi.assignRunPeriod();
+      TString period = lumi->assignRunPeriod();
       
       //////////////////
       // CORRECTIONS //
       ////////////////      
-      jec.smearJetEnergies(ev);
+      jec->smearJetEnergies(ev);
              
       ///////////////////////////
       // RECO LEVEL SELECTION //
       /////////////////////////
-      TString chTag = selector.flagFinalState(ev);
-      std::vector<Particle> &photons     = selector.getSelPhotons(); 
-      std::vector<Particle> &leptons     = selector.getSelLeptons(); 
-      std::vector<Jet>      &alljets     = selector.getJets();  
+      TString chTag = selector->flagFinalState(ev);
+      std::vector<Particle> &photons     = selector->getSelPhotons(); 
+      std::vector<Particle> &leptons     = selector->getSelLeptons(); 
+      std::vector<Jet>      &alljets     = selector->getJets();  
       std::vector<Jet> jets;
+
+      //Pileup jet id
       for(auto j : alljets) {
         int idx=j.getJetIndex();
         int jid=ev.j_id[idx];
@@ -171,20 +65,26 @@ void RunVBFVectorBoson(TString filename,
         if(!passLoosePu) continue;
         jets.push_back(j);
       }
+      //Category selection
       if(chTag!="A" && chTag!="MM") continue;
+      float cat[6] = {0,0,0,0,0,0};
+      if(chTag == "A") cat[1] = 1;
+      if(chTag == "MM") cat[0] = 1;
 
       //jet related variables and selection
-      float mjj(jets.size()>=2 ?  (jets[0]+jets[1]).M() : 0.);
-      float detajj(jets.size()>=2 ? fabs(jets[0].Eta()-jets[1].Eta()) : -1.);
-      float dphijj(jets.size()>=2 ? jets[0].DeltaPhi(jets[1]) : -1.);
-      float jjpt=(jets.size()>=2 ? (jets[0]+jets[1]).Pt() : 0.);
-      float scalarht(0.);
+      mjj = (jets.size()>=2 ?  (jets[0]+jets[1]).M() : 0.);
+      detajj = (jets.size()>=2 ? fabs(jets[0].Eta()-jets[1].Eta()) : -1.);
+      dphijj = (jets.size()>=2 ? jets[0].DeltaPhi(jets[1]) : -1.);
+      jjpt = (jets.size()>=2 ? (jets[0]+jets[1]).Pt() : 0.);
+
+      scalarht = 0.;
       TLorentzVector mhtP4(0,0,0,0);
+      mht = 0;
       for(auto j : jets) {
         scalarht += j.Pt();
         mhtP4 += j;
       }
-      float mht(mhtP4.Pt());
+      mht = mhtP4.Pt();
       bool passJets(jets.size()>=2 && mjj>400);
       bool passVBFJetsTrigger(passJets && detajj>3.0);
 
@@ -192,19 +92,18 @@ void RunVBFVectorBoson(TString filename,
       //for the photon refine also the category according to the trigger  bit
       TLorentzVector boson(0,0,0,0);     
       bool isHighPt(false),isVBF(false),isHighPtAndVBF(false),isBosonPlusOneJet(false);
-      float sihih(0),chiso(0),r9(0),hoe(0);
-      if(chTag=="A") {
-        
+      sihih = 0, chiso = 0 ,r9 = 0, hoe = 0;
+      if(chTag=="A") {        
         boson += photons[0];
         sihih = ev.gamma_sieie[photons[0].originalReference()];
         chiso = ev.gamma_chargedHadronIso[photons[0].originalReference()];
         r9    = ev.gamma_r9[photons[0].originalReference()];
         hoe   = ev.gamma_hoe[photons[0].originalReference()];
-        isVBF    = (selector.hasTriggerBit("HLT_Photon75_R9Id90_HE10_IsoM_EBOnly_PFJetsMJJ300DEta3_v", ev.triggerBits) 
+        isVBF    = (selector->hasTriggerBit("HLT_Photon75_R9Id90_HE10_IsoM_EBOnly_PFJetsMJJ300DEta3_v", ev.triggerBits) 
                     && photons[0].Pt()>75 
                     && fabs(photons[0].Eta())<1.442
                     && passVBFJetsTrigger);
-        isHighPt = ( selector.hasTriggerBit("HLT_Photon200_v", ev.triggerBits) 
+        isHighPt = ( selector->hasTriggerBit("HLT_Photon200_v", ev.triggerBits) 
                      && photons[0].Pt()>200 );       
         isHighPtAndVBF = (isHighPt && isVBF);
         isBosonPlusOneJet=(isHighPt && alljets.size()==1);
@@ -216,8 +115,7 @@ void RunVBFVectorBoson(TString filename,
           isHighPtAndVBF = false;
         }
           
-      }
-      else {
+      } else {
         boson   += leptons[0];
         boson   += leptons[1];
         isVBF    = boson.Pt()>75 && fabs(boson.Rapidity())<1.442 && passVBFJetsTrigger;
@@ -225,20 +123,22 @@ void RunVBFVectorBoson(TString filename,
         isHighPtAndVBF = (isHighPt && isVBF);
         isBosonPlusOneJet=(isHighPt && alljets.size()==1);
       }
+
       if(!isVBF && !isHighPt && !isBosonPlusOneJet) continue;      
 
       std::vector<TString> chTags;      
-      if(isVBF)             chTags.push_back("VBF"+chTag);
-      if(isHighPt)          chTags.push_back("HighPt"+chTag);
-      if(isHighPtAndVBF)    chTags.push_back("HighPtVBF"+chTag);
-      if(isBosonPlusOneJet) chTags.push_back("V1J"+chTag);
+      if(isVBF)             {chTags.push_back("VBF"+chTag); cat[2]=1;}
+      if(isHighPt)	    {chTags.push_back("HighPt"+chTag); cat[3]=1;}
+      if(isHighPtAndVBF)    {chTags.push_back("HighPtVBF"+chTag); cat[4]=1;}
+      if(isBosonPlusOneJet) {chTags.push_back("V1J"+chTag); cat[5]=1;}
 
+      category.set(cat);
       //leptons and boson
-      double mindrl(9999.);
+      mindrl = 9999.;
       for(auto &l: leptons) mindrl=min(l.DeltaR(boson),mindrl);
 
       //system variables and event shapes
-      float ystar(0),balance(0);
+      ystar = 0,balance = 0;
       if(passJets) {
         ystar=boson.Rapidity()-0.5*(jets[0].Rapidity()+jets[1].Rapidity());
         balance=(boson+jets[0]+jets[1]).Pt();
@@ -248,6 +148,13 @@ void RunVBFVectorBoson(TString filename,
       for(size_t ij=0; ij<min(size_t(2),jets.size());ij++) {
         inputVectors.push_back( math::XYZVector(jets[ij].Px(),jets[ij].Py(),jets[ij].Pz()) );
       EventShapeVariables esv(inputVectors);
+      isotropy =  esv.isotropy();
+      circularity =  esv.circularity();
+      sphericity =    esv.sphericity(1.);
+      aplanarity =   esv.aplanarity(1.);
+      C =            esv.C(1.);
+      D =            esv.D(1.);
+
 
       ////////////////////
       // EVENT WEIGHTS //
@@ -259,23 +166,23 @@ void RunVBFVectorBoson(TString filename,
         wgt  = (normH? normH->GetBinContent(1) : 1.0);
             
         // pu weight
-        ht.fill("puwgtctr",0,plotwgts);
-        double puWgt(lumi.pileupWeight(ev.g_pu,period)[0]);
+        ht->fill("puwgtctr",0,plotwgts);
+        double puWgt(lumi->pileupWeight(ev.g_pu,period)[0]);
         std::vector<double>puPlotWgts(1,puWgt);
-        ht.fill("puwgtctr",1,puPlotWgts);
+        ht->fill("puwgtctr",1,puPlotWgts);
         
         // photon trigger*selection weights
         float trigSF(1.0), selSF(1.0);
         if(chTag=="A")
           {
-            trigSF *= gammaEffWR.getTriggerCorrection({},photons,{}, period).first;
-            selSF  *= gammaEffWR.getOfflineCorrection(photons[0], period).first;
+            trigSF *= gammaEffWR->getTriggerCorrection({},photons,{}, period).first;
+            selSF  *= gammaEffWR->getOfflineCorrection(photons[0], period).first;
           }
         else
           {
-            trigSF *=gammaEffWR.getTriggerCorrection(leptons,{},{}, period).first;
-            selSF  *=gammaEffWR.getOfflineCorrection(leptons[0], period).first;
-            selSF  *=gammaEffWR.getOfflineCorrection(leptons[1], period).first;
+            trigSF *=gammaEffWR->getTriggerCorrection(leptons,{},{}, period).first;
+            selSF  *=gammaEffWR->getOfflineCorrection(leptons[0], period).first;
+            selSF  *=gammaEffWR->getOfflineCorrection(leptons[1], period).first;
           }
         wgt *= puWgt*trigSF*selSF;
         
@@ -289,7 +196,6 @@ void RunVBFVectorBoson(TString filename,
 
       //control histograms
       for( auto c : chTags) {
-
         std::vector<double> cplotwgts(plotwgts);
 
         //photon pT weighting
@@ -305,48 +211,11 @@ void RunVBFVectorBoson(TString filename,
           cplotwgts[0]*=photonPtWgt;
         } 
 
-        ht.fill("nvtx",   ev.nvtx,          cplotwgts,c);        
-
-        //boson histos
-        ht.fill("vpt",    boson.Pt(),       cplotwgts,c);
-        ht.fill("vy",     boson.Rapidity(), cplotwgts,c);   
-        ht.fill("sihih",  sihih,            cplotwgts,c);   
-        ht.fill("r9",     r9,               cplotwgts,c);   
-        ht.fill("hoe",    hoe,             cplotwgts,c);   
-        ht.fill("chiso",  chiso,            cplotwgts,c);   
-        ht.fill("mindrl", mindrl,           cplotwgts,c);   
-
-        //jet histos
-        double minEta(9999),maxEta(-9999);
-        for(size_t ij=0; ij<min(size_t(2),jets.size());ij++) {
-          TString jtype(ij==0?"lead":"sublead");
-          ht.fill(jtype+"pt",       jets[ij].Pt(),        cplotwgts,c);          
-          ht.fill(jtype+"pumva",    jets[ij].getPUMVA(),  cplotwgts,c);
-          ht.fill("dr"+jtype+"b",   jets[ij].DeltaR(boson),  cplotwgts,c);
-          minEta=min(minEta,fabs(jets[ij].Eta()));
-          maxEta=max(maxEta,fabs(jets[ij].Eta()));
-        }
-        ht.fill("njets",        jets.size(), cplotwgts,c);
-        ht.fill("ht",           scalarht,    cplotwgts,c);
-        ht.fill("mht",          mht,         cplotwgts,c);
-        ht.fill("centraleta",   minEta,      cplotwgts,c);
-        ht.fill("forwardeta",   maxEta,      cplotwgts,c);
-        ht.fill("dijetpt",      jjpt,        cplotwgts,c);
-        ht.fill("detajj",       detajj,      cplotwgts,c);
-        ht.fill("dphijj",       dphijj,      cplotwgts,c);
-        ht.fill("mjj", 	        mjj,         cplotwgts,c);
-
-        //visible system histos
-        ht.fill("vystar",       ystar,              cplotwgts,c);        
-        ht.fill("balance",      balance,            cplotwgts,c);
-        ht.fill("isotropy",     esv.isotropy(),     cplotwgts,c);
-        ht.fill("circularity",  esv.circularity(),  cplotwgts,c);
-        ht.fill("sphericity",   esv.sphericity(1.), cplotwgts,c);
-        ht.fill("aplanarity",   esv.aplanarity(1.), cplotwgts,c);
-        ht.fill("C",            esv.C(1.),          cplotwgts,c);
-        ht.fill("D",            esv.D(1.),          cplotwgts,c);
-
-        }
+	//What is the final weight? 0 or 1 in the array?
+	evtWeight = cplotwgts[0]*xsec;
+	training = useForTraining(); 
+	fill( ev,  boson,  jets,  cplotwgts, c);
+       }
       }
     }
   
@@ -361,23 +230,250 @@ void RunVBFVectorBoson(TString filename,
   }
   
   //save histos to file  
-  fOut->cd();
-  for (auto& it : ht.getPlots())  { 
-    for(auto &wit : photonPtWgtCtr){
-      if(!it.first.Contains(wit.first)) continue;
-      cout << "Scaling " << it.first << " by "<< wit.second.first <<endl;
-      it.second->Scale(wit.second.first);
-      break;
-    }
-    it.second->SetDirectory(fOut); it.second->Write(); 
-  }
-  for (auto& it : ht.get2dPlots())  { 
-    for(auto &wit : photonPtWgtCtr){
-      if(!it.first.Contains(wit.first)) continue;
-      it.second->Scale(wit.second.first);
-      break;
-    }
-    it.second->SetDirectory(fOut); it.second->Write(); 
+  saveHistos();
+  if(skimtree){
+  	fMVATree->cd();
+  	newTree->Write();
+  	fMVATree->Close();
   }
-  fOut->Close();
+}
+void VBFVectorBoson::saveHistos(){
+	fOut->cd();
+	for (auto& it : ht->getPlots())  { 
+		for(auto &wit : photonPtWgtCtr){
+			if(!it.first.Contains(wit.first)) continue;
+			cout << "Scaling " << it.first << " by "<< wit.second.first <<endl;
+			it.second->Scale(wit.second.first);
+			break;
+		}
+		it.second->SetDirectory(fOut); it.second->Write(); 
+	}
+	for (auto& it : ht->get2dPlots())  { 
+		for(auto &wit : photonPtWgtCtr){
+			if(!it.first.Contains(wit.first)) continue;
+			it.second->Scale(wit.second.first);
+			break;
+		}
+		it.second->SetDirectory(fOut); it.second->Write(); 
+	}
+  	fOut->Close();		
+}
+
+void VBFVectorBoson::readTree(){
+	f = TFile::Open(filename);
+	triggerList=(TH1 *)f->Get("analysis/triggerList");
+	t = (TTree*)f->Get("analysis/data");
+	attachToMiniEventTree(t,ev,true);
+	nentries = t->GetEntriesFast();
+	if (debug) nentries = 10000; //restrict number of entries for testing
+	t->GetEntry(0);
+	isQCDEMEnriched = filename.Contains("MC13TeV_QCDEM");
+}
+
+void VBFVectorBoson::prepareOutput(){
+	baseName=gSystem->BaseName(outname); 
+	TString dirName=gSystem->DirName(outname);
+
+	if(skimtree){
+		fMVATree=TFile::Open(dirName+"/MVATree_"+baseName,"RECREATE");
+	  	newTree = t->CloneTree(0);
+	}
+	fOut=TFile::Open(dirName+"/"+baseName,"RECREATE");
+	fOut->cd();
+}
+
+void VBFVectorBoson::bookHistograms(){
+	ht = new HistTool(0);
+	ht->addHist("puwgtctr",      new TH1F("puwgtctr",         ";Weight sums;Events",2,0,2));  
+	ht->addHist("qscale",        new TH1F("qscale",           ";Q^{2} scale;Events",100,0,2000));  
+	ht->addHist("nvtx",          new TH1F("nvtx",             ";Vertex multiplicity;Events",100,-0.5,99.5));  
+	ht->addHist("vpt", 	     new TH1F("vectorbosonPt",    ";Boson p_{T}[GeV];Events",25,0,500));  
+	ht->addHist("vy", 	     new TH1F("vectorbosony",     ";Boson rapidity;Events",25,0,3));  
+	ht->addHist("mindrl", 	     new TH1F("mindrl",           ";min #Delta R(boson,lepton);Events",25,0,6));  
+	ht->addHist("sihih", 	     new TH1F("sihih",            ";#sigma(i#eta,i#eta);Events",50,0,0.1));  
+	ht->addHist("hoe", 	     new TH1F("hoe",              ";h/e;Events",25,0,0.1));  
+	ht->addHist("r9", 	     new TH1F("r9",               ";r9;Events",25,0,1.0));  
+	ht->addHist("chiso", 	     new TH1F("chiso",            ";Charged isolation [GeV];Events",25,0,0.10));  
+	ht->addHist("vystar",        new TH1F("vectorbosonystar", ";y-(1/2)(y_{j1}+y_{j2});Events",25,0,5));  
+	ht->addHist("njets",         new TH1F("njets",            ";Jet multiplicity;Events",10,-0.5,9.5));  
+	ht->addHist("mjj", 	     new TH1F("mjj",              ";Dijet invariant mass [GeV];Events",40,0,4000));  
+	ht->addHist("leadpt",        new TH1F("leadpt",           ";Leading jet p_{T} [GeV];Events",25,0,500));  
+	ht->addHist("subleadpt",     new TH1F("subleadpt"   ,     ";Sub-leading jet p_{T} [GeV];Events",25,0,500));  
+	ht->addHist("drj1b",         new TH1F("drj1b",            ";#DeltaR(j_{1},boson);Events",25,0,8));  
+	ht->addHist("drj2b",         new TH1F("drj2b"   ,         ";#DeltaR(j_{2},boson);Events",25,0,8));  
+	ht->addHist("leadpumva",     new TH1F("leadpumva",        ";Pileup MVA;Events",25,-1,1));  
+	ht->addHist("subleadpumva",  new TH1F("subleadpumva"   ,  ";Pileup MVA;Events",25,-1,1));  
+	ht->addHist("centraleta",    new TH1F("centraleta",       ";Most central jet |#eta|;Events",25,0,5));  
+	ht->addHist("forwardeta",    new TH1F("forwardeta",       ";Most forward jet |#eta|;Events",25,0,5));  
+	ht->addHist("dijetpt",       new TH1F("dijetpt",          ";Dijet p_{T} [GeV];Events",20,0,1000));  
+	ht->addHist("detajj",        new TH1F("detajj" ,          ";#Delta#eta(J,J);Events",20,0,8));  
+	ht->addHist("dphijj",        new TH1F("dphijj" ,          ";#Delta#phi(J,J) [rad];Events",20,-3.15,3.15));  
+	ht->addHist("ht",            new TH1F("ht",               ";H_{T} [GeV];Events",20,0,4000));  
+	ht->addHist("mht",           new TH1F("mht",              ";Missing H_{T} [GeV];Events",20,0,500));  
+	ht->addHist("balance",       new TH1F("balance",          ";System p_{T} balance [GeV];Events",20,0,300));  
+	ht->addHist("sphericity",    new TH1F("sphericity",       ";Sphericity;Events",20,0,1.0));  
+	ht->addHist("aplanarity",    new TH1F("aplanarity",       ";Aplanarity;Events",20,0,1.0));  
+	ht->addHist("C",             new TH1F("C",                ";C;Events",20,0,1.0));  
+	ht->addHist("D",             new TH1F("D",                ";D;Events",20,0,1.0));  
+	ht->addHist("isotropy",      new TH1F("isotropy",         ";Isotropy;Events",20,0,1.0));  
+	ht->addHist("circularity",   new TH1F("circularity",      "Circularity;;Events",20,0,1.0));
+	// Study of jet variables
+	ht->addHist("jet_c1_00", 	  new TH1F("jet_c1_00",          ";Jet shape var. c1_00;Jets",100,-1,1));  
+	ht->addHist("jet_c1_02", 	  new TH1F("jet_c1_02",          ";Jet shape var. c1_02;Jets",100,-1,1));  
+	ht->addHist("jet_c1_05", 	  new TH1F("jet_c1_05",          ";Jet shape var. c1_05;Jets",100,-1,1));  
+	ht->addHist("jet_c2_00", 	  new TH1F("jet_c2_00",          ";Jet shape var. c2_00;Jets",100,-1,1));  
+	ht->addHist("jet_c2_02", 	  new TH1F("jet_c2_02",          ";Jet shape var. c2_02;Jets",100,-1,1));  
+	ht->addHist("jet_c2_05", 	  new TH1F("jet_c2_05",          ";Jet shape var. c2_05;Jets",100,-1,1));  
+	ht->addHist("jet_c3_00", 	  new TH1F("jet_c3_00",          ";Jet shape var. c3_00;Jets",100,-1,1));  
+	ht->addHist("jet_c3_02", 	  new TH1F("jet_c3_02",          ";Jet shape var. c3_02;Jets",100,-1,1));  
+	ht->addHist("jet_c3_05", 	  new TH1F("jet_c3_05",          ";Jet shape var. c3_05;Jets",100,-1,1));  
+	ht->addHist("jet_zg", 	  	  new TH1F("jet_zg",          ";Jet shape var. zg;Jets",100,-1,1));  
+	ht->addHist("jet_gaptd", 	  new TH1F("jet_gaptd",          ";Jet shape var. gaptd;Jets",100,-1,1));  
+	ht->addHist("jet_gawidth",        new TH1F("jet_gawidth",          ";Jet shape var. gawidth;Jets",100,-1,1));
+	//additional variables from https://link.springer.com/content/pdf/10.1140/epjc/s10052-017-5315-6.pdf
+	ht->addHist("jjetas", 	  	  new TH1F("jjetas",          ";#eta_{j1}#eta_{j2};Events",200,-25,25));  
+	ht->addHist("centjy",		  new TH1F("centjy",          ";Central jet rapidity;Jets",25,0,3));  
+	ht->addHist("ncentj", 	  	  new TH1F("ncentjj",          ";Number of central jets;Events",10,-0.5,9.5));  
+	ht->addHist("dphivj0", 	  new TH1F("dphivj0",          ";#Delta#phi(V,j0);Jets",20,0,4));  
+	ht->addHist("dphivj1", 	  new TH1F("dphivj1",          ";#Delta#phi(V,j1);Jets",20,0,4));  
+	ht->addHist("dphivj2", 	  new TH1F("dphivj2",          ";#Delta#phi(V,j2);Jets",20,0,4));  
+	ht->addHist("dphivj3", 	  new TH1F("dphivj3",          ";#Delta#phi(V,j3);Jets",20,0,4));
+}
+void VBFVectorBoson::setGammaZPtWeights(){
+	TString wgtUrl("${CMSSW_BASE}/src/TopLJets2015/TopAnalysis/test/analysis/VBFVectorBoson/raw/plots/ratio_plotter.root");
+	gSystem->ExpandPathName(wgtUrl);
+	TFile *wgtF=TFile::Open(wgtUrl);
+	if(wgtF) {
+		cout << "Reading photon/Z pT weights" << endl;
+		TString pfix(baseName.Contains("Data13TeV_") ? "" : "_mc_MC");
+		photonPtWgts["VBFA"]      = new TGraph((TH1* )wgtF->Get("VBFA_vectorbosonPt_ratio/VBFA_vectorbosonPt"+pfix));
+		photonPtWgts["HighPtA"]   = new TGraph((TH1* )wgtF->Get("HighPtA_vectorbosonPt_ratio/HighPtA_vectorbosonPt"+pfix));
+		photonPtWgtCtr["VBFA"]    = std::pair<double,double>(0.0,0.0);
+		photonPtWgtCtr["HighPtA"] = std::pair<double,double>(0.0,0.0);
+		wgtF->Close();
+	} else {
+		cout << "Requested to reweight photon spectrum but could not find " << wgtUrl << endl
+           		<< "Proceeding without" << endl;
+	}
+}
+void VBFVectorBoson::loadCorrections(){
+	lumi = new LumiTools(era,genPU);
+	gammaEffWR = new EfficiencyScaleFactorsWrapper(filename.Contains("Data13TeV"),era);
+	jec = new JECTools(era);
+	if(anFlag>0) this->setGammaZPtWeights();
+}
+void VBFVectorBoson::addMVAvars(){
+	newTree->Branch("centralEta", &centraleta);
+	newTree->Branch("mjj", &mjj);
+	newTree->Branch("detajj", &detajj);
+	newTree->Branch("jjpt", &jjpt);
+	newTree->Branch("dphijj", &dphijj);
+	newTree->Branch("forwardeta", &forwardeta);
+	newTree->Branch("jjetas", &jjetas);
+	newTree->Branch("centjy", &centjy);
+	newTree->Branch("ncentjj", &ncentjj);
+	newTree->Branch("dphivj0", &dphivj0);
+	newTree->Branch("dphivj1", &dphivj1);
+	newTree->Branch("dphivj2", &dphivj2);
+	newTree->Branch("dphivj3", &dphivj3);
+	newTree->Branch("evtWeight", &evtWeight);
+        newTree->Branch("mht", &mht);
+        newTree->Branch("balance", &balance);
+        newTree->Branch("ht", &scalarht);
+	newTree->Branch("isotropy", &isotropy);
+	newTree->Branch("circularity",&circularity);
+	newTree->Branch("sphericity",&sphericity);
+	newTree->Branch("aplanarity",&aplanarity);
+	newTree->Branch("C",&C);
+	newTree->Branch("D",&D);
+	newTree->Branch("training",&training);
+	newTree->Branch("category", &category, "MM:A:VBF:HighPt:HighPtVBF:V1J");
+}
+
+void VBFVectorBoson::fill(MiniEvent_t ev, TLorentzVector boson, std::vector<Jet> jets, std::vector<double> cplotwgts, TString c){
+	ht->fill("nvtx",   ev.nvtx,          cplotwgts,c);        
+
+	//boson histos
+	ht->fill("vpt",    boson.Pt(),       cplotwgts,c);
+	ht->fill("vy",     boson.Rapidity(), cplotwgts,c);   
+	ht->fill("sihih",  sihih,            cplotwgts,c);   
+	ht->fill("r9",     r9,               cplotwgts,c);   
+	ht->fill("hoe",    hoe,             cplotwgts,c);   
+	ht->fill("chiso",  chiso,            cplotwgts,c);   
+	ht->fill("mindrl", mindrl,           cplotwgts,c);   
+
+        //jet histos
+        centraleta = 9999;
+	forwardeta = -9999;
+        for(size_t ij=0; ij<min(size_t(2),jets.size());ij++) {
+	  TString jtype(ij==0?"lead":"sublead");
+	  ht->fill(jtype+"pt",       jets[ij].Pt(),        cplotwgts,c);          
+	  ht->fill(jtype+"pumva",    jets[ij].getPUMVA(),  cplotwgts,c);
+	  ht->fill("dr"+jtype+"b",   jets[ij].DeltaR(boson),  cplotwgts,c);
+	  ht->fill("jet_c1_00", 	ev.j_c1_00[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c1_02", 	ev.j_c1_02[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c1_05", 	ev.j_c1_05[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c2_00", 	ev.j_c2_00[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c2_02", 	ev.j_c2_02[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c2_05",	ev.j_c2_05[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c3_00", 	ev.j_c3_00[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c3_02",	ev.j_c3_02[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_c3_05",	ev.j_c3_05[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_zg", 		ev.j_zg[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_gaptd", 	ev.j_gaptd[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  ht->fill("jet_gawidth", ev.j_gawidth[jets[ij].getJetIndex()]	  ,  cplotwgts,c);
+	  centraleta=min(centraleta,float(fabs(jets[ij].Eta())));
+	  forwardeta=max(forwardeta,float(fabs(jets[ij].Eta())));
+	}
+	jjetas = 9999;
+	dphivj0 = 9999; dphivj1 = 9999;
+	if(jets.size() >= 2){
+	  jjetas = jets[0].Eta()*jets[1].Eta();
+	  dphivj0 = fabs(jets[0].DeltaPhi(boson));
+	  dphivj1 = fabs(jets[1].DeltaPhi(boson));
+	  ht->fill("jjetas",  jjetas,   cplotwgts,c);
+	  ht->fill("dphivj0", dphivj0 ,  cplotwgts,c);
+	  ht->fill("dphivj1", dphivj1 ,  cplotwgts,c);
+	}
+	dphivj2 = 9999; dphivj3 = 9999;
+	centjy = 9999; ncentjj = 0;  
+	if(jets.size() > 2){
+		dphivj2 = fabs(jets[2].DeltaPhi(boson));
+		ht->fill("dphivj2", dphivj2 ,  cplotwgts,c);
+		for(unsigned int iJet = 2; iJet < jets.size(); iJet++){	
+			float dy = fabs(jets[0].Rapidity() - jets[1].Rapidity())/2;
+			float sumy = (jets[0].Rapidity() + jets[1].Rapidity())/2;
+			if(fabs(jets[iJet].Rapidity() - sumy) < dy){
+				centjy =  jets[iJet].Rapidity();
+				ht->fill("centjy",centjy,  cplotwgts,c);
+				ncentjj++;
+			}
+		}
+		ht->fill("ncentj", ncentjj, cplotwgts, c);
+	}
+	if(jets.size() > 3){
+		dphivj3 =  fabs(jets[3].DeltaPhi(boson));
+		ht->fill("dphivj3", dphivj3 , cplotwgts,c);
+	}
+	ht->fill("njets",        jets.size(), cplotwgts,c);
+	ht->fill("ht",           scalarht,    cplotwgts,c);
+	ht->fill("mht",          mht,         cplotwgts,c);
+	ht->fill("centraleta",   centraleta,      cplotwgts,c);
+	ht->fill("forwardeta",   forwardeta,      cplotwgts,c);
+	ht->fill("dijetpt",      jjpt,        cplotwgts,c);
+	ht->fill("detajj",       detajj,      cplotwgts,c);
+	ht->fill("dphijj",       dphijj,      cplotwgts,c);
+	ht->fill("mjj", 	        mjj,         cplotwgts,c);
+	
+	//visible system histos
+	ht->fill("vystar",       ystar,              cplotwgts,c);        
+	ht->fill("balance",      balance,            cplotwgts,c);
+	ht->fill("isotropy",     isotropy,     cplotwgts,c);
+	ht->fill("circularity",  circularity,  cplotwgts,c);
+	ht->fill("sphericity",   sphericity, cplotwgts,c);
+	ht->fill("aplanarity",   aplanarity, cplotwgts,c);
+	ht->fill("C",            C,          cplotwgts,c);
+	ht->fill("D",            D,          cplotwgts,c);
+	cout<< "hist filling is done ... "<<endl;
+	if(skimtree) newTree->Fill();
 }
diff --git a/TopAnalysis/test/analysis/steerVBFVectorBoson.sh b/TopAnalysis/test/analysis/steerVBFVectorBoson.sh
index b9188671a3536ddd26d788b23681ea61f2cb8925..5a1425667360ceb6b94887fee13a9960f5d704fd 100644
--- a/TopAnalysis/test/analysis/steerVBFVectorBoson.sh
+++ b/TopAnalysis/test/analysis/steerVBFVectorBoson.sh
@@ -27,24 +27,25 @@ NC='\e[0m'
 case $WHAT in
 
     TESTSEL )
-        #input=${eosdir}/MC13TeV_DY50toInf/MergedMiniEvents_0_ext0.root
-        #output=MC13TeV_DY.root
-        #tag="--tag MC13TeV_DY50toInf"
-        input=${eosdir}/MC13TeV_GJets_HT100to200/MergedMiniEvents_0_ext0.root
-        #output=MC13TeV_GJets_HT100to200.root
-        #tag="--tag MC13TeV_GJets_HT100to200 --flag 1"
-        output=MC13TeV_GJets_HT100to200_raw.root
-        tag="--tag MC13TeV_GJets_HT100to200"
+        input=${eosdir}/MC13TeV_DY50toInf/MergedMiniEvents_4_ext0.root
+        output=MC13TeV_DY4Jets50toInf.root
+        tag="--tag MC13TeV_DY50toInf"
+        #input=${eosdir}/Data13TeV_SingleMuon_2017D/MergedMiniEvents_0_ext0.root
+        #output=Data13TeV_SingleMuon_2017D.root
+        #input=${eosdir}/MC13TeV_GJets_HT100to200_DR04/MergedMiniEvents_0_ext0.root
+        #output=MC13TeV_GJets_HT100to200_DR04.root \
         #input=${eosdir}/Data13TeV_SingleMuon_2017C/MergedMiniEvents_0_ext0.root
         #output=Data13TeV_SingleMuon_2017C.root
 	python scripts/runLocalAnalysis.py \
             -i ${input} -o ${output} ${tag} \
-            --njobs 1 -q local --debug \
+            --njobs 1 -q local --debug --mvatree \
             --era era2017 -m VBFVectorBoson::RunVBFVectorBoson --ch 0 --runSysts;
         ;;
 
     SEL )
 	python scripts/runLocalAnalysis.py -i ${eosdir} \
+#			--only DY,EWKZJJ,GJets_HT,QCDEM \
+#            -o ${outdir} \
             -o ${outdir}/raw \
             -q ${queue} \
             --era era2017 -m VBFVectorBoson::RunVBFVectorBoson --ch 0 --runSysts;
@@ -62,10 +63,12 @@ case $WHAT in
 	;;
 
     PLOT )
+#	commonOpts="-i ${outdir} --puNormSF puwgtctr -l 1  --saveLog"
+#	python scripts/plotter.py ${commonOpts} -j data/era2017/vbf_samples.json  --noStack --skip TT,ZZ,WW,WZ,Single,QCD,GJets,DY1Jets,DY2Jets,DY3Jets,DY4Jets,DY50toInf_HT -O ${outdir}/plots_DY ;
 	commonOpts="-i ${outdir}/${EXTRA} --puNormSF puwgtctr -l ${fulllumi}  --saveLog --mcUnc ${lumiUnc} --lumiSpecs VBFA:${vbflumi}"
 	python scripts/plotter.py ${commonOpts} -j data/era2017/vbf_samples.json; 
-        #python test/analysis/computeVBFRatios.py -i ${outdir}/${EXTRA}/plots/plotter.root -o ${outdir}/${EXTRA}/plots/ratio_plotter.root
         python test/analysis/computeVBFRatios.py -t -i ${outdir}/${EXTRA}/plots/plotter.root -o ${outdir}/${EXTRA}/plots/trigger_ratio_plotter.root
+
 	;;
 
     WWW )