//////////////////////////////////////////////////////////////
//  File: PlotVBF.cpp
//  Created: July 30, 2007 based on Plot() in wz_cafreco
///////////////////////////////////////////////////////////////


#include "vbf_cafe/PlotVBF.hpp"
#include "vbf_cafe/topology.cpp"
#include "cafe/Event.hpp"
#include <algorithm>

using namespace std;

namespace vbf_cafe {

  PlotVBF::PlotVBF(const char *name) : cafe::Processor(name)
  {
  cafe::Config config(name);

  _debug_level = config.get("debug_level",0);
  _root_file_name = config.get("root_file_name", "VBFPlot.root");

  _muonBranch = config.get("muonBranch", "selectedMuons") ;
  _jetBranch = config.get("jetBranch", "selectedJets") ;

  _zmumu_mass_low = config.get("Zmumu_mass_low",65);
  _zmumu_mass_high = config.get("Zmumu_mass_high",115);
  _zee_mass_low = config.get("Zee_mass_low",65);
  _zee_mass_high = config.get("Zee_mass_high",115);

  _usealljets     =config.get("UseAllJetsForTopo",   true);

  _jet_pt_cut = config.get("leading_jet_pt_cut", 20);
  _oppjet_cut = config.get("opp_jet_phi_cut", TMath::Pi()/2);

  nerr=0;
  nevt=0;
  }

  PlotVBF::~PlotVBF()
  {
  }


