DiTauSelectionCuts.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
 
// local include(s)
#include "TauAnalysisTools/DiTauSelectionCuts.h"
#include "TauAnalysisTools/DiTauSelectionTool.h"
 
using namespace TauAnalysisTools;
 
//______________________________________________________________________________
DiTauSelectionCut::DiTauSelectionCut(const std::string& sName, TauAnalysisTools::DiTauSelectionTool* tDTST)
  : m_sName(sName)
  , m_hHistCutPre(nullptr)
  , m_hHistCut(nullptr)
  , m_tDTST(tDTST)
{}
 
//______________________________________________________________________________
DiTauSelectionCut::~DiTauSelectionCut()
{
  // FIXME: could use unique_ptr
  delete m_hHistCutPre;
  delete m_hHistCut;
}
 
//______________________________________________________________________________
void DiTauSelectionCut::writeControlHistograms()
{
  m_hHistCutPre->Write();
  m_hHistCut->Write();
}
 
//______________________________________________________________________________
TH1F* DiTauSelectionCut::CreateControlPlot(const char* sName, const char* sTitle, int iBins, double dXLow, double dXUp)
{
  if (m_tDTST->m_bCreateControlPlots)
  {
    TH1F* hHist = new TH1F(sName, sTitle, iBins, dXLow, dXUp);
    hHist->SetDirectory(0);
    return hHist;
  }
 
  return nullptr;
}
 
//______________________________________________________________________________
void DiTauSelectionCut::fillHistogramCutPre(const xAOD::DiTauJet& xDiTau)
{
  fillHistogram(xDiTau, *m_hHistCutPre);
}
 
//______________________________________________________________________________
void DiTauSelectionCut::fillHistogramCut(const xAOD::DiTauJet& xDiTau)
{
  fillHistogram(xDiTau, *m_hHistCut);
}
 
 
//______________________________________________________________________________
void DiTauSelectionCut::setProperty(const std::string& name, const std::string& value)
{
  std::map<std::string, std::string&>::iterator it = m_mProperties.find(name);
  if(it == m_mProperties.end() )
    throw std::runtime_error (("Undeclared property: " + name + "\n").c_str());
  it->second = value;
}
 
//______________________________________________________________________________
void DiTauSelectionCut::declareProperty(const std::string& name, std::string& loc)
{
  std::pair<std::string, std::string&> p(name, loc);
  m_mProperties.insert(p);
}
 
//______________________________________________________________________________
std::string DiTauSelectionCut::getProperty(const std::string& name)
{
  std::map<std::string, std::string&>::iterator it = m_mProperties.find(name);
  if(it == m_mProperties.end() )
    throw std::runtime_error (("Undeclared property: " + name + "\n").c_str());
 
  return it->second;
}
 
 
//_______________________________SelectionCutPt_________________________________
//______________________________________________________________________________
DiTauSelectionCutPt::DiTauSelectionCutPt(DiTauSelectionTool* tDTST)
  : DiTauSelectionCut("CutPt", tDTST)
{
  m_hHistCutPre = CreateControlPlot("hPt_pre","Pt_pre;di-#tau-p_{T} [GeV]; events",1000,0,1000);
  m_hHistCut = CreateControlPlot("hPt_cut","Pt_cut;di-#tau-p_{T} [GeV]; events",1000,0,1000);
}
 
//______________________________________________________________________________
void DiTauSelectionCutPt::fillHistogram(const xAOD::DiTauJet& xDiTau, TH1F& hHist) const
{
  hHist.Fill(xDiTau.pt()/1000.);
}
 
