#include "ZtautauAnalysis.h"
#include "ZtautauTree.h"
#include "EtTree.h"
#include "TauTree.h"
#include "ElectronTree.h"
#include "MuonTree.h"
#include "WeightTree.h"

#include "TTree.h"
#include "TH1F.h"
#include "TFile.h"
#include "TString.h"

#include <iostream>
#include <iomanip>
#include <cmath>

int nEventMET(0);
int nEventCand(0);

int nTauMET(0);
int nTauCand(0);

ZtautauAnalysis::ZtautauAnalysis(std::string name, float xsect, float lumi, TTree* ztautauTree, TTree* etTree, TTree* tauTree, TTree* elecTree, TTree* muonTree, TTree* weightTree): 
  m_cutFlow(10, 0),
  m_cutFlowSq(10, 0),
  m_cutNames(10, ""),
  m_name(name),
  m_xsect(xsect),
  m_lumi(lumi),
  m_lumiScale(1.),
  m_doSafe(false)
{
  m_ztautauTree = (ztautauTree ? new ZtautauTree(ztautauTree): 0);
  m_etTree      = (etTree      ? new EtTree(etTree): 0);
  m_tauTree     = (tauTree     ? new TauTree(tauTree): 0);
  m_elecTree    = (elecTree    ? new ElectronTree(elecTree): 0);
  m_muonTree    = (muonTree    ? new MuonTree(muonTree): 0);
  m_weightTree  = (weightTree  ? new WeightTree(weightTree): 0);

  if (etTree->GetEntries() != tauTree->GetEntries()   || elecTree->GetEntries() != tauTree->GetEntries() ||
      muonTree->GetEntries() != tauTree->GetEntries() || weightTree->GetEntries() != tauTree->GetEntries()) {
    std::cout << "Control trees have different numbers of entries" << std::endl;
    exit(1);
  }
}

ZtautauAnalysis::~ZtautauAnalysis(){}

void ZtautauAnalysis::Initialise() {

  // Cut names
  m_cutNames[0] = "lepton-tau candidates";
  m_cutNames[1] = "DeltaPhi";
  m_cutNames[2] = "Tau ID";
  m_cutNames[3] = "OS events";
  m_cutNames[4] = "Mass window";
  m_cutNames[5] = "SS events";
  m_cutNames[6] = "Mass Window";
  
  // Scale
  SetLumiScale();

  // Hists
  BookHists();
}