  bool PlotVBF::processEvent(cafe::Event& event)  {
  nevt++;
  int NWmu=0, NZmu=0, NWe=0, NZee=0;
  bool good_zmumu, good_zee, good_wenu, good_wmunu;
  float VB_phi, VB_pt, VB_rap, lep_high_pt, lep_low_pt;
  float VB2mu_phi, VB2mu_pt, VB2e_phi, VB2e_pt;
  float lep_eta_p_1, lep_eta_p_2, lep_eta_d_1, lep_eta_d_2;
  

  ///////////////// Get the event weight: ////////////////////
  event_weight=1;
  event.get("StatPointer", statPointer);
  event_weight=statPointer.pointer()->eventWeight();

  float trigger_weight ;
  event.get("TriggerProbability", trigger_weight) ;
  if(_debug_level>2) {
    cout << "Event Weight " << event_weight << endl;;
    cout << "Trigger Weight " << trigger_weight << endl;
  }


  //////////////// Debug muon info: ////////////////////////////////////
  if(_debug_level>2) {
    muons = event.getCollection<TMBMuon>(_muonBranch.c_str());

    
    cout<<"Found: "<<muons.size()<<" Muons."<<endl;
    if(_debug_level>5) {
      cout<<"Found: "<<muons.size()<<" Muons."<<endl;
      for (int i=0; i<muons.size(); i++ ){
	cout <<i<<" "<<muons[i].isMedium()<<" "<<muons[i].hasCentral();
	if(muons[i].isMedium() && muons[i].hasCentral() ) {
	  TMBTrack * track = muons[i].GetChargedTrack();
	  if(track==0){
	    cout <<" Error track!"<<endl;
	    continue;
	  }
	  TString name = track->IsA()->GetName();
	  if (name == "TMBCellContainer") {
	    std::cout << "track is not a TMBTrack --> skipping it" << std::endl;
	    continue;
	  }
	  cout <<" "<<muons[i].GetChargedTrack()->Pt();
	}//has central
	cout<<endl;
      }
    }
  }

  //////////////// Z->mumu basic plots: ////////////////////////////////
  Zmumu_vec.clear();
  good_zmumu=false;
  if(event.get("Zmumu candidates", Zmumu_vec)) {
    
    NZmu = Zmumu_vec.size();
    if(_debug_level>2)  cout <<"Found "<<NZmu<<" Z->mumu candidates"<<endl;
    
    zmu_num_plot->Fill(NZmu, event_weight);
    
    if(_debug_level>0) {
      if(NZmu>1) {
	for(int i=0; i<NZmu; i++) {
	  cout <<"Z->mumu_vec "<<i<<" "<<Zmumu_vec[i]->InvMass()<<" "<<Zmumu_vec[i]->LepPhi(0)<<" "<<Zmumu_vec[i]->LepPhi(1)<<endl;
	}
      }
    }
    
    if(NZmu>0 && Zmumu_vec[0]->CutsFailed()==0){
      zmu_mass_plot->Fill(Zmumu_vec[0]->InvMass(), event_weight);
      
      if(Zmumu_vec[0]->InvMass() > _zmumu_mass_low && Zmumu_vec[0]->InvMass() < _zmumu_mass_high) {
	good_zmumu=true;
	VB_rap = Zmumu_vec[0]->ZRap();
	zmu_rap_plot->Fill(VB_rap, event_weight);
	VB_pt = Zmumu_vec[0]->ZPt();
	zmu_pt_plot->Fill(VB_pt, event_weight);
	VB_phi = Zmumu_vec[0]->ZPhi();
	zmu_phi_plot->Fill(VB_phi, event_weight);
	
	if(Zmumu_vec[0]->LepPt(1) > Zmumu_vec[0]->LepPt(0)) {
	  cout <<"Wrong pT's!!"<<Zmumu_vec[0]->LepPt(0)<< " " << Zmumu_vec[0]->LepPt(1)<<endl;
	  cout <<"# SMT hits:"<<Zmumu_vec[0]->nsmt(0)<< " " << Zmumu_vec[0]->nsmt(1)<<endl;
	}
	lep_high_pt = Zmumu_vec[0]->LepPt(0);
	lep_low_pt = Zmumu_vec[0]->LepPt(1);
	lep_eta_d_1 = Zmumu_vec[0]->LepDetEta(0);
	lep_eta_d_2 = Zmumu_vec[0]->LepDetEta(1);
	lep_eta_p_1 = Zmumu_vec[0]->LepEta(0);
	lep_eta_p_2 = Zmumu_vec[0]->LepEta(1);
	zmuon_highpt_plot->Fill( lep_high_pt, event_weight );
	zmuon_lowpt_plot->Fill( lep_low_pt, event_weight );
	zmuon_eta_p_plot->Fill( lep_eta_p_1, event_weight );
	zmuon_eta_p_plot->Fill( lep_eta_p_2, event_weight );
	zmuon_eta_d_plot->Fill( lep_eta_d_1, event_weight );
	zmuon_eta_d_plot->Fill( lep_eta_d_2, event_weight );
	float dphi = fabs(Zmumu_vec[0]->LepPhi(0) - Zmumu_vec[0]->LepPhi(1));
	if(dphi>TMath::Pi()) dphi = 2*TMath::Pi() - dphi;
	zmuon_dphi_plot->Fill(dphi, event_weight);
	
	zmuon_dz0_plot->Fill(Zmumu_vec[0]->LepZ0(0) - Zmumu_vec[0]->LepZ0(1), event_weight);
	zmuon_z0_plot->Fill(0.5*(Zmumu_vec[0]->LepZ0(0) + Zmumu_vec[0]->LepZ0(1) ), event_weight );
	zmuon_trackchi_plot->Fill(Zmumu_vec[0]->TrackChiSq(0), event_weight );
	zmuon_trackchi_plot->Fill(Zmumu_vec[0]->TrackChiSq(1), event_weight );
      }
    }
  }
    
  

  
  //////////////// Z->ee basic plots: //////////////////////////////////
  Zee_vec.clear();
  good_zee=false;
  if(event.get("Zee candidates", Zee_vec)) {
    
    NZee = Zee_vec.size();
    if(_debug_level>2)  cout <<"Found "<<NZee<<" Z->ee candidates"<<endl;    
    zee_num_plot->Fill(NZee, event_weight);   
    
    if(NZee>0) {
      if(_debug_level>2) {  
	cout <<"Z[0] fails "<<Zee_vec[0]->CutsFailed()<<" cuts"<<endl;
      }
      
      if(Zee_vec[0]->CutsFailed()==0){
	zee_mass_plot->Fill(Zee_vec[0]->InvMass(), event_weight);
	
	if(Zee_vec[0]->InvMass() > _zee_mass_low && Zee_vec[0]->InvMass() < _zee_mass_high) {
	  good_zee=true;
	  VB_rap = Zee_vec[0]->ZRap();
	  zee_rap_plot->Fill(VB_rap, event_weight);
	  VB_pt = Zee_vec[0]->ZPt();
	  zee_pt_plot->Fill(VB_pt, event_weight);
	  VB_phi = Zee_vec[0]->ZPhi();
	  zee_phi_plot->Fill(VB_phi, event_weight);
	  
	  if(Zee_vec[0]->LepPt(1) > Zee_vec[0]->LepPt(0)) cout <<"Wrong pT's!!"<<Zee_vec[0]->LepPt(1)<< " " << Zee_vec[0]->LepPt(0)<<endl;

	  lep_high_pt = Zee_vec[0]->LepPt(0);
	  lep_low_pt = Zee_vec[0]->LepPt(1);
	  lep_eta_d_1 = Zee_vec[0]->LepDetEta(0);
	  lep_eta_d_2 = Zee_vec[0]->LepDetEta(1);
	  lep_eta_p_1 = Zee_vec[0]->LepEta(0);
	  lep_eta_p_2 = Zee_vec[0]->LepEta(1);
	  
	  zelec_highpt_plot->Fill( lep_high_pt, event_weight );
	  zelec_lowpt_plot->Fill( lep_low_pt, event_weight );
	  zelec_eta_d_plot->Fill( lep_eta_d_1, event_weight );
	  zelec_eta_d_plot->Fill( lep_eta_d_2, event_weight );
	  zelec_eta_p_plot->Fill( lep_eta_p_1, event_weight );
	  zelec_eta_p_plot->Fill( lep_eta_p_2, event_weight );
	  float dphi = fabs(Zee_vec[0]->LepPhi(0) - Zee_vec[0]->LepPhi(1));
	  if(dphi>TMath::Pi()) dphi = 2*TMath::Pi() - dphi;
	  zelec_dphi_plot->Fill(dphi, event_weight);
	  
	  zelec_dz0_plot->Fill(Zee_vec[0]->LepZ0(0) - Zee_vec[0]->LepZ0(1), event_weight);
	  zelec_z0_plot->Fill(0.5*(Zee_vec[0]->LepZ0(0) + Zee_vec[0]->LepZ0(1) ), event_weight );
	  zelec_trackchisq_plot->Fill(Zee_vec[0]->track_match_chi2prob(0), event_weight );
	  zelec_trackchisq_plot->Fill(Zee_vec[0]->track_match_chi2prob(1), event_weight );
	}
      }
    } 
  }

 

  
  //////////////// W->munu basic plots: //////////////////////////////////
  jets = event.getCollection<TMBJet>(_jetBranch.c_str());
 
  Wmunu_vec.clear();
  good_wmunu=false;
  double myMETPhi;
  if(event.get("Wmunu candidates", Wmunu_vec)) {
   
    NWmu = Wmunu_vec.size();
    if(_debug_level>2)  cout <<"Found "<<NWmu<<" W candidates"<<endl;
    
    wmu_num_plot->Fill(NWmu, event_weight);

    if(_debug_level>2) {  
      for(int i=0; i<NWmu; i++) {
	cout <<"W "<<i<<" fails "<<Wmunu_vec[i]->CutsFailed()<<" cuts, mu pT: "<<Wmunu_vec[i]->LepPt()<<endl;
      }
    }
    if(NWmu>0 && Wmunu_vec[0]->CutsFailed()==0) {
      good_wmunu=true;

      double dphi_mujet, dphi_jetmet, dphi_mumet;
      dphi_jetmet = fabs(Wmunu_vec[0]->METPhi()-jets[0].Phi());
      if (dphi_jetmet > M_PI ) dphi_jetmet = 2*M_PI - dphi_jetmet;
      dphi_mumet = fabs(Wmunu_vec[0]->METPhi()-Wmunu_vec[0]->LepPhi());
      if (dphi_mumet > M_PI ) dphi_mumet = 2*M_PI - dphi_mumet;

      // start filling plots
      dphi_Wmunu_jetmet_plot->Fill(Wmunu_vec[0]->MET(),dphi_jetmet);
      dphi_mumet_plot->Fill(Wmunu_vec[0]->MET(),dphi_mumet);
      wmu_mass_plot->Fill(Wmunu_vec[0]->TrMass(), event_weight);
      VB_pt = Wmunu_vec[0]->WPt();
      wmu_pt_plot->Fill(VB_pt, event_weight);
      VB_phi = Wmunu_vec[0]->WPhi();
      wmu_phi_plot->Fill(VB_phi, event_weight);
      VB_rap = -999;//Wmunu_vec[0]->WEta();
      wmu_rap_plot->Fill(VB_rap, event_weight);
      wmuon_pt_plot->Fill(Wmunu_vec[0]->LepPt(), event_weight);
      wmuon_eta_plot->Fill(Wmunu_vec[0]->LepEta(), event_weight);
      wmuon_deteta_plot->Fill(Wmunu_vec[0]->LepDetEta(), event_weight);
      wmuon_phi_plot->Fill(Wmunu_vec[0]->LepPhi(), event_weight);
      wmuon_dca_plot->Fill(Wmunu_vec[0]->LepDCA(), event_weight);
      wmunu_met_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
      myMETPhi = atan2(Wmunu_vec[0]->METy(),Wmunu_vec[0]->METx());
      wmunu_metphi_plot->Fill(Wmunu_vec[0]->METPhi(), event_weight);
      wmuon_dphi_plot->Fill(Wmunu_vec[0]->DPhi(), event_weight);
      wmuon_z0_plot->Fill(Wmunu_vec[0]->LepZ0(), event_weight);
      wmuon_trackchi_plot->Fill(Wmunu_vec[0]->TrackChiSq(), event_weight);
      wmuon_ncft_plot->Fill(Wmunu_vec[0]->ncft(), event_weight);
      wmuon_nsmt_plot->Fill(Wmunu_vec[0]->nsmt(), event_weight);
      wmuon_triso_plot->Fill(Wmunu_vec[0]->TrackIso(), event_weight);
      wmuon_caliso_plot->Fill(Wmunu_vec[0]->CalIso(), event_weight);
      bool backtoback = false;
      for( int i = 0; i<jets.size(); i++){
        // compute delta-phi between muon and nearest jet
        dphi_mujet = fabs(Wmunu_vec[0]->LepPhi()-jets[i].Phi());
        if (dphi_mujet > M_PI ) dphi_mujet = 2*M_PI - dphi_mujet;
        if (dphi_mujet > M_PI - 0.2 ) backtoback = true;
      }
      if( backtoback){
        if( Wmunu_vec[0]->TrackIso()<=2.5 && Wmunu_vec[0]->CalIso()<=2.5 ){
          wmuon_qcdiso1_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
        }
        if ( jets.size()==2 ){
          wmunu_met2jet_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
          if( Wmunu_vec[0]->TrackIso()<=2.5 && Wmunu_vec[0]->CalIso()<=2.5 ){
            wmuon_qcdiso2_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
          }
        }
        else if( jets.size()==3 ){
          wmunu_met3jet_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
          if( Wmunu_vec[0]->TrackIso()<=2.5 && Wmunu_vec[0]->CalIso()<=2.5 ){
            wmuon_qcdiso3_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
          }
        }
        else if( jets.size()==4 ){
          wmunu_met4jet_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
          if( Wmunu_vec[0]->TrackIso()<=2.5 && Wmunu_vec[0]->CalIso()<=2.5 ){
            wmuon_qcdiso4_plot->Fill(Wmunu_vec[0]->MET(), event_weight);
          }
        }
      }
      // plot for distance between muon and leading/subleading jet
      float dR = 0;
      if( jets.size()>0 ){
        float dR1 = sqrt((jets[0].Eta()-Wmunu_vec[0]->LepEta())*(jets[0].Eta()-Wmunu_vec[0]->LepEta())
                       + (jets[0].Phi()-Wmunu_vec[0]->LepPhi())*(jets[0].Phi()-Wmunu_vec[0]->LepPhi()));
        wmuon_dRjet1_plot->Fill(dR1,event_weight);
        dR = dR1;
        if( jets.size()>1 ){
          float dR2 = sqrt((jets[1].Eta()-Wmunu_vec[0]->LepEta())*(jets[1].Eta()-Wmunu_vec[0]->LepEta())
                         + (jets[1].Phi()-Wmunu_vec[0]->LepPhi())*(jets[1].Phi()-Wmunu_vec[0]->LepPhi()));
          wmuon_dRjet2_plot->Fill(dR2,event_weight);
          if( dR2<dR1 ) dR = dR2; 
          if( jets.size()>2 ){
            float dR3 = sqrt((jets[2].Eta()-Wmunu_vec[0]->LepEta())*(jets[2].Eta()-Wmunu_vec[0]->LepEta())
                           + (jets[2].Phi()-Wmunu_vec[0]->LepPhi())*(jets[2].Phi()-Wmunu_vec[0]->LepPhi()));
            if( dR3<dR2 || dR3<dR1 ) dR = dR3;
          }
        }
        wmuon_dRjet_plot->Fill(dR,event_weight);
        // plot for distance between muon and nearest jet
        double deltaR;
        wmu_deltaR_nearestJet = 100.0;
        for (unsigned int i=0; i < jets.size(); i++) {
          deltaR = kinem::delta_R(jets[i].Eta(), jets[i].Phi(), Wmunu_vec[0]->LepEta(), Wmunu_vec[0]->LepPhi());
          if (deltaR < wmu_deltaR_nearestJet) wmu_deltaR_nearestJet = deltaR;
        } // Loop over jets
        wmuon_dRjetPI_plot->Fill(wmu_deltaR_nearestJet,event_weight);
      }
    }

    if(NWmu>1 && Wmunu_vec[1]->CutsFailed()==0) {
      wmu2_mass_plot->Fill(Wmunu_vec[1]->TrMass(), event_weight);
      VB2mu_pt = Wmunu_vec[1]->WPt();
      wmu2_pt_plot->Fill(VB2mu_pt, event_weight);
      VB2mu_phi = Wmunu_vec[1]->WPhi();
      wmu2_phi_plot->Fill(VB2mu_phi, event_weight);
      wmuon2_pt_plot->Fill(Wmunu_vec[1]->LepPt(), event_weight);
      wmuon2_eta_plot->Fill(Wmunu_vec[1]->LepEta(), event_weight);
      wmuon2_deteta_plot->Fill(Wmunu_vec[1]->LepDetEta(), event_weight);
      wmuon2_dca_plot->Fill(Wmunu_vec[1]->LepDCA(), event_weight);
      wmunu2_met_plot->Fill(Wmunu_vec[1]->MET(), event_weight);
      double myMETPhi2 = atan2(Wmunu_vec[1]->METx(),Wmunu_vec[1]->METy());
      wmunu2_metphi_plot->Fill(Wmunu_vec[1]->METPhi(), event_weight);
    }
    
  }
  
  // DPDP
  //////////////// W->enu basic plots: //////////////////////////////////
  jets = event.getCollection<TMBJet>(_jetBranch.c_str());

  //  cout << "DPDP:: " << _jetBranch.c_str() << endl;

  Wenu_vec.clear();
  good_wenu=false;
  if(event.get("Wenu candidates", Wenu_vec)) {
   
    NWe = Wenu_vec.size();
    if(_debug_level>2)  cout <<"Found "<<NWe<<" W candidates"<<endl;
    
    we_num_plot->Fill(NWe, event_weight);
    
    for(int i=0; i<NWe; i++) {
      if(_debug_level>2) {  
	cout <<"W "<<i<<" fails "<<Wenu_vec[i]->CutsFailed()<<" cuts, e pT: "<<Wenu_vec[i]->LepPt()<<endl;
      }
    }
    
    //    std::cout << "Jets size: " << jets.size() << std::endl;
    if(NWe>0 && Wenu_vec[0]->CutsFailed()==0) {
      good_wenu=true;

      if(_debug_level>2) cout <<"Good event is run/event: "<<event.getGlobal()->runno() << " " << event.getGlobal()->evtno() << endl;
      
      double dphi_ejet, dphi_jetmet, dphi_emet;
      dphi_jetmet = -999;
      if (jets.size()>0) {
	dphi_jetmet = fabs(Wenu_vec[0]->METPhi()-jets[0].Phi());
	if (dphi_jetmet > M_PI ) dphi_jetmet = 2*M_PI - dphi_jetmet;
      }
      dphi_emet = fabs(Wenu_vec[0]->METPhi()-Wenu_vec[0]->LepPhi());
      if (dphi_emet > M_PI ) dphi_emet = 2*M_PI - dphi_emet;

      dphi_Wenu_jetmet_plot->Fill(Wenu_vec[0]->MET(),dphi_jetmet);
      dphi_emet_plot->Fill(Wenu_vec[0]->MET(),dphi_emet);
      we_mass_plot->Fill(Wenu_vec[0]->TrMass(), event_weight);
      VB_pt = Wenu_vec[0]->WPt();
      we_pt_plot->Fill(VB_pt, event_weight);
      VB_phi = Wenu_vec[0]->WPhi();
      we_phi_plot->Fill(VB_phi, event_weight);
      
      welec_jetnum_plot->Fill(jets.size(), event_weight);
      if (jets.size()>0) {
	welec_jet1eta_plot->Fill(jets[0].Eta(), event_weight);
	if (jets.size()>1) {
	  welec_jet2eta_plot->Fill(jets[1].Eta(), event_weight);
	  if (jets.size()>2) {
	    welec_jet3eta_plot->Fill(jets[2].Eta(), event_weight);
	  }
	}	
      }

      VB_rap = -999;//Wenu_vec[0]->WEta();
      we_rap_plot->Fill(VB_rap, event_weight);
      welec_pt_plot->Fill(Wenu_vec[0]->LepPt(), event_weight);
      welec_eta_plot->Fill(Wenu_vec[0]->LepEta(), event_weight);
      welec_deteta_plot->Fill(Wenu_vec[0]->LepDetEta(), event_weight);
      welec_phi_plot->Fill(Wenu_vec[0]->LepPhi(), event_weight);
      welec_dca_plot->Fill(Wenu_vec[0]->LepDCA(), event_weight);
      wenu_met_plot->Fill(Wenu_vec[0]->MET(), event_weight);
      myMETPhi = atan2(Wenu_vec[0]->METy(),Wenu_vec[0]->METx());
      wenu_metphi_plot->Fill(Wenu_vec[0]->METPhi(), event_weight);
      welec_dphi_plot->Fill(Wenu_vec[0]->DPhi(), event_weight);
      welec_z0_plot->Fill(Wenu_vec[0]->LepZ0(), event_weight);
      welec_trackchi_plot->Fill(Wenu_vec[0]->TrackChiSq(), event_weight);
      welec_ncft_plot->Fill(Wenu_vec[0]->ncft(), event_weight);
      welec_nsmt_plot->Fill(Wenu_vec[0]->nsmt(), event_weight);
      welec_triso_plot->Fill(-999.0,event_weight);//Wenu_vec[0]->TrackIso(), event_weight);
      welec_caliso_plot->Fill(-999.0,event_weight);//Wenu_vec[0]->CalIso(), event_weight);
      bool backtoback = false;
      for( int i = 0; i<jets.size(); i++){
        // compute delta-phi between elec and nearest jet
        dphi_ejet = fabs(Wenu_vec[0]->LepPhi()-jets[i].Phi());
        if (dphi_ejet > M_PI ) dphi_ejet = 2*M_PI - dphi_ejet;
        if (dphi_ejet > M_PI - 0.2 ) backtoback = true;
      }
      bool electronIsIsolated =  true;//Wenu_vec[0]->TrackIso()<=2.5 && Wenu_vec[0]->CalIso()<=2.5;
      if( backtoback ){
        if( electronIsIsolated ){
          welec_qcdiso1_plot->Fill(Wenu_vec[0]->MET(), event_weight);
        }
        if ( jets.size()==2 ){
          wenu_met2jet_plot->Fill(Wenu_vec[0]->MET(), event_weight);
          if( electronIsIsolated ){
            welec_qcdiso2_plot->Fill(Wenu_vec[0]->MET(), event_weight);
          }
        }
        else if( jets.size()==3 ){
          wenu_met3jet_plot->Fill(Wenu_vec[0]->MET(), event_weight);
          if( electronIsIsolated ){
            welec_qcdiso3_plot->Fill(Wenu_vec[0]->MET(), event_weight);
          }
        }
        else if( jets.size()==4 ){
          wenu_met4jet_plot->Fill(Wenu_vec[0]->MET(), event_weight);
          if( electronIsIsolated ){
            welec_qcdiso4_plot->Fill(Wenu_vec[0]->MET(), event_weight);
          }
        }
      }
      // plot for distance between elec and leading/subleading jet
      float dR = 0;
      if( jets.size()>0 ){
        float dR1 = sqrt((jets[0].Eta()-Wenu_vec[0]->LepEta())*(jets[0].Eta()-Wenu_vec[0]->LepEta())
                       + (jets[0].Phi()-Wenu_vec[0]->LepPhi())*(jets[0].Phi()-Wenu_vec[0]->LepPhi()));
        welec_dRjet1_plot->Fill(dR1,event_weight);
        dR = dR1;
        if( jets.size()>1 ){
          float dR2 = sqrt((jets[1].Eta()-Wenu_vec[0]->LepEta())*(jets[1].Eta()-Wenu_vec[0]->LepEta())
                         + (jets[1].Phi()-Wenu_vec[0]->LepPhi())*(jets[1].Phi()-Wenu_vec[0]->LepPhi()));
          welec_dRjet2_plot->Fill(dR2,event_weight);
          if( dR2<dR1 ) dR = dR2; 
          if( jets.size()>2 ){
            float dR3 = sqrt((jets[2].Eta()-Wenu_vec[0]->LepEta())*(jets[2].Eta()-Wenu_vec[0]->LepEta())
                           + (jets[2].Phi()-Wenu_vec[0]->LepPhi())*(jets[2].Phi()-Wenu_vec[0]->LepPhi()));
            if( dR3<dR2 || dR3<dR1 ) dR = dR3;
          }
        }
        welec_dRjet_plot->Fill(dR,event_weight);
        // plot for distance between elec and nearest jet
        double deltaR;
        we_deltaR_nearestJet = 100.0;
        for (unsigned int i=0; i < jets.size(); i++) {
          deltaR = kinem::delta_R(jets[i].Eta(), jets[i].Phi(), Wenu_vec[0]->LepEta(), Wenu_vec[0]->LepPhi());
          if (deltaR < we_deltaR_nearestJet) we_deltaR_nearestJet = deltaR;
        } // Loop over jets
        welec_dRjetPI_plot->Fill(we_deltaR_nearestJet,event_weight);
      }
    }

    if(NWe>1 && Wenu_vec[1]->CutsFailed()==0) {
      we2_mass_plot->Fill(Wenu_vec[1]->TrMass(), event_weight);
      VB2e_pt = Wenu_vec[1]->WPt();
      we2_pt_plot->Fill(VB2e_pt, event_weight);
      VB2e_phi = Wenu_vec[1]->WPhi();
      we2_phi_plot->Fill(VB2e_phi, event_weight);
      welec2_pt_plot->Fill(Wenu_vec[1]->LepPt(), event_weight);
      welec2_eta_plot->Fill(Wenu_vec[1]->LepEta(), event_weight);
      welec2_deteta_plot->Fill(Wenu_vec[1]->LepDetEta(), event_weight);
      welec2_dca_plot->Fill(Wenu_vec[1]->LepDCA(), event_weight);
      wenu2_met_plot->Fill(Wenu_vec[1]->MET(), event_weight);
      double myMETPhi2 = atan2(Wenu_vec[1]->METx(),Wenu_vec[1]->METy());
      wenu2_metphi_plot->Fill(Wenu_vec[1]->METPhi(), event_weight);
    }
 
  }
  
  ////////////////////////////////////////////////////////////////////////


  //if we dont have any candidates, return false (so this event is not skimmed)
  //if(!NWmu && !NZmu && !NZee && !NWe) return false;
  //if(!NWmu) return false;
  //if we have candidates, but not in the desired mass window, return true, but don't look at the jets:
  if(!good_wmunu && !good_wenu && !good_zmumu && !good_zee ) return true;
  

  
  //////////////// Jets basic plots: ////////////////////////////////////
  Njet = jets.size();
  if(_debug_level>2){
    cout <<"Found "<<Njet<<" jets"<<endl;
    if(_debug_level>3) {
      for(int i=0; i<Njet; i++) {
	cout <<jets[i].Pt()<<endl;
      }
    }
  }
  jet_num_plot->Fill(Njet, event_weight);
    
  if(Njet) {
    jet_leadingpt_plot->Fill(jets[0].Pt(), event_weight);
    jet_leadingphi_plot->Fill(jets[0].Phi(), event_weight);
    jet_leadingrap_plot->Fill(jets[0].Eta(), event_weight);
    if(NZmu) {
      zmu_rap_1j_plot->Fill(VB_rap, event_weight);
      zmu_pt_1j_plot->Fill(VB_pt, event_weight);
      zmuon_highpt_1j_plot->Fill(lep_high_pt, event_weight);
      zmuon_lowpt_1j_plot->Fill(lep_low_pt, event_weight);
      zmuon_eta_p_1j_plot->Fill(lep_eta_p_1, event_weight);
      zmuon_eta_p_1j_plot->Fill(lep_eta_p_2, event_weight);
      zmuon_eta_d_1j_plot->Fill(lep_eta_d_1, event_weight);
      zmuon_eta_d_1j_plot->Fill(lep_eta_d_2, event_weight);
    }
    if(NZee) {
      zee_rap_1j_plot->Fill(VB_rap, event_weight);
      zee_pt_1j_plot->Fill(VB_pt, event_weight);
      zelec_highpt_1j_plot->Fill(lep_high_pt, event_weight);
      zelec_lowpt_1j_plot->Fill(lep_low_pt, event_weight);
      zelec_eta_p_1j_plot->Fill(lep_eta_p_1, event_weight);
      zelec_eta_p_1j_plot->Fill(lep_eta_p_2, event_weight);
      zelec_eta_d_1j_plot->Fill(lep_eta_d_1, event_weight);
      zelec_eta_d_1j_plot->Fill(lep_eta_d_2, event_weight);
    }
  }
  double dphi;
  if(Njet>1){
    jet_leadingpt_2j_plot->Fill(jets[0].Pt(), event_weight);
    jet_leadingrap_2j_plot->Fill(jets[0].Eta(), event_weight);
    jet_firstsecondrap_plot->Fill(jets[0].Eta(), event_weight);
    jet_leadingphi_2j_plot->Fill(jets[0].Phi(), event_weight);
    jet_firstsecondphi_plot->Fill(jets[0].Phi(), event_weight);
    jet_secondpt_plot->Fill(jets[1].Pt(), event_weight);
    jet_secondrap_plot->Fill(jets[1].Eta(), event_weight);
    jet_firstsecondrap_plot->Fill(jets[1].Eta(), event_weight);
    jet_secondphi_plot->Fill(jets[1].Phi(), event_weight);
    jet_firstsecondphi_plot->Fill(jets[1].Phi(), event_weight);
    jet_deta_plot->Fill(fabs(jets[0].Eta()-jets[1].Eta()), event_weight);
    dphi = fabs(jets[0].Phi()-jets[1].Phi());
    if (dphi > TMath::Pi()) dphi = 2*TMath::Pi() - dphi;
    jet_dphi_plot->Fill(dphi, event_weight);
    invmass = sqrt( (jets[0].E()+jets[1].E())*(jets[0].E()+jets[1].E()) -
                    (jets[0].Px()+jets[1].Px())*(jets[0].Px()+jets[1].Px()) - 
                    (jets[0].Py()+jets[1].Py())*(jets[0].Py()+jets[1].Py()) - 
                    (jets[0].Pz()+jets[1].Pz())*(jets[0].Pz()+jets[1].Pz()) );
    jet_invmass_plot->Fill(invmass, event_weight);
    if(NZmu) {
      zmu_rap_2j_plot->Fill(VB_rap, event_weight);
      zmu_pt_2j_plot->Fill(VB_pt, event_weight);
      zmuon_highpt_2j_plot->Fill(lep_high_pt, event_weight);
      zmuon_lowpt_2j_plot->Fill(lep_low_pt, event_weight);
      zmuon_eta_p_2j_plot->Fill(lep_eta_p_1, event_weight);
      zmuon_eta_p_2j_plot->Fill(lep_eta_p_2, event_weight);
      zmuon_eta_d_2j_plot->Fill(lep_eta_d_1, event_weight);
      zmuon_eta_d_2j_plot->Fill(lep_eta_d_2, event_weight);
    }
    if(NZee) {
      zee_rap_2j_plot->Fill(VB_rap, event_weight);
      zee_pt_2j_plot->Fill(VB_pt, event_weight);
      zelec_highpt_2j_plot->Fill(lep_high_pt, event_weight);
      zelec_lowpt_2j_plot->Fill(lep_low_pt, event_weight);
      zelec_eta_p_2j_plot->Fill(lep_eta_p_1, event_weight);
      zelec_eta_p_2j_plot->Fill(lep_eta_p_2, event_weight);
      zelec_eta_d_2j_plot->Fill(lep_eta_d_1, event_weight);
      zelec_eta_d_2j_plot->Fill(lep_eta_d_2, event_weight);
    }
  }
  if(Njet>2) {
    jet_thirdpt_plot->Fill(jets[2].Pt(), event_weight);
    jet_thirdrap_plot->Fill(jets[2].Eta(), event_weight);
    jet_thirdphi_plot->Fill(jets[2].Phi(), event_weight);
  }
  
  for(int i=0; i<Njet+1; i++)
    jet_num_inclusive_plot->Fill(i, event_weight);
  
 

  //////////////// Jets Probability plots: ////////////////////////////////////

  int jet_there=0;
  float jet_eta[Njet];
  for(int i=0; i<Njet; i++) { jet_eta[i] = 999.; jeteta_ind[i] = -99; }
  for(int i=0; i<Njet; i++) {
    
    jet_allpt_plot->Fill(jets[i].Pt(), event_weight);
    jet_allrap_plot->Fill(jets[i].Eta(), event_weight);
    jet_allphi_plot->Fill(jets[i].Phi(), event_weight);
    if(jets[i].Pt() < _jet_pt_cut) continue; 
    
    float  dphij = jets[i].Phi() -  VB_phi;
    if(dphij>TMath::Pi()) dphij-=2*TMath::Pi();
    if(dphij<-TMath::Pi()) dphij+=2*TMath::Pi();
    dphij = fabs(dphij);
    jet_dphiW_plot->Fill(dphij, event_weight);
    if(dphij > _oppjet_cut) jet_there = 1;

    //fill jet array for sorting
    jet_eta[i] = jets[i].Eta();
    
  }

  //sort jets in eta
  sort(jet_eta,jet_eta+Njet);
  for(int i=0; i<Njet; i++){
    for(int j=0; j<Njet; j++){
      if(jet_eta[i]==(float)jets[j].Eta()){jeteta_ind[i]=j;}
    }
  }
  
  if(jet_there) jet_yes_plot->Fill(VB_pt, event_weight);
  else jet_no_plot->Fill(VB_pt, event_weight);
 
  if (NWmu>0) {
  //////////////// Ntuple for NN /////////////////
    wmu_mupt = Wmunu_vec[0]->LepPt();
    wmu_mueta = Wmunu_vec[0]->LepEta();
    wmu_muphi = Wmunu_vec[0]->LepPhi();
    wmu_met = Wmunu_vec[0]->MET();
    wmu_metphi = Wmunu_vec[0]->METPhi();
    wmu_wpt = Wmunu_vec[0]->WPt();
    wmu_wtrmass = Wmunu_vec[0]->TrMass();
    wmu_wdphi = Wmunu_vec[0]->DPhi();
    wmu_lepz0 = Wmunu_vec[0]->LepZ0();
    wmu_caliso = Wmunu_vec[0]->CalIso();
    wmu_trackiso = Wmunu_vec[0]->TrackIso();
    for(int i=0; i<Njet; i++){
      wmu_jetpt[i] = jets[i].Pt();
      wmu_jeteta[i] = jets[i].Eta();
      wmu_jetphi[i] = jets[i].Phi();
      wmu_jetpx[i] = jets[i].Px();
      wmu_jetpy[i] = jets[i].Py();
      wmu_jetpz[i] = jets[i].Pz();
      wmu_jetE[i] = jets[i].E();
    }
  }
  if (NWe>0) {
  //////////////// Ntuple for NN /////////////////
    we_ept = Wenu_vec[0]->LepPt();
    we_eeta = Wenu_vec[0]->LepEta();
    we_ephi = Wenu_vec[0]->LepPhi();
    we_met = Wenu_vec[0]->MET();
    we_metphi = Wenu_vec[0]->METPhi();
    we_wpt = Wenu_vec[0]->WPt();
    we_wtrmass = Wenu_vec[0]->TrMass();
    we_wdphi = Wenu_vec[0]->DPhi();
    we_lepz0 = Wenu_vec[0]->LepZ0();
    we_caliso = -999;//Wenu_vec[0]->CalIso();
    we_trackiso = -999;//Wenu_vec[0]->TrackIso();
    for(int i=0; i<Njet; i++){
      we_jetpt[i] = jets[i].Pt();
      we_jeteta[i] = jets[i].Eta();
      we_jetphi[i] = jets[i].Phi();
      we_jetpx[i] = jets[i].Px();
      we_jetpy[i] = jets[i].Py();
      we_jetpz[i] = jets[i].Pz();
      we_jetE[i] = jets[i].E();
    }
  }

  // trace the JES corrections
  if(NWmu>0 && Njet>1){
    jets[0].ActAsUnCorrected(); jet_leadingpt_uncorr->Fill(jets[0].Pt());
    jets[0].ActAsJESCorrected(); jet_leadingpt_jes->Fill(jets[0].Pt());
    jets[0].ActAsJESMUCorrected(); jet_leadingpt_jesmu->Fill(jets[0].Pt());
    jets[0].ActAsSmearedMU(); jet_leadingpt_jesmusmear->Fill(jets[0].Pt());
    jets[1].ActAsUnCorrected(); jet_secondpt_uncorr->Fill(jets[1].Pt());
    jets[1].ActAsJESCorrected(); jet_secondpt_jes->Fill(jets[1].Pt());
    jets[1].ActAsJESMUCorrected(); jet_secondpt_jesmu->Fill(jets[1].Pt());
    jets[1].ActAsSmearedMU(); jet_secondpt_jesmusmear->Fill(jets[1].Pt());
  }

  //----------------------------------------
  // topological variables
  //----------------------------------------
  if (NWe>0) {
    topology topo(_usealljets);
    
    int _njt           = jets.size();
    double _met        = Wenu_vec[0]->MET();
    double _metx       = Wenu_vec[0]->METx();
    double _mety       = Wenu_vec[0]->METy();
    double _jetpt[Njet];
    double _jeteta[Njet];
    double _jetphi[Njet];
    double _jetE[Njet];
    for(int i=0; i<Njet; i++){
      _jetpt[i] = jets[i].Pt();
      _jeteta[i] = jets[i].Eta();
      _jetphi[i] = jets[i].Phi();
      _jetE[i] = jets[i].E();
    }

    topo.eval(_njt,
	      _jetpt,_jeteta,_jetphi,_jetE,
	      we_ept,we_eeta,we_ephi,Wenu_vec[0]->LepE(),
	      _met,_metx,_mety );

    wmt      =topo.wmt();
    wet      =topo.wet();
    weta     =topo.weta();
    ht25     =topo.ht25();
    ht20     =topo.ht20();
    ht       =topo.ht();
    apla     =topo.aplanarity();
    cent     =topo.centrality();
    plan     =topo.planarity();
    spher    =topo.spheresity();
    ktminp   =topo.ktminp();
    ktminpme =topo.ktminpme();
    
    jetmt    =topo.jetmt();
    mevent   =topo.mevent();
    jetamax  =topo.etamax();
    jetasqsum=topo.etasqsum();
    drmin    =topo.drmin();
    dijetmass  = topo.dijetmass();
    dijetptrel = topo.dijetptrel();
    ht2      =topo.ht2();
    ht3      =topo.ht3();
    htall    =topo.htall();
    htl      =topo.htl();
    htp      =topo.htp();
    ht2p     =topo.ht2p();
    ht3p     =topo.ht3p();
    htpp     =topo.htpp();
    ht2pp    =topo.ht2pp();
    ht3pp    =topo.ht3pp();
    wpt      =topo.wpT();

    /*
    cout << "wmt " << wmt << endl;
    cout << "wet " << wet << endl;
    cout << "weta " << weta << endl;
    cout << "ht25 " << ht25 << endl;
    cout << "ht20 " << ht20 << endl;
    cout << "ht " << ht << endl;
    cout << "apla " << apla << endl;
    cout << "cent " << cent << endl;
    cout << "plan " << plan << endl;
    cout << "spher " << spher << endl;
    cout << "ktminp " << ktminp << endl;
    cout << "ktminpme " << ktminpme << endl;
    cout << "jetmt " << jetmt << endl;
    cout << "mevent " << mevent << endl;
    cout << "jetamax " << jetamax << endl;
    cout << "jetasqsum " << jetasqsum << endl;
    cout << "drmin " << drmin << endl;
    */

    // end of topological variables
  }

  if (NWmu>0 || NWe>0) { 
    nntree->Fill();
  }


 
  /////////////////////// Finish this event: ////////////////////////////
  Zmumu_vec.clear();
  Zee_vec.clear();
  Wmunu_vec.clear();
  Wenu_vec.clear();
  
  // ntuple vars
  invmass=-9999;
  for(int i=0; i<Njet; i++){
    jeteta_ind[i]=-9999;
  }
  wmu_mupt=-9999;
  wmu_mueta=-9999;
  wmu_muphi=-9999;
  wmu_met=-9999;
  wmu_metphi=-9999;
  wmu_wpt=-9999;
  wmu_wtrmass=-9999;
  wmu_wdphi=-9999;
  wmu_lepz0=-9999;
  wmu_trackiso=-9999;
  wmu_caliso=-9999;
  event_weight=-9999;
  if (NWmu>0) {
    for(int i=0; i<Njet; i++){
      wmu_jetpt[i]=-9999;
      wmu_jeteta[i]=-9999;
      wmu_jetphi[i]=-9999;
      wmu_jetpx[i]=-9999;
      wmu_jetpy[i]=-9999;
      wmu_jetpz[i]=-9999;
      wmu_jetE[i]=-9999;
    }
  }
  wmu_deltaR_nearestJet=-9999;
  we_ept=-9999;
  we_eeta=-9999;
  we_ephi=-9999;
  we_met=-9999;
  we_metphi=-9999;
  we_wpt=-9999;
  we_wtrmass=-9999;
  we_wdphi=-9999;
  we_lepz0=-9999;
  we_trackiso=-9999;
  we_caliso=-9999;
  if (NWe>0) {
    for(int i=0; i<Njet; i++){
      we_jetpt[i]=-9999;
      we_jeteta[i]=-9999;
      we_jetphi[i]=-9999;
      we_jetpx[i]=-9999;
      we_jetpy[i]=-9999;
      we_jetpz[i]=-9999;
      we_jetE[i]=-9999;
    }
  }
  we_deltaR_nearestJet=-9999;
  wmt=-9999;
  wet=-9999;
  weta=-9999;
  ht25=-9999;
  ht20=-9999;
  ht=-9999;
  apla=-9999;
  cent=-9999;
  plan=-9999;
  spher=-9999;
  ht2p=-9999;
  ktminp=-9999;
  ktminpme=-9999;
  jetmt=-9999;
  mevent=-9999;
  jetamax=-9999;
  jetasqsum=-9999;
  drmin=-9999;
  ht2=-9999;
  ht3=-9999;
  htall=-9999;
  htl=-9999;
  dijetmass=-9999;
  dijetptrel=-9999;
  htp=-9999;
  //     ht2p=-9999;
  ht3p=-9999;
  htpp=-9999;
  ht2pp=-9999;
  ht3pp=-9999;
  wpt=-9999;
  ////
  Njet=-9999; // must be after njet loops!! O_o
  
  
  return true;
  
  }