//______________________________________________________________________________
void DiTauSelectionCutPt::setAcceptInfo(asg::AcceptInfo& info) const
{
  info.addCut( "Pt",
               "Selection of ditaus according to their transverse momentum" );
}
//______________________________________________________________________________
bool DiTauSelectionCutPt::accept(const xAOD::DiTauJet& xDiTau,
                            asg::AcceptData& acceptData)
{
  // save ditau pt in GeV
  double pt = xDiTau.pt() / 1000.;
  // in case of only one entry in vector, run for lower limits
  if (m_tDTST->m_vPtRegion.size() == 1)
  {
    if ( pt >= m_tDTST->m_vPtRegion.at(0) )
    {
      acceptData.setCutResult( "Pt", true );
      return true;
    }
  }
  unsigned int iNumPtRegion = m_tDTST->m_vPtRegion.size()/2;
  for( unsigned int iPtRegion = 0; iPtRegion < iNumPtRegion; iPtRegion++ )
  {
    if ( pt >= m_tDTST->m_vPtRegion.at(iPtRegion*2) and pt <= m_tDTST->m_vPtRegion.at(iPtRegion*2+1))
    {
      acceptData.setCutResult( "Pt", true );
      return true;
    }
  }
  m_tDTST->msg() << MSG::VERBOSE << "DiTau failed pt requirement, ditau pt [GeV]: " << pt << endmsg;
  return false;
}
 
//_____________________________SelectionCutAbsEta_______________________________
//______________________________________________________________________________
DiTauSelectionCutAbsEta::DiTauSelectionCutAbsEta(DiTauSelectionTool* tDTST)
  : DiTauSelectionCut("CutAbsEta", tDTST)
{
  m_hHistCutPre = CreateControlPlot("hEta_pre","Eta_pre;di-#tau-#eta; events",100,-3,3);
  m_hHistCut = CreateControlPlot("hEta_cut","Eta_cut;di-#tau-#eta; events",100,-3,3);
}
 
//______________________________________________________________________________
void DiTauSelectionCutAbsEta::fillHistogram(const xAOD::DiTauJet& xDiTau, TH1F& hHist) const
{
  hHist.Fill(xDiTau.eta());
}
 
//______________________________________________________________________________
void DiTauSelectionCutAbsEta::setAcceptInfo(asg::AcceptInfo& info) const
{
  info.addCut( "AbsEta",
               "Selection of ditaus according to their absolute pseudorapidity" );
}
//______________________________________________________________________________
bool DiTauSelectionCutAbsEta::accept(const xAOD::DiTauJet& xDiTau,
                                asg::AcceptData& acceptData)
{
  // check regions of eta, if ditau is in one region then return true; false otherwise
  unsigned int iNumEtaRegion = m_tDTST->m_vAbsEtaRegion.size()/2;
  for( unsigned int iEtaRegion = 0; iEtaRegion < iNumEtaRegion; iEtaRegion++ )
  {
    if ( std::abs( xDiTau.eta() ) >= m_tDTST->m_vAbsEtaRegion.at(iEtaRegion*2) and std::abs( xDiTau.eta() ) <= m_tDTST->m_vAbsEtaRegion.at(iEtaRegion*2+1))
    {
      acceptData.setCutResult( "AbsEta", true );
      return true;
    }
  }
  m_tDTST->msg() << MSG::VERBOSE << "DiTau failed eta requirement, ditau eta: " << xDiTau.eta() << endmsg;
  return false;
}
 
//_______________________________SelectionCutNSubjets_________________________________
//______________________________________________________________________________
DiTauSelectionCutNSubjets::DiTauSelectionCutNSubjets(DiTauSelectionTool* tDTST)
  : DiTauSelectionCut("CutNSubjets", tDTST)
{
  m_hHistCutPre = CreateControlPlot("hNSubjets_pre","NSubjets_pre;di-#tau NSubjets; events",30,0,30);
  m_hHistCut = CreateControlPlot("hNSubjets_cut","NSubjets_cut;di-#tau NSubjets; events",30,0,30);
}
 
//______________________________________________________________________________
void DiTauSelectionCutNSubjets::fillHistogram(const xAOD::DiTauJet& xDiTau, TH1F& hHist) const
{
  hHist.Fill(xDiTau.nSubjets());
}
 