void ZtautauAnalysis::BookHists() {
  m_hist_vismass_os = new TH1F("VisMassOS","VisMassSS;m_{#tau,lep};# Events", 40, 0, 200);  m_hist_vismass_os->Sumw2();
  m_hist_vismass_ss = new TH1F("VisMassSS","VisMassOS;m_{#tau,lep};# Events", 40, 0, 200);  m_hist_vismass_ss->Sumw2();
  m_hist_invmass_os = new TH1F("InvMassOS","InvMassSS;m_{#tau#tau};# Events", 60, 0, 300);  m_hist_invmass_os->Sumw2();
  m_hist_invmass_ss = new TH1F("InvMassSS","InvMassOS;m_{#tau#tau};# Events", 60, 0, 300);  m_hist_invmass_ss->Sumw2();

  m_hist_invpt      = new TH1F("InvPt",    "InvPt;P_{t}^{z};# Events",        25, 0, 100);    m_hist_invpt->Sumw2();
  m_hist_invphi     = new TH1F("InvPhi",   "InvPhi;#phi^{z};# Events",        18, -180, 180); m_hist_invphi->Sumw2();
  m_hist_inveta     = new TH1F("InvEta",   "InvEta;#eta^{z};# Events",        25, 0, 2.5);    m_hist_inveta->Sumw2();

  m_hist_vispt      = new TH1F("VisPt",    "VisPt;P_{t}^{z} (vis);# Events",  25, 0, 100);    m_hist_vispt->Sumw2();
  m_hist_visphi     = new TH1F("VisPhi",   "VisPhi;#phi^{z} (vis);# Events",  18, -180, 180); m_hist_visphi->Sumw2();
  m_hist_viseta     = new TH1F("VisEta",   "VisEta;#eta^{z} (vis);# Events",  25, 0, 2.5);    m_hist_viseta->Sumw2();

  m_hist_visdphi    = new TH1F("VisDPhi",   "VisDPhi;#delta#phi;# Events",    36, 0, 180);  m_hist_visdphi->Sumw2();

  m_hist_etmiss     = new TH1F("EtMiss",   "EtMiss;MET;# Events",             50, 0, 100);  m_hist_etmiss->Sumw2();
  m_hist_etsum      = new TH1F("EtSum",    "EtSum; #Sigma Et;# Events",       60, 0, 600);  m_hist_etsum->Sumw2();
  m_hist_mt         = new TH1F("Mt",       "Mt;m_{t}^{lep.MET};# Events",     50, 0, 150);  m_hist_mt->Sumw2();
  m_hist_mt_etmiss  = new TH1F("MtEtMiss", "Mt(MET);m_{t}^{lep.MET};# Events",50, 0, 150);  m_hist_mt_etmiss->Sumw2();
  m_hist_etmiss_mt  = new TH1F("EtMissMt", "EtMiss(Mt);MET;# Events",         50, 0, 100);  m_hist_etmiss_mt->Sumw2();
  m_hist_etsum_mtmet= new TH1F("EtSumMtEtMiss", "EtSum(Mt MET); #Sigma Et;# Events", 30, 0, 600);  m_hist_etsum_mtmet->Sumw2();

  m_hist_nprong       = new TH1F("Nprong",       "NProng;Track Mult.;# Taus",         10, 0, 10);  m_hist_nprong->Sumw2();
  m_hist_nprong_met   = new TH1F("NprongMET",    "NProng(MET);Track Mult.;# Taus",    10, 0, 10);  m_hist_nprong_met->Sumw2();
  m_hist_nprong_mt    = new TH1F("NprongMt",     "NProng(Mt);Track Mult.;# Taus",     10, 0, 10);  m_hist_nprong_mt->Sumw2();
  m_hist_nprong_sumet = new TH1F("NprongSumEt",  "NProng(SumEt);Track Mult.;# Taus",  10, 0, 10);  m_hist_nprong_sumet->Sumw2();
  m_hist_nprong_dphi  = new TH1F("NprongDphi",   "NProng(Dphi);Track Mult.;# Taus",   10, 0, 10);  m_hist_nprong_dphi->Sumw2();
  m_hist_nprong_cand  = new TH1F("NprongCand",   "NProng(Cand);Track Mult.;# Taus",   10, 0, 10);  m_hist_nprong_cand->Sumw2();
  m_hist_nprong_tau   = new TH1F("NprongTau",    "NProng(Tau);Track Mult.;# Taus",    10, 0, 10);  m_hist_nprong_tau->Sumw2();
  m_hist_nprong_chrg  = new TH1F("NprongChrg",   "NProng(Chrg);Track Mult.;# Taus",   10, 0, 10);  m_hist_nprong_chrg->Sumw2();
  m_hist_nprong_vism  = new TH1F("NprongVisM",   "NProng(VisM);Track Mult.;# Taus",   10, 0, 10);  m_hist_nprong_vism->Sumw2();

  m_hist_nprong_chrg_os = new TH1F("NprongChrgOS", "NProng(ChrgOS);Track Mult.;# Taus", 10, 0, 10); m_hist_nprong_chrg_os->Sumw2();
  m_hist_nprong_vism_os = new TH1F("NprongVisMOS", "NProng(VisMOS);Track Mult.;# Taus", 10, 0, 10); m_hist_nprong_vism_os->Sumw2();

  m_hist_charge_tau  = new TH1F("ChargeTau",     "Charge(Tau);#tau charge;# Cands",       10, -5, 5);  m_hist_charge_tau->Sumw2();
  m_hist_charge_lep  = new TH1F("ChargeLep",     "Charge(Lep);lep charge; # Cands",       10, -5, 5);  m_hist_charge_lep->Sumw2();
  m_hist_charge_prod = new TH1F("ChargeProd",    "Charge(Prod);lep*#tau charge;# Cands",  10, -5, 5);  m_hist_charge_prod->Sumw2();

  m_hist_emradius    = new TH1F("EmRadius",      "EMRadius;EM Radius;# Taus",             40, 0, 0.4 );  m_hist_emradius->Sumw2();
  m_hist_width2      = new TH1F("Width2",        "Width2;StripWidth2;# Taus",              50, 0, 0.05);  m_hist_width2->Sumw2();
  m_hist_isolfrac    = new TH1F("IsolFrac",      "IsolFrac;Isolation fraction;# Taus",    50, 0, 1   );  m_hist_isolfrac->Sumw2();
  m_hist_etratio     = new TH1F("EtRatio",       "EtRatio;E_{t}(EM)/E_{t};# Taus",        55, 0, 1.1 );  m_hist_etratio->Sumw2();


  int nCuts(m_cutFlow.size());
  m_hist_cuts       = new TH1F("CutHist","CutHist;Cut Number;NEvents",nCuts, 0, nCuts);             m_hist_cuts->Sumw2();
  m_hist_cutsscaled = new TH1F("CutHistScaled","CutHistScaled;Cut Number;NEvents",nCuts, 0, nCuts); m_hist_cutsscaled->Sumw2();
}

