// file:        VisMassCycle.cxx
// author:      Ryan Reece       <ryan.reece@cern.ch>
//              Stan Lai         <stan.lai@cern.ch>
//              Justin Griffiths <griffith@cern.ch>
// created:     February 2010
//------------------------------------------------------------------------------

// std
#include <cmath> // fabs, sqrt, cosh
#include <iomanip> // setw
#include <algorithm> // sort
#include <sstream> // stringstream
#include <fstream>
#include <utility> // pair

// ROOT
#include "TH1F.h"
#include "TParameter.h"
#include "TLorentzVector.h"
#include "TFile.h"

// Local
#include "../include/VisMassCycle.h"
#include "../include/GoodRunsListTool.h"
#include "../include/TriggerTool.h"
#include "../include/TrigNavTool.h"
#include "../include/TrigMatchTool.h"
#include "../include/VertexTool.h"
#include "../include/CleaningTool.h"
#include "../include/MuonSelectorTool.h"
#include "../include/ElectronSelectorTool.h"
#include "../include/TauSelectorTool.h"
#include "../include/JetSelectorTool.h"
#include "../include/OverlapRemovalTool.h"
#include "../include/MuonIsolationTool.h"
#include "../include/ElectronIsolationTool.h"
#include "../include/ElectronScaleFactorsTool.h"
#include "../include/ElectronEnergyRescalerTool.h"
#include "../include/MuonSmearingTool.h"
#include "../include/METTool.h"
#include "../include/Candidate.h"
#include "../include/WriteTool.h"
#include "../include/DuplicateFilterTool.h"
#include "../include/TrueLepHadFilterTool.h"
#include "../include/TauMETEnergyScaleCorrectionTool.h"

// AtlasSFrameUtils
#include "CycleMacros.h"
#include "Muon.h"
#include "Electron.h"
#include "Tau.h"
#include "Jet.h"


const float GeV = 1000.0;

ClassImp( VisMassCycle );

//______________________________________________________________________________
VisMassCycle::VisMassCycle()
{
    DeclareProperty("channel", c_channel = 2);
    DeclareProperty("tree_name", c_tree_name = "tauPerf");
    DeclareProperty("muon_container", c_muon_container = "mu_staco_");
    DeclareProperty("jet_container", c_jet_container = "jet_");
    DeclareProperty("control_regions", c_control_regions = "");
    DeclareProperty("per_obj_trigger", c_per_obj_trigger = false );
    DeclareProperty("apply_electron_sf", c_apply_electron_sf = false );
    DeclareProperty("apply_muon_sf", c_apply_muon_sf = false );
    DeclareProperty("sys_sigma_muon_sf", c_sys_sigma_muon_sf = 0);
    DeclareProperty("use_met_locHadTopo_forEl", c_use_met_locHadTopo_forEl = false );
    DeclareProperty("use_met_locHadTopo_forMu", c_use_met_locHadTopo_forMu = false );
    DeclareProperty( "apply_jetCleaning_forMC", c_apply_jetCleaning_forMC = true );
    DeclareProperty("apply_tau_jet_electron_sf", c_apply_tau_jet_electron_sf = false );
    DeclareProperty("apply_tau_jet_electron_sf_sigma", c_apply_tau_jet_electron_sf_sigma = 0 );

    DECLARE_TOOL(GoodRunsListTool, "GoodRunsList");
    DECLARE_TOOL(TriggerTool, "ElectronTrigger");
    DECLARE_TOOL(TriggerTool, "MuonTrigger");
    DECLARE_TOOL(TriggerTool, "TauTrigger");
    DECLARE_TOOL(TrigNavTool, "TrigNav");
    DECLARE_TOOL(TrigMatchTool, "TrigMatch");
    DECLARE_TOOL(VertexTool, "Vertex");
    DECLARE_TOOL(CleaningTool, "Cleaning");
    DECLARE_TOOL(MuonSelectorTool, "MuonPreselector");
    DECLARE_TOOL(ElectronSelectorTool, "ElectronPreselector");
    DECLARE_TOOL(TauSelectorTool, "TauPreselector");
    DECLARE_TOOL(JetSelectorTool, "JetPreselector");
    DECLARE_TOOL(MuonSelectorTool, "MuonSelector");
    DECLARE_TOOL(ElectronSelectorTool, "ElectronSelector");
    DECLARE_TOOL(TauSelectorTool, "TauSelector");
    DECLARE_TOOL(TauSelectorTool, "LooseTauSelector");
    DECLARE_TOOL(JetSelectorTool, "JetSelector");
    DECLARE_TOOL(OverlapRemovalTool, "OverlapRemoval");
    DECLARE_TOOL(MuonIsolationTool, "MuonIsolation");
    DECLARE_TOOL(ElectronIsolationTool, "ElectronIsolation");
    DECLARE_TOOL(ElectronScaleFactorsTool, "ElectronScaleFactors");
    DECLARE_TOOL(ElectronEnergyRescalerTool, "ElectronEnergyRescaler");
    DECLARE_TOOL(MuonSmearingTool, "MuonSmearing");
    DECLARE_TOOL(METTool, "METTool");
    DECLARE_TOOL(WriteTool, "MiniTree");
    DECLARE_TOOL(DuplicateFilterTool, "DuplicateFilter");
    DECLARE_TOOL(TrueLepHadFilterTool, "TruthFilter");
    DECLARE_TOOL(TauMETEnergyScaleCorrectionTool, "TauMETReScaler");

}


//______________________________________________________________________________
VisMassCycle::~VisMassCycle()
{ m_logger << DEBUG  << "Destroying Cycle "                         << SLogger::endmsg; }


//______________________________________________________________________________
void VisMassCycle::BeginCycleImp() throw( SError )
{m_logger << DEBUG  << "BeginCycleImp "                         << SLogger::endmsg; }


//______________________________________________________________________________

void VisMassCycle::BeginMasterInputDataImp( const SInputData& id) throw( SError )
{m_logger << DEBUG  << "BeginMasterInputDataImp "                         << SLogger::endmsg; }


//______________________________________________________________________________
void VisMassCycle::BeginInputDataImp( const SInputData& id ) throw( SError )
{
    parse_control_regions();
    m_event_d3pdobject = new D3PDReader::EventInfoD3PDObject( &m_entry_number );
    TriggerTool* trigger_decisions = 0; 
    if(c_channel==1) 
      trigger_decisions = dynamic_cast<TriggerTool*> (GetTool("ElectronTrigger"));
    if(c_channel==2)
      trigger_decisions = dynamic_cast<TriggerTool*> (GetTool("MuonTrigger"));
    if(trigger_decisions == 0) {
      m_logger << FATAL << "No TriggerTool Found" << SLogger::endmsg;
      throw SError( SError::StopExecution );
    }

    m_trigger = new Trigger( m_entry_number, trigger_decisions->getCurrentTrigger().c_str(),  D3PDReader::TriggerD3PDObject( m_entry_number ) );
    m_trigdecision_d3pdobject = new D3PDReader::TrigDecisionD3PDObject( m_entry_number );
    std::string trigelstr("trig_EF_el_EF_e15_medium");//make this configurable? does it matter
    m_trigefel = new ElectronROIObject( m_entry_number, trigelstr.c_str(), D3PDReader::TrigEFelD3PDObject( m_entry_number ) );
    m_trigefmu = new D3PDReader::TrigEFmuD3PDObject( m_entry_number );
    m_trigeftau= new D3PDReader::TrigEFtauD3PDObject( m_entry_number );
    m_vxp_d3pdobject = new D3PDReader::VxpD3PDObject( m_entry_number );
    m_muon_d3pdobject = new D3PDReader::MuonD3PDObject( m_entry_number, c_muon_container.c_str() );
    m_electron_d3pdobject = new D3PDReader::ElectronD3PDObject( m_entry_number, "el_" );
    m_tau_d3pdobject = new D3PDReader::TauD3PDObject( m_entry_number, "tau_" );
    m_jet_d3pdobject = new D3PDReader::JetD3PDObject( m_entry_number, c_jet_container.c_str() );
    m_met_d3pdobject = new D3PDReader::METD3PDObject( m_entry_number );
    m_truth_d3pdobject = new D3PDReader::TruthD3PDObject( m_entry_number );
    m_truetau_d3pdobject = new D3PDReader::TrueTauD3PDObject( m_entry_number );
    m_mcevent_d3pdobject = new D3PDReader::MCEventD3PDObject( m_entry_number );
    m_cluster_d3pdobject = new D3PDReader::ClusterD3PDObject( m_entry_number );

    GET_TOOL(write_tool, WriteTool, "MiniTree");
    if(write_tool->IsOn())
    {
        gROOT->ProcessLine("#include <vector>"); //needed for root528a
        TTree* out_tree = GetOutputTree("VisMassCycleTree");
        m_event_d3pdobject->SetActive();
        m_event_d3pdobject->detmask0.SetActive(0);
        m_event_d3pdobject->detmask1.SetActive(0);
        m_event_d3pdobject->WriteTo(out_tree);
    }

    TString version = id.GetVersion();
    version.ToLower();
    if(!(version.Contains("zee") ||
	 //	 version.Contains("ttbar") ||
	 (version.Contains("dy") && version.Contains("ee")))
       ) c_apply_tau_jet_electron_sf = false;
    m_logger << INFO << "apply_tau_jet_electron_sf " << c_apply_tau_jet_electron_sf << " Sample: " << id.GetVersion() << SLogger::endmsg;

}


//______________________________________________________________________________
void VisMassCycle::BeginInputFileImp( const SInputData& id) throw( SError ) 
{ 
    
    m_logger << DEBUG  << "BeginInputFileImp "                         << SLogger::endmsg;
    m_event_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    //    m_trigger_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_trigdecision_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
//    m_trigger->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_trigger->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_trigefel->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_trigefmu->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_trigeftau->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_vxp_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_muon_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_electron_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_tau_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_jet_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_met_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    //m_collcand_d3pdobject->ReadFrom( GetInputTree("tauPerf") );
    m_truth_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_truetau_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_mcevent_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
    m_cluster_d3pdobject->ReadFrom( GetInputTree(c_tree_name.c_str()) );
				    
    m_current_file = GetInputTree(c_tree_name.c_str())->GetDirectory()->GetFile()->GetName();
}

