#define protos_anal_cxx
#include "protos_anal.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <iostream>
#include <fstream>
#include "TLorentzVector.h"

void protos_anal::Loop()
{

  double LEP_PTCUT = 25.;
  double LEP_ETACUT = 2.5;
  
  double JET_PTCUT = 25.;
  double JET_ETACUT = 2.5;

  // HERWIG    sc = 155   (POWHEG & MC@NLO)
  // HERWIG++  sc =  11   (POWHEG)
  // PYTHIA    sc =   2,3 (POWHEG & MG, PROTOS)
  // PYTHIA++  sc =  52   (POWHEG)
  // SHERPA
  // ALPGEN


  // ACER

  bool doHERWIG   = false;
  bool doHERWIGPP = false;
  bool doPYTHIA6  = true;
  bool doPYTHIA8  = false;
  
  int  WSTATCODE;
  if(doHERWIG)   WSTATCODE = 155;
  if(doHERWIGPP) WSTATCODE =  11;
  if(doPYTHIA6)  WSTATCODE =   2; // 3?
  if(doPYTHIA8)  WSTATCODE =  52;

  int HQUARK = 6;

  TFile* f;
  f = new TFile("histos_mB900.root","RECREATE");
  
  TH1F *h_weight     = new TH1F("h_weight","Event Weight",200,0-10,10);

  TH1F *h_Q_pt     = new TH1F("h_Q_pT","Heavy Quark pt",100,0,2000);
  TH1F *h_Q_eta    = new TH1F("h_Q_eta","Heavy Quark eta",50,-5,5);
  TH1F *h_Q_mass   = new TH1F("h_Q_mass","Heavy Quark mass",200,500,1500);
  TH1F *h_Q_ID     = new TH1F("h_Q_ID","Heavy Quark eta",20,-10,10);

  TH1F *h_Q_dID    = new TH1F("h_Q_dID","Heavy Quark daughter ID",60,-30,30);

  TH1F *h_W_pt     = new TH1F("h_W_pT","Decay W pt",100,0,2000);
  TH1F *h_W_eta    = new TH1F("h_W_eta","Decay W eta",50,-5,5);
  TH1F *h_Z_pt     = new TH1F("h_Z_pT","Decay Z pt",100,0,2000);
  TH1F *h_Z_eta    = new TH1F("h_Z_eta","Decay Z eta",50,-5,5);
  TH1F *h_H_pt     = new TH1F("h_H_pT","Decay H pt",100,0,2000);
  TH1F *h_H_eta    = new TH1F("h_H_eta","Decay H eta",50,-5,5);

  TH1F *h_qW_pt     = new TH1F("h_qW_pT","Decay quark (w/W) pt",100,0,2000);
  TH1F *h_qW_eta    = new TH1F("h_qW_eta","Decay quark (w/W) eta",50,-5,5);
  TH1F *h_qZ_pt     = new TH1F("h_qZ_pT","Decay quark (w/Z) pt",100,0,2000);
  TH1F *h_qZ_eta    = new TH1F("h_qZ_eta","Decay quark (w/Z) eta",50,-5,5);
  TH1F *h_qH_pt     = new TH1F("h_qH_pT","Decay quark (w/H) pt",100,0,2000);
  TH1F *h_qH_eta    = new TH1F("h_qH_eta","Decay quark (w/H) eta",50,-5,5);

  TH1F *h_qW_dR     = new TH1F("h_qW_dR","dR(Wq)",60,0,6);
  TH1F *h_qZ_dR     = new TH1F("h_qZ_dR","dR(Zq)",60,0,6);
  TH1F *h_qH_dR     = new TH1F("h_qH_dR","dR(Hq)",60,0,6);

  // ----  
  
  TH1F *h_lep_n     = new TH1F("h_lep_n","N Leptons",5,0,5);
  TH1F *h_lep_n_WW  = new TH1F("h_lep_n_WW","N Leptons (WW)",5,0,5);
  TH1F *h_lep_n_ZZ  = new TH1F("h_lep_n_ZZ","N Leptons (ZZ)",5,0,5);
  TH1F *h_lep_n_HH  = new TH1F("h_lep_n_HH","N Leptons (HH)",5,0,5);
  TH1F *h_lep_n_WZ  = new TH1F("h_lep_n_WZ","N Leptons (WZ)",5,0,5);
  TH1F *h_lep_n_WH  = new TH1F("h_lep_n_WH","N Leptons (WH)",5,0,5);
  TH1F *h_lep_n_ZH  = new TH1F("h_lep_n_ZH","N Leptons (ZH)",5,0,5);

  TH1F *h_el_n     = new TH1F("h_el_n","N Electrons",5,0,5);
  TH1F *h_el_pt    = new TH1F("h_el_pt","Electron pT",75,0,750); 
  TH1F *h_el_eta   = new TH1F("h_el_eta","Electron eta",50,-5,5);
  TH1F *h_el_PID   = new TH1F("h_el_PID","Electron parent ID",30,0,30);
  
  TH1F *h_mu_n     = new TH1F("h_mu_n","N Muon",5,0,5);
  TH1F *h_mu_pt    = new TH1F("h_mu_pt","Muon pT",75,0,750);
  TH1F *h_mu_eta   = new TH1F("h_mu_eta","Muon eta",50,-5,5);
  TH1F *h_mu_PID   = new TH1F("h_mu_PID","Electron parent ID",30,0,30);

  TH1F *h_MET     = new TH1F("h_MET","MET",75,0,750);
  TH1F *h_MET_WW  = new TH1F("h_MET_WW","MET (WW)",75,0,750);
  TH1F *h_MET_ZZ  = new TH1F("h_MET_ZZ","MET (ZZ)",75,0,750);
  TH1F *h_MET_HH  = new TH1F("h_MET_HH","MET (HH)",75,0,750);
  TH1F *h_MET_WZ  = new TH1F("h_MET_WZ","MET (WZ)",75,0,750);
  TH1F *h_MET_WH  = new TH1F("h_MET_WH","MET (WH)",75,0,750);
  TH1F *h_MET_ZH  = new TH1F("h_MET_ZH","MET (ZH)",75,0,750);
  
  TH1F *h_jet4pre_n  = new TH1F("h_jet4pre_n","Pre Jet4 N",20,0,20); 
  TH1F *h_jet4_n     = new TH1F("h_jet4_n","Jet4 N",20,0,20); 
  TH1F *h_jet4_pt    = new TH1F("h_jet4_pt","Jet4 pT",75,0,750); 
  TH1F *h_jet4_eta   = new TH1F("h_jet4_eta","Jet4 eta",50,-5.0,5.0);   
  TH1F *h_jet4_mass  = new TH1F("h_jet4_mass","Jet4 mass",40,0,200);   
  
  TH1F *h_jet6pre_n  = new TH1F("h_jet6pre_n","Pre Jet6 N",20,0,20); 
  TH1F *h_jet6_n     = new TH1F("h_jet6_n","Jet6 N",20,0,20); 
  TH1F *h_jet6_pt    = new TH1F("h_jet6_pt","Jet6 pT",75,0,750); 
  TH1F *h_jet6_eta   = new TH1F("h_jet6_eta","Jet6 eta",50,-5.0,5.0);   
  TH1F *h_jet6_mass  = new TH1F("h_jet6_mass","Jet6 mass",40,0,200);   
  
  TH1F *h_dR_eljet   = new TH1F("h_dR_eljet","dR(el,jet)",60,0,6); 
  TH1F *h_dR_mujet   = new TH1F("h_dR_mujet","dR(mu,jet)",60,0,6); 
  
  //--
  //TH1F *hem_el_pt    = new TH1F("hem_el_pt","Electron pT, em",30,0,300); 
  //TH1F *hem_el_eta   = new TH1F("hem_el_eta","Electron eta, em",50,-2.5,2.5); 
  //TH1F *hem_mu_pt    = new TH1F("hem_mu_pt","Muon pT, em",30,0,300); 
  //TH1F *hem_mu_eta   = new TH1F("hem_mu_eta","Muon eta, em",50,-2.5,2.5); 
  //TH1F *hem_jet_n      = new TH1F("hem_jet_n","Jet N, em",10,0,10); 
  //TH1F *hem_jet_pt     = new TH1F("hem_jet_pt","Jet pT, em",50,0,250); 
  //TH1F *hem_jet_eta    = new TH1F("hem_jet_eta","Jet eta, em",50,-5.0,5.0); 
  //TH1F *hem_dR_eljet   = new TH1F("hem_dR_eljet","dR(el,jet)",60,0,6); 
  //TH1F *hem_dR_mujet   = new TH1F("hem_dR_mujet","dR(mu,jet)",60,0,6); 
  //--
  
  TH1F *h_top_pt      = new TH1F("h_top_pt","Truth Top pT",50,0,1000); // 1000 
  TH1F *h_top_y       = new TH1F("h_top_y","Truth Top Rapidity",80,-4,4); 
  TH1F *h_ttsys_pt    = new TH1F("h_ttsys_pt","Truth pT(tt)",50,0,500); 
  TH1F *h_ttsys_dphi  = new TH1F("h_ttsys_dphi","Truth Dphi(tt)",64,0,3.2); 
  
  TH1F *h_bottom_pt    = new TH1F("h_bottom_pt","Truth Bottom pT",50,0,1000); 
  TH1F *h_bottom_y     = new TH1F("h_bottom_y","Truth Bottom Rapidity",80,-4,4); 
  
  //

  TH1F *h_N_truthW     = new TH1F("h_N_truthW","N W's passing truth",10,0,10); 
  TH1F *h_N_truthWlep  = new TH1F("h_N_truthWlep","N W's passing truth (lep)",10,0,10); 
  TH1F *h_N_truthWhad  = new TH1F("h_N_truthWhad","N W's passing truth (had)",10,0,10);  

  //TH1F *h_W_pt       = new TH1F("h_w_pt","Truth W pT",50,0,1000); 
  //TH1F *h_W_y        = new TH1F("h_w_y","Truth W Rapidity",80,-4,4); 
  TH1F *h_Wdecays     = new TH1F("h_Wdecays","W Decays: 0=leplep, 1=lephad",2,0,2);
  TH1F *hem_Wdecays   = new TH1F("hem_Wdecays","W Decays: 0=leplep, 1=lephad",2,0,2); 
  TH1F *h_Wchild0    = new TH1F("h_Wchild0","W Child 0",60,-30,30); 
  TH1F *h_Wchild1    = new TH1F("h_Wchild1","W Child 1",60,-30,30); 
  
  TH1F *h_evtype    = new TH1F("h_evtype","Event Type",6,0,6); 
 
  int nevs    = 0;
  int nevs_em = 0;
  int nevs_ee = 0;
  int nevs_mm = 0;
  
  TLorentzVector v_jetpre[20];
  TLorentzVector v_jet[20];
  TLorentzVector v_el[20];
  TLorentzVector v_mu[20];

  TLorentzVector v_top[2];
  TLorentzVector v_ttsys;
  TLorentzVector v_bottom[2];
  
  if (fChain == 0) return;
  Long64_t nentries = fChain->GetEntriesFast();
  Long64_t nbytes = 0, nb = 0;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    
    nevs++;
  
    h_weight->Fill(mc_event_weight);
  
    //if(jentry>5) break;
    cout << "Event:\t" << nevs << endl;
    
    int ntops = 0;
    int nbottoms = 0;
    
    bool haveWlep = false;
    bool haveWhad  = false; 

    int ntruthw = 0;
    int ntruthwlep = 0;
    int ntruthwhad = 0;

    bool haveWq = false;
    bool haveZq = false;
    bool haveHq = false;

    for(int i=0;i<mc_n;i++){

      // -------------------
      // Heavy Quarks (7,8)
      // ------------------
      if(fabs(mc_pdgId->at(i))==7){
	//cout << "PDGID:\t" << mc_pdgId->at(i) << "\tM:\t" << mc_m->at(i)/1000. << "\tpT:\t" << mc_pt->at(i)/1000. << endl;
	
	h_Q_pt->Fill(mc_pt->at(i)/1000.);
	h_Q_eta->Fill(mc_eta->at(i));
	h_Q_mass->Fill(mc_m->at(i)/1000.);
	h_Q_ID->Fill(mc_pdgId->at(i));

	// quark first, boson second
	//cout << mc_child_index->at(i).size() << "\t" << mc_pdgId->at(mc_child_index->at(i).at(0)) << "\t" << mc_pdgId->at(mc_child_index->at(i).at(1)) << endl;

	if( fabs(  mc_pdgId->at(mc_child_index->at(i).at(1)) )==24){
	  haveWq = true;
	  h_W_pt->Fill(mc_pt->at(mc_child_index->at(i).at(1))/1.e3);
	  h_W_eta->Fill(mc_eta->at(mc_child_index->at(i).at(1)));
	  h_qW_pt->Fill(mc_pt->at(mc_child_index->at(i).at(0))/1.e3);
	  h_qW_eta->Fill(mc_eta->at(mc_child_index->at(i).at(0)));
	  h_qW_dR->Fill(dR(mc_eta->at(mc_child_index->at(i).at(1)),mc_phi->at(mc_child_index->at(i).at(1)),mc_eta->at(mc_child_index->at(i).at(0)),mc_phi->at(mc_child_index->at(i).at(0))));
	}	
	if( fabs(  mc_pdgId->at(mc_child_index->at(i).at(1)) )==23){
	  haveZq = true;
	  h_Z_pt->Fill(mc_pt->at(mc_child_index->at(i).at(1))/1.e3);
	  h_Z_eta->Fill(mc_eta->at(mc_child_index->at(i).at(1)));
	  h_qZ_pt->Fill(mc_pt->at(mc_child_index->at(i).at(1))/1.e3);
	  h_qZ_eta->Fill(mc_eta->at(mc_child_index->at(i).at(1)));
	  h_qZ_dR->Fill(dR(mc_eta->at(mc_child_index->at(i).at(1)),mc_phi->at(mc_child_index->at(i).at(1)),mc_eta->at(mc_child_index->at(i).at(0)),mc_phi->at(mc_child_index->at(i).at(0))));
	}
	if( fabs(  mc_pdgId->at(mc_child_index->at(i).at(1)) )==25){
	  haveHq = true;
	  h_H_pt->Fill(mc_pt->at(mc_child_index->at(i).at(1))/1.e3);
	  h_H_eta->Fill(mc_eta->at(mc_child_index->at(i).at(1)));
	  h_qH_pt->Fill(mc_pt->at(mc_child_index->at(i).at(0))/1.e3);
	  h_qH_eta->Fill(mc_eta->at(mc_child_index->at(i).at(0)));
	  h_qH_dR->Fill(dR(mc_eta->at(mc_child_index->at(i).at(1)),mc_phi->at(mc_child_index->at(i).at(1)),mc_eta->at(mc_child_index->at(i).at(0)),mc_phi->at(mc_child_index->at(i).at(0))));
	}

	h_Q_dID->Fill( mc_pdgId->at(mc_child_index->at(i).at(0)) );
	h_Q_dID->Fill( mc_pdgId->at(mc_child_index->at(i).at(1)) );
	
      }
      
      //if(mc_status->at(i)==22) 
      //cout << "\t Have top \t" << mc_pdgId->at(i) << "\t" << mc_barcode->at(i) << "\t" << mc_status->at(i) << "\t" <<  mc_pt->at(i) << "\t" << mc_m->at(i) << endl;

      // ----------
      // Truth Tops
      // ----------
      // Find tops that decay to W's, W can be at(0) or at(1) in Pythia6
      
      //if( fabs( mc_pdgId->at(i) )==HQUARK){ 
      if( fabs( mc_pdgId->at(i) )==HQUARK &&
	  mc_child_index->at(i).size() && fabs(mc_pdgId->at(mc_child_index->at(i).at(1)))==24){
	
	//&& mc_pdgId->at(mc_parent_index->at(i).size()) && fabs(mc_pdgId->at(mc_parent_index->at(i).at(0)))!=6){
	//cout << "\t Have top \t" << mc_pdgId->at(i) << "\t" << mc_barcode->at(i) << "\t" << mc_status->at(i) << "\t" <<  mc_pt->at(i) << endl;
	//if(mc_pdgId->at(mc_parent_index->at(i).size())) cout << "\t Parent:\t" << mc_pdgId->at(mc_parent_index->at(i).at(0)) << endl;
	//if(mc_pdgId->at(mc_child_index->at(i).size())) cout << "\tChild:\t" << mc_pdgId->at(mc_child_index->at(i).at(1)) << endl;

	h_top_pt->Fill(mc_pt->at(i)/1.e3);

	if(ntops<2){
	  v_top[ntops].SetPtEtaPhiM(mc_pt->at(i)/1.e3,
				    mc_eta->at(i),
				    mc_phi->at(i),
				    mc_m->at(i)
				    );
	}
	ntops++;
      } // tops
      
      // ----------
      // Truth W
      // ----------
      // W's with specific status codes
      // HERWIG    sc = 155   (POWHEG & MC@NLO)
      // HERWIG++  sc =  11
      // PYTHIA    sc =   2,3 (POWHEG & MG)
      // PYTHIA++  sc =  52

      // powheg 3
      // no status code 3 in mc@nlo / jimmy (use 155)
      //if( fabs( mc_pdgId->at(i) )==24){
      if( fabs( mc_pdgId->at(i) )==24 && mc_status->at(i)==WSTATCODE){
	//cout << "\t Have W \t" << mc_pdgId->at(i) << "\t" << mc_barcode->at(i) << "\t" << mc_status->at(i) << "\t" <<  mc_pt->at(i) << endl;
	//cout << mc_child_index->at(i).size() << endl;
	//for(int j=0;j<mc_child_index->at(i).size();j++) cout << "\t" << mc_pdgId->at(mc_child_index->at(i).at(j)) << endl;
	//if(mc_parent_index->at(i).size() ) cout << "Parent:\t" << mc_pdgId->at(mc_parent_index->at(i).at(0)) << endl;
	if(mc_child_index->at(i).size() > 1){
	  //cout << "\t Have W \t" << mc_pdgId->at(i) << "\t" << mc_barcode->at(i) << "\t" << mc_status->at(i) << "\t" <<  mc_pt->at(i) //<< endl; 
	  //     << "\t" <<  mc_eta->at(i) <<  "\t" << mc_phi->at(i) << endl;
	  //cout << "\tNchild:\t" << mc_child_index->at(i).size() << endl;
	  //for(int j=0;j<mc_child_index->at(i).size();j++) cout << "Child" << j << "\t" << mc_pdgId->at(mc_child_index->at(i).at(j)) << endl;
	  //cout << "\tChild0:\t" << mc_pdgId->at(mc_child_index->at(i).at(0)) << endl;
	  //cout << "\tChild1:\t" << mc_pdgId->at(mc_child_index->at(i).at(1)) << endl;
	
	  ntruthw++;
	  //ntruthwlep = 0;
	  //ntruthwhad = 0;

	  //h_W_pt->Fill(mc_pt->at(i)/1.e3);
	  //h_W_y->Fill(mc_eta->at(i));
	  
	  if( fabs(mc_pdgId->at(mc_child_index->at(i).at(0))) > 4 ){
	    ntruthwlep++;
	    haveWlep = true;
	  }else{ 
	    ntruthwhad++;
	    haveWhad = true;
	  }
	  
	  h_Wchild0->Fill( mc_pdgId->at(mc_child_index->at(i).at(0)) ); // 11, 13, 15 ; 1, 3
	  h_Wchild1->Fill( mc_pdgId->at(mc_child_index->at(i).at(1)) ); // 12, 14, 16 ; 2, 4
	}   
      }   
      
      // ------------------
      // Truth b's from top
      // ------------------
      // b's with specific status codes
      // powheg 3
      // no status code 3 in mc@nlo (use 124)
      //if( fabs( mc_pdgId->at(i) )==5  && fabs( mc_status->at(i) )==3 ){
      if( fabs( mc_pdgId->at(i) )==5 && mc_parent_index->at(i).size()>0 && fabs(mc_pdgId->at(mc_parent_index->at(i).at(0)))==6 ){ // herwig++
	//cout << "\t Have b \t" << mc_pdgId->at(i) << "\t" << mc_barcode->at(i) << "\t" << mc_status->at(i) << "\t" <<  mc_pt->at(i) << endl;
	if(mc_parent_index->at(i).size() >0){
	  //cout << "\tParent:\t" << mc_pdgId->at(mc_parent_index->at(i).at(0)) << endl;
	}
	
	if(nbottoms<2){
	  v_bottom[nbottoms].SetPtEtaPhiM(mc_pt->at(i)/1.e3,
					  mc_eta->at(i),
					  mc_phi->at(i),
					  mc_m->at(i)
					  );
	}
	nbottoms++;
	
	h_bottom_pt->Fill(mc_pt->at(i)/1.e3);
	h_bottom_y->Fill(mc_eta->at(i));

      }

    } // particle loop

    int evtype;
    //cout << haveWq << "\t" << haveZq << "\t" << haveHq << endl;
    if(haveWq==true  && haveZq == false && haveHq == false){ cout << "WqWq event" << endl; evtype=0; h_evtype->Fill(evtype);}
    if(haveWq==false && haveZq == true  && haveHq == false){ cout << "ZqZq event" << endl; evtype=1; h_evtype->Fill(evtype);}
    if(haveWq==false && haveZq == false && haveHq == true ){ cout << "HqHq event" << endl; evtype=2; h_evtype->Fill(evtype);}
    if(haveWq==true  && haveZq == true  && haveHq == false){ cout << "WqZq event" << endl; evtype=3; h_evtype->Fill(evtype);}
    if(haveWq==true  && haveZq == false && haveHq == true ){ cout << "WqHq event" << endl; evtype=4; h_evtype->Fill(evtype);}
    if(haveWq==false && haveZq == true  && haveHq == true ){ cout << "ZqHq event" << endl; evtype=5; h_evtype->Fill(evtype);}

    v_ttsys = v_top[0] + v_top[1];

    h_top_y->Fill(v_top[0].PseudoRapidity());
    h_top_y->Fill(v_top[1].PseudoRapidity());

    h_ttsys_pt->Fill(v_ttsys.Pt());
    h_ttsys_dphi->Fill(fabs(v_top[0].DeltaPhi(v_top[1])));    

    //cout << haveWlep << "\t" << haveWhad << endl;
    h_N_truthW->Fill(ntruthw);
    h_N_truthWlep->Fill(ntruthwlep);
    h_N_truthWhad->Fill(ntruthwhad);

    if(haveWhad) h_Wdecays->Fill(1.);
    else  h_Wdecays->Fill(0.);

    // ---------------
    // Truth Electrons
    // ---------------
    
    int nsel_el = 0;
    for(int i=0;i<el_n;i++){
      if( fabs(mc_pdgId->at(el_parent_index->at(i).at(0)))==24 ||  fabs(mc_pdgId->at(el_parent_index->at(i).at(0)))==23){
	//cout << "ele:\t" << el_dressed_pt->at(i)/1.e3 << "\t" <<  el_dressed_pt->at(i)/1.e3 << "\t" <<mc_pdgId->at(el_parent_index->at(i).at(0)) << endl;
        
	// if(el_truth_pt->at(i)/1.e3 > 25. && fabs(el_truth_eta->at(i))<2.5 ){
	nsel_el++;
	v_el[i].SetPtEtaPhiM(el_dressed_pt->at(i)/1.e3,
	 		     el_dressed_eta->at(i),
	 		     el_dressed_phi->at(i),
	 		     el_dressed_m->at(i)
	 		     );
	
	h_el_pt->Fill(el_dressed_pt->at(i)/1.e3);
	h_el_eta->Fill(el_dressed_eta->at(i));
	h_el_PID->Fill(  fabs(mc_pdgId->at(el_parent_index->at(i).at(0))) );
	
      }
    }
    h_el_n->Fill(nsel_el);
    
    
    // -----------
    // Truth Muons
    // -----------
    
    int nsel_mu = 0;
    for(int i=0;i<mu_n;i++){
      if( fabs(mc_pdgId->at(mu_parent_index->at(i).at(0)))==24 ||  fabs(mc_pdgId->at(mu_parent_index->at(i).at(0)))==23){

	// if(mu_staco_truth_pt->at(i)/1.e3 > 25. && fabs(mu_staco_truth_eta->at(i))<2.5 ){
	nsel_mu++;
	v_mu[i].SetPtEtaPhiM(mu_dressed_pt->at(i)/1.e3,
			     mu_dressed_eta->at(i),
			     mu_dressed_phi->at(i),
			     mu_dressed_m->at(i)
			     );

	h_mu_pt->Fill(mu_dressed_pt->at(i)/1.e3);
	h_mu_eta->Fill(mu_dressed_eta->at(i));
	h_mu_PID->Fill(  fabs(mc_pdgId->at(mu_parent_index->at(i).at(0))) );
	
      }
    }
    h_mu_n->Fill(nsel_mu);
    
    h_lep_n->Fill(nsel_mu+nsel_el);  
    if(evtype==0) h_lep_n_WW->Fill(nsel_mu+nsel_el); 
    if(evtype==1) h_lep_n_ZZ->Fill(nsel_mu+nsel_el); 
    if(evtype==2) h_lep_n_HH->Fill(nsel_mu+nsel_el); 
    if(evtype==3) h_lep_n_WZ->Fill(nsel_mu+nsel_el); 
    if(evtype==4) h_lep_n_WH->Fill(nsel_mu+nsel_el); 
    if(evtype==5) h_lep_n_ZH->Fill(nsel_mu+nsel_el); 

    // ---
    // MET
    // ---

    double MET = TMath::Sqrt(MET_Truth_NonInt_etx*MET_Truth_NonInt_etx+MET_Truth_NonInt_ety*MET_Truth_NonInt_ety)/1.e3;
    
    h_MET->Fill(MET);  
    if(evtype==0) h_MET_WW->Fill(MET);
    if(evtype==1) h_MET_ZZ->Fill(MET);
    if(evtype==2) h_MET_HH->Fill(MET);
    if(evtype==3) h_MET_WZ->Fill(MET);
    if(evtype==4) h_MET_WH->Fill(MET);
    if(evtype==5) h_MET_ZH->Fill(MET);


    // ---------------
    // Truth Jet (Pre)
    // ---------------
    int nsel_jetpre = 0;
    for(int i=0;i<jet_AntiKt4TruthJets_WZ_n;i++){
      if(jet_AntiKt4TruthJets_WZ_pt->at(i)/1.e3 > JET_PTCUT && fabs(jet_AntiKt4TruthJets_WZ_eta->at(i))<JET_ETACUT ){
     	nsel_jetpre++;
	// v_jetpre[i].SetPtEtaPhiM(jet_AntiKt4TruthJets_pt->at(i)/1.e3,
	// 			 jet_AntiKt4TruthJets_eta->at(i),
	// 			 jet_AntiKt4TruthJets_phi->at(i),
	// 			 jet_AntiKt4TruthJets_m->at(i)
	// 			 );
	
	h_jet4_pt->Fill(jet_AntiKt4TruthJets_WZ_pt->at(i)/1.e3);
    	h_jet4_eta->Fill(jet_AntiKt4TruthJets_WZ_eta->at(i));
	h_jet4_mass->Fill(jet_AntiKt4TruthJets_WZ_m->at(i)/1.e3);
     	//cout << "\t Have jet" << endl;
      }
    }
    h_jet4pre_n->Fill(nsel_jetpre);


    nsel_jetpre = 0;
    for(int i=0;i<jet_AntiKt6TruthJets_WZ_n;i++){
      if(jet_AntiKt6TruthJets_WZ_pt->at(i)/1.e3 > JET_PTCUT && fabs(jet_AntiKt6TruthJets_WZ_eta->at(i))<JET_ETACUT ){
     	nsel_jetpre++;
	// v_jetpre[i].SetPtEtaPhiM(jet_AntiKt4TruthJets_pt->at(i)/1.e3,
	// 			 jet_AntiKt4TruthJets_eta->at(i),
	// 			 jet_AntiKt4TruthJets_phi->at(i),
	// 			 jet_AntiKt4TruthJets_m->at(i)
	// 			 );
	
	h_jet6_pt->Fill(jet_AntiKt6TruthJets_WZ_pt->at(i)/1.e3);
    	h_jet6_eta->Fill(jet_AntiKt6TruthJets_WZ_eta->at(i));
	h_jet6_mass->Fill(jet_AntiKt6TruthJets_WZ_m->at(i)/1.e3);
     	//cout << "\t Have jet" << endl;
      }
    }
    h_jet6pre_n->Fill(nsel_jetpre);



    // // ---------
    // // Truth Jet
    // // --> exclude jets close to electrons and muons
    // // ---------
    // int nsel_jet = 0;
    
    // for(int i=0;i<nsel_jetpre;i++){
      
    //   bool isOverlap = false;
    //   for(int j=0;j<nsel_el;j++) 
    // 	if( v_jetpre[i].DeltaR(v_el[j]) < 0.1 ) isOverlap=true; 
      
    //   for(int j=0;j<nsel_mu;j++) 
    // 	if( v_jetpre[i].DeltaR(v_mu[j]) < 0.1 ) isOverlap=true; 
      
    //   if(!isOverlap){
    // 	v_jet[nsel_jet].SetPtEtaPhiM(v_jetpre[i].Pt(),
    // 				     v_jetpre[i].Eta(),
    // 				     v_jetpre[i].Phi(),
    // 				     v_jetpre[i].M()
    // 				     );
    // 	nsel_jet++;
	
    // 	h_jet_pt->Fill(jet_AntiKt4TruthJets_pt->at(i)/1.e3);
    // 	h_jet_eta->Fill(jet_AntiKt4TruthJets_eta->at(i));
    //   }
    // }
    
    // h_jet_n->Fill(nsel_jet);
    
    // // xxxxxxxxxxxxxxxxxxxxxxxxxxx
    
    // //if(nsel_el==1 && nsel_mu==1 && !haveWhad){
    // if(nsel_el==1 && nsel_mu==1){
    //   nevs_em++;
      
    //   //cout << "Have em event" << endl;
    //   //cout << haveWlep << "\t" << haveWhad << endl;

    //   if(haveWhad) hem_Wdecays->Fill(1.);
    //   else  hem_Wdecays->Fill(0.);

    //   hem_jet_n->Fill(nsel_jet);
    //   for(int j=0;j<nsel_jet;j++){
    // 	hem_dR_eljet->Fill( v_jet[j].DeltaR(v_el[0]) );	
    // 	hem_dR_mujet->Fill( v_jet[j].DeltaR(v_mu[0]) );	
    //   	hem_jet_pt->Fill(v_jet[j].Pt());
    // 	hem_jet_eta->Fill(v_jet[j].Eta());
    //   }

    //   hem_el_pt->Fill(v_el[0].Pt());
    //   hem_el_eta->Fill(v_el[0].Eta());
      
    //   hem_mu_pt->Fill(v_mu[0].Pt());
    //   hem_mu_eta->Fill(v_mu[0].Eta());

    //   //cout << "Electron:\t" << v_el[0].Pt() << "\t" << v_el[0].Eta() << "\t" << v_el[0].Phi() << endl;
    //   //cout << "\t" << v_el[0].DeltaR(v_bottom[0]) << "\t" << v_el[0].DeltaR(v_bottom[1]) << endl;
    //   //cout << "Muon    :\t" << v_mu[0].Pt() << "\t" << v_mu[0].Eta() << "\t" << v_mu[0].Phi() << endl;
    //   //cout << "\t" << v_mu[0].DeltaR(v_bottom[0]) << "\t" << v_mu[0].DeltaR(v_bottom[1]) << endl;
    //   //cout << "Bottom 1:\t" << v_bottom[0].Pt() << "\t" << v_bottom[0].Eta() << "\t" << v_bottom[0].Phi() << endl;
    //   //cout << "Bottom 2:\t" << v_bottom[1].Pt() << "\t" << v_bottom[1].Eta() << "\t" << v_bottom[1].Phi() << endl;
    //   //cout << "Jets" << endl;
    //   //for(int j=0;j<nsel_jet;j++){
    //   //cout << v_jet[j].Pt() << "\t" << v_jet[j].Eta() << "\t" << v_jet[j].Phi() << endl;
    //   //}
    
    // }
    
    //if(nsel_el==2) nevs_ee++;
    //if(nsel_mu==2) nevs_mm++;
  }
  
  f->Write();
  
  cout << "Total Events    :\t" << nevs    << endl;
  cout << "Total Events, em:\t" << nevs_em << endl;
  cout << "Total Events, ee:\t" << nevs_ee << endl;
  cout << "Total Events, mm:\t" << nevs_mm << endl;
}

Double_t dR( double eta1 , double phi1, double eta2 , double phi2 ){
  
  Double_t dR, deta, dphi;
  
  deta = fabs(eta1-eta2);
  dphi = fabs(phi1-phi2);
  if( dphi > 3.14159 ){
    dphi = 2*3.14159 - dphi;
  }
  
  dR = sqrt( TMath::Power(deta,2) + TMath::Power(dphi,2) ); 
  return dR;
}