void ZtautauAnalysis::EtMiss(int nevents)
{
  float sumw(0.), sumws(0.);
  int ntot(0), n1lep(0);

  for (int ievent(0); ievent < nevents; ievent++) {
    m_etTree->LoadTree(ievent);
    m_etTree->GetEntry(ievent);

    m_weightTree->LoadTree(ievent);
    m_weightTree->b_weight_weight->GetEntry(ievent);
    float weight = m_weightTree->weight_weight;

    //float weight = m_etTree->et_weight;
    sumw += weight;
    sumws += weight;
    ntot++;

    m_elecTree->LoadTree(ievent);
    m_elecTree->b_electron_ismedium->GetEntry(ievent);

    m_muonTree->LoadTree(ievent);
    m_muonTree->b_muon_pt->GetEntry(ievent);

    m_tauTree->LoadTree(ievent);
    m_tauTree->b_tau_id_TauCutTight->GetEntry(ievent);
    m_tauTree->b_tau_id_TauCutSafeCaloTight->GetEntry(ievent);

    int nElec(0), nMuon(m_muonTree->muon_pt->size()), nLeptons(nMuon);
    for (int i(0); i < m_elecTree->electron_ismedium->size(); i++) {
      if (m_elecTree->electron_ismedium->at(i) == 1) {nElec++, nLeptons++;}
    }

    //cout << "Ne, Ne(tight), Nmu, Nl = (" << m_elecTree->electron_ismedium->size() << ") " << nElec  << " " << nMuon << " " << nLeptons << endl;
    //cout << nElec  << " " << nMuon << " " << ((nLeptons==1) ? "passed" : "failed") << endl;

    if (nLeptons!=1) continue;

    n1lep++;

    bool hasTightTau(false);    
    for (int i(0); i < m_tauTree->tau_id_TauCutTight->size(); i++) {
      bool tauId = (m_doSafe ? (m_tauTree->tau_id_TauCutSafeCaloTight->at(i) == 1) : (m_tauTree->tau_id_TauCutTight->at(i) == 1));
      //if (m_tauTree->tau_id_TauCutTight->at(i) == 1) {hasTightTau = true; break;}
      if (tauId) {hasTightTau = true; break;}
    }

    if (!hasTightTau) continue;

    m_hist_etmiss->Fill(m_etTree->et_et, weight);
    m_hist_etsum ->Fill(m_etTree->et_sumet, weight);
    m_hist_mt    ->Fill(m_etTree->et_lepton_mt, weight);    

    if (m_etTree->et_lepton_mt < 30) m_hist_etmiss_mt->Fill(m_etTree->et_et, weight);    
    if (m_etTree->et_et > 20) m_hist_mt_etmiss->Fill(m_etTree->et_lepton_mt, weight);    
    if (m_etTree->et_et > 20 && m_etTree->et_lepton_mt < 30) m_hist_etsum_mtmet->Fill(m_etTree->et_sumet, weight);    
    
  }

  std::cout << "MET Events = " << nevents << " " << ntot << " " << sumw << " " << " " << sumws << " " << sumw * m_lumiScale << " " << n1lep << std::endl;
}