//______________________________________________________________________________
void DiTauSelectionCutNSubjets::setAcceptInfo(asg::AcceptInfo& info) const
{
  info.addCut( "NSubjets",
               "Selection of ditaus according to their number of subjets" );
}
//______________________________________________________________________________
bool DiTauSelectionCutNSubjets::accept(const xAOD::DiTauJet& xDiTau,
                            asg::AcceptData& acceptData)
{
  float nsubjets = xDiTau.nSubjets();
  // in case of only one entry in vector, run for lower limits
  if (m_tDTST->m_vNSubjetsRegion.size() == 1)
  {
    if ( nsubjets >= m_tDTST->m_vNSubjetsRegion.at(0) )
    {
      acceptData.setCutResult( "NSubjets", true );
      return true;
    }
  }
  unsigned int iNumNSubjetsRegion = m_tDTST->m_vNSubjetsRegion.size()/2;
  for( unsigned int iNSubjetsRegion = 0; iNSubjetsRegion < iNumNSubjetsRegion; iNSubjetsRegion++ )
  {
    if ( nsubjets >= m_tDTST->m_vNSubjetsRegion.at(iNSubjetsRegion*2) and nsubjets <= m_tDTST->m_vNSubjetsRegion.at(iNSubjetsRegion*2+1))
    {
      acceptData.setCutResult( "NSubjets", true );
      return true;
    }
  }
  m_tDTST->msg() << MSG::VERBOSE << "DiTau failed NSubjets requirement, ditau number of subjets: " << nsubjets << endmsg;
  return false;
}
 
//____________________________SelectionCutAbsCharge_____________________________
//______________________________________________________________________________
DiTauSelectionCutAbsCharge::DiTauSelectionCutAbsCharge(DiTauSelectionTool* tDTST)
  : DiTauSelectionCut("CutAbsCharge", tDTST),
    m_bDiTauCharge(-1234) 
{
  m_hHistCutPre = CreateControlPlot("hCharge_pre","Charge_pre;charge; events",7,-3.5,3.5);
  m_hHistCut = CreateControlPlot("hCharge_cut","Charge_cut;charge; events",7,-3.5,3.5);
}
 
//______________________________________________________________________________
void DiTauSelectionCutAbsCharge::fillHistogram(const xAOD::DiTauJet& /*xTau*/, TH1F& hHist) const
{
  hHist.Fill(m_bDiTauCharge); 
}
 
//______________________________________________________________________________
void DiTauSelectionCutAbsCharge::setAcceptInfo(asg::AcceptInfo& info) const
{
  info.addCut( "AbsCharge",
               "Selection of taus according to their absolute charge" );
}
//______________________________________________________________________________
bool DiTauSelectionCutAbsCharge::accept(const xAOD::DiTauJet& xTau,
                            asg::AcceptData& acceptData)
{
  m_bDiTauCharge = 0;
  static const SG::ConstAccessor<float> acc_charge ("charge");
  if ( acc_charge.isAvailable(xTau) ) {
     m_bDiTauCharge = acc_charge(xTau);
  } else {
     for (const auto& xTrack : xTau.trackLinks()) {
        if (!xTrack.isValid())
           continue;
 
        if(xTau.nSubjets() >= 2){
           for (int i = 0; i < 2; ++i) { // loop over two leading subjets 
              TLorentzVector tlvSubjet = TLorentzVector();
              tlvSubjet.SetPtEtaPhiE(xTau.subjetPt(i), xTau.subjetEta(i),
                                     xTau.subjetPhi(i), xTau.subjetE(i));
              double dR = tlvSubjet.DeltaR((*xTrack)->p4());
              if (dR < 0.1) {
                 m_bDiTauCharge += (*xTrack)->charge();
                 break; //prevents double counting of tracks
              }
           }  // loop over subjets
        }      
     } // loop over tracks
  }
 
  // check charge, if ditau has one of the charges requiered then return true; false otherwise
  for( unsigned int iCharge = 0; iCharge < m_tDTST->m_vAbsCharges.size(); iCharge++ )
  {
    if ( std::abs( m_bDiTauCharge ) == m_tDTST->m_vAbsCharges.at(iCharge) )
    {
      acceptData.setCutResult( "AbsCharge", true );
      return true;
    }
  }
  m_tDTST->msg() << MSG::VERBOSE << "DiTau failed charge requirement, ditau charge: " << m_bDiTauCharge << endmsg;
  return false;
 
}
 