  void PlotVBF::begin()
  {
  zee_num_plot = new TH1F("zee_num","Number of Z->ee candidates", 10,0,10);
  zmu_num_plot = new TH1F("zmumu_num","Number of Z->mumu candidates", 10,0,10);
  wmu_num_plot = new TH1F("wmu_num","Number of W->munu candidates", 10,0,10);
  we_num_plot = new TH1F("we_num","Number of W->enu candidates", 10,0,10);

  zee_num_plot->Sumw2();
  zmu_num_plot->Sumw2();
  wmu_num_plot->Sumw2();
  we_num_plot->Sumw2();

  zee_mass_plot = new TH1F("zee_mass","Di-EM Mass", 300,0,300);
  zmu_mass_plot = new TH1F("zmumu_mass","Di-Muon Mass", 300,0,300);
  wmu_mass_plot = new TH1F("wmu_trmass","W->munu Transverse Mass", 120,0,300);
  wmu2_mass_plot = new TH1F("wmu2_trmass","W2->munu Transverse Mass", 120,0,300);
  we_mass_plot = new TH1F("we_trmass","W->enu Transverse Mass", 120,0,300);
  we2_mass_plot = new TH1F("we2_trmass","W2->enu Transverse Mass", 120,0,300);

  zee_mass_plot->Sumw2();
  zmu_mass_plot->Sumw2();
  wmu_mass_plot->Sumw2();
  wmu2_mass_plot->Sumw2();
  we_mass_plot->Sumw2();
  we2_mass_plot->Sumw2();

  zee_rap_plot = new TH1F("zee_rap","Di-EM Rapidity", 60,-3,3);
  zmu_rap_plot = new TH1F("zmumu_rap","Di-Muon Rapidity", 60,-3,3);
  zee_rap_1j_plot = new TH1F("zee_rap_1j","Di-EM Rapidity, >= 1 jet", 30,-3,3);
  zmu_rap_1j_plot = new TH1F("zmumu_rap_1j","Di-Muon Rapidity, >= 1 jet", 30,-3,3);
  zee_rap_2j_plot = new TH1F("zee_rap_2j","Di-EM Rapidity, >= 2 jets", 20,-3,3);
  zmu_rap_2j_plot = new TH1F("zmumu_rap_2j","Di-Muon Rapidity, >= 2 jets", 20,-3,3);

  zee_rap_plot->Sumw2();
  zmu_rap_plot->Sumw2();
  zee_rap_1j_plot->Sumw2();
  zmu_rap_1j_plot->Sumw2();
  zee_rap_2j_plot->Sumw2();
  zmu_rap_2j_plot->Sumw2();

  zee_pt_plot = new TH1F("zee_pt","Di-EM pT", 100,0,200);
  zmu_pt_plot = new TH1F("zmumu_pt","Di-Muon pT", 100,0,200);
  zee_pt_1j_plot = new TH1F("zee_pt_1j","Di-EM pT, >= 1 jet", 50,0,200);
  zmu_pt_1j_plot = new TH1F("zmumu_pt_1j","Di-Muon pT, >= 1 jet", 50,0,200);
  zee_pt_2j_plot = new TH1F("zee_pt_2j","Di-EM pT, >= 2 jet", 25,0,200);
  zmu_pt_2j_plot = new TH1F("zmumu_pt_2j","Di-Muon pT, >= 2 jet", 25,0,200);

  zee_pt_plot->Sumw2();
  zmu_pt_plot->Sumw2();
  zee_pt_1j_plot->Sumw2();
  zmu_pt_1j_plot->Sumw2();
  zee_pt_2j_plot->Sumw2();
  zmu_pt_2j_plot->Sumw2();

  zmuon_highpt_plot = new TH1F("muon_high_pt_Z","Leading Muon pT, Z Events", 100,0,100);
  zmuon_lowpt_plot = new TH1F("muon_low_pt_Z","Second Muon pT, Z Events", 100,0,100);
  zmuon_highpt_1j_plot = new TH1F("muon_high_pt_1j_Z","Leading Muon pT, Z Events, >= 1 jet", 50,0,100);
  zmuon_lowpt_1j_plot = new TH1F("muon_low_pt_1j_Z","Second Muon pT, Z Events, >= 1 jet", 50,0,100);
  zmuon_highpt_2j_plot = new TH1F("muon_high_pt_2j_Z","Leading Muon pT, Z Events, >= 2 jet", 20,0,100);
  zmuon_lowpt_2j_plot = new TH1F("muon_low_pt_2j_Z","Second Muon pT, Z Events, >= 2 jet", 20,0,100);

  zmuon_highpt_plot->Sumw2();
  zmuon_lowpt_plot->Sumw2();
  zmuon_highpt_1j_plot->Sumw2();
  zmuon_lowpt_1j_plot->Sumw2();
  zmuon_highpt_2j_plot->Sumw2();
  zmuon_lowpt_2j_plot->Sumw2();

  zelec_highpt_plot = new TH1F("elec_high_pt_Z","Leading Electron pT, Z Events", 100,0,100);
  zelec_lowpt_plot = new TH1F("elec_low_pt_Z","Second Electron pT, Z Events", 100,0,100);
  zelec_highpt_1j_plot = new TH1F("elec_high_pt_1j_Z","Leading Electron pT, Z Events, >= 1 jet", 100,0,100);
  zelec_lowpt_1j_plot = new TH1F("elec_low_pt_1j_Z","Second Electron pT, Z Events, >= 1 jet", 100,0,100);
  zelec_highpt_2j_plot = new TH1F("elec_high_2j_pt_Z","Leading Electron pT, Z Events, >= 2 jet", 100,0,100);
  zelec_lowpt_2j_plot = new TH1F("elec_low_pt_2j_Z","Second Electron pT, Z Events, >= 2 jet", 100,0,100);

  zelec_highpt_plot->Sumw2();
  zelec_lowpt_plot->Sumw2();
  zelec_highpt_1j_plot->Sumw2();
  zelec_lowpt_1j_plot->Sumw2();
  zelec_highpt_2j_plot->Sumw2();
  zelec_lowpt_2j_plot->Sumw2();

  zmuon_eta_p_plot = new TH1F("muon_eta_p_Z","Muon (physics) #eta, Z Events", 100,-3,3);
  zmuon_eta_d_plot = new TH1F("muon_eta_d_Z","Muon (detector) #eta, Z Events", 100,-3,3);
  zmuon_eta_p_1j_plot = new TH1F("muon_eta_p_1j_Z","Muon (physics) #eta, Z Events, >= 1 jet", 50,-3,3);
  zmuon_eta_d_1j_plot = new TH1F("muon_eta_d_1j_Z","Muon (detector) #eta, Z Events, >= 1 jet", 50,-3,3);
  zmuon_eta_p_2j_plot = new TH1F("muon_eta_p_2j_Z","Muon (physics) #eta, Z Events, >= 2 jets", 20,-3,3);
  zmuon_eta_d_2j_plot = new TH1F("muon_eta_d_2j_Z","Muon (detector) #eta, Z Events, >= 2 jets", 20,-3,3);

  zmuon_eta_p_plot->Sumw2();
  zmuon_eta_d_plot->Sumw2();
  zmuon_eta_p_1j_plot->Sumw2();
  zmuon_eta_d_1j_plot->Sumw2();
  zmuon_eta_p_2j_plot->Sumw2();
  zmuon_eta_d_2j_plot->Sumw2();

  zelec_eta_p_plot = new TH1F("elec_eta_p_Z","Electron (physics) #eta, Z Events", 100,-3,3);
  zelec_eta_d_plot = new TH1F("elec_eta_d_Z","Electron (detector) #eta, Z Events", 100,-3,3);
  zelec_eta_p_1j_plot = new TH1F("elec_eta_p_1j_Z","Electron (physics) #eta, Z Events, >= 1 jet", 50,-3,3);
  zelec_eta_d_1j_plot = new TH1F("elec_eta_d_1j_Z","Electron (detector) #eta, Z Events, >= 1 jet", 50,-3,3);
  zelec_eta_p_2j_plot = new TH1F("elec_eta_p_2j_Z","Electron (physics) #eta, Z Events, >= 2 jet", 20,-3,3);
  zelec_eta_d_2j_plot = new TH1F("elec_eta_d_2j_Z","Electron (detector) #eta, Z Events, >= 2 jet", 20,-3,3);

  zelec_eta_p_plot->Sumw2();
  zelec_eta_d_plot->Sumw2();
  zelec_eta_p_1j_plot->Sumw2();
  zelec_eta_d_1j_plot->Sumw2();
  zelec_eta_p_2j_plot->Sumw2();
  zelec_eta_d_2j_plot->Sumw2();

  zee_phi_plot = new TH1F("zee_phi","Di-EM #phi", 100,0,6.4);
  zmu_phi_plot = new TH1F("zmumu_phi","Di-Muon #phi", 100,0,6.4);
  wmu_phi_plot = new TH1F("wmu_phi","W #phi (W->munu)", 100,0,6.4);
  wmu_rap_plot = new TH1F("wmu_eta","W #rap (W->munu)", 60,-3,3);
  wmu2_phi_plot = new TH1F("wmu2_phi","W #phi (W->munu)", 100,0,6.4);
  we_phi_plot = new TH1F("we_phi","W #phi (W->enu)", 100,0,6.4);
  we_rap_plot = new TH1F("we_eta","W #rap (W->enu)", 60,-3,3);
  we2_phi_plot = new TH1F("we2_phi","W #phi (W->enu)", 100,0,6.4);
  
  welec_jetnum_plot = new TH1F("welec_jetnum_plot","Jet multiplicity", 15, 0,15);
  welec_jet1eta_plot = new TH1F("welec_jet1eta_plot","Rapidity of leading jet", 25, -5,5);
  welec_jet2eta_plot = new TH1F("welec_jet2eta_plot","Rapidity of second jet", 25, -5,5);
  welec_jet3eta_plot = new TH1F("welec_jet3eta_plot","Rapidity of third jet", 25, -5,5);

  zee_phi_plot->Sumw2();
  zmu_phi_plot->Sumw2();
  wmu_phi_plot->Sumw2();
  wmu_rap_plot->Sumw2();
  wmu2_phi_plot->Sumw2();
  we_phi_plot->Sumw2();
  we2_phi_plot->Sumw2();

  welec_jetnum_plot->Sumw2();
  welec_jet1eta_plot->Sumw2();
  welec_jet2eta_plot->Sumw2();
  welec_jet3eta_plot->Sumw2();
  

  zmuon_dphi_plot = new TH1F("muon_dphi_Z","Muon #Delta #phi, Z Events", 100,0,3.2);
  zmuon_z0_plot = new TH1F("muon_z0_Z","Muon z0, Z Events", 100,-100,100);
  zmuon_pv_plot = new TH1F("muon_pv_Z","Primary Vertex z0, Z->#mu#mu Events", 100,-100,100);
  zmuon_dz0_plot = new TH1F("muon_dz0_Z","Muon #Delta z0, Z Events", 100,-10,10);
  zmuon_trackchi_plot = new TH1F("muon_track_chisq","Central track chi^2, Z Events", 100,0,10);

  zmuon_dphi_plot->Sumw2();
  zmuon_z0_plot->Sumw2();
  zmuon_pv_plot->Sumw2();
  zmuon_dz0_plot->Sumw2();
  zmuon_trackchi_plot->Sumw2();

  zelec_dphi_plot = new TH1F("elec_dphi_Z","Electron #Delta #phi, Z Events", 100,0,3.2);
  zelec_z0_plot = new TH1F("elec_z0_Z","Electron z0, Z Events", 100,-100,100);
  zelec_pv_plot = new TH1F("elec_pv_Z","Primary Vertex z0, Z->ee Events", 100,-100,100);
  zelec_dz0_plot = new TH1F("elec_dz0_Z","Electron #Delta z0, Z Events", 100,-10,10);
  zelec_trackchisq_plot = new TH1F("em_track_chisq","EM track match chisq, Z Events", 100,0,1);

  zelec_dphi_plot->Sumw2();
  zelec_z0_plot->Sumw2();
  zelec_pv_plot->Sumw2();
  zelec_dz0_plot->Sumw2();
  zelec_trackchisq_plot->Sumw2();

  wmuon_pt_plot = new TH1F("muon_pt_W","Leading Muon PT, W Events", 100,0,200);
  wmuon_eta_plot = new TH1F("muon_eta_W","Muon #eta, W Events", 60,-3,3);
  wmuon_deteta_plot = new TH1F("muon_deteta_W","Muon #eta, W Events", 60,-3,3);
  wmuon2_pt_plot = new TH1F("muon_pt_W2","Second Muon PT, W Events", 100,0,200);
  wmuon2_eta_plot = new TH1F("muon_eta_W2","Second Muon eta, W Events", 60,-3,3);
  wmuon2_deteta_plot = new TH1F("muon_deteta_W2","Second Muon eta, W Events", 60,-3,3);
  wmuon_dphi_plot = new TH1F("muon_dphi_W","Muon #Delta #phi (mu,met #phi), W Events", 64,0,3.15);
  wmuon_z0_plot = new TH1F("muon_z0_W","Muon z0, W Events", 100,-100,100);
  wmuon_trackchi_plot = new TH1F("muon_track_chisq_W","Central track chi^2, W Events", 20,0,5);
  wmuon_dca_plot = new TH1F("muon_dca_W","Muon DCA, W Events",25,0,0.5);
  wmuon2_dca_plot = new TH1F("muon_dca_W2","Second Muon DCA, W Events",25,0,0.5);
  wmuon_phi_plot = new TH1F("muon_phi_W","Muon #phi, W Events",32,0,6.3);
  wmuon_ncft_plot = new TH1F("muon_ncft_W","Muon cft tracks, W Events",20,0,20);
  wmuon_nsmt_plot = new TH1F("muon_nsmt_W","Muon smt tracks, W Events",20,0,20);
  wmuon_triso_plot = new TH1F("muon_triso_W","Muon track isolation, W Events",20,0,10);
  wmuon_caliso_plot = new TH1F("muon_caliso_W","Muon cal isolation, W Events",20,0,10);
  wmuon_qcdiso1_plot = new TH1F("muon_qcdiso1_W","Tight/Loose Muon vs MET, W Events",50,0,100);
  wmuon_qcdiso2_plot = new TH1F("muon_qcdiso2_W","Tight/Loose Muon vs MET, W Events",50,0,100);
  wmuon_qcdiso3_plot = new TH1F("muon_qcdiso3_W","Tight/Loose Muon vs MET, W Events",50,0,100);
  wmuon_qcdiso4_plot = new TH1F("muon_qcdiso4_W","Tight/Loose Muon vs MET, W Events",50,0,100);
  wmuon_dRjet_plot = new TH1F("muon_dRjet_W","#Delta R muon-jet",64,0,6.3);
  wmuon_dRjetPI_plot = new TH1F("muon_dRjetPI_W","#Delta R muon-jet",64,0,6.3);
  wmuon_dRjet1_plot = new TH1F("muon_dRjet1_W","#Delta R muon-jet1",64,0,6.3);
  wmuon_dRjet2_plot = new TH1F("muon_dRjet2_W","#Delta R muon-jet2",64,0,6.3);
  wmunu_met_plot = new TH1F("munu_met_W","MET, W Events",50,0,100);
  wmunu_met2jet_plot = new TH1F("munu_met2jet_W","MET, 2-jet bin, W Events",50,0,100);
  wmunu_met3jet_plot = new TH1F("munu_met3jet_W","MET, 3-jet bin, W Events",50,0,100);
  wmunu_met4jet_plot = new TH1F("munu_met4jet_W","MET, 4-jet bin, W Events",50,0,100);
  wmunu_metphi_plot = new TH1F("munu_metphi_W","MET #phi, W Events",64,0,6.3);
  wmunu2_met_plot = new TH1F("munu_met_W2","Second MET, W Events",50,0,100);
  wmunu2_metphi_plot = new TH1F("munu_metphi_W2","Second MET #phi, W Events",64,0,6.3);
 
  dphi_mumet_plot = new TH2F("dphi_mumet_plot","Triangle #Delta #phi(#mu-MET) vs. MET",100,0,200,32,0,3.2);
  dphi_Wmunu_jetmet_plot = new TH2F("dphi_Wmunu_jetmet_plot","Triangle #Delta #phi(jet-MET) vs. MET",100,0,200,32,0,3.2);


  welec_pt_plot = new TH1F("elec_pt_W","Leading Electron PT, W Events", 100,0,200);
  welec_eta_plot = new TH1F("elec_eta_W","Electron #eta, W Events", 60,-3,3);
  welec_deteta_plot = new TH1F("elec_deteta_W","Electron #eta, W Events", 60,-3,3);
  welec2_pt_plot = new TH1F("elec_pt_W2","Second Electron PT, W Events", 100,0,200);
  welec2_eta_plot = new TH1F("elec_eta_W2","Second Electron eta, W Events", 60,-3,3);
  welec2_deteta_plot = new TH1F("elec_deteta_W2","Second Electron eta, W Events", 60,-3,3);
  welec_dphi_plot = new TH1F("elec_dphi_W","Electron #Delta #phi (e,met #phi), W Events", 64,0,3.15);
  welec_z0_plot = new TH1F("elec_z0_W","Electron z0, W Events", 100,-100,100);
  welec_trackchi_plot = new TH1F("elec_track_chisq_W","Central track chi^2, W Events", 20,0,5);
  welec_dca_plot = new TH1F("elec_dca_W","Electron DCA, W Events",25,0,0.5);
  welec2_dca_plot = new TH1F("elec_dca_W2","Second Electron DCA, W Events",25,0,0.5);
  welec_phi_plot = new TH1F("elec_phi_W","Electron #phi, W Events",32,0,6.3);
  welec_ncft_plot = new TH1F("elec_ncft_W","Electron cft tracks, W Events",20,0,20);
  welec_nsmt_plot = new TH1F("elec_nsmt_W","Electron smt tracks, W Events",20,0,20);
  welec_triso_plot = new TH1F("elec_triso_W","Electron track isolation, W Events",20,0,10);
  welec_caliso_plot = new TH1F("elec_caliso_W","Electron cal isolation, W Events",20,0,10);
  welec_qcdiso1_plot = new TH1F("elec_qcdiso1_W","Tight/Loose Electron vs MET, W Events",50,0,100);
  welec_qcdiso2_plot = new TH1F("elec_qcdiso2_W","Tight/Loose Electron vs MET, W Events",50,0,100);
  welec_qcdiso3_plot = new TH1F("elec_qcdiso3_W","Tight/Loose Electron vs MET, W Events",50,0,100);
  welec_qcdiso4_plot = new TH1F("elec_qcdiso4_W","Tight/Loose Electron vs MET, W Events",50,0,100);
  welec_dRjet_plot = new TH1F("elec_dRjet_W","#Delta R elec-jet",64,0,6.3);
  welec_dRjetPI_plot = new TH1F("elec_dRjetPI_W","#Delta R elec-jet",64,0,6.3);
  welec_dRjet1_plot = new TH1F("elec_dRjet1_W","#Delta R elec-jet1",64,0,6.3);
  welec_dRjet2_plot = new TH1F("elec_dRjet2_W","#Delta R elec-jet2",64,0,6.3);
  wenu_met_plot = new TH1F("enu_met_W","MET, W Events",50,0,100);
  wenu_met2jet_plot = new TH1F("enu_met2jet_W","MET, 2-jet bin, W Events",50,0,100);
  wenu_met3jet_plot = new TH1F("enu_met3jet_W","MET, 3-jet bin, W Events",50,0,100);
  wenu_met4jet_plot = new TH1F("enu_met4jet_W","MET, 4-jet bin, W Events",50,0,100);
  wenu_metphi_plot = new TH1F("enu_metphi_W","MET #phi, W Events",64,0,6.3);
  wenu2_met_plot = new TH1F("enu_met_W2","Second MET, W Events",50,0,100);
  wenu2_metphi_plot = new TH1F("enu_metphi_W2","Second MET #phi, W Events",64,0,6.3);
 
  dphi_emet_plot = new TH2F("dphi_emet_plot","Triangle #Delta #phi(#e-MET) vs. MET",100,0,200,32,0,3.2);
  dphi_Wenu_jetmet_plot = new TH2F("dphi_Wenu_jetmet_plot","Triangle #Delta #phi(jet-MET) vs. MET",100,0,200,32,0,3.2);

  wmuon_pt_plot->Sumw2();
  wmuon_eta_plot->Sumw2();
  wmuon_deteta_plot->Sumw2();
  wmuon2_pt_plot->Sumw2();
  wmuon2_eta_plot->Sumw2();
  wmuon2_deteta_plot->Sumw2();
  wmuon_dphi_plot->Sumw2();
  wmuon_z0_plot->Sumw2();
  wmuon_trackchi_plot->Sumw2();
  wmuon_dca_plot->Sumw2();
  wmuon2_dca_plot->Sumw2();
  wmuon_phi_plot->Sumw2();
  wmuon_ncft_plot->Sumw2();
  wmuon_nsmt_plot->Sumw2();
  wmuon_triso_plot->Sumw2();
  wmuon_caliso_plot->Sumw2();
  wmuon_qcdiso1_plot->Sumw2();
  wmuon_qcdiso2_plot->Sumw2();
  wmuon_qcdiso3_plot->Sumw2();
  wmuon_qcdiso4_plot->Sumw2();
  wmuon_dRjet_plot->Sumw2();
  wmuon_dRjetPI_plot->Sumw2();
  wmuon_dRjet1_plot->Sumw2();
  wmuon_dRjet2_plot->Sumw2();
  wmunu_met_plot->Sumw2();
  wmunu_met2jet_plot->Sumw2();
  wmunu_met3jet_plot->Sumw2();
  wmunu_met4jet_plot->Sumw2();
  wmunu_metphi_plot->Sumw2();
  wmunu2_met_plot->Sumw2();
  wmunu2_metphi_plot->Sumw2();

  dphi_mumet_plot->Sumw2();
  dphi_Wmunu_jetmet_plot->Sumw2();


  welec_pt_plot->Sumw2();
  welec_eta_plot->Sumw2();
  welec_deteta_plot->Sumw2();
  welec2_pt_plot->Sumw2();
  welec2_eta_plot->Sumw2();
  welec2_deteta_plot->Sumw2();
  welec_dphi_plot->Sumw2();
  welec_z0_plot->Sumw2();
  welec_trackchi_plot->Sumw2();
  welec_dca_plot->Sumw2();
  welec2_dca_plot->Sumw2();
  welec_phi_plot->Sumw2();
  welec_ncft_plot->Sumw2();
  welec_nsmt_plot->Sumw2();
  welec_triso_plot->Sumw2();
  welec_caliso_plot->Sumw2();
  welec_qcdiso1_plot->Sumw2();
  welec_qcdiso2_plot->Sumw2();
  welec_qcdiso3_plot->Sumw2();
  welec_qcdiso4_plot->Sumw2();
  welec_dRjet_plot->Sumw2();
  welec_dRjetPI_plot->Sumw2();
  welec_dRjet1_plot->Sumw2();
  welec_dRjet2_plot->Sumw2();
  wenu_met_plot->Sumw2();
  wenu_met2jet_plot->Sumw2();
  wenu_met3jet_plot->Sumw2();
  wenu_met4jet_plot->Sumw2();
  wenu_metphi_plot->Sumw2();
  wenu2_met_plot->Sumw2();
  wenu2_metphi_plot->Sumw2();

  dphi_emet_plot->Sumw2();
  dphi_Wenu_jetmet_plot->Sumw2();


  jet_num_plot = new TH1F("jet_num","Jet Multiplicity", 15,0,15);
  jet_num_inclusive_plot = new TH1F("jet_num_inclusive","Jet Multiplicity (Inclusive)", 15,0,15);

  jet_num_plot->Sumw2();
  jet_num_inclusive_plot->Sumw2();

  jet_allpt_plot = new TH1F("jet_pt_all","pT of all jets", 100,0,200);
  jet_leadingpt_plot = new TH1F("jet_pt_leading","pT of leading jet", 40,0,200);
  jet_leadingpt_uncorr = new TH1F("jet_pt_leading_uncorr","uncorrected pT of leading jet", 40,0,200);
  jet_leadingpt_jes = new TH1F("jet_pt_leading_jes","jes pT of leading jet", 40,0,200);
  jet_leadingpt_jesmu = new TH1F("jet_pt_leading_jesmu","jesmu pT of leading jet", 40,0,200);
  jet_leadingpt_jesmusmear = new TH1F("jet_pt_leading_jesmusmear","jesmusmear pT of leading jet", 40,0,200);
  jet_leadingpt_2j_plot = new TH1F("jet_pt_leading_2j","pT of leading jet", 40,0,200);
  jet_secondpt_plot = new TH1F("jet_pt_second","pT of second leading jet", 40,0,200);
  jet_secondpt_uncorr = new TH1F("jet_pt_second_uncorr","uncorrected pT of second jet", 40,0,200);
  jet_secondpt_jes = new TH1F("jet_pt_second_jes","jes pT of second jet", 40,0,200);
  jet_secondpt_jesmu = new TH1F("jet_pt_second_jesmu","jesmu pT of second jet", 40,0,200);
  jet_secondpt_jesmusmear = new TH1F("jet_pt_second_jesmusmear","jesmusmear pT of second jet", 40,0,200);
  jet_thirdpt_plot = new TH1F("jet_pt_third","pT of third leading jet", 40,0,200);

  jet_allpt_plot->Sumw2();
  jet_leadingpt_plot->Sumw2();
  jet_leadingpt_uncorr->Sumw2();
  jet_leadingpt_jes->Sumw2();
  jet_leadingpt_jesmu->Sumw2();
  jet_leadingpt_jesmusmear->Sumw2();
  jet_leadingpt_2j_plot->Sumw2();
  jet_secondpt_plot->Sumw2();
  jet_secondpt_uncorr->Sumw2();
  jet_secondpt_jes->Sumw2();
  jet_secondpt_jesmu->Sumw2();
  jet_secondpt_jesmusmear->Sumw2();
  jet_thirdpt_plot->Sumw2();

  jet_allrap_plot = new TH1F("jet_eta_all","#eta of all jets", 50,-5,5);
  jet_leadingrap_plot = new TH1F("jet_eta_leading","Rapidity of leading jet, >= 1 jet", 25, -5,5);
  jet_leadingrap_2j_plot = new TH1F("jet_eta_leading_2j","Rapidity of leading jet, >= 2 jet", 25,-5,5);
  jet_secondrap_plot = new TH1F("jet_eta_second_eta","Rapidity of second jet, >= 2 jet", 12,-5,5);
  jet_firstsecondrap_plot = new TH1F("jet_eta_first_second","Rapidity of first and second jets, >= 2 jet", 12,-5,5);
  jet_thirdrap_plot = new TH1F("jet_eta_third_eta","Rapidity of third jet, >= 3 jet", 12,-5,5);

  jet_allrap_plot->Sumw2();
  jet_leadingrap_plot->Sumw2();
  jet_leadingrap_2j_plot->Sumw2();
  jet_secondrap_plot->Sumw2();
  jet_firstsecondrap_plot->Sumw2();
  jet_thirdrap_plot->Sumw2();

  jet_allphi_plot = new TH1F("jet_phi_all","#phi of all jets", 16,0,6.3);
  jet_leadingphi_plot = new TH1F("jet_phi_leading","#phi of leading jet, >= 1 jet", 16,0,6.3);
  jet_leadingphi_2j_plot = new TH1F("jet_phi_leading_2j","#phi of leading jet, >= 2 jet", 16,0,6.3);
  jet_secondphi_plot = new TH1F("jet_phi_second","#phi of second jet, >= 2 jet", 16,0,6.3);
  jet_firstsecondphi_plot = new TH1F("jet_phi_first_second","#phi of first and second jets, >= 2 jet", 16,0,6.3);
  jet_thirdphi_plot = new TH1F("jet_phi_third","#phi of third jet, >= 3 jet", 16,0,6.3);

  jet_dphiW_plot = new TH1F("jet_dphiW","#Delta #phi, boson-jet", 32,0,3.2);
  jet_dphi_plot = new TH1F("jet_dphi","#Delta #phi, jet-jet", 32,0,3.2);
  jet_deta_plot = new TH1F("jet_deta","#Delta #eta, jet-jet", 60,0,6.0);
  jet_invmass_plot = new TH1F("jet_dijetinvmass","Invariant mass of two leading jets", 120,0,300);

  jet_allphi_plot->Sumw2();
  jet_leadingphi_plot->Sumw2();
  jet_leadingphi_2j_plot->Sumw2();
  jet_secondphi_plot->Sumw2();
  jet_firstsecondphi_plot->Sumw2();
  jet_thirdphi_plot->Sumw2();

  jet_dphiW_plot->Sumw2();
  jet_dphi_plot->Sumw2();
  jet_deta_plot->Sumw2();
  jet_invmass_plot->Sumw2();


  _bins = 100;
  _min = 0;
  _max = 200;
  wmu_pt_plot = new TH1F("wmu_pt","W->munu PT", _bins,_min,_max);
  wmu2_pt_plot = new TH1F("wmu2_pt","W2->munu PT", _bins,_min,_max);
  we_pt_plot = new TH1F("we_pt","W->enu PT", _bins,_min,_max);
  we2_pt_plot = new TH1F("we2_pt","W2->enu PT", _bins,_min,_max);
  jet_yes_plot = new TH1F("jet_yes","Boson pT, opposite jet", _bins,_min,_max);
  jet_no_plot = new TH1F("jet_no","Boson pT, no opposite jet", _bins,_min,_max);
  jet_prob_plot = new TH1F("jet_prob", "Jet Probability vs Boson pT",_bins,_min,_max);
  jet_prob_plot->SetStats(0);

  wmu_pt_plot->Sumw2();
  wmu2_pt_plot->Sumw2();
  we_pt_plot->Sumw2();
  we2_pt_plot->Sumw2();
  jet_yes_plot->Sumw2();
  jet_no_plot->Sumw2();
  jet_prob_plot->Sumw2();

  ///////////////  Ntuple for NN /////////////////
  nntree = new TTree("nntree","Tree for VBF NN");

  invmass = 0;
  wmu_deltaR_nearestJet = 0;
  we_deltaR_nearestJet = 0;
  event_weight = 0;

  nntree->Branch("Njet",&Njet,"Njet/I");
  nntree->Branch("invmass",&invmass,"invmass/D");
  nntree->Branch("jeteta_ind",&jeteta_ind,"jeteta_ind[Njet]/I");
 
  nntree->Branch("wmu_mupt",&wmu_mupt,"wmu_mupt/F");
  nntree->Branch("wmu_mueta",&wmu_mueta,"wmu_mueta/F");
  nntree->Branch("wmu_muphi",&wmu_muphi,"wmu_muphi/F");
  nntree->Branch("wmu_met",&wmu_met,"wmu_met/F");
  nntree->Branch("wmu_metphi",&wmu_metphi,"wmu_metphi/F");
  nntree->Branch("wmu_wpt",&wmu_wpt,"wmu_wpt/F");
  nntree->Branch("wmu_wtrmass",&wmu_wtrmass,"wmu_wtrmass/F");
  nntree->Branch("wmu_wdphi",&wmu_wdphi,"wmu_wdphi/F");
  nntree->Branch("wmu_lepz0",&wmu_lepz0,"wmu_lepz0/F");
  nntree->Branch("wmu_trackiso",&wmu_trackiso,"wmu_trackiso/F");
  nntree->Branch("wmu_caliso",&wmu_caliso,"wmu_caliso/F");
  nntree->Branch("wmu_weight",&event_weight,"wmu_weight/D");
  nntree->Branch("wmu_jetpt",&wmu_jetpt,"wmu_jetpt[Njet]/F");
  nntree->Branch("wmu_jeteta",&wmu_jeteta,"wmu_jeteta[Njet]/F");
  nntree->Branch("wmu_jetphi",&wmu_jetphi,"wmu_jetphi[Njet]/F");
  nntree->Branch("wmu_jetpx",&wmu_jetpx,"wmu_jetpx[Njet]/F");
  nntree->Branch("wmu_jetpy",&wmu_jetpy,"wmu_jetpy[Njet]/F");
  nntree->Branch("wmu_jetpz",&wmu_jetpz,"wmu_jetpz[Njet]/F");
  nntree->Branch("wmu_jetE",&wmu_jetE,"wmu_jetE[Njet]/F");
  nntree->Branch("wmu_deltaR_nearestJet",&wmu_deltaR_nearestJet,"wmu_deltaR_nearestJet/D");

  nntree->Branch("we_ept",&we_ept,"we_ept/F");
  nntree->Branch("we_eeta",&we_eeta,"we_eeta/F");
  nntree->Branch("we_ephi",&we_ephi,"we_ephi/F");
  nntree->Branch("we_met",&we_met,"we_met/F");
  nntree->Branch("we_metphi",&we_metphi,"we_metphi/F");
  nntree->Branch("we_wpt",&we_wpt,"we_wpt/F");
  nntree->Branch("we_wtrmass",&we_wtrmass,"we_wtrmass/F");
  nntree->Branch("we_wdphi",&we_wdphi,"we_wdphi/F");
  nntree->Branch("we_lepz0",&we_lepz0,"we_lepz0/F");
  nntree->Branch("we_trackiso",&we_trackiso,"we_trackiso/F");
  nntree->Branch("we_caliso",&we_caliso,"we_caliso/F");
  nntree->Branch("we_weight",&event_weight,"we_weight/D");
  nntree->Branch("we_jetpt",&we_jetpt,"we_jetpt[Njet]/F");
  nntree->Branch("we_jeteta",&we_jeteta,"we_jeteta[Njet]/F");
  nntree->Branch("we_jetphi",&we_jetphi,"we_jetphi[Njet]/F");
  nntree->Branch("we_jetpx",&we_jetpx,"we_jetpx[Njet]/F");
  nntree->Branch("we_jetpy",&we_jetpy,"we_jetpy[Njet]/F");
  nntree->Branch("we_jetpz",&we_jetpz,"we_jetpz[Njet]/F");
  nntree->Branch("we_jetE",&we_jetE,"we_jetE[Njet]/F");
  nntree->Branch("we_deltaR_nearestJet",&we_deltaR_nearestJet,"we_deltaR_nearestJet/D");
  
  // topological variables
  nntree->Branch("wmt",     &wmt,   "wmt/D");
  nntree->Branch("wet",     &wet,   "wet/D");
  nntree->Branch("weta",    &weta,  "weta/D");
  nntree->Branch("ht25",    &ht25,  "ht25/D");
  nntree->Branch("ht20",    &ht20,  "ht20/D");
  nntree->Branch("ht",      &ht,    "ht/D");
  nntree->Branch("apla",    &apla,  "apla/D");
  nntree->Branch("cent",    &cent,  "cent/D");
  nntree->Branch("plan",    &plan,  "plan/D");
  nntree->Branch("spher",   &spher, "spher/D");
  nntree->Branch("ht2p",    &ht2p,  "ht2p/D");
  nntree->Branch("ktminp",  &ktminp,"ktminp/D");
  nntree->Branch("ktminpme",&ktminpme,"ktminpme/D");

  nntree->Branch("jetmt",&jetmt,"jetmt/D");
  nntree->Branch("mevent",&mevent,"mevent/D");
  nntree->Branch("jetamax",&jetamax,"jetamax/D");
  nntree->Branch("jetasqsum",&jetasqsum,"jetasqsum/D");
  nntree->Branch("drmin",&drmin,"drmin/D");
  nntree->Branch("ht2",&ht2,"ht2/D");
  nntree->Branch("ht3",&ht3,"ht3/D");
  nntree->Branch("htall",&htall,"htall/D");
  nntree->Branch("htl",&htl,"htl/D");
  nntree->Branch("dijetmass",&dijetmass,"dijetmass/D");
  nntree->Branch("dijetptrel",&dijetptrel,"dijetptrel/D");
  nntree->Branch("htp",&htp,"htp/D");
  //     nntree->Branch("ht2p",&ht2p,"ht2p/D");
  nntree->Branch("ht3p",&ht3p,"ht3p/D");
  nntree->Branch("htpp",&htpp,"htpp/D");
  nntree->Branch("ht2pp",&ht2pp,"ht2pp/D");
  nntree->Branch("ht3pp",&ht3pp,"ht3pp/D");
  nntree->Branch("wpt",&wpt,"wpt/D");
  // end topological vars
  
  }