void ZtautauAnalysis::Nprong(int nevents) 
{
   for (int ievent(0); ievent < nevents; ievent++) {
     m_etTree->LoadTree(ievent);
     m_etTree->GetEntry(ievent);
     
     m_weightTree->LoadTree(ievent);
     m_weightTree->b_weight_weight->GetEntry(ievent);
     float weight = m_weightTree->weight_weight;

     //float weight = m_etTree->et_weight;
     
     m_elecTree->LoadTree(ievent);
     m_elecTree->b_electron_ismedium->GetEntry(ievent);
     
     m_muonTree->LoadTree(ievent);
     m_muonTree->b_muon_pt->GetEntry(ievent);
     
     m_tauTree->LoadTree(ievent);
     m_tauTree->b_tau_id_TauCutTight->GetEntry(ievent);
     m_tauTree->b_tau_id_TauCutSafeCaloTight->GetEntry(ievent);
     m_tauTree->b_tau_numTracks->GetEntry(ievent);
     
     int nElec(0), nMuon(m_muonTree->muon_pt->size()), nLeptons(nMuon);
     for (int i(0); i < m_elecTree->electron_ismedium->size(); i++) {
       if (m_elecTree->electron_ismedium->at(i) == 1) {nElec++, nLeptons++;}
     } 
     
     if (nLeptons!=1) continue;
     
     bool cutMET   = (m_etTree->et_et > 20);
     bool cutMt    = (m_etTree->et_lepton_mt < 30);
     bool cutSumEt = (m_etTree->et_sumet < 400);

     if (cutMET && cutMt && cutSumEt && m_tauTree->tau_numTracks->size()>0) ++nEventMET;

     int nTau(0);
     for (int i(0); i < m_tauTree->tau_id_TauCutTight->size(); i++) {
       bool tauId = (m_doSafe ? (m_tauTree->tau_id_TauCutSafeCaloTight->at(i) == 1) : (m_tauTree->tau_id_TauCutTight->at(i) == 1));
       //if (m_tauTree->tau_id_TauCutTight->at(i) != 1) continue;
       if (!tauId) continue;
       
       nTau++;

       int nTrk = m_tauTree->tau_numTracks->at(i);

       m_hist_nprong->Fill(nTrk, weight);

       if (!cutMET) continue;   
       m_hist_nprong_met->Fill(nTrk, weight);
       if (!cutMt) continue;
       m_hist_nprong_mt->Fill(nTrk, weight);
       if (!cutSumEt) continue;
       m_hist_nprong_sumet->Fill(nTrk, weight);

       nTauMET++;
     }

     //std::cout << "Ntau (Event)" << ievent << " = " << nTau << std::endl;
   }
}

void ZtautauAnalysis::Tau(int nevents) 
{

  float n1lep(0), nMet(0), nMt(0), nSum(0), nAll(0);
  for (int ievent(0); ievent < nevents; ievent++) {
    m_etTree->LoadTree(ievent);
    m_etTree->GetEntry(ievent);
    
    m_weightTree->LoadTree(ievent);
    m_weightTree->b_weight_weight->GetEntry(ievent);
    float weight = m_weightTree->weight_weight;

    //float weight = m_etTree->et_weight;

    m_elecTree->LoadTree(ievent);
    m_elecTree->b_electron_ismedium->GetEntry(ievent);
    
    m_muonTree->LoadTree(ievent);
    m_muonTree->b_muon_pt->GetEntry(ievent);
    
    int nElec(0), nMuon(m_muonTree->muon_pt->size()), nLeptons(nMuon);
    for (int i(0); i < m_elecTree->electron_ismedium->size(); i++) {
      if (m_elecTree->electron_ismedium->at(i) == 1) {nElec++, nLeptons++;}
    } 
    
    if (nLeptons!=1) continue;
    
    n1lep += weight;

    bool cutMET   = (m_etTree->et_et > 20);
    bool cutMt    = (m_etTree->et_lepton_mt < 30);
    bool cutSumEt = (m_etTree->et_sumet < 400);
    
    if (cutMET) nMet += weight;
    if (cutMET && cutMt) nMt += weight;
    if (cutMET && cutMt && cutSumEt) nSum += weight;

    if (!(cutMET && cutMt && cutSumEt)) continue;
    nAll += weight; 

    m_tauTree->LoadTree(ievent);
    m_tauTree->b_tau_pt->GetEntry(ievent);
    m_tauTree->b_tau_isTauRec->GetEntry(ievent);
    m_tauTree->b_tau_isTau1P3P->GetEntry(ievent);
    m_tauTree->b_tau_taurec_emRadius->GetEntry(ievent);
    m_tauTree->b_tau_taurec_stripWidth2->GetEntry(ievent);
    m_tauTree->b_tau_taurec_isolationFraction->GetEntry(ievent);
    m_tauTree->b_tau_taurec_etEMCalib->GetEntry(ievent);
    m_tauTree->b_tau_taurec_etHadCalib->GetEntry(ievent);
    m_tauTree->b_tau_tau1p3p_emRadius->GetEntry(ievent);
    m_tauTree->b_tau_tau1p3p_stripWidth2->GetEntry(ievent);
    m_tauTree->b_tau_tau1p3p_isolationFraction->GetEntry(ievent);
    m_tauTree->b_tau_pt->GetEntry(ievent);

    for (int i(0); i < m_tauTree->tau_isTauRec->size(); i++) {

      float pt = m_tauTree->tau_pt->at(i);
      if (pt < 25 || pt >= 45) continue;

      bool isRec = (m_tauTree->tau_isTauRec->at(i)==1);
      bool is1P3P = (m_tauTree->tau_isTau1P3P->at(i)==1);

      float emRadius(0), isolFrac(0), width2(0), etEM(0), etHad(0);

      if (isRec) {
	emRadius = m_tauTree->tau_taurec_emRadius->at(i);
	width2   = m_tauTree->tau_taurec_stripWidth2->at(i);
	isolFrac = m_tauTree->tau_taurec_isolationFraction->at(i);
	etEM     = m_tauTree->tau_taurec_etEMCalib->at(i);
	etHad    = m_tauTree->tau_taurec_etHadCalib->at(i);
      } else if (is1P3P) {
	emRadius = m_tauTree->tau_tau1p3p_emRadius->at(i);
	width2   = m_tauTree->tau_tau1p3p_stripWidth2->at(i);
	isolFrac = m_tauTree->tau_tau1p3p_isolationFraction->at(i);	
      }

      m_hist_emradius->Fill(emRadius, weight);
      m_hist_width2->Fill(width2, weight);
      m_hist_isolfrac->Fill(isolFrac, weight);
      if (isRec)  m_hist_etratio->Fill(etEM/(etEM+etHad), weight);
    }
  }
  std::cout << "N Events after Et selction = " << n1lep << " " << nMet << " " << nMt << " " << nSum << " " << nAll << std::endl;
}