//______________________________________________________________________________
void VisMassCycle::ExecuteEventImp( const SInputData& id, Double_t weight ) throw( SError )
{
    if(!is_data())
    {
        std::map< Long64_t, std::pair<std::set<std::string>, int> >::iterator map_itr = m_event_map.find((*m_event_d3pdobject).EventNumber());
	    if( map_itr != m_event_map.end() )
        {
            (*map_itr).second.second++;
	    (*map_itr).second.first.insert(m_current_file);
            throw SError(SError::SkipEvent);
        }
    	else
        {
	  std::set<std::string> file_names;
	  file_names.insert(m_current_file);
	        m_event_map.insert( std::make_pair( (*m_event_d3pdobject).EventNumber() , std::make_pair( file_names, 1 ) ) );
        }
    }
    GET_TOOL(duplicate_filter, DuplicateFilterTool, "DuplicateFilter");
    duplicate_filter->FilterDuplicate((*m_event_d3pdobject).EventNumber());    
 
    GET_TOOL(truth_filter, TrueLepHadFilterTool, "TruthFilter");
    bool true_lephad_event = true;
    if(truth_filter->GetChannel()!=TrueLepHadFilterTool::NO_FILTER && !is_data()){
      //Only build these heavy objects if requesting filter
      vector<TruthParticle> all_truth; 
      ParticleElementBuilder::build(all_truth, *(m_truth_d3pdobject));
      vector<ParticleElement<D3PDReader::TrueTauD3PDObjectElement> > all_true_taus;
      //      ParticleElementBuilder::build(all_true_taus, (*m_truetau_d3pdobject));
      true_lephad_event = truth_filter->FilterEvent(all_truth, all_true_taus);
    }
    GET_TOOL(write_tool, WriteTool, "MiniTree");
    write_tool->Write("true_lephad_event", true_lephad_event);
    //==========================================================================
    // 1 Initialize event info
    //==========================================================================
    
    float hist_weight = 1.0;
    const int RunNumber = is_data() ? (*m_event_d3pdobject).RunNumber() : 167521;
    const int lbn = is_data() ? (*m_event_d3pdobject).lbn() : 1;

    //--------------------------------------------------------------------------
    //Apply Vertex Weight before initial cut
    //--------------------------------------------------------------------------
    GET_TOOL(vertex_tool, VertexTool, "Vertex");
    float vxp_weight = vertex_tool->vertex_weight( *m_vxp_d3pdobject );
    hist_weight *= vxp_weight;
    write_tool->Write("vxp_weight", vxp_weight);

    //--------------------------------------------------------------------------
    //Apply mcatnlo event weight
    //--------------------------------------------------------------------------
    
    if(!is_data())
    {
        const float mcevt_weight = (*m_mcevent_d3pdobject)[0].weight();
	write_tool->Write("mcevt_weight", mcevt_weight);
	if(mcevt_weight == 0.0)
	{
	    m_logger << DEBUG << "mcevt_weight = 0, using 1 (not affecting event weight)." << SLogger::endmsg;
	}
	else
	{
	    if(fabs(mcevt_weight) != 1.0)
	    {
	        m_logger << WARNING << "fabs(mcevt_weight) != 1.0.  mcevt_weight = " << mcevt_weight << SLogger::endmsg;
	    }	    
	    hist_weight *= mcevt_weight;
	}
    }

    //--------------------------------------------------------------------------
    //First Cut(is weighted)
    //--------------------------------------------------------------------------

    CutFlow("ALL", true, true, hist_weight);
    Book("h_RunNumber_all", ";RunNumber;entries", 50000, 150000-0.5, 200000-0.5)->Fill(RunNumber);

    //==========================================================================
    // 2 Event preselection
    //==========================================================================

    //--------------------------------------------------------------------------
    // 2.1 Check GRL
    //--------------------------------------------------------------------------
 
    GET_TOOL(grl, GoodRunsListTool, "GoodRunsList");
    const bool passed_grl = grl->Select(*m_event_d3pdobject);

    //    CutFlow("GRL", passed_grl, true, hist_weight); // <--- default
//    CutFlow("GRL", passed_grl || true, true, hist_weight);
    CutFlow("GRL", passed_grl, false, hist_weight); // <--- check the GRL but do not reject the event

    //--------------------------------------------------------------------------
    // 2.2 Check for trigger
    //--------------------------------------------------------------------------
    GET_TOOL(trigger_decisions, TriggerTool, "ElectronTrigger"); //reads (EF)trig decision from chainId vecs/trigger metadata maps 
    bool passed_trigger = false;

    if (c_channel == 1)
    {
        m_logger << DEBUG << "electron trigger set" << SLogger::endmsg;
    }
    else if(c_channel == 2)
    {
        //Macro will not keep the correct triggerTool in scope
        trigger_decisions = dynamic_cast<TriggerTool*> (GetTool("MuonTrigger"));
        m_logger<<DEBUG<<"muon trigger set"<<SLogger::endmsg;
    }
    else 
    {
        m_logger<<WARNING<<"Trigger set to fail. Invalid channel: "<<c_channel<<"  allowed channels: 1=electron channel, 2=muon channel"<<SLogger::endmsg;
    }

    if( trigger_decisions->goodMetadata() )
    {
        passed_trigger = trigger_decisions->SelectTrigger(*m_event_d3pdobject, *m_trigdecision_d3pdobject);
    }
    else
    {
        passed_trigger = trigger_decisions->SelectTrigger( *m_event_d3pdobject, *m_trigger );
    }
    //    Cut("EF_mu40_MSonly", trigger_decisions->mu40_MSonly(), false, hist_weight);
    trigger_decisions->PlotCurrentTrigger( hist_weight, "all" );
    CutFlow("TRIGGER", passed_trigger, true, hist_weight);
    trigger_decisions->PlotCurrentTrigger( hist_weight, "passed" );

    //--------------------------------------------------------------------------
    // 2.3 Check for vertex
    //--------------------------------------------------------------------------
    D3PDReader::VxpD3PDObject selected_vxp_d3pdobject =  D3PDReader::VxpD3PDObject();
    bool passed_vertex = vertex_tool->Select(*m_vxp_d3pdobject, selected_vxp_d3pdobject);
    vertex_tool->PlotVertexInfo(*m_vxp_d3pdobject, hist_weight);
    vertex_tool->PlotVertexProfile(*m_vxp_d3pdobject, RunNumber, lbn, hist_weight);
    CutFlow("VERTEX", passed_vertex, true, hist_weight);
    write_tool->Write("vxp_n_selected", selected_vxp_d3pdobject.n());
    // 2.4 Jet cleaning -- done after we have the vector of jet / tau objets

    //==========================================================================
    // 3 Read particles from the D3PD
    //==========================================================================

    //--------------------------------------------------------------------------
    // 3.1 Get muons
    //--------------------------------------------------------------------------
    // build
    std::vector<Muon> all_muons;
    //MyParticleElementBuilder::build( all_muons, (*m_muon_d3pdobject) );
    ParticleElementBuilder::build( all_muons, (*m_muon_d3pdobject) );
    std::sort(all_muons.begin(), all_muons.end(), Muon::pt_cmp);
     
    // plot
    GET_TOOL(muon_selector, MuonSelectorTool, "MuonSelector"); // use for making plots
    //corrected muon pt
    std::vector<float> mu_pt_corrected;
    GET_TOOL(muon_smearing, MuonSmearingTool, "MuonSmearing");
    //mu_pt_corrected is simply muon_pt if tool is switched off
    muon_smearing->ApplySmearing(all_muons, mu_pt_corrected, (*m_event_d3pdobject).EventNumber());
    muon_selector->RecalcTLV(all_muons, mu_pt_corrected);
    std::sort(all_muons.begin(), all_muons.end(), Muon::pt_cmp);

    muon_selector->PlotKinematics(all_muons, hist_weight, "all");
    muon_selector->PlotIdentification(all_muons, hist_weight, "all");

    //--------------------------------------------------------------------------
    // 3.2 Get electrons
    //--------------------------------------------------------------------------
    GET_TOOL(electron_selector, ElectronSelectorTool, "ElectronSelector"); // use for making plots
    // build
    std::vector<Electron> all_electrons;
    ParticleElementBuilder::build( all_electrons, (*m_electron_d3pdobject) );

    //corrected electron energy
    std::vector<float> el_E_corrected;    
    GET_TOOL(electron_energy_rescaler, ElectronEnergyRescalerTool, "ElectronEnergyRescaler"); 
    m_logger<<DEBUG<< "Event number: " << (*m_event_d3pdobject).EventNumber() << SLogger::endmsg;
    //If this tool is turned off el_E_corrected is cl_E
    electron_energy_rescaler->ApplyEnergyRescaler(all_electrons, el_E_corrected, (*m_event_d3pdobject).EventNumber());
    electron_selector->RecalcTLV(all_electrons, el_E_corrected);

    std::sort(all_electrons.begin(), all_electrons.end(), Electron::pt_cmp);

    // plot
    electron_selector->PlotKinematics(all_electrons, hist_weight, "all");
    electron_selector->PlotIdentification(all_electrons, hist_weight, "all");

    //--------------------------------------------------------------------------
    // 3.3 Get taus
    //--------------------------------------------------------------------------
    // build
    std::vector<Tau> all_taus;
    ParticleElementBuilder::build( all_taus, (*m_tau_d3pdobject) );    
    GET_TOOL(tau_met_rescaler, TauMETEnergyScaleCorrectionTool, "TauMETReScaler");
    tau_met_rescaler->scale_tau(all_taus);
    std::sort(all_taus.begin(), all_taus.end(), Tau::pt_cmp);

    // plot
    GET_TOOL(tau_selector, TauSelectorTool, "TauSelector"); // use for making plots
    tau_selector->PlotKinematics(all_taus, hist_weight, "all");
    tau_selector->PlotIdentification(all_taus, hist_weight, "all");
    tau_selector->PlotIDVariables(all_taus, hist_weight, "all");

    //--------------------------------------------------------------------------
    // 3.4 Get jets
    //--------------------------------------------------------------------------
    // build
    std::vector<Jet> all_jets;
    ParticleElementBuilder::build( all_jets, (*m_jet_d3pdobject) );
    std::sort(all_jets.begin(), all_jets.end(), Jet::pt_cmp);

    // plot
    GET_TOOL(jet_selector, JetSelectorTool, "JetSelector"); // use for making plots
    jet_selector->PlotKinematics(all_jets, hist_weight, "all");
    jet_selector->PlotIdentification(all_jets, hist_weight, "all");
    jet_selector->PlotMultiplicity(all_jets, hist_weight, "all");
    


    //==========================================================================
    // 4 Object preselection
    //==========================================================================

    //--------------------------------------------------------------------------
    // 4.1 Muon preselection
    //--------------------------------------------------------------------------
    GET_TOOL(muon_preselector, MuonSelectorTool, "MuonPreselector");
    std::vector<Muon> preselected_muons;
    muon_preselector->Select(all_muons, preselected_muons);
    muon_selector->PlotKinematics(preselected_muons, hist_weight, "preselected");
    muon_selector->PlotIdentification(preselected_muons, hist_weight, "preselected");

    //--------------------------------------------------------------------------
    // 4.2 Electron preselection
    //--------------------------------------------------------------------------
    GET_TOOL(electron_preselector, ElectronSelectorTool, "ElectronPreselector");
    std::vector<Electron> preselected_electrons;
    electron_preselector->Select(all_electrons, preselected_electrons, RunNumber);
    electron_selector->PlotKinematics(preselected_electrons, hist_weight, "preselected");
    electron_selector->PlotIdentification(preselected_electrons, hist_weight, "preselected");

    //--------------------------------------------------------------------------
    // 4.3 Tau preselection
    //--------------------------------------------------------------------------
    GET_TOOL(tau_preselector, TauSelectorTool, "TauPreselector");
    std::vector<Tau> preselected_taus;
    tau_preselector->Select(all_taus, preselected_taus);
    tau_selector->PlotKinematics(preselected_taus, hist_weight, "preselected");
    tau_selector->PlotIdentification(preselected_taus, hist_weight, "preselected");
    tau_selector->PlotIDVariables(preselected_taus, hist_weight, "preselected");
    
    //--------------------------------------------------------------------------
    // 4.4 Jet preselection
    //--------------------------------------------------------------------------
    GET_TOOL(jet_preselector, JetSelectorTool, "JetPreselector");
    std::vector<Jet> preselected_jets;
    jet_preselector->Select(all_jets, preselected_jets);
    jet_selector->PlotKinematics(preselected_jets, hist_weight, "preselected");
    jet_selector->PlotIdentification(preselected_jets, hist_weight, "preselected");
    jet_selector->PlotMultiplicity(preselected_jets, hist_weight, "preselected");

    int n_preselected_channel_leptons = 0;
    switch(c_channel)
    {
        case 1: n_preselected_channel_leptons = preselected_electrons.size(); break;
        case 2: n_preselected_channel_leptons = preselected_muons.size(); break;
        default:
                m_logger << FATAL << "Invalid c_channel = " << c_channel << SLogger::endmsg;
                throw SError( SError::StopExecution );
    }
    CutFlow("PRESELECTED_LEPTON", n_preselected_channel_leptons > 0, true, hist_weight);
    vertex_tool->PlotVertexInfo(*m_vxp_d3pdobject, hist_weight, "PRESELECTED_LEPTON");

    CutFlow("PRESELECTED_TAU", preselected_taus.size(), true, hist_weight);

    //==========================================================================
    // 5 Overlap removal
    //==========================================================================
    GET_TOOL(overlap_remover, OverlapRemovalTool, "OverlapRemoval");
    std::vector<Muon> olr_muons(preselected_muons);
    std::vector<Electron> olr_electrons(preselected_electrons);
    std::vector<Tau> olr_taus(preselected_taus);
    std::vector<Jet> olr_jets(preselected_jets);
    std::vector<Jet> olr_jets_for_cleaning(all_jets);
    std::vector<Tau> olr_taus_for_cleaning;
    overlap_remover->DoOverlapRemoval(olr_muons, 0.2);
    overlap_remover->DoOverlapRemoval(olr_muons, olr_electrons, 0.2);
    overlap_remover->DoOverlapRemoval(olr_electrons, 0.2);
    overlap_remover->DoOverlapRemoval(olr_muons, olr_taus, 0.4);
    overlap_remover->DoOverlapRemoval(olr_electrons, olr_taus, 0.4);
    overlap_remover->DoOverlapRemoval(olr_taus, 0.4);
    overlap_remover->DoOverlapRemoval(olr_taus, olr_jets, 0.4);
    overlap_remover->DoOverlapRemoval(olr_electrons,olr_jets_for_cleaning, 0.4);
    overlap_remover->DoOverlapRemoval(olr_muons,olr_jets_for_cleaning, 0.4);

    const int n_olr_leptons = olr_electrons.size()*(c_channel==1) + olr_muons.size() * (c_channel==2);
    CutFlow("OVERLAP_REMOVED_LEPTON", n_olr_leptons > 0, true, hist_weight);

    //--------------------------------------------------------------------------
    // 2.4 Jet cleaning
    //--------------------------------------------------------------------------
    if ((is_data())||(c_apply_jetCleaning_forMC))
    {
        GET_TOOL(cleaning, CleaningTool, "Cleaning");
        bool passed_jet_cleaning = !(cleaning->Select(olr_jets_for_cleaning,olr_taus_for_cleaning)) ;
        CutFlow("JET_CLEANING", passed_jet_cleaning, true, hist_weight);
    }
    else
    {
        CutFlow("JET_CLEANING", true, true, hist_weight);
    }
    vertex_tool->PlotVertexInfo(*m_vxp_d3pdobject, hist_weight, "JET_CLEANING");


    //==========================================================================
    // 6 Trigger-object matching
    //==========================================================================
    std::vector<Electron> trigmatch_electrons( olr_electrons );
    std::vector<Muon>     trigmatch_muons( olr_muons );
    if( c_per_obj_trigger )
    {
        //GET_TOOL(trigger_decisions, TriggerTool, "Trigger"); //reads (EF)trig decision from chainId vecs/trigger metadata maps
        GET_TOOL(trigger_nav, TrigNavTool, "TrigNav");         //reads per event trig nav maps, chain->obj
        GET_TOOL(trigger_match, TrigMatchTool, "TrigMatch");   //matches trig obj -> offline obj
        if( c_channel==1 )
        {
            std::string trigelstr("trig_EF_el_");
            trigelstr+=trigger_decisions->getCurrentTrigger();
            //m_logger<<DEBUG<<"TrigMatch Cut: Setting current trigger to "<<trigelstr<<SLogger::endmsg;
            m_trigefel->CurrentTriggerROI.SetName( trigelstr.c_str() ); //set internal var of m_trigefel to correct trigger 
            std::vector<int> trigefel_ptrs = trigger_nav->GetTrigObjIndxs( *m_trigefel );
            trigger_match->Select( trigmatch_electrons, trigefel_ptrs, *m_trigefel );
        }
        else if( c_channel==2 )
        {
            //m_logger<<DEBUG<<"TrigMatch Cut: Setting current trig chain id to "<<trigger_decisions->getCurrentChainId( *m_trigdecision_d3pdobject )<<" for current trigger "<<trigger_decisions->getCurrentTrigger()<<SLogger::endmsg;
            int trig_chainId = trigger_decisions->getCurrentChainId( *m_trigdecision_d3pdobject );
            std::vector<int> trigefmu_ptrs = trigger_nav->GetTrigObjIndxs( trig_chainId, "mu", *m_trigdecision_d3pdobject );
            trigger_match->Select( trigmatch_muons, trigefmu_ptrs, *m_trigefmu );
        } 
    }

    int n_trigger_matched_leptons = 0;
    switch(c_channel)
    {
        case 1: n_trigger_matched_leptons = trigmatch_electrons.size(); break;
        case 2: n_trigger_matched_leptons = trigmatch_muons.size(); break;
        default:
                m_logger << WARNING << "Invalid c_channel = " << c_channel << "no per object trigger cut was applied"<< SLogger::endmsg;
    }
    CutFlow("TRIGGER_MATCHED_LEPTON", n_trigger_matched_leptons>0, true, hist_weight);

    const int n_olr_taus = olr_taus.size();
    CutFlow("OVERLAP_REMOVED_TAU", n_olr_taus > 0, true, hist_weight);
    write_tool->Write("n_olr_taus", n_olr_taus);

    //==========================================================================
    // 7 Object selection
    //==========================================================================

    //--------------------------------------------------------------------------
    // 7.1 Muon selection
    //--------------------------------------------------------------------------
    std::vector<Muon> selected_muons;
    muon_selector->Select(trigmatch_muons, selected_muons);
    
    muon_selector->PlotKinematics(selected_muons, hist_weight, "selected");
    muon_selector->PlotIdentification(selected_muons, hist_weight, "selected");

    //--------------------------------------------------------------------------
    // 7.2 Electron selection
    //--------------------------------------------------------------------------
    std::vector<Electron> selected_electrons;
    electron_selector->Select(trigmatch_electrons, selected_electrons, RunNumber);
    electron_selector->PlotKinematics(selected_electrons, hist_weight, "selected");
    electron_selector->PlotIdentification(selected_electrons, hist_weight, "selected");

    //--------------------------------------------------------------------------
    // 7.3 Tau selection
    //--------------------------------------------------------------------------
    std::vector<Tau> selected_taus;
    tau_selector->Select(olr_taus, selected_taus);
    tau_selector->PlotKinematics(selected_taus, hist_weight, "selected");
    tau_selector->PlotIdentification(selected_taus, hist_weight, "selected");
    tau_selector->PlotIDVariables(selected_taus, hist_weight, "selected");

    std::vector<Tau> loose_taus;
    GET_TOOL(loose_tau_selector, TauSelectorTool, "LooseTauSelector");
    loose_tau_selector->Select(olr_taus, loose_taus);
    tau_selector->PlotKinematics(loose_taus, hist_weight, "loose");
    tau_selector->PlotIdentification(loose_taus, hist_weight, "loose");
    tau_selector->PlotIDVariables(loose_taus, hist_weight, "loose");
    
    //--------------------------------------------------------------------------
    // 7.4 Jet selection
    //--------------------------------------------------------------------------
    std::vector<Jet> selected_jets;
    jet_selector->Select(olr_jets, selected_jets);
    jet_selector->PlotKinematics(selected_jets, hist_weight, "selected");
    jet_selector->PlotIdentification(selected_jets, hist_weight, "selected");
    jet_selector->PlotMultiplicity(selected_jets, hist_weight, "selected");

    int n_selected_channel_leptons = 0;
    switch(c_channel)
    {
        case 1: n_selected_channel_leptons = selected_electrons.size(); break;
        case 2: n_selected_channel_leptons = selected_muons.size(); break;
    }
    CutFlow("SELECTED_LEPTON", n_selected_channel_leptons > 0, true, hist_weight);

    //==========================================================================
    // 8 Lepton isolation
    //==========================================================================

    //--------------------------------------------------------------------------
    // 8.1 Muon isolation
    //--------------------------------------------------------------------------
    GET_TOOL(muon_isolation, MuonIsolationTool, "MuonIsolation");
    muon_isolation->PlotIsolation(selected_muons, hist_weight, "selected");
    muon_isolation->PlotIsolation2D(selected_muons, hist_weight, "selected");
    std::vector<Muon> isolated_muons;
    muon_isolation->Select(selected_muons, isolated_muons);
    muon_isolation->PlotIsolation(isolated_muons, hist_weight, "isolated");
    muon_selector->PlotKinematics(isolated_muons, hist_weight, "isolated");

    //--------------------------------------------------------------------------
    // 8.2 Electron isolation
    //--------------------------------------------------------------------------
    GET_TOOL(electron_isolation, ElectronIsolationTool, "ElectronIsolation");
    electron_isolation->PlotIsolation(selected_electrons, hist_weight, "selected");
    electron_isolation->PlotIsolation2D(selected_electrons, hist_weight, "selected");
    std::vector<Electron> isolated_electrons;
    electron_isolation->Select(selected_electrons, isolated_electrons);
    electron_isolation->PlotIsolation(isolated_electrons, hist_weight, "isolated");
    electron_selector->PlotKinematics(isolated_electrons, hist_weight, "isolated");


    int n_isolated_channel_leptons = 0;
    switch(c_channel)
    {
        case 1: n_isolated_channel_leptons = isolated_electrons.size(); break;
        case 2: n_isolated_channel_leptons = isolated_muons.size(); break;
    }

    
    // prepare for a possible application of the electron SF in MC
    const float hist_weight_orig = hist_weight;
    
    GET_TOOL(electron_sf, ElectronScaleFactorsTool, "ElectronScaleFactors");
    std::pair<float,float> m_electron_id_sf, m_electron_iso_sf, m_electron_id_sf_no_et_correction;

    if(c_channel == 1) // electron + tau channel
    {
        // apply electron scale factors in MC if required (efficiency + trigger SF)
        // m_electron_id_sf.first - SF, m_electron_id_sf.second - system. uncertainty on SF
      if ((!is_data())&&(c_apply_electron_sf))
        {
	    //All electron scale factors applied (trigger, ID, isolation)
            Electron* selected_el = 0;
	    if(isolated_electrons.size()) selected_el = &(isolated_electrons.front());
	    else selected_el = &(selected_electrons.front());
            m_electron_id_sf = electron_sf->Trigger_ID_ScaleFactors(*selected_el);
	    m_electron_id_sf_no_et_correction = electron_sf->Trigger_ID_ScaleFactors(*selected_el, false);
	    //Ensure garbage not set for isolation scale factors for non-iso leptons
	    m_electron_iso_sf.first = 1.0;
	    m_electron_iso_sf.second = 0.0;    
            if (isolated_electrons.size()) 
	      {
		m_electron_iso_sf = electron_sf->IsolationScaleFactor(*selected_el);
	      }
	    //Systematic variation applied if requested by ElectronScaleFactorsTool
	    hist_weight *= electron_sf->VaryScaleFactor(m_electron_id_sf, m_electron_iso_sf);
        }
      write_tool->Write("isol_scale_factor", m_electron_iso_sf.first);
      write_tool->Write("id_scale_factor", m_electron_id_sf.first);            
      write_tool->Write("id_scale_factor_no_et", m_electron_id_sf_no_et_correction.first);
    }
    else if(c_channel == 2) // electron + tau channel
    {
        // apply muon scale factors in MC if required
        if (!is_data() && c_apply_muon_sf )
        {
            const Muon* mu = 0;
	    if(isolated_muons.size()) mu = &(isolated_muons.front());
	    else mu = &(selected_muons.front());
	    TLorentzVector def_muon_tlv = mu->tlv;
	    def_muon_tlv.SetPtEtaPhiM( mu->pt(), mu->eta(), mu->phi(), 0);
            float trig_scale_factor = m_staco_trig_scale_fac.scaleFactor( mu->tlv, trigger_decisions->getCurrentTrigger() );
            float trig_scale_factor_stat_err = m_staco_trig_scale_fac.scaleFactorStatUncertainty( mu->tlv, trigger_decisions->getCurrentTrigger() );
            float trig_scale_factor_sys_err = m_staco_trig_scale_fac.scaleFactorSysUncertainty( mu->tlv, trigger_decisions->getCurrentTrigger() );
	    float reco_scale_factor = m_staco_cb_scale_fac.scaleFactor( mu->tlv );
	    float reco_scale_factor_err = m_staco_cb_scale_fac.scaleFactorUncertainty( mu->tlv );
	    float isol_scale_factor = 1.0;
	    float isol_scale_factor_stat_err = 0.0;
	    float isol_scale_factor_sys_err = 0.0;
	    if(isolated_muons.size()) {
	       isol_scale_factor = m_staco_isol_scale_fac.scaleFactor( def_muon_tlv );
	       isol_scale_factor_stat_err = m_staco_isol_scale_fac.scaleFactorStatUncertainty( def_muon_tlv );
	       isol_scale_factor_sys_err = m_staco_isol_scale_fac.scaleFactorSysUncertainty( def_muon_tlv );
	    }
	    write_tool->Write("isol_scale_factor", isol_scale_factor);
	    write_tool->Write("id_scale_factor", reco_scale_factor*trig_scale_factor);
	      
	    float scale_factor = trig_scale_factor*reco_scale_factor*isol_scale_factor *
	      ( 1 +
		c_sys_sigma_muon_sf * std::sqrt(  std::pow(trig_scale_factor_stat_err/trig_scale_factor,2) + 
						  std::pow(isol_scale_factor_stat_err/isol_scale_factor,2) +
						  std::pow(trig_scale_factor_sys_err/trig_scale_factor,2) +
						  std::pow(isol_scale_factor_sys_err/isol_scale_factor,2) +
						  std::pow(reco_scale_factor_err/reco_scale_factor,2) )
		);
	    hist_weight = hist_weight_orig*scale_factor;
	}
    }


    //==========================================================================
    // 9 Define MET
    //==========================================================================

    // TODO: make configurable
    GET_TOOL(met_tool, METTool, "METTool");
    TLorentzVector met_tlv = met_tool->DefineMET(m_met_d3pdobject,olr_muons , olr_electrons, electron_energy_rescaler, muon_smearing, c_channel);
    if(c_channel==1) {
      if(isolated_electrons.size()) tau_met_rescaler->scale_met(m_cluster_d3pdobject, met_tlv, &(isolated_electrons.front()));
      else if(selected_electrons.size()) tau_met_rescaler->scale_met(m_cluster_d3pdobject, met_tlv, &(selected_electrons.front()));
    }
    else {
      if(isolated_muons.size()) tau_met_rescaler->scale_met(m_cluster_d3pdobject, met_tlv, 0, &(isolated_muons.front()));
      else if(selected_muons.size()) tau_met_rescaler->scale_met(m_cluster_d3pdobject, met_tlv, 0, &(selected_muons.front()));
    }
    float met_x = met_tlv.Px();
    float met_y = met_tlv.Py();
    float met_sumet =  met_tool->sumet();
    // met_tool->PlotMETInfo(*m_vxp_d3pdobject, selected_vxp_d3pdobject, hist_weight , "test");

    CutFlow("ISOLATED_LEPTON", n_isolated_channel_leptons > 0, false, hist_weight);
    write_tool->Write("n_olr_leptons", n_olr_leptons);
    write_tool->Write("n_leptons_dlv", (int) (olr_muons.size() + olr_electrons.size()));
    write_tool->Write("n_preselected_electrons", (int) preselected_electrons.size());
    write_tool->Write("n_preselected_muons", (int) preselected_muons.size());
    write_tool->Write("n_selected_electrons", (int) selected_electrons.size());
    write_tool->Write("n_selected_muons", (int) selected_muons.size());
    write_tool->Write("n_isolated_electrons", (int) isolated_electrons.size());
    write_tool->Write("n_isolated_muons", (int) isolated_muons.size());

    const int n_selected_taus = selected_taus.size();
    CutFlow("SELECTED_TAU", n_selected_taus > 0, false, hist_weight);

    //==========================================================================
    // 10 Build and loop over candidates
    //==========================================================================

    //Choose Tau
    Tau* selected_tau = 0;
    if(selected_taus.size()) selected_tau = &(selected_taus.front());
    else if(loose_taus.size()) selected_tau = &(loose_taus.front());
    else if(olr_taus.size()) selected_tau = &(olr_taus.front());

    //--------------------------------------------------------------------------
    // 10.1 Build candidates
    //--------------------------------------------------------------------------

    std::vector<Candidate> candidates;
    if(c_channel == 1) // electron + tau channel
    {
        Electron* selected_electron = 0;
        if(isolated_electrons.size()) selected_electron = &(isolated_electrons.front());
        else if(selected_electrons.size()) selected_electron = &(selected_electrons.front());

        if( selected_electron && selected_tau ) 
        {
            candidates.push_back( Candidate( selected_electron, 0, selected_tau ) );
        }
    }
    else if(c_channel == 2) // muon + tau channel
    {
        Muon* selected_muon = 0;
        if(isolated_muons.size()) selected_muon = &(isolated_muons.front());
        else if(selected_muons.size()) selected_muon = &(selected_muons.front());

        if( selected_muon && selected_tau ) 
        {
            candidates.push_back( Candidate( 0, selected_muon, selected_tau ) );
        }

    }

    bool tau_is_jet = true;
    if(selected_tau && c_apply_tau_jet_electron_sf){
	//Truth Pariticles
      vector<TruthParticle> all_truth; 
      ParticleElementBuilder::build(all_truth, *(m_truth_d3pdobject));
      for(std::vector<TruthParticle>::const_iterator itr = all_truth.begin(); itr != all_truth.end(); ++itr){
	TruthParticle tp = (*itr);
	if(std::abs(tp.pdgId()) != 11) continue;
	if(selected_tau->tlv.DeltaR(tp.tlv) > 0.2) continue;
	TruthParticle mother = tp.GetMother();
	int tp_barcode = tp.barcode();
	while(tp_barcode != mother.barcode()){
	  tp_barcode = mother.barcode();
	  if(mother.pdgId() == 23 || std::abs(mother.pdgId()) == 24) {
	    tau_is_jet = false;
	    Book("h_tau_electron_match_mother",";;", 60, -30, 30)->Fill(mother.pdgId(), hist_weight);
	    break;
	  }
	  else if (std::abs(mother.pdgId()) == 11 || std::abs(mother.pdgId()) == 15 ) continue;
	  else {
	    //This electron is junk continue;
	    Book("h_tau_electron_match_mother_junk",";;", 60, -30, 30)->Fill(mother.pdgId(), hist_weight);
	    break;
	  }
	}//while
	if(tau_is_jet==0) break;
      }
    }      

    float tau_jet_electron_sf = 1.0;
    write_tool->Write("tau_is_jet", tau_is_jet);
    if( c_apply_tau_jet_electron_sf && !tau_is_jet && selected_tau && 
	selected_tau->numTrack()==1 && selected_tau->tauCutSafeMedium()){
      if(fabs(selected_tau->tlv.Eta()) < 1.37) {
	tau_jet_electron_sf = (1.23 + c_apply_tau_jet_electron_sf_sigma*0.36);
      }
      else if(fabs(selected_tau->tlv.Eta()) < 1.52) {
	tau_jet_electron_sf = (0.22 + c_apply_tau_jet_electron_sf_sigma*0.14);
      }
      else {
	tau_jet_electron_sf = (1.19 + c_apply_tau_jet_electron_sf_sigma*0.47); 
      }
    }
    write_tool->Write("tau_jet_electron_sf", tau_jet_electron_sf);		    
    hist_weight *= tau_jet_electron_sf;
    
    


    //--------------------------------------------------------------------------
    // 10.2 Loop over candidates
    //--------------------------------------------------------------------------

    //weight to be applied to the histograms entry (added because of electron SF)
    hist_weight = hist_weight_orig;

    for(std::vector<Candidate>::iterator cand_it = candidates.begin(); cand_it != candidates.end(); ++cand_it)
    {
        Tau& tau = *(cand_it->tau);
        TLorentzVector tau_tlv = cand_it->tau->tlv;
        TLorentzVector lepton_tlv;
        std::string lepton_name;
        if(c_channel == 1)
        {
            lepton_tlv = cand_it->electron->tlv;
            float el_eta, el_phi;
            if ((cand_it->electron->nPixHits()+cand_it->electron->nSCTHits())<4)
            {
                el_eta = cand_it->electron->cl_eta();
                el_phi = cand_it->electron->cl_phi();
            }
            else 
            {
                el_eta = cand_it->electron->tracketa();
                el_phi = cand_it->electron->trackphi();
            }
            //call Pt from tlv where the rescaled energy is saved
            const float el_pt = cand_it->electron->tlv.Pt();
            m_logger<<DEBUG<<"Candidate Pt " << el_pt/GeV<< SLogger::endmsg;
            lepton_tlv.SetPtEtaPhiM(el_pt, el_eta, el_phi, 0.);
            lepton_name = "e";

            //Apply electron scale factors in MC if required in the histograms
            //m_electron_sf.first - SF, m_electron_sf.second - system. uncertainty on SF
            if ((!is_data())&&(c_apply_electron_sf))
            {
                m_electron_id_sf.first = 0;
                m_electron_id_sf.second = 0;
		//Ensure non-isolated leptons are scaled to 1
                m_electron_iso_sf.first = 1.0;
                m_electron_iso_sf.second = 0;
                m_electron_id_sf = electron_sf->Trigger_ID_ScaleFactors(*cand_it->electron);
		if (isolated_electrons.size()) 
		  {
		    m_electron_iso_sf = electron_sf->IsolationScaleFactor(*cand_it->electron);
		  }
		//Systematic variation applied if requested by ElectronScaleFactorsTool
		hist_weight *= electron_sf->VaryScaleFactor(m_electron_id_sf, m_electron_iso_sf);
            }	    
            write_tool->Write("el_ptcone20", cand_it->electron->mvaptcone20());
            write_tool->Write("el_ptcone30", cand_it->electron->mvaptcone30());
            write_tool->Write("el_ptcone40", cand_it->electron->mvaptcone40());
            write_tool->Write("el_Etcone20", cand_it->electron->Etcone20());
            write_tool->Write("el_Etcone30", cand_it->electron->Etcone30());
            write_tool->Write("el_Etcone40", cand_it->electron->Etcone40());
        }
        else
        {
            lepton_tlv = cand_it->muon->tlv;
            lepton_name = "#mu";
            if(c_channel == 2) // electron + tau channel
            {
                //Apply muon scale factors in MC if required
                if (!is_data() && c_apply_muon_sf ) 
                {
		  const Muon* mu = cand_it->muon;
		  TLorentzVector def_muon_tlv = mu->tlv;
		  def_muon_tlv.SetPtEtaPhiM( mu->pt(), mu->eta(), mu->phi(), 0);
		  float trig_scale_factor = m_staco_trig_scale_fac.scaleFactor( mu->tlv, trigger_decisions->getCurrentTrigger() );
		  float trig_scale_factor_stat_err = m_staco_trig_scale_fac.scaleFactorStatUncertainty( mu->tlv, trigger_decisions->getCurrentTrigger() );
		  float trig_scale_factor_sys_err = m_staco_trig_scale_fac.scaleFactorSysUncertainty( mu->tlv, trigger_decisions->getCurrentTrigger() );
		  float reco_scale_factor = m_staco_cb_scale_fac.scaleFactor( mu->tlv );
		  float reco_scale_factor_err = m_staco_cb_scale_fac.scaleFactorUncertainty( mu->tlv );
		  float isol_scale_factor = 1.0;
		  float isol_scale_factor_stat_err = 0.0;
		  float isol_scale_factor_sys_err = 0.0;
		  if(isolated_muons.size()) {
		    isol_scale_factor = m_staco_isol_scale_fac.scaleFactor( def_muon_tlv );
		    isol_scale_factor_stat_err = m_staco_isol_scale_fac.scaleFactorStatUncertainty( def_muon_tlv );
		    isol_scale_factor_sys_err = m_staco_isol_scale_fac.scaleFactorSysUncertainty( def_muon_tlv );
		  }
	      

		  float scale_factor = trig_scale_factor*reco_scale_factor*isol_scale_factor *
		    ( 1 +
		      c_sys_sigma_muon_sf * std::sqrt( std::pow(trig_scale_factor_stat_err/trig_scale_factor,2) +
						       std::pow(isol_scale_factor_stat_err/isol_scale_factor,2) +
						       std::pow(trig_scale_factor_sys_err/trig_scale_factor,2) +
						       std::pow(isol_scale_factor_sys_err/isol_scale_factor,2) +
						       std::pow(reco_scale_factor_err/reco_scale_factor,2) )
		      );

		  hist_weight = hist_weight_orig*scale_factor;
                }
            }	    
            write_tool->Write("mu_ptcone20", cand_it->muon->ptcone20());
            write_tool->Write("mu_ptcone30", cand_it->muon->ptcone30());
            write_tool->Write("mu_ptcone40", cand_it->muon->ptcone40());
            write_tool->Write("mu_etcone20", cand_it->muon->etcone20());
            write_tool->Write("mu_etcone30", cand_it->muon->etcone30());
            write_tool->Write("mu_etcone40", cand_it->muon->etcone40());
        }
        write_tool->Write("lepton_tlv_Pt", (float) lepton_tlv.Pt());
        write_tool->Write("lepton_tlv_Eta", (float) lepton_tlv.Eta());
        write_tool->Write("lepton_tlv_Phi", (float) lepton_tlv.Phi());
        write_tool->Write("tau_tlv_Pt", (float) tau_tlv.Pt());
        write_tool->Write("tau_tlv_Eta", (float) tau_tlv.Eta());
        write_tool->Write("tau_tlv_Phi", (float) tau_tlv.Phi());
	// ---- some extra properties -------------------------
	write_tool->Write("tau_leadTrkPt", (float) cand_it->tau->leadTrkPt());
	write_tool->Write("tau_EF_dr",     (float) cand_it->tau->EF_dr());
	write_tool->Write("tau_EF_E",      (float) cand_it->tau->EF_E());
	write_tool->Write("tau_EF_eta",    (float) cand_it->tau->EF_eta());
	write_tool->Write("tau_EF_phi",    (float) cand_it->tau->EF_phi());
	write_tool->Write("tau_EF_matched",  (int) cand_it->tau->EF_matched());
	// list of tau triggers:..............................
	// Event Level quantities
	//write_tool->Write("twotau16_loose", (bool) m_trigger->twotau16_loose()); //--> this will crash
// 	write_tool->Write("twotau20_loose", (bool) m_trigger->twotau20_loose());
// 	write_tool->Write("twotau29_loose", (bool) m_trigger->twotau29_loose());
// 	write_tool->Write("twotau29_loose1", (bool) m_trigger->twotau29_loose1());
// 	write_tool->Write("twotau29_loose1_EFxe15_noMu", (bool) m_trigger->twotau29_loose1_EFxe15_noMu());
// 	write_tool->Write("twotau29_loose1_EFxe40_noMu", (bool) m_trigger->twotau29_loose1_EFxe40_noMu());
// 	write_tool->Write("twotau29_medium1", (bool) m_trigger->twotau29_medium1());
// 	write_tool->Write("twotau38_loose", (bool) m_trigger->twotau38_loose());
	write_tool->Write("j35_L1TAU_HV", (bool) m_trigger->j35_L1TAU_HV());
	write_tool->Write("j35_L1TAU_HV_jetNoEF", (bool) m_trigger->j35_L1TAU_HV_jetNoEF());
	write_tool->Write("j35_L1TAU_HV_jetNoEF_cosmic", (bool) m_trigger->j35_L1TAU_HV_jetNoEF_cosmic());
	write_tool->Write("j35_L1TAU_HV_jetNoEF_firstempty", (bool) m_trigger->j35_L1TAU_HV_jetNoEF_firstempty());
	write_tool->Write("j35_L1TAU_HV_jetNoEF_unpaired_iso", (bool) m_trigger->j35_L1TAU_HV_jetNoEF_unpaired_iso());
	write_tool->Write("j35_L1TAU_HV_jetNoEF_unpaired_noniso", (bool) m_trigger->j35_L1TAU_HV_jetNoEF_unpaired_noniso());
	write_tool->Write("tau100_loose", (bool) m_trigger->tau100_loose());
	write_tool->Write("tau100_medium", (bool) m_trigger->tau100_medium());
	write_tool->Write("tau125_loose", (bool) m_trigger->tau125_loose());
	write_tool->Write("tau125_medium", (bool) m_trigger->tau125_medium());
	write_tool->Write("tau16_IDTrkNoCut", (bool) m_trigger->tau16_IDTrkNoCut());
	write_tool->Write("tau16_j75_a4_EFFS_xe40_loose_noMu", (bool) m_trigger->tau16_j75_a4_EFFS_xe40_loose_noMu());
	write_tool->Write("tau16_loose", (bool) m_trigger->tau16_loose());
	write_tool->Write("tau16_loose_e10_loose", (bool) m_trigger->tau16_loose_e10_loose());
	write_tool->Write("tau16_loose_e15_loose", (bool) m_trigger->tau16_loose_e15_loose());
	write_tool->Write("tau16_loose_e15_tight", (bool) m_trigger->tau16_loose_e15_tight());
	write_tool->Write("tau16_loose_mu10", (bool) m_trigger->tau16_loose_mu10());
	write_tool->Write("tau16_loose_mu15", (bool) m_trigger->tau16_loose_mu15());
	write_tool->Write("tau16_loose_xe20_noMu", (bool) m_trigger->tau16_loose_xe20_noMu());
	write_tool->Write("tau16_loose_xe25_noMu", (bool) m_trigger->tau16_loose_xe25_noMu());
	write_tool->Write("tau16_loose_xe30_noMu", (bool) m_trigger->tau16_loose_xe30_noMu());
	write_tool->Write("tau16_medium", (bool) m_trigger->tau16_medium());
	write_tool->Write("tau16_medium_e15_tight", (bool) m_trigger->tau16_medium_e15_tight());
	write_tool->Write("tau16_medium_mu10", (bool) m_trigger->tau16_medium_mu10());
	write_tool->Write("tau16_medium_mu15", (bool) m_trigger->tau16_medium_mu15());
	write_tool->Write("tau16_medium_xe22_noMu", (bool) m_trigger->tau16_medium_xe22_noMu());
	write_tool->Write("tau16_medium_xe25_noMu", (bool) m_trigger->tau16_medium_xe25_noMu());
	write_tool->Write("tau16_medium_xe25_tight_noMu", (bool) m_trigger->tau16_medium_xe25_tight_noMu());
	write_tool->Write("tau20_loose", (bool) m_trigger->tau20_loose());
	write_tool->Write("tau20_loose_e15_tight", (bool) m_trigger->tau20_loose_e15_tight());
	write_tool->Write("tau20_loose_mu15", (bool) m_trigger->tau20_loose_mu15());
	write_tool->Write("tau20_loose_tau29_loose", (bool) m_trigger->tau20_loose_tau29_loose());
	write_tool->Write("tau20_loose_xe25_noMu", (bool) m_trigger->tau20_loose_xe25_noMu());
	write_tool->Write("tau20_medium1", (bool) m_trigger->tau20_medium1());
	write_tool->Write("tau20_medium_tau29_medium_j30", (bool) m_trigger->tau20_medium_tau29_medium_j30());
	write_tool->Write("tau29_loose", (bool) m_trigger->tau29_loose());
	write_tool->Write("tau29_loose1", (bool) m_trigger->tau29_loose1());
	write_tool->Write("tau29_loose_e15_tight", (bool) m_trigger->tau29_loose_e15_tight());
	write_tool->Write("tau29_loose_mu15", (bool) m_trigger->tau29_loose_mu15());
	write_tool->Write("tau29_loose_tau38_loose", (bool) m_trigger->tau29_loose_tau38_loose());
	write_tool->Write("tau29_loose_tau38_loose_j30", (bool) m_trigger->tau29_loose_tau38_loose_j30());
	write_tool->Write("tau29_loose_tau50_loose", (bool) m_trigger->tau29_loose_tau50_loose());
	write_tool->Write("tau29_loose_xs45_loose_noMu_3L1J10", (bool) m_trigger->tau29_loose_xs45_loose_noMu_3L1J10());
	write_tool->Write("tau29_loose_xs80_loose_noMu", (bool) m_trigger->tau29_loose_xs80_loose_noMu());
	write_tool->Write("tau29_medium", (bool) m_trigger->tau29_medium());
	write_tool->Write("tau29_medium1", (bool) m_trigger->tau29_medium1());
	write_tool->Write("tau29_medium1_tau20_medium1", (bool) m_trigger->tau29_medium1_tau20_medium1());
	write_tool->Write("tau29_medium_tau38_loose", (bool) m_trigger->tau29_medium_tau38_loose());
	write_tool->Write("tau29_medium_tau38_medium", (bool) m_trigger->tau29_medium_tau38_medium());
	write_tool->Write("tau29_medium_xe35_noMu", (bool) m_trigger->tau29_medium_xe35_noMu());
	write_tool->Write("tau29_medium_xe40_loose_noMu", (bool) m_trigger->tau29_medium_xe40_loose_noMu());
	write_tool->Write("tau29_medium_xs45_loose_noMu_3L1J10", (bool) m_trigger->tau29_medium_xs45_loose_noMu_3L1J10());
	write_tool->Write("tau29_medium_xs80_loose_noMu", (bool) m_trigger->tau29_medium_xs80_loose_noMu());
	write_tool->Write("tau38_loose", (bool) m_trigger->tau38_loose());
	write_tool->Write("tau38_medium", (bool) m_trigger->tau38_medium());
	write_tool->Write("tau50_IDTrkNoCut", (bool) m_trigger->tau50_IDTrkNoCut());
	write_tool->Write("tau50_loose", (bool) m_trigger->tau50_loose());
	write_tool->Write("tau50_loose_IdScan", (bool) m_trigger->tau50_loose_IdScan());
	write_tool->Write("tau50_medium", (bool) m_trigger->tau50_medium());
	write_tool->Write("tau84_loose", (bool) m_trigger->tau84_loose());
	write_tool->Write("tau84_medium", (bool) m_trigger->tau84_medium());
	write_tool->Write("tauNoCut", (bool) m_trigger->tauNoCut());
	write_tool->Write("tauNoCut_L1TAU50", (bool) m_trigger->tauNoCut_L1TAU50());
	write_tool->Write("tauNoCut_cosmic", (bool) m_trigger->tauNoCut_cosmic());
	write_tool->Write("tauNoCut_firstempty", (bool) m_trigger->tauNoCut_firstempty());
	write_tool->Write("tauNoCut_unpaired_iso", (bool) m_trigger->tauNoCut_unpaired_iso());
	write_tool->Write("tauNoCut_unpaired_noniso", (bool) m_trigger->tauNoCut_unpaired_noniso());
	//
	// match to trigger ..............................................................................
	GET_TOOL(trigger_nav, TrigNavTool, "TrigNav");         //reads per event trig nav maps, chain->obj
	GET_TOOL(trigger_match, TrigMatchTool, "TrigMatch");   //matches trig obj -> offline obj
	
	/* ****/ // example how to access the trigger
	/**********
	int chainId = 619;//trigger_decisions->chainId("trig_tau_EF_tau20_loose", *m_trigdecision_d3pdobject);
	std::vector<int> tauTrigIndxs = 
	  trigger_nav->GetTrigObjIndxs(chainId, std::string("tau"), *m_trigdecision_d3pdobject);

	vector<float>* etaTaus = (*m_trigeftau).eta();
	vector<float>* phiTaus = (*m_trigeftau).phi();
	m_logger<<ALWAYS<<" TETS trig objs: #" << tauTrigIndxs.size() << " and etaTaus.size: " 
		<< etaTaus->size() <<SLogger::endmsg;
	m_logger<<ALWAYS<<" my tau:      eta:" << tau.tlv.Eta() << " phi: " << tau.tlv.Phi() <<SLogger::endmsg;
	m_logger<<ALWAYS<<" printing my trigger indices: " <<SLogger::endmsg;
	for (int i=0; i< tauTrigIndxs.size(); ++i) {
	  m_logger<<ALWAYS<<"  Ind: "<< i << " : "<< tauTrigIndxs[i] <<SLogger::endmsg;
	}
	m_logger<<ALWAYS<<"Trigger Objs: " <<SLogger::endmsg;
	for (int i=0; i<etaTaus->size(); ++i) {
	  m_logger<<ALWAYS<<" -->" << i <<" eta/phi: " << (*etaTaus)[i] << "/" << (*phiTaus)[i] 
		  << " DR from tau: " << dr(tau.tlv.Eta(), (*etaTaus)[i], tau.tlv.Phi(), (*phiTaus)[i])
		  << SLogger::endmsg;
	}
	
	bool decision = trigger_match->Select(tau, tauTrigIndxs, *m_trigeftau);
	if (decision) {
	  m_logger<<ALWAYS<<" TEST tau trigger match: EF_tau16_loose " << chainId << " passed: " << decision
		  <<SLogger::endmsg;
	} else {
	  m_logger<<ALWAYS<<"NO MATCH" <<SLogger::endmsg;
	}
	********************/
	// -- DEBUG: test why tau16_loose may fire the event, but the tau is not matched to trigger
	/********
	if (m_trigger->tau16_loose()) {
	  int chainId = 619;
	  std::vector<int> tauTrigIndxs = 
	    trigger_nav->GetTrigObjIndxs(chainId, std::string("tau"), *m_trigdecision_d3pdobject);
	  bool decision = trigger_match->Select(tau, tauTrigIndxs, *m_trigeftau);
	  if (not decision) {
	    vector<float>* etaTaus = (*m_trigeftau).eta();
	    vector<float>* phiTaus = (*m_trigeftau).phi();
	    m_logger<<ALWAYS<<" TETS trig objs: #" << tauTrigIndxs.size() << " and etaTaus.size: " 
		    << etaTaus->size() <<SLogger::endmsg;
	    m_logger<<ALWAYS<<" my tau:      eta:" << tau.tlv.Eta() << " phi: " << tau.tlv.Phi() <<SLogger::endmsg;
	    m_logger<<ALWAYS<<" printing my trigger indices: " <<SLogger::endmsg;
	    for (int i=0; i< tauTrigIndxs.size(); ++i) {
	      m_logger<<ALWAYS<<"  Ind: "<< i << " : "<< tauTrigIndxs[i] <<SLogger::endmsg;
	    }
	    m_logger<<ALWAYS<<"Trigger Objs: " <<SLogger::endmsg;
	    for (int i=0; i<etaTaus->size(); ++i) {
	      m_logger<<ALWAYS<<" -->" << i <<" eta/phi: " << (*etaTaus)[i] << "/" << (*phiTaus)[i] 
		      << " DR from tau: " << dr(tau.tlv.Eta(), (*etaTaus)[i], tau.tlv.Phi(), (*phiTaus)[i])
		      << SLogger::endmsg;
	    }
	  }
	}
	****************/
	
	// -- DEBUG
	write_tool->Write("matched2_EF_tau100_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(2360,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau100_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(977,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau125_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(360,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau125_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(94,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_IDTrkNoCut", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(470,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_j75_a4_EFFS_xe40_loose_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4063,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(619,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_e10_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(213,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_e15_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(203,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_e15_tight", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(945,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_mu10", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(207,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_mu15", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(238,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_xe20_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(226,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_xe25_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(700,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_loose_xe30_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(702,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(720,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium_e15_tight", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4307,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium_mu10", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(201,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium_mu15", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4311,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium_xe22_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(723,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium_xe25_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(724,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau16_medium_xe25_tight_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(725,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(620,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_loose_e15_tight", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4306,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_loose_mu15", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4310,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_loose_tau29_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4304,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_loose_xe25_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(701,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_medium1", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(976,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau20_medium_tau29_medium_j30", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4315,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(621,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose1", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(721,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_e15_tight", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4305,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_mu15", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4309,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_tau38_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4303,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_tau38_loose_j30", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4313,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_tau50_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4302,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_xs45_loose_noMu_3L1J10", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(5173,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_loose_xs80_loose_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(5171,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(92,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium1", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(974,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium1_tau20_medium1", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(986,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium_tau38_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4301,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium_tau38_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4300,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium_xe35_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(70,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium_xe40_loose_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(71,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium_xs45_loose_noMu_3L1J10", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(5172,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau29_medium_xs80_loose_noMu", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(5170,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau38_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(88,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau38_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(364,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau50_IDTrkNoCut", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(784,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau50_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(90,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau50_loose_IdScan", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(630,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau50_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(366,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau84_loose", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(91,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tau84_medium", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(4000,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tauNoCut", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(93,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tauNoCut_L1TAU50", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(333,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tauNoCut_cosmic", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(739,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tauNoCut_firstempty", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(376,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tauNoCut_unpaired_iso", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(368,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));
	write_tool->Write("matched2_EF_tauNoCut_unpaired_noniso", (bool) trigger_match->Select(tau, trigger_nav->GetTrigObjIndxs(369,std::string("tau"), *m_trigdecision_d3pdobject), *m_trigeftau));



	// ................................................................................................
	//
	// ....................................................
	//
	hist_weight*=tau_jet_electron_sf;
        write_tool->Write("hist_weight", hist_weight);

        bool passcaloisol = 0;
        bool passtrackerisol = 0;

        if(c_channel == 1) // electron
        {
            passcaloisol = electron_isolation->PassCaloIsolation(*(cand_it->electron));
            passtrackerisol = electron_isolation->PassTrackerIsolation(*cand_it->electron);
        }
        else if(c_channel == 2) // muon 
        {
            passcaloisol = muon_isolation->PassCaloIsolation(*cand_it->muon);
            passtrackerisol = muon_isolation->PassTrackerIsolation(*cand_it->muon);
        }

        Cut("TRACK_ISOLATION", passtrackerisol, false );
        Cut("CALO_ISOLATION", passcaloisol, false );
        Cut("LEP_20", lepton_tlv.Pt() > 20*GeV, false );	
        Cut("LEP_40", lepton_tlv.Pt() > 40*GeV, false );	

	const bool is_loose_tau = loose_tau_selector->Select(tau);
	const bool is_tight_tau = tau_selector->Select(tau);      

        Cut("LOOSE_TAU", is_loose_tau, false );
        Cut("TIGHT_TAU", is_tight_tau, false );
        Cut("TAU_20", tau_tlv.Pt() > 20*GeV, false );


        //======================================================================
        // 11 Event selection
        //======================================================================

        //----------------------------------------------------------------------
        // 11.1 Dilpeton veto
        //----------------------------------------------------------------------
        const int n_leptons_dlv = olr_electrons.size() + olr_muons.size();
        CutFlow("DILEPTON_VETO", n_leptons_dlv == 1, false, hist_weight);
	
        //----------------------------------------------------------------------
        // 11.2 Sum cos dphi
        //----------------------------------------------------------------------
        const float sum_cos_dphi = cos( tau_tlv.DeltaPhi(met_tlv) ) + cos( lepton_tlv.DeltaPhi(met_tlv) );
        write_tool->Write("sum_cos_dphi", sum_cos_dphi);
        CutFlow("SUM_COS_DPHI", sum_cos_dphi > -0.15, false, hist_weight);

        //----------------------------------------------------------------------
        // 11.3 Transverse mass
        //----------------------------------------------------------------------
        const float trans_mass = sqrt( 2.0 * lepton_tlv.Pt() * met_tlv.Pt() * (1.0 - cos( lepton_tlv.DeltaPhi(met_tlv) ) ) );
        write_tool->Write("trans_mass", trans_mass);
        CutFlow("TRANS_MASS", trans_mass < 50.0*GeV, false, hist_weight);

        //----------------------------------------------------------------------
        // 11.4 Visible mass
        //----------------------------------------------------------------------
        const float vis_mass = ( lepton_tlv + tau_tlv ).M();
        write_tool->Write("vis_mass", vis_mass);
        CutFlow("VIS_MASS_1", (0.0 <= vis_mass)&&(vis_mass < 200.0*GeV), false, hist_weight);
        CutFlow("VIS_MASS_2", (35.0*GeV <= vis_mass)&&(vis_mass < 75.0*GeV), false, hist_weight);
        Cut("VIS_MASS_3", (80.0*GeV <= vis_mass)&&(vis_mass < 100.0*GeV), false);

        //----------------------------------------------------------------------
        // 11.5 Tau numTrack and charge
        //----------------------------------------------------------------------
	const int tau_numTrack = tau.seedCalo_numTrack();
	const float tau_charge = tau.charge();	   
	
        CutFlow("TAU_NUMTRACK_1_OR_3", (tau_numTrack == 1)||(tau_numTrack == 3), false, hist_weight);
        CutFlow("TAU_UNIT_CHARGE", fabs(tau_charge) == 1, false, hist_weight);
        Cut("TAU_NUMTRACK_1", tau_numTrack == 1, false);
        Cut("TAU_NUMTRACK_3", tau_numTrack == 3, false);
        write_tool->Write("tau_numTrack", tau_numTrack);
        write_tool->Write("tau_charge", tau_charge);
	write_tool->Write("tauCutSafeLoose", tau.tauCutSafeLoose());
	write_tool->Write("tauCutSafeMedium", tau.tauCutSafeMedium());
	write_tool->Write("tauCutSafeTight", tau.tauCutSafeTight());
	write_tool->Write("tau_seedCalo_EMRadius", tau.seedCalo_EMRadius());

        //----------------------------------------------------------------------
        // 11.6 OS
        //----------------------------------------------------------------------
        float lepton_charge = 0.0;
        if(c_channel == 1)
            lepton_charge = cand_it->electron->charge();
        else
            lepton_charge = cand_it->muon->charge();
        const float charge_product = lepton_charge * tau_charge;
	write_tool->Write("charge_product", lepton_charge * tau_charge);

        CutFlow("OS", charge_product < 0.0, false, hist_weight);
        Cut("SS", charge_product > 0.0, false);

        //----------------------------------------------------------------------
        // 11.7 MET
        //----------------------------------------------------------------------
        CutFlow("MET", met_tlv.Pt() > 20.0*GeV, false, hist_weight);
	write_tool->Write("met_tlv_Pt", (float) met_tlv.Pt());
	write_tool->Write("met_tlv_Phi", (float) met_tlv.Phi());

        //----------------------------------------------------------------------
        // 11.8 dphi
        //----------------------------------------------------------------------
        const float tau_lep_dphi = tau_tlv.DeltaPhi(lepton_tlv);
        const float fabs_tau_lep_dphi = fabs(tau_lep_dphi);
        CutFlow("TAU_LEP_DPHI", (1 <= fabs_tau_lep_dphi)&&(fabs_tau_lep_dphi < 2.9), false, hist_weight);
	write_tool->Write("tau_lep_dphi", tau_lep_dphi);

        //----------------------------------------------------------------------
        // 11.9 Invariant mass
        //----------------------------------------------------------------------
        const float inv_mass = calc_inv_mass(lepton_tlv,tau_tlv,met_tlv); // TODO
        //float test_mass = ((lepton_tlv+tau_tlv).M()) / ((tau_tlv.Vect().Cross(lepton_tlv.Vect())).Z() * sqrt(((tau_tlv+met_tlv).Vect().Cross(lepton_tlv.Vect())).Z()   * ((tau_tlv.Vect()).Cross((lepton_tlv+met_tlv).Vect())).Z())    ) ;
        //std::cout << inv_mass << "   " << test_mass << std::endl; 

        CutFlow("INV_MASS_1", (0. < inv_mass)&&(inv_mass < 200.*GeV), false, hist_weight);
        CutFlow("INV_MASS_2", (60.*GeV < inv_mass)&&(inv_mass < 150.*GeV), false, hist_weight);

        const float plus_mass = sqrt( 2.0 * lepton_tlv.Pt() * tau_tlv.Pt() * (1.0 - cos( lepton_tlv.DeltaPhi(tau_tlv) ) ) );
        const float minus_mass = sqrt( 2.0 * lepton_tlv.Pt() * tau_tlv.Pt() * (1.0 + cos( lepton_tlv.DeltaPhi(tau_tlv) ) ) );

        //======================================================================
        // 12 Loop over regions
        //======================================================================
        evaluate_region_results();

        //sanity check
        if( m_region_results.size() != m_control_region_names.size())
            m_logger << WARNING << "Duplicate Control Regions Defined " << SLogger::endmsg;
        
        std::map<std::string, bool>::const_iterator cr_itr = m_region_results.begin();
        std::map<std::string, bool>::const_iterator cr_itr_end = m_region_results.end();	
        for( ; cr_itr != cr_itr_end; ++cr_itr)
        {
            const std::string region_name = (*cr_itr).first;
            const std::string region_dir = std::string("regions/") + region_name;
            const std::string slash_region_dir = std::string("/") + region_dir;

            if((*cr_itr).second)// passed region
            {
                // Event-level histograms
                Book("h_sum_cos_dphi", ";#scale[0.7]{#sum}cos#Delta#phi;Events / 0.05", 80, -2.0, 2.0, region_dir)
                    ->Fill( sum_cos_dphi, hist_weight );
                Book("h_trans_mass", ";m_{T}(" + lepton_name + ", E_{T}^{miss})  [GeV];Events / 2.5 GeV", 60, 0.0, 150.0, region_dir)
                    ->Fill( trans_mass/GeV, hist_weight );
                Book("h_n_leptons_dlv", ";# Leptons;Events", 11, -0.5, 10.5, region_dir)
                    ->Fill( n_leptons_dlv, hist_weight );
                Book("h_charge_product", ";charge(" + lepton_name + ") #times charge(#tau_{h});Events", 11, -5.5, 5.5, region_dir)
                    ->Fill( charge_product, hist_weight );
                Book("h_vis_mass", ";m_{vis}(" + lepton_name + ", #tau_{h})  [GeV];Events / 2.5 GeV", 80, 0.0, 200.0, region_dir)
                    ->Fill( vis_mass/GeV, hist_weight );
                Book("h_inv_mass", ";m_{inv}(" + lepton_name + ", #tau_{h})  [GeV];Events / 2.5 GeV", 80, 0.0, 200.0, region_dir)
                    ->Fill( inv_mass/GeV, hist_weight );
                Book("h_met", ";E_{T}^{miss}  [GeV];Events / 2.5 GeV", 60, 0.0, 150.0, region_dir)
                    ->Fill( met_tlv.Pt()/GeV, hist_weight );
                Book("h_sumet", ";#scale[0.7]{#sum}E_{T}  [GeV];Events / 20 GeV", 50, 0.0, 1000.0, region_dir)
                    ->Fill( met_sumet/GeV, hist_weight );
                Book("h_tau_lep_dphi", ";#Delta#phi(#tau_{h}, " + lepton_name + ");Events / (#pi/30)", 30, 0.0, 3.1416, region_dir)
                    ->Fill( fabs_tau_lep_dphi, hist_weight );
                Book("h_plus_mass", ";m_{T}(" + lepton_name + " + #tau_{h})  [GeV];Events / 2.5 GeV", 80, 0.0, 200.0, region_dir)
                    ->Fill( plus_mass/GeV, hist_weight );
                Book("h_minus_mass", ";m_{T}(" + lepton_name + " - #tau_{h})  [GeV];Events / 2.5 GeV", 80, 0.0, 200.0, region_dir)
                    ->Fill( minus_mass/GeV, hist_weight );
                Book("h_RunNumber", ";RunNumber;entries", 50000, 150000-0.5, 200000-0.5, region_dir)
                    ->Fill( (*m_event_d3pdobject).RunNumber());
                // Vertex Info
                vertex_tool->PlotVertexInfo(*m_vxp_d3pdobject, hist_weight, slash_region_dir);

                // tau histograms
		tau_selector->PlotKinematics(tau, hist_weight, slash_region_dir);
		tau_selector->PlotIdentification(tau, hist_weight, slash_region_dir);
		tau_selector->PlotIDVariables(tau, hist_weight, slash_region_dir);

                // lepton histograms
                if(c_channel == 1)
                {
                    electron_selector->PlotKinematics(*(cand_it->electron), hist_weight,  slash_region_dir);
                    electron_selector->PlotIdentification(*(cand_it->electron), hist_weight,  slash_region_dir);
                    electron_isolation->PlotIsolation(*(cand_it->electron), hist_weight,  slash_region_dir);
                }
                else
                {
                    muon_selector->PlotKinematics(*(cand_it->muon), hist_weight,  slash_region_dir);
                    muon_selector->PlotIdentification(*(cand_it->muon), hist_weight,  slash_region_dir);
                    muon_isolation->PlotIsolation(*(cand_it->muon), hist_weight,  slash_region_dir);
                }
            }
        }
    }

    //==========================================================================
    // 13 Write output ntuple
    //==========================================================================
    // TODO
    write_tool->Write( m_cut_results );
    write_tool->Write( m_region_results );
    if(write_tool->IsOn()){
      m_event_d3pdobject->ReadAllActive();
 //      m_tau_d3pdobject->ReadAllActive();
//       m_tau_d3pdobject->Clear();
//       if( selected_tau ) m_tau_d3pdobject->Add( *(dynamic_cast<D3PDReader::TauD3PDObjectElement*> (selected_tau)) );
    }
}


//______________________________________________________________________________
void VisMassCycle::EndInputDataImp( const SInputData& id ) throw( SError )
{
    delete m_event_d3pdobject;
    //    delete m_trigger_d3pdobject;
    delete m_trigger;
    delete m_trigdecision_d3pdobject;
    delete m_trigefel;
    delete m_trigefmu;
    delete m_trigeftau;
    delete m_vxp_d3pdobject;
    delete m_muon_d3pdobject;
    delete m_electron_d3pdobject;
    delete m_tau_d3pdobject;
    delete m_jet_d3pdobject;
    delete m_met_d3pdobject;
    delete m_truth_d3pdobject;
    delete m_truetau_d3pdobject;
    delete m_mcevent_d3pdobject;
    delete m_cluster_d3pdobject;

    // clean up control regions
    m_control_region_names.clear();
    m_control_region_reqs.clear();
    m_control_region_vetoes.clear();
    m_region_results.clear();

}

//______________________________________________________________________________
void VisMassCycle::EndMasterInputDataImp( const SInputData& id) throw( SError )
{
    PrintCutFlow(id);
    m_event_map.clear();
}

//______________________________________________________________________________
void VisMassCycle::EndCycleImp() throw( SError )
{ }

//______________________________________________________________________________
void VisMassCycle::parse_control_regions() throw( SError )
{
    // split colons to get region:CUT1,CUT2,... region2:CUT3,CUT4,...
    std::vector<std::string> semicolon_split;
    StringHelper::split(semicolon_split, c_control_regions, ';');

    for(unsigned int i=0; i < semicolon_split.size(); i++)
    {
        std::string line = semicolon_split[i];
        line = StringHelper::strip(line);
        if(line.empty()) continue;

        std::vector<std::string> colon_split;
        StringHelper::split(colon_split, line, ':');
        if(colon_split.size() != 2)
        {
            // must be even (have words on either side of the colons)
            m_logger << FATAL << "Error in parse_control_regions(): colon_split.size() != 2\n";
            m_logger << FATAL << "    line = " << line << SLogger::endmsg;
            throw SError( SError::StopExecution );
        }

        std::string region_name = colon_split[0];
        region_name = StringHelper::strip(region_name);
        if(region_name.empty())
        {
            // must not be empty
            m_logger << FATAL << "Error in parse_control_regions(): region_name.empty()" << SLogger::endmsg;
            throw SError( SError::StopExecution );
        }
        // split commas to get cuts for a region
        std::vector<std::string> comma_split;
        StringHelper::split(comma_split, colon_split[1], ',');
        std::vector<std::string> reqs;
        std::vector<std::string> vetoes;
        for(unsigned int j=0; j < comma_split.size(); j++)
        {
            std::string cut_name = StringHelper::strip(comma_split[j]);
            if(cut_name[0] == '!')
            {
                cut_name = cut_name.substr(1, cut_name.size()-1);
                vetoes.push_back(cut_name);
            }
            else
            {
                reqs.push_back(cut_name);
            }
            if(cut_name.empty())
            {
                // must not be empty
                m_logger << FATAL << "Error in parse_control_regions(): cut_name.empty()" << SLogger::endmsg;
                throw SError( SError::StopExecution );
            }
        }
        // set regions
        m_control_region_names.push_back(region_name);
        m_control_region_reqs.push_back(reqs);
        m_control_region_vetoes.push_back(vetoes);
    }
}

//______________________________________________________________________________
void VisMassCycle::evaluate_region_results()
{
    for(unsigned int i_region=0; i_region < m_control_region_names.size(); i_region++)
    {
        const std::vector<std::string>& region_reqs = m_control_region_reqs[i_region];
        const std::vector<std::string>& region_vetoes = m_control_region_vetoes[i_region];
        bool passed_region = true;

        for(unsigned int i_cut = 0; i_cut < region_reqs.size(); i_cut++)
        {
            const std::string& cut_name = region_reqs[i_cut];
            std::map<std::string,bool>::iterator it = m_cut_results.find(cut_name);
            if(it != m_cut_results.end())
            {
                passed_region = passed_region && it->second;
            }
            else 
            {
                m_logger << ERROR << "The cut " << cut_name << " has not yet been applied -> result is undefined!"
                    << SLogger::endmsg;
            }

            if(!passed_region) break;
        }

        for(unsigned int i_cut = 0; i_cut < region_vetoes.size(); i_cut++)
        {
            const std::string& cut_name = region_vetoes[i_cut];
            std::map<std::string,bool>::iterator it = m_cut_results.find(cut_name);
            if(it != m_cut_results.end())
            {
                passed_region = passed_region && !(it->second);
            }
            else
            {
                m_logger << ERROR << "The cut " << cut_name << " has not yet been applied -> result is undefined!"
                    << SLogger::endmsg;
            }

            if(!passed_region) break;
        }

        m_region_results[ m_control_region_names[i_region] ] = passed_region;
    }
}

//______________________________________________________________________________
float VisMassCycle::calc_inv_mass( const TLorentzVector& lep, const TLorentzVector& tau, const TLorentzVector& met ){
    //https://svnweb.cern.ch/trac/atlasoff/browser/PhysicsAnalysis/StandardModelPhys/WZBenchmarks/ZToTwoTaus/trunk/src/InvMassCalculatorTool.cxx

    float rexy        = tau.Py()/tau.Px();
    float rex         = tau.Px()/tau.P();
    float rcx         = lep.Px()/lep.P();
    float rcy         = lep.Py()/lep.P();
    float ta          = met.Py() - rexy*met.Px();
    float tb          = rcy-rexy*rcx;
    float eNu1        = ta/tb;
    float eNu2        = (met.Px()-rcx*eNu1)/rex;

    if(eNu1 > 0 && eNu2 > 0) {
        float eTau1 = lep.P()+eNu1;
        float eTau2 = tau.P()+eNu2;
        float cosTheta = (lep.Px()*tau.Px() + lep.Py()*tau.Py() + lep.Pz()*tau.Pz()) /  (lep.P()*tau.P());

        float mass2 = 2 * eTau1 * eTau2 * (1. - cosTheta);
        return std::sqrt(mass2);
    }

    return -100.;
}

// EOF