  void PlotVBF::finish()
  {

  std::cout <<"PlotVBF: processed "<<nevt<<" events with "<<nerr<<" errors";
  if(nevt>0) std::cout <<" = "<<(float)nerr/nevt<<" errors per event";
  std::cout<<std::endl;


  char filename[50];  
  sprintf(filename,"%s", _root_file_name.c_str() );
  TFile *file = new TFile(filename,"RECREATE");
  
  file->mkdir("Zmumu"); 
  file->cd("Zmumu"); 

  zmu_num_plot->Write(); 
  zmu_mass_plot->Write(); 
  zmu_pt_plot->Write(); 
  zmu_pt_1j_plot->Write(); 
  zmu_pt_2j_plot->Write(); 
  zmu_rap_plot->Write(); 
  zmu_rap_1j_plot->Write(); 
  zmu_rap_2j_plot->Write(); 
  zmu_phi_plot->Write(); 
  zmuon_highpt_plot->Write();
  zmuon_lowpt_plot->Write();
  zmuon_highpt_1j_plot->Write();
  zmuon_lowpt_1j_plot->Write();
  zmuon_highpt_2j_plot->Write();
  zmuon_lowpt_2j_plot->Write();
  zmuon_eta_p_plot->Write();
  zmuon_eta_p_1j_plot->Write();
  zmuon_eta_p_2j_plot->Write();
  zmuon_eta_d_plot->Write();
  zmuon_eta_d_1j_plot->Write();
  zmuon_eta_d_2j_plot->Write();
  zmuon_dphi_plot->Write();
  zmuon_dz0_plot->Write();
  zmuon_z0_plot->Write();
  zmuon_trackchi_plot->Write();

  file->cd(); 
  file->mkdir("Zee"); 
  file->cd("Zee"); 

  zee_num_plot->Write(); 
  zee_mass_plot->Write(); 
  zee_pt_plot->Write(); 
  zee_pt_1j_plot->Write(); 
  zee_pt_2j_plot->Write(); 
  zee_rap_plot->Write(); 
  zee_rap_1j_plot->Write(); 
  zee_rap_2j_plot->Write(); 
  zee_phi_plot->Write(); 
  zelec_highpt_plot->Write();
  zelec_lowpt_plot->Write();
  zelec_highpt_1j_plot->Write();
  zelec_lowpt_1j_plot->Write();
  zelec_highpt_2j_plot->Write();
  zelec_lowpt_2j_plot->Write();
  zelec_eta_p_plot->Write();
  zelec_eta_p_1j_plot->Write();
  zelec_eta_p_2j_plot->Write();
  zelec_eta_d_plot->Write();
  zelec_eta_d_1j_plot->Write();
  zelec_eta_d_2j_plot->Write();
  zelec_dphi_plot->Write();
  zelec_dz0_plot->Write();
  zelec_z0_plot->Write();
  zelec_trackchisq_plot->Write();

  file->cd(); 
  file->mkdir("Wmunu"); 
  file->cd("Wmunu"); 
  wmu_num_plot->Write();
  wmu_mass_plot->Write();
  wmu2_mass_plot->Write();
  wmu_pt_plot->Write();
  wmu2_pt_plot->Write();
  wmu_phi_plot->Write();
  wmu_rap_plot->Write();
  wmu2_phi_plot->Write();
  wmuon_pt_plot->Write();
  wmuon_eta_plot->Write();
  wmuon_deteta_plot->Write();
  wmuon2_pt_plot->Write();
  wmuon2_eta_plot->Write();
  wmuon2_deteta_plot->Write();
  wmuon_phi_plot->Write();
  wmuon_dphi_plot->Write();
  wmuon_z0_plot->Write();
  wmuon_trackchi_plot->Write();
  wmuon_dca_plot->Write();
  wmuon2_dca_plot->Write();
  wmuon_ncft_plot->Write();
  wmuon_nsmt_plot->Write();
  wmuon_triso_plot->Write();
  wmuon_caliso_plot->Write();
  wmuon_qcdiso1_plot->Write();
  wmuon_qcdiso2_plot->Write();
  wmuon_qcdiso3_plot->Write();
  wmuon_qcdiso4_plot->Write();
  wmuon_dRjet_plot->Write();
  wmuon_dRjetPI_plot->Write();
  wmuon_dRjet1_plot->Write();
  wmuon_dRjet2_plot->Write();
  wmunu_met_plot->Write();
  wmunu_met2jet_plot->Write();
  wmunu_met3jet_plot->Write();
  wmunu_met4jet_plot->Write();
  wmunu_metphi_plot->Write();
  wmunu2_met_plot->Write();
  wmunu2_metphi_plot->Write();
  dphi_mumet_plot->Write();
  dphi_Wmunu_jetmet_plot->Write();
  
  file->cd(); 
  file->mkdir("Wenu"); 
  file->cd("Wenu"); 
  we_num_plot->Write();
  we_mass_plot->Write();
  we_pt_plot->Write();
  we_phi_plot->Write();
  we_rap_plot->Write();
  we2_mass_plot->Write();
  we2_pt_plot->Write();
  we2_phi_plot->Write();
  welec_pt_plot->Write();
  welec_eta_plot->Write();
  welec_deteta_plot->Write();
  welec2_pt_plot->Write();
  welec2_eta_plot->Write();
  welec2_deteta_plot->Write();
  welec_dphi_plot->Write();
  welec_z0_plot->Write();
  welec_trackchi_plot->Write();
  welec_dca_plot->Write();
  welec2_dca_plot->Write();
  welec_phi_plot->Write();
  welec_ncft_plot->Write();
  welec_nsmt_plot->Write();
  welec_triso_plot->Write();
  welec_caliso_plot->Write();
  welec_qcdiso1_plot->Write();
  welec_qcdiso2_plot->Write();
  welec_qcdiso3_plot->Write();
  welec_qcdiso4_plot->Write();
  welec_dRjet_plot->Write();
  welec_dRjetPI_plot->Write();
  welec_dRjet1_plot->Write();
  welec_dRjet2_plot->Write();
  wenu_met_plot->Write();
  wenu_met2jet_plot->Write();
  wenu_met3jet_plot->Write();
  wenu_met4jet_plot->Write();
  wenu_metphi_plot->Write();
  wenu2_met_plot->Write();
  wenu2_metphi_plot->Write();

  dphi_emet_plot->Write();
  dphi_Wenu_jetmet_plot->Write();

  welec_jetnum_plot->Write();
  welec_jet1eta_plot->Write();
  welec_jet2eta_plot->Write();
  welec_jet3eta_plot->Write();

  
  file->cd(); 
  file->mkdir("Jets"); 
  file->cd("Jets"); 
  

  jet_num_plot->Write(); 
  jet_num_inclusive_plot->Write(); 

  jet_allpt_plot->Write();
  jet_leadingpt_plot->Write(); 
  jet_leadingpt_uncorr->Write(); 
  jet_leadingpt_jes->Write(); 
  jet_leadingpt_jesmu->Write(); 
  jet_leadingpt_jesmusmear->Write(); 
  jet_leadingpt_2j_plot->Write(); 
  jet_secondpt_plot->Write(); 
  jet_secondpt_uncorr->Write(); 
  jet_secondpt_jes->Write(); 
  jet_secondpt_jesmu->Write(); 
  jet_secondpt_jesmusmear->Write(); 
  jet_thirdpt_plot->Write(); 

  jet_allrap_plot->Write();
  jet_leadingrap_plot->Write(); 
  jet_leadingrap_2j_plot->Write(); 
  jet_secondrap_plot->Write(); 
  jet_firstsecondrap_plot->Write(); 
  jet_thirdrap_plot->Write(); 

  jet_allphi_plot->Write();
  jet_leadingphi_plot->Write(); 
  jet_leadingphi_2j_plot->Write(); 
  jet_secondphi_plot->Write(); 
  jet_firstsecondphi_plot->Write(); 
  jet_thirdphi_plot->Write(); 

  
  jet_dphiW_plot->Write(); 
  jet_dphi_plot->Write(); 
  jet_deta_plot->Write();
  jet_invmass_plot->Write();
  jet_yes_plot->Write(); 
  jet_no_plot->Write(); 

  jet_yes_plot->Sumw2();
  jet_no_plot->Sumw2();
  TH1F *den = new TH1F("den", "",_bins,_min,_max);
  den->Add(jet_yes_plot, jet_no_plot);
  jet_prob_plot->Divide(jet_yes_plot,den,1,1,"B");

  jet_prob_plot->Write();

  file->cd();
  nntree->Write();

  file->Write();

  
  std::cout <<"Written root file: "<<filename<<std::endl;

  }




} // namespace vbf_cafe

ClassImp(vbf_cafe::PlotVBF);