void ZtautauAnalysis::Ztautau(int nevents) 
{
  std::cout << "Name = " << m_ztautauTree->fChain->GetCurrentFile()->GetName() << std::endl;

  for (int ievent(0); ievent < nevents; ievent++) {
    m_ztautauTree->LoadTree(ievent);
    m_ztautauTree->GetEntry(ievent);

    float weight = m_ztautauTree->Ztautau_weight; 
    int   nCand  = m_ztautauTree->Ztautau_m->size();
    int nCandPass(0);

    ++ nEventCand;

    for (int icand(0); icand < nCand; icand++) {
      
      int icut(0);

      m_cutFlow[icut] += weight;
      m_cutFlowSq[icut++] += weight*weight;

      int nTrk = m_ztautauTree->Ztautau_element_numTracks->at(icand).at(tau);
      float dPhi = m_ztautauTree->Ztautau_visible_deltaphi->at(icand); 

      m_hist_visdphi->Fill(dPhi * 180/M_PI, weight);

      bool tauId = (m_doSafe ? (m_ztautauTree->Ztautau_element_tau_id_TauCutSafeCaloTight->at(icand).at(tau) == 1) : 
      		    (m_ztautauTree->Ztautau_element_tau_id_TauCutTight->at(icand).at(tau) == 1));

      if (tauId) {
        m_hist_nprong_cand->Fill(nTrk, weight);
	nTauCand++;
      }

      //bool tightElec = (abs(m_ztautauTree->Ztautau_element_pdgid->at(icand).at(lepton)) != 11 || m_ztautauTree->Ztautau_element_elec_isTight->at(icand).at(lepton));
      //if (!tightElec) continue;

      if (!(dPhi > 1. && dPhi < 3.1)) continue;
      m_cutFlow[icut] += weight;
      m_cutFlowSq[icut++] += weight*weight;

      int charge     = m_ztautauTree->Ztautau_charge->at(icand);
      float mass     = m_ztautauTree->Ztautau_visible_m->at(icand);
      float invmass  = m_ztautauTree->Ztautau_m->at(icand);
      int chargeprod = m_ztautauTree->Ztautau_element_charge->at(icand).at(tau) 
	             * m_ztautauTree->Ztautau_element_charge->at(icand).at(lepton);

      //bool tauId = (m_doSafe ? (m_ztautauTree->Ztautau_element_tau_id_TauCutSafeCaloTight->at(icand).at(tau) == 1) : 
      //		    (m_ztautauTree->Ztautau_element_tau_id_TauCutTight->at(icand).at(tau) == 1));
      bool tauId1or3 = (tauId && (m_ztautauTree->Ztautau_element_numTracks->at(icand).at(tau)==1 || 
      				  m_ztautauTree->Ztautau_element_numTracks->at(icand).at(tau) ==3));

//       bool tauId = (m_ztautauTree->Ztautau_element_tau_id_TauCutTight->at(icand).at(tau)==1 &&
// 		    //bool tauId = (m_ztautauTree->Ztautau_element_tau_id_TauCutSafeCaloTight->at(icand).at(tau)==1 &&
// 		    (m_ztautauTree->Ztautau_element_numTracks->at(icand).at(tau)==1 || 
// 		     m_ztautauTree->Ztautau_element_numTracks->at(icand).at(tau) ==3));            
	
      //if (m_ztautauTree->Ztautau_element_tau_id_TauCutTight->at(icand).at(tau)==1) {
      if (tauId) {
	m_hist_nprong_dphi->Fill(nTrk, weight);

	m_hist_charge_tau->Fill(m_ztautauTree->Ztautau_element_charge->at(icand).at(tau), weight); 
	m_hist_charge_lep->Fill(m_ztautauTree->Ztautau_element_charge->at(icand).at(lepton), weight); 
	m_hist_charge_prod->Fill(m_ztautauTree->Ztautau_element_charge->at(icand).at(tau) *
				 m_ztautauTree->Ztautau_element_charge->at(icand).at(lepton), weight); 

	//if (charge == 0) {
	if (chargeprod < 0) {
	  m_hist_nprong_chrg->Fill(nTrk, weight);
	  if (mass > 37 && mass < 75) {
	    m_hist_nprong_vism->Fill(nTrk, weight);
	    if (tauId1or3) {
	      m_hist_nprong_tau->Fill(nTrk, weight);
	    }
	  }
	}

	if (chargeprod > 0) {
	  m_hist_nprong_chrg_os->Fill(nTrk, weight);
	  if (mass > 37 && mass < 75) {
	    m_hist_nprong_vism_os->Fill(nTrk, weight);
	  }
	}       
	
      }

      //if (!tauId) continue;
      if (!tauId1or3) continue;
      m_cutFlow[icut] += weight;      
      m_cutFlowSq[icut++] += weight*weight;

      //if (charge == 0) {
      if (chargeprod < 0) {
	m_cutFlow[icut] += weight;
	m_cutFlowSq[icut++] += weight*weight;

	m_hist_vismass_os->Fill(mass, weight);
	m_hist_invmass_os->Fill(invmass, weight);

	m_hist_invpt->Fill(m_ztautauTree->Ztautau_pt->at(icand), weight);
	m_hist_invphi->Fill(m_ztautauTree->Ztautau_phi->at(icand)*180/M_PI, weight);
	m_hist_inveta->Fill(m_ztautauTree->Ztautau_eta->at(icand), weight);

	m_hist_vispt->Fill(m_ztautauTree->Ztautau_visible_pt->at(icand), weight);
	m_hist_visphi->Fill(m_ztautauTree->Ztautau_visible_phi->at(icand)*180/M_PI, weight);
	m_hist_viseta->Fill(m_ztautauTree->Ztautau_visible_eta->at(icand), weight);

	if (mass > 37 && mass < 75) {
	  m_cutFlow[icut] += weight;
	  m_cutFlowSq[icut++] += weight*weight;
	  
	  if ((m_name == "Ztautau" && (ievent==566 || ievent == 1119)) ||
	      (m_name == "Zee" && (ievent==502)) ||
	      (m_name == "Ttbar" && (ievent==1135))
	      ) {
	    std::cout << "Tau (Pt, Eta, Phi) = " 
		      << m_ztautauTree->Ztautau_element_pt->at(icand).at(tau)  << " " 
		      << m_ztautauTree->Ztautau_element_eta->at(icand).at(tau) << " "
		      << m_ztautauTree->Ztautau_element_phi->at(icand).at(tau) << " "
		      << m_ztautauTree->Ztautau_element_tau_id_ElectronVeto->at(icand).at(tau) << " "
		    << m_ztautauTree->Ztautau_element_tau_id_MuonVeto->at(icand).at(tau) << " "
		      << std::endl;
	  }
	}

	nCandPass++;
      } else icut += 2;

      //if (fabs(charge) == 2) {
      if (chargeprod > 0) {
	m_cutFlow[icut] += weight;
	m_cutFlowSq[icut++] += weight*weight;

	m_hist_vismass_ss->Fill(mass, weight);
	m_hist_invmass_ss->Fill(invmass, weight);

	if (mass > 37 && mass < 75) {
	  m_cutFlow[icut] += weight;
	  m_cutFlowSq[icut++] += weight*weight;
	}
      } else icut+=2;

    }
    
    if (nCandPass>1) {
      std::cout << "N Candidates passed = " << nCandPass << " for event "  << ievent << std::endl;
    }
  }

  std::cout << "NEvents = " << nEventMET << ", " << nEventCand << std::endl;
  std::cout << "NTaus = " << nTauMET << ", " << nTauCand << std::endl;
  CutFlow();
}