//___________________________SelectionCutOmniScore_____________________________
//______________________________________________________________________________
DiTauSelectionCutOmniScore::DiTauSelectionCutOmniScore(DiTauSelectionTool* tDTST)
  : DiTauSelectionCut("CutOmniScore", tDTST)
{
  m_hHistCutPre = CreateControlPlot("hOmniScore_pre","OmniScore_pre;OmniScore; events",100,0,1);
  m_hHistCut = CreateControlPlot("hOmniScore_cut","OmniScore_cut;OmniScore; events",100,0,1);
}
//______________________________________________________________________________
void DiTauSelectionCutOmniScore::fillHistogram(const xAOD::DiTauJet& xTau, TH1F& hHist) const
{
  static const SG::ConstAccessor<float> acc_OmniScore("omni_score");    
  hHist.Fill(acc_OmniScore(xTau));
}
//______________________________________________________________________________
void DiTauSelectionCutOmniScore::setAcceptInfo(asg::AcceptInfo& info) const
{
  info.addCut( "OmniScore",
               "Selection of taus according to their OmniScore" );
}
//______________________________________________________________________________
bool DiTauSelectionCutOmniScore::accept(const xAOD::DiTauJet& xTau,
                                              asg::AcceptData& acceptData)
{
  // check OmniScore score, if tau has a OmniScore score in one of the regions requiered then return true; false otherwise
  static const SG::ConstAccessor<float> acc ("omni_score");
  float dOmniScore = acc(xTau); 
  unsigned int iNumOmniScoreRegion = m_tDTST->m_vOmniScoreRegion.size()/2;
  for( unsigned int iOmniScoreRegion = 0; iOmniScoreRegion < iNumOmniScoreRegion; iOmniScoreRegion++ )
  {
    if ( dOmniScore >= m_tDTST->m_vOmniScoreRegion.at(iOmniScoreRegion*2) and dOmniScore <= m_tDTST->m_vOmniScoreRegion.at(iOmniScoreRegion*2+1))
    {
      acceptData.setCutResult( "OmniScore", true );
      return true;
    }
  }
  m_tDTST->msg() << MSG::VERBOSE << "Tau failed OmniScore requirement, tau OmniScore: " << dOmniScore << endmsg;
  return false;
}
 
//_____________________________SelectionCutOmniIDWP______________________________
//_______________________________________________________________________________
DiTauSelectionCutOmniIDWP::DiTauSelectionCutOmniIDWP(DiTauSelectionTool* tDTST)
  : DiTauSelectionCut("CutOmniIDWP", tDTST)
{
  m_hHistCutPre = CreateControlPlot("hOmniIDWP_pre","OmniIDWP_pre;; events",8,-.5,7.5);
  m_hHistCut = CreateControlPlot("hOmniIDWP_cut","OmniIDWP_cut;; events",8,-.5,7.5);
  // only proceed if histograms are defined
  if (!m_hHistCutPre or !m_hHistCut)
    return;
 
  m_hHistCutPre->GetXaxis()->SetBinLabel(1,"!VeryLoose");
  m_hHistCutPre->GetXaxis()->SetBinLabel(2,"VeryLoose");
  m_hHistCutPre->GetXaxis()->SetBinLabel(3,"!Loose");
  m_hHistCutPre->GetXaxis()->SetBinLabel(4,"Loose");
  m_hHistCutPre->GetXaxis()->SetBinLabel(5,"!Medium");
  m_hHistCutPre->GetXaxis()->SetBinLabel(6,"Medium");
  m_hHistCutPre->GetXaxis()->SetBinLabel(7,"!Tight");
  m_hHistCutPre->GetXaxis()->SetBinLabel(8,"Tight");
  m_hHistCut->GetXaxis()->SetBinLabel(1,"!VeryLoose");
  m_hHistCut->GetXaxis()->SetBinLabel(2,"VeryLoose");
  m_hHistCut->GetXaxis()->SetBinLabel(3,"!Loose");
  m_hHistCut->GetXaxis()->SetBinLabel(4,"Loose");
  m_hHistCut->GetXaxis()->SetBinLabel(5,"!Medium");
  m_hHistCut->GetXaxis()->SetBinLabel(6,"Medium");
  m_hHistCut->GetXaxis()->SetBinLabel(7,"!Tight");
  m_hHistCut->GetXaxis()->SetBinLabel(8,"Tight");
}
 