void ZtautauAnalysis::Finalise() 
{
  ScaleHists();
  WriteHists();
}

void ZtautauAnalysis::CutFlow() 
{
  // Cut flow
  std::cout << "==========================================================================================================================" << std::endl; 
  std::cout << " Dumping cuts for " << m_name << ", Scale = " << m_lumiScale << std::endl;
  std::cout << "==========================================================================================================================" << std::endl; 
  std::cout << "--------------------------------------------------------------------------------------------------------------------------" << std::endl;
  std::cout << std::setw(50) << "Cut" << " | " << std::endl; 
  std::cout << "--------------------------------------------------------------------------------------------------------------------------" << std::endl;
  std::vector<float>::const_iterator cutItr(m_cutFlow.begin());
  std::vector<float>::const_iterator cutEnd(m_cutFlow.end());
  for (int icut(0); cutItr != cutEnd; ++cutItr, icut++) {
    std::cout << std::setw(50) << m_cutNames.at(icut) << " | " 
	      << std::setw(15) << (*cutItr) << " +- " << std::setw(15) << sqrt(m_cutFlowSq.at(icut)) << " | " 
	      << std::setw(15) << (*cutItr) * m_lumiScale  << " +- " << std::setw(15) << m_lumiScale * sqrt(m_cutFlowSq.at(icut)) << " | "
	      << std::endl;

    m_hist_cuts->SetBinContent(icut+1, *cutItr);
    m_hist_cuts->SetBinError(icut+1, sqrt(m_cutFlowSq.at(icut)));

    m_hist_cutsscaled->SetBinContent(icut+1, (*cutItr)); // * m_lumiScale);
    m_hist_cutsscaled->SetBinError(icut+1, sqrt(m_cutFlowSq.at(icut)));
   std::cout << "--------------------------------------------------------------------------------------------------------------------------" << std::endl;
  }

  m_hist_cutsscaled->Scale(m_lumiScale);

  std::cout << "==========================================================================================================================" << std::endl; 

}

void ZtautauAnalysis::ScaleHists()
{
  m_hist_vismass_os->Scale(m_lumiScale);  
  m_hist_vismass_ss->Scale(m_lumiScale);  
  m_hist_invmass_os->Scale(m_lumiScale);  
  m_hist_invmass_ss->Scale(m_lumiScale);  

  m_hist_invpt->Scale(m_lumiScale); 
  m_hist_invphi->Scale(m_lumiScale);
  m_hist_inveta->Scale(m_lumiScale);
  m_hist_vispt->Scale(m_lumiScale); 
  m_hist_visphi->Scale(m_lumiScale);
  m_hist_viseta->Scale(m_lumiScale);
  m_hist_visdphi->Scale(m_lumiScale);
  
  m_hist_etmiss->Scale(m_lumiScale);      
  m_hist_etsum->Scale(m_lumiScale);       
  m_hist_mt->Scale(m_lumiScale);          
  m_hist_mt_etmiss->Scale(m_lumiScale);          
  m_hist_etmiss_mt->Scale(m_lumiScale);   
  m_hist_etsum_mtmet->Scale(m_lumiScale);              

  m_hist_nprong->Scale(m_lumiScale);       
  m_hist_nprong_met->Scale(m_lumiScale);   
  m_hist_nprong_mt->Scale(m_lumiScale);    
  m_hist_nprong_sumet->Scale(m_lumiScale); 
  m_hist_nprong_dphi->Scale(m_lumiScale);
  m_hist_nprong_cand->Scale(m_lumiScale);
  m_hist_nprong_tau->Scale(m_lumiScale);  
  m_hist_nprong_chrg->Scale(m_lumiScale); 
  m_hist_nprong_vism->Scale(m_lumiScale); 
  m_hist_nprong_chrg_os->Scale(m_lumiScale); 
  m_hist_nprong_vism_os->Scale(m_lumiScale); 

  m_hist_charge_tau->Scale(m_lumiScale); 
  m_hist_charge_lep->Scale(m_lumiScale); 
  m_hist_charge_prod->Scale(m_lumiScale); 

  m_hist_emradius->Scale(m_lumiScale);
  m_hist_width2->Scale(m_lumiScale);  
  m_hist_isolfrac->Scale(m_lumiScale);
  m_hist_etratio->Scale(m_lumiScale);   
} 