//______________________________________________________________________________
void DiTauSelectionCutOmniIDWP::fillHistogram(const xAOD::DiTauJet& xTau, TH1F& hHist) const
{
  if(m_tDTST->m_useOmniScore){ 
     static const SG::ConstAccessor<char> acc_OmniVeryLoose("omni_score_VL");
     static const SG::ConstAccessor<char> acc_OmniLoose("omni_score_L");
     static const SG::ConstAccessor<char> acc_OmniMedium("omni_score_M");
     static const SG::ConstAccessor<char> acc_OmniTight("omni_score_T");
     hHist.Fill(acc_OmniVeryLoose(xTau));
     hHist.Fill(acc_OmniLoose(xTau)+2);
     hHist.Fill(acc_OmniMedium(xTau)+4);
     hHist.Fill(acc_OmniTight(xTau)+6);
  } 
}
 
//______________________________________________________________________________
void DiTauSelectionCutOmniIDWP::setAcceptInfo(asg::AcceptInfo& info) const
{
  info.addCut( "OmniIDWP",
               "Selection of ditaus according to their OmniIDScore" );
}
//______________________________________________________________________________
bool DiTauSelectionCutOmniIDWP::accept(const xAOD::DiTauJet& xTau,
                                 asg::AcceptData& acceptData)
{
  // check Omni ID working point, if ditau passes OmniID working point then return true; false otherwise
  bool bPass = false;
  switch (m_tDTST->m_iOmniIDWP)
  {
  case OMNIIDNONE:
    bPass = true;
    break;
  case OMNIIDVERYLOOSE:
    static const SG::ConstAccessor<char> acc_OmniVeryLoose("omni_score_VL");
    if (!acc_OmniVeryLoose.isAvailable(xTau)) m_tDTST->msg() << MSG::WARNING << "Omni VeryLoose WP not available" << endmsg;
    else bPass = acc_OmniVeryLoose(xTau);
    break;
  case OMNIIDLOOSE:
    static const SG::ConstAccessor<char> acc_OmniLoose("omni_score_L");
    if (!acc_OmniLoose.isAvailable(xTau)) m_tDTST->msg() << MSG::WARNING << "Omni Loose WP not available" << endmsg;
    else bPass = acc_OmniLoose(xTau);
    break;
  case OMNIIDMEDIUM:
    static const SG::ConstAccessor<char> acc_OmniMedium("omni_score_M");
    if (!acc_OmniMedium.isAvailable(xTau)) m_tDTST->msg() << MSG::WARNING << "Omni Medium WP not available" << endmsg;
    else bPass = acc_OmniMedium(xTau);
    break;
  case OMNIIDTIGHT:
    static const SG::ConstAccessor<char> acc_OmniTight("omni_score_T");
    if (!acc_OmniTight.isAvailable(xTau)) m_tDTST->msg() << MSG::WARNING << "Omni Tight WP not available" << endmsg;
    else bPass = acc_OmniTight(xTau);
    break;
  default:
    m_tDTST->msg() << MSG::WARNING << "The Omni ID working point with the enum " << m_tDTST->m_iOmniIDWP << " is not available" << endmsg;
    break;
  }
  if (bPass)
  {
    acceptData.setCutResult( "OmniIDWP", true );
    return true;
  }
  m_tDTST->msg() << MSG::VERBOSE << "DiTau failed OmniIDWP requirement" << endmsg;  
  return false;
}