void ZtautauAnalysis::WriteHists() 
{
  TFile* histFile = new TFile(TString(m_name) +"Hists.root", "recreate");

  m_hist_vismass_os->Write();          delete  m_hist_vismass_os;
  m_hist_vismass_ss->Write();          delete  m_hist_vismass_ss;  
  m_hist_invmass_os->Write();          delete  m_hist_invmass_os;
  m_hist_invmass_ss->Write();          delete  m_hist_invmass_ss;  

  m_hist_invpt->Write();               delete m_hist_invpt; 
  m_hist_invphi->Write();	       delete m_hist_invphi;
  m_hist_inveta->Write();	       delete m_hist_inveta;
  m_hist_vispt->Write(); 	       delete m_hist_vispt; 
  m_hist_visphi->Write();	       delete m_hist_visphi;
  m_hist_viseta->Write();	       delete m_hist_viseta;
  m_hist_visdphi->Write();	       delete m_hist_visdphi;

  m_hist_etmiss->Write();              delete m_hist_etmiss;
  m_hist_etsum->Write();               delete m_hist_etsum;
  m_hist_mt->Write();                  delete m_hist_mt;
  m_hist_mt_etmiss->Write();           delete m_hist_mt_etmiss;
  m_hist_etmiss_mt->Write();           delete m_hist_etmiss_mt;
  m_hist_etsum_mtmet->Write();         delete m_hist_etsum_mtmet;

  m_hist_nprong->Write();              delete m_hist_nprong;           
  m_hist_nprong_met->Write();          delete m_hist_nprong_met;   
  m_hist_nprong_mt->Write();           delete m_hist_nprong_mt;    
  m_hist_nprong_sumet->Write();        delete m_hist_nprong_sumet; 
  m_hist_nprong_cand->Write();         delete m_hist_nprong_cand;
  m_hist_nprong_dphi->Write();         delete m_hist_nprong_dphi;
  m_hist_nprong_tau->Write();          delete m_hist_nprong_tau; 
  m_hist_nprong_chrg->Write();         delete m_hist_nprong_chrg;
  m_hist_nprong_vism->Write();         delete m_hist_nprong_vism;
  m_hist_nprong_chrg_os->Write();      delete m_hist_nprong_chrg_os;
  m_hist_nprong_vism_os->Write();      delete m_hist_nprong_vism_os;

  m_hist_charge_tau->Write();          delete  m_hist_charge_tau;
  m_hist_charge_lep->Write();          delete  m_hist_charge_lep;
  m_hist_charge_prod->Write();         delete  m_hist_charge_prod;

  m_hist_emradius->Write();            delete m_hist_emradius;
  m_hist_width2->Write();  	       delete m_hist_width2;  
  m_hist_isolfrac->Write();	       delete m_hist_isolfrac;
  m_hist_etratio->Write(); 	       delete m_hist_etratio; 

  // Cuts
  m_hist_cuts->Write();                delete m_hist_cuts;
  m_hist_cutsscaled->Write();          delete m_hist_cutsscaled;
}

void ZtautauAnalysis::SetLumiScale() 
{

  TFile file(TString("/hepstore/store2/gwilliam/ZtautauNtuples/Ntuple.") + m_name +  ".All.root");

  TH1F* cutHist = dynamic_cast<TH1F*>(file.Get("Cuts/CutHist"));

  m_lumiScale = m_lumi * m_xsect / cutHist->GetBinContent(1);

  std::cout << "Scale = " << m_lumiScale << std::endl;

  delete cutHist;
}