diff --git a/Reconstruction/tauRecTools/CMakeLists.txt b/Reconstruction/tauRecTools/CMakeLists.txt index cba8060dbeb003f9a1d0165b4653c6e4bbcc6a3f..8d2bbe9c02de0016fab00a79cf4dec5d71123c7a 100644 --- a/Reconstruction/tauRecTools/CMakeLists.txt +++ b/Reconstruction/tauRecTools/CMakeLists.txt @@ -24,7 +24,7 @@ atlas_depends_on_subdirs( PUBLIC Event/xAOD/xAODPFlow GaudiKernel InnerDetector/InDetRecTools/InDetRecToolInterfaces - PhysicsAnalysis/TauID/TauAnalysisTools + InnerDetector/InDetRecTools/InDetTrackSelectionTool Reconstruction/Jet/JetEDM Reconstruction/Particle Reconstruction/RecoTools/ITrackToVertex diff --git a/Reconstruction/tauRecTools/Root/CombinedP4FromRecoTaus.cxx b/Reconstruction/tauRecTools/Root/CombinedP4FromRecoTaus.cxx index 28cfabd82a5aad6895a83a75fc00767465d11222..2c4fb4cac4e78329578fdde9886e8f3b69426d55 100644 --- a/Reconstruction/tauRecTools/Root/CombinedP4FromRecoTaus.cxx +++ b/Reconstruction/tauRecTools/Root/CombinedP4FromRecoTaus.cxx @@ -5,7 +5,7 @@ // Framework include(s) #include "PathResolver/PathResolver.h" -#include "TauAnalysisTools/HelperFunctions.h" +//#include "TauAnalysisTools/HelperFunctions.h" // local include(s) #include "tauRecTools/CombinedP4FromRecoTaus.h" @@ -139,7 +139,7 @@ StatusCode CombinedP4FromRecoTaus::initialize() { StatusCode CombinedP4FromRecoTaus::execute(xAOD::TauJet& xTau) { xAOD::TauJet* Tau = &xTau; int tmpDecayModeProto; - xTau.panTauDetail(xAOD::TauJetParameters::PanTauDetails::pantau_CellBasedInput_DecayModeProto, tmpDecayModeProto); + xTau.panTauDetail(xAOD::TauJetParameters::PanTauDetails::PanTau_DecayModeProto, tmpDecayModeProto); if(tmpDecayModeProto>xAOD::TauJetParameters::Mode_3pXn){ xTau.auxdecor<float>("pt_combined") = 1.; return StatusCode::SUCCESS; @@ -554,7 +554,7 @@ TLorentzVector CombinedP4FromRecoTaus::getCombinedP4(const xAOD::TauJet* tau) { ATH_MSG_DEBUG( "TauRecET: " << tauRecP4.Et() ); xAOD::TauJetParameters::DecayMode decayMode = xAOD::TauJetParameters::DecayMode::Mode_Error; int tmpDecayMode; - if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::pantau_CellBasedInput_DecayMode, tmpDecayMode)) { + if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::PanTau_DecayMode, tmpDecayMode)) { decayMode = static_cast< xAOD::TauJetParameters::DecayMode>(tmpDecayMode); } ATH_MSG_DEBUG( "Decaymode is: " << decayMode ); @@ -598,7 +598,7 @@ TLorentzVector CombinedP4FromRecoTaus::getCalibratedConstituentP4(const xAOD::Ta xAOD::TauJetParameters::DecayMode decayMode = xAOD::TauJetParameters::DecayMode::Mode_Error; int tmpDecayMode; - if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::pantau_CellBasedInput_DecayMode, tmpDecayMode)) { + if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::PanTau_DecayMode, tmpDecayMode)) { decayMode = static_cast< xAOD::TauJetParameters::DecayMode>(tmpDecayMode); } ATH_MSG_DEBUG( "Decaymode is: " << decayMode ); @@ -638,7 +638,7 @@ TLorentzVector CombinedP4FromRecoTaus::getCalibratedTauRecP4(const xAOD::TauJet* ATH_MSG_DEBUG( "TauRecET: " << tauRecP4.Et() ); xAOD::TauJetParameters::DecayMode decayMode = xAOD::TauJetParameters::DecayMode::Mode_Error; int tmpDecayMode; - if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::pantau_CellBasedInput_DecayMode, tmpDecayMode)) { + if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::PanTau_DecayMode, tmpDecayMode)) { decayMode = static_cast< xAOD::TauJetParameters::DecayMode>(tmpDecayMode); } ATH_MSG_DEBUG( "Decaymode is: " << decayMode ); @@ -679,7 +679,7 @@ TLorentzVector CombinedP4FromRecoTaus::getWeightedP4(const xAOD::TauJet* tau) { ATH_MSG_DEBUG( "TauRecET: " << tauRecP4.Et() ); xAOD::TauJetParameters::DecayMode decayMode = xAOD::TauJetParameters::DecayMode::Mode_Error; int tmpDecayMode; - if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::pantau_CellBasedInput_DecayMode, tmpDecayMode)) { + if (tau->panTauDetail(xAOD::TauJetParameters::PanTauDetails::PanTau_DecayMode, tmpDecayMode)) { decayMode = static_cast< xAOD::TauJetParameters::DecayMode>(tmpDecayMode); } ATH_MSG_DEBUG( "Decaymode is: " << decayMode ); diff --git a/Reconstruction/tauRecTools/Root/LinkDef.h b/Reconstruction/tauRecTools/Root/LinkDef.h index c13ca4fb2dbc36b659aaab4d27185b83e7e8d05a..6f73af0eb33889a53e5df85281d75de8c42cd768 100644 --- a/Reconstruction/tauRecTools/Root/LinkDef.h +++ b/Reconstruction/tauRecTools/Root/LinkDef.h @@ -12,6 +12,7 @@ #include "tauRecTools/ITauToolBase.h" #include "tauRecTools/MvaTESVariableDecorator.h" #include "tauRecTools/MvaTESEvaluator.h" +#include "tauRecTools/TauTrackClassifier.h" #include "tauRecTools/CombinedP4FromRecoTaus.h" #ifdef __CINT__ @@ -23,7 +24,6 @@ #endif - #ifdef __CINT__ #pragma link C++ class TauCalibrateLC+; @@ -37,6 +37,8 @@ #pragma link C++ class ITauToolBase+; #pragma link C++ class MvaTESVariableDecorator+; #pragma link C++ class MvaTESEvaluator+; +#pragma link C++ class TauTrackClassifier+; +#pragma link C++ class TrackMVABDT+; #pragma link C++ class CombinedP4FromRecoTaus+; #endif diff --git a/Reconstruction/tauRecTools/Root/MvaTESEvaluator.cxx b/Reconstruction/tauRecTools/Root/MvaTESEvaluator.cxx index 8c15a3bc1bcab88fdc720892709ea105203726ef..3a92c114731253ee635ec98199653c495b04beb5 100644 --- a/Reconstruction/tauRecTools/Root/MvaTESEvaluator.cxx +++ b/Reconstruction/tauRecTools/Root/MvaTESEvaluator.cxx @@ -7,7 +7,7 @@ // tools include(s) -#include "TauAnalysisTools/HelperFunctions.h" +//#include "TauAnalysisTools/HelperFunctions.h" //_____________________________________________________________________________ MvaTESEvaluator::MvaTESEvaluator(const std::string& name) @@ -75,36 +75,34 @@ StatusCode MvaTESEvaluator::execute(xAOD::TauJet& xTau){ nVtx = xTau.auxdata<int>("nVtx"); // Retrieve seed jet info - center_lambda = xTau.auxdata<double>("center_lambda"); - first_eng_dens = xTau.auxdata<double>("first_eng_dens"); - em_probability = xTau.auxdata<double>("em_probability"); - second_lambda = xTau.auxdata<double>("second_lambda"); - presampler_frac = xTau.auxdata<double>("presampler_frac"); - nMuSeg = (float)xTau.auxdata<int>("GhostMuonSegmentCount"); - + xTau.detail(xAOD::TauJetParameters::ClustersMeanCenterLambda, center_lambda); + xTau.detail(xAOD::TauJetParameters::ClustersMeanFirstEngDens, first_eng_dens); + xTau.detail(xAOD::TauJetParameters::ClustersMeanEMProbability,em_probability); + xTau.detail(xAOD::TauJetParameters::ClustersMeanSecondLambda, second_lambda); + xTau.detail(xAOD::TauJetParameters::ClustersMeanPresamplerFrac, presampler_frac); + int nMuSeg_i=0; + xTau.detail(xAOD::TauJetParameters::GhostMuonSegmentCount, nMuSeg_i); + nMuSeg=nMuSeg_i; + // Retrieve pantau and LC-precalib TES - seedCalo_eta = xTau.auxdata<float>("seedCalo_eta"); - float pT_LC = xTau.auxdata<float>("LC_TES_precalib"); + seedCalo_eta = xTau.etaDetectorAxis(); + float pT_LC = xTau.ptDetectorAxis(); float pT_pantau = xTau.ptPanTauCellBased(); - interpolPt = xTau.auxdata<double>("LC_pantau_interpolPt"); + xTau.detail(xAOD::TauJetParameters::LC_pantau_interpolPt, interpolPt); LC_D_interpolPt = pT_LC / interpolPt; pantau_D_interpolPt = pT_pantau / interpolPt; // Retrieve substructures info nTracks = (float)xTau.nTracks(); - nPi0PFOs = (float)xTau.auxdata<int>("nPi0PFOs"); - PFOEngRelDiff = xTau.auxdata<double>("PFOEngRelDiff"); + nPi0PFOs = (float) xTau.nPi0s(); + xTau.detail(xAOD::TauJetParameters::PFOEngRelDiff, PFOEngRelDiff); // "Retrieve" spectator variables truthPtVis = 0.; - // xTau.auxdecor<float>("ptMvaTES") = float( interpolPt * reader->EvaluateRegression( 0, "BDTG" ) ); float ptMVA = float( interpolPt * reader->EvaluateRegression( 0, "BDTG" ) ); if(ptMVA<1) ptMVA=1; - xTau.auxdecor<float>("ptFinalCalib") = ptMVA; - xTau.auxdecor<float>("etaFinalCalib") = xTau.auxdata<float>("etaPanTauCellBased"); - xTau.auxdecor<float>("phiFinalCalib") = xTau.auxdata<float>("phiPanTauCellBased"); - xTau.auxdecor<float>("mFinalCalib") = 0; + xTau.setP4(xAOD::TauJetParameters::FinalCalib, ptMVA, xTau.etaPanTauCellBased(), xTau.phiPanTauCellBased(), 0); return StatusCode::SUCCESS; diff --git a/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx b/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx index bfaca23b3f9a5d7dd9d061a8f4400208f0a6e258..173bbb35c2b6c0efc35b6eee46bf5349ea74f8b2 100644 --- a/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx +++ b/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx @@ -6,7 +6,7 @@ #include "tauRecTools/MvaTESVariableDecorator.h" // tools include(s) -#include "TauAnalysisTools/HelperFunctions.h" +//#include "TauAnalysisTools/HelperFunctions.h" //_____________________________________________________________________________ MvaTESVariableDecorator::MvaTESVariableDecorator(const std::string& name) @@ -30,7 +30,7 @@ StatusCode MvaTESVariableDecorator::eventInitialize() m_mu = m_xEventInfo->averageInteractionsPerCrossing(); if(evtStore()->contains<xAOD::VertexContainer>("PrimaryVertices")){ - ATH_CHECK(evtStore()->retrieve(m_xVertexContainer, "PrimaryVertices")); + ATH_CHECK(evtStore()->retrieve(m_xVertexContainer, "PrimaryVertices")); m_nVtx = (int)m_xVertexContainer->size(); } else { @@ -45,6 +45,7 @@ StatusCode MvaTESVariableDecorator::eventInitialize() StatusCode MvaTESVariableDecorator::execute(xAOD::TauJet& xTau) { // Decorate event info + xTau.auxdata<double>("mu") = m_mu; xTau.auxdata<int>("nVtx") = m_nVtx; @@ -57,10 +58,10 @@ StatusCode MvaTESVariableDecorator::execute(xAOD::TauJet& xTau) { double clE=0., Etot=0.; TLorentzVector LC_P4; - LC_P4.SetPtEtaPhiM(xTau.auxdata<float>("LC_TES_precalib"), - xTau.auxdata<float>("seedCalo_eta"), - xTau.auxdata<float>("seedCalo_phi"), - xTau.m()); + LC_P4.SetPtEtaPhiM(xTau.ptDetectorAxis(), + xTau.etaDetectorAxis(), + xTau.phiDetectorAxis(), + xTau.m()); // ----loop over jet seed constituents xAOD::JetConstituentVector vec = jet_seed->getConstituents(); @@ -100,30 +101,27 @@ StatusCode MvaTESVariableDecorator::execute(xAOD::TauJet& xTau) { } // ----retrieve Ghost Muon Segment Count (for punch-through studies) - const int nMuSeg = jet_seed->getAttribute<int>("GhostMuonSegmentCount"); + int nMuSeg=0; + if(!jet_seed->getAttribute<int>("GhostMuonSegmentCount", nMuSeg)) nMuSeg=0; // ----decorating jet seed information to tau - xTau.auxdecor<double>("center_lambda") = mean_center_lambda; - xTau.auxdecor<double>("first_eng_dens") = mean_first_eng_dens; - xTau.auxdecor<double>("em_probability") = mean_em_probability; - xTau.auxdecor<double>("second_lambda") = mean_second_lambda; - xTau.auxdecor<double>("presampler_frac") = mean_presampler_frac; - xTau.auxdecor<int>("GhostMuonSegmentCount") = nMuSeg; + xTau.setDetail(xAOD::TauJetParameters::ClustersMeanCenterLambda, (float) mean_center_lambda); + xTau.setDetail(xAOD::TauJetParameters::ClustersMeanFirstEngDens, (float) mean_first_eng_dens); + xTau.setDetail(xAOD::TauJetParameters::ClustersMeanEMProbability, (float) mean_em_probability); + xTau.setDetail(xAOD::TauJetParameters::ClustersMeanSecondLambda, (float) mean_second_lambda); + xTau.setDetail(xAOD::TauJetParameters::ClustersMeanPresamplerFrac, (float) mean_presampler_frac); + xTau.setDetail(xAOD::TauJetParameters::GhostMuonSegmentCount, nMuSeg); // calculate PFO energy relative difference // ----summing corrected Pi0 PFO energies TLorentzVector Pi0_totalP4; Pi0_totalP4.SetPtEtaPhiM(0,0,0,0); - std::vector<TLorentzVector> Pi0PFOs; - TauAnalysisTools::createPi0Vectors(&xTau,Pi0PFOs); - - for(size_t i=0; i<Pi0PFOs.size(); i++){ - Pi0_totalP4 += Pi0PFOs.at(i); - }; + //This should be available in EDM as of TauJet_v3 + // TauAnalysisTools::createPi0Vectors(&xTau,Pi0PFOs); + for( size_t i=0; i != xTau.nPi0s(); ++i ) Pi0_totalP4+= xTau.pi0(i)->p4(); double Pi0_totalE = Pi0_totalP4.E(); - int nPi0PFOs = (int)Pi0PFOs.size(); // ----summing charged PFO energies TLorentzVector charged_totalP4; @@ -134,27 +132,24 @@ StatusCode MvaTESVariableDecorator::execute(xAOD::TauJet& xTau) { }; double charged_totalE = charged_totalP4.E(); - int nChargedPFOs = (int)xTau.nChargedPFOs(); // ----calculate relative difference and decorate to tau double relDiff=0.; if(Pi0_totalE+charged_totalE){ relDiff = (charged_totalE - Pi0_totalE) / (charged_totalE + Pi0_totalE) ; } - xTau.auxdecor<int>("nPi0PFOs") = nPi0PFOs; - xTau.auxdecor<int>("nChargedPFOs") = nChargedPFOs; - xTau.auxdecor<double>("PFOEngRelDiff") = relDiff; + xTau.setDetail(xAOD::TauJetParameters::PFOEngRelDiff, (float) relDiff); // calculate interpolated pT double GeV = 1000.; double pt_pantau = xTau.ptPanTauCellBased(); - double pt_LC = xTau.auxdata<float>("LC_TES_precalib"); + double pt_LC = xTau.ptDetectorAxis(); double interpolWeight; interpolWeight = 0.5 * ( 1. + TMath::TanH( ( pt_LC/GeV - 250. ) / 20. ) ); double LC_pantau_interpolPt = interpolWeight*pt_LC + (1.-interpolWeight)*pt_pantau; - xTau.auxdecor<double>("LC_pantau_interpolPt") = LC_pantau_interpolPt; + xTau.setDetail(xAOD::TauJetParameters::LC_pantau_interpolPt, (float) LC_pantau_interpolPt); return StatusCode::SUCCESS; diff --git a/Reconstruction/tauRecTools/Root/TauBuilderTool.cxx b/Reconstruction/tauRecTools/Root/TauBuilderTool.cxx index 239735e88a813f92568e64df768723f5ccd466e3..eb210a64ee6e1e60ce888933c53fae71f32a9acc 100644 --- a/Reconstruction/tauRecTools/Root/TauBuilderTool.cxx +++ b/Reconstruction/tauRecTools/Root/TauBuilderTool.cxx @@ -10,15 +10,18 @@ #include "xAODTau/TauJetContainer.h" #include "xAODTau/TauJetAuxContainer.h" #include "xAODTau/TauDefs.h" +#include "xAODTau/TauTrackContainer.h" +#include "xAODTau/TauTrackAuxContainer.h" #include "tauRecTools/TauEventData.h" - //________________________________________ TauBuilderTool::TauBuilderTool(const std::string& type) : asg::AsgTool(type), m_tauContainerName("TauJets"), m_tauAuxContainerName("TauJetsAux."), + m_tauTrackContainerName("TauTracks"), + m_tauTrackAuxContainerName("TauTracksAux."), m_seedContainerName(""), m_maxEta(2.5), m_minPt(10000), @@ -26,6 +29,8 @@ TauBuilderTool::TauBuilderTool(const std::string& type) : { declareProperty("TauContainer", m_tauContainerName); declareProperty("TauAuxContainer", m_tauAuxContainerName); + declareProperty("TauTrackContainer", m_tauTrackContainerName); + declareProperty("TauTrackAuxContainer", m_tauTrackAuxContainerName); declareProperty("SeedContainer", m_seedContainerName); declareProperty("MaxEta", m_maxEta); declareProperty("MinPt", m_minPt); @@ -95,7 +100,8 @@ StatusCode TauBuilderTool::execute(){ xAOD::TauJetContainer * pContainer = 0; xAOD::TauJetAuxContainer* pAuxContainer = 0; - + xAOD::TauTrackContainer* pTracks = 0; + xAOD::TauTrackAuxContainer* pAuxTracks = 0; if (m_doCreateTauContainers) { @@ -123,6 +129,13 @@ StatusCode TauBuilderTool::execute(){ ATH_MSG_DEBUG( "Recorded xAOD tau jets with key: " << m_tauAuxContainerName ); + + pTracks = new xAOD::TauTrackContainer(); + ATH_CHECK( evtStore()->record( pTracks, m_tauTrackContainerName) ); + pAuxTracks = new xAOD::TauTrackAuxContainer(); + ATH_CHECK( evtStore()->record( pAuxTracks, m_tauTrackAuxContainerName)); + pTracks->setStore( pAuxTracks ); + } else { //------------------------------------------------------------------------- // retrieve Tau Containers from StoreGate @@ -192,7 +205,8 @@ StatusCode TauBuilderTool::execute(){ ATH_MSG_VERBOSE("Number of seeds in the container: " << pSeedContainer->size()); - for (; itS != itSE; ++itS) { + for (; itS != itSE; ++itS) { + const xAOD::Jet *pSeed = (*itS); ATH_MSG_VERBOSE("Seeds eta:" << pSeed->eta() << ", pt:" << pSeed->pt()); @@ -261,11 +275,24 @@ StatusCode TauBuilderTool::execute(){ (*p_itT1)->cleanup(&m_data); (*p_itT1)->cleanup(&m_data); */ + + //remove orphaned tracks before tau is deleted via pop_back + xAOD::TauJet* bad_tau = pContainer->back(); + ATH_MSG_DEBUG("Deleting " << bad_tau->nAllTracks() << "Tracks associated with tau: "); + pTracks->erase(pTracks->end()-bad_tau->nAllTracks(), pTracks->end()); + //m_data.xAODTauContainer->pop_back(); pContainer->pop_back(); - } else + } else{ + + //remove orphaned tracks before tau is deleted via pop_back + xAOD::TauJet* bad_tau = pContainer->back(); + ATH_MSG_DEBUG("Deleting " << bad_tau->nAllTracks() << "Tracks associated with tau: "); + pTracks->erase(pTracks->end()-bad_tau->nAllTracks(), pTracks->end()); + //m_data.xAODTauContainer->pop_back(); pContainer->pop_back(); + } } @@ -282,6 +309,7 @@ StatusCode TauBuilderTool::execute(){ return StatusCode::FAILURE; } + //keep this here for future use (in case more than one seeding algo exist) /* p_itET = m_endTools.begin(); diff --git a/Reconstruction/tauRecTools/Root/TauCalibrateLC.cxx b/Reconstruction/tauRecTools/Root/TauCalibrateLC.cxx index e0dbb09b514ece12848abe924a551fe45e48ea00..adb81f2b6db00d6992e5ae8d884523baef18b92a 100644 --- a/Reconstruction/tauRecTools/Root/TauCalibrateLC.cxx +++ b/Reconstruction/tauRecTools/Root/TauCalibrateLC.cxx @@ -43,12 +43,6 @@ m_clusterCone(0.2) //not used declareProperty("isCaloOnly", m_isCaloOnly); } -TauCalibrateLC::TauCalibrateLC() : - TauCalibrateLC("TauCalibrateLC")//c++11 -{ -} - - /********************************************************************/ TauCalibrateLC::~TauCalibrateLC() { } @@ -73,7 +67,7 @@ StatusCode TauCalibrateLC::initialize() { etaBinHist = dynamic_cast<TH1 *> (obj); } if (etaBinHist) { - TH1 * tmp = const_cast<TH1*> (etaBinHist); + TH1 * tmp = dynamic_cast<TH1*> (obj); tmp->SetDirectory(0); } else { ATH_MSG_FATAL("Failed to get an object with key " << key); @@ -99,7 +93,7 @@ StatusCode TauCalibrateLC::initialize() { etaCorrectionHist = dynamic_cast<TH1 *> (obj); } if (etaCorrectionHist) { - TH1 * tmp = const_cast<TH1*> (etaCorrectionHist); + TH1 * tmp = dynamic_cast<TH1*> (obj); tmp->SetDirectory(0); } else { ATH_MSG_FATAL("Failed to get an object with key " << key); @@ -116,7 +110,7 @@ StatusCode TauCalibrateLC::initialize() { slopeNPVHist[i] = dynamic_cast<TH1 *> (obj); } if (slopeNPVHist[i]) { - TH1 * tmp = const_cast<TH1*> (slopeNPVHist[i]); + TH1 * tmp = dynamic_cast<TH1*> (obj); tmp->SetDirectory(0); } else { ATH_MSG_FATAL("Failed to get an object with key " << tmpSlopKey[i]); diff --git a/Reconstruction/tauRecTools/Root/TauCommonCalcVars.cxx b/Reconstruction/tauRecTools/Root/TauCommonCalcVars.cxx index c5f3d986560b36cad90cf491b397c97007d1024d..a8617e7f773ba4122a0643727dc5b3d71fa77c55 100644 --- a/Reconstruction/tauRecTools/Root/TauCommonCalcVars.cxx +++ b/Reconstruction/tauRecTools/Root/TauCommonCalcVars.cxx @@ -23,7 +23,7 @@ #include "tauRecTools/TauCommonCalcVars.h" #include "tauRecTools/KineUtils.h" - +#include <vector> //----------------------------------------------------------------------------- // Constructor //----------------------------------------------------------------------------- @@ -106,54 +106,45 @@ StatusCode TauCommonCalcVars::execute(xAOD::TauJet& pTau) { // } // invariant mass of track system - if ((pTau.nTracks() + pTau.nWideTracks()) > 0) { + std::vector<const xAOD::TauTrack*> tauTracks = pTau.tracks(xAOD::TauJetParameters::TauTrackFlag::classifiedCharged); + for( const xAOD::TauTrack* trk : pTau.tracks(xAOD::TauJetParameters::TauTrackFlag::classifiedIsolation) ) tauTracks.push_back(trk); + // if ((pTau.nTracks() + pTau.nWideTracks()) > 0) { + if (tauTracks.size()> 0) { TLorentzVector sumOfTrackVector; TLorentzVector tempTrackVector; - for (unsigned int i = 0; i != pTau.nTracks(); ++i) - { - tempTrackVector.SetPtEtaPhiM( pTau.track(i)->pt(), pTau.track(i)->eta(), pTau.track(i)->phi(), pTau.track(i)->m()); - sumOfTrackVector += tempTrackVector; - } - for (unsigned int i = 0; i != pTau.nWideTracks(); ++i) + for (const xAOD::TauTrack* tauTrk : tauTracks) { - tempTrackVector.SetPtEtaPhiM( pTau.wideTrack(i)->pt(), pTau.wideTrack(i)->eta(), pTau.wideTrack(i)->phi(), pTau.wideTrack(i)->m()); + tempTrackVector=tauTrk->p4(); sumOfTrackVector += tempTrackVector; } - pTau.setDetail( xAOD::TauJetParameters::massTrkSys, static_cast<float>( sumOfTrackVector.M() ) ); float tempfloat = 0; if ( pTau.detail( xAOD::TauJetParameters::massTrkSys, tempfloat ) ) - ATH_MSG_VERBOSE("set massTrkSys " << tempfloat); + ATH_MSG_VERBOSE("set massTrkSys " << tempfloat); } // width of track system squared (trkWidth2) - if (pTau.nTracks() > 1) { + // if (pTau.nTracks() > 1) { // originally looped over core+wide tracks below + if (tauTracks.size()> 0 && pTau.nTracks()>1) { - double leadTrkPhi = pTau.track(0)->phi(); - double leadTrkEta = pTau.track(0)->eta(); + double leadTrkPhi = pTau.track(0)->phi();//fix this depending on how we choose to define this + double leadTrkEta = pTau.track(0)->eta();//dito. new TauJet_v3 track sorting doesn't guarantee this will be the same track double ptSum = 0; double sumWeightedDR = 0; double sumWeightedDR2 = 0; - for (unsigned int i = 0; i != pTau.nTracks(); ++i) { - double deltaR = Tau1P3PKineUtils::deltaR(leadTrkEta, leadTrkPhi, pTau.track(i)->eta(), pTau.track(i)->phi() ); + for (const xAOD::TauTrack* tauTrk : tauTracks){ - ptSum += pTau.track(i)->pt(); - sumWeightedDR += deltaR * (pTau.track(i)->pt()); - sumWeightedDR2 += deltaR * deltaR * (pTau.track(i)->pt()); - } + double deltaR = Tau1P3PKineUtils::deltaR(leadTrkEta, leadTrkPhi, tauTrk->eta(), tauTrk->phi() ); - for (unsigned int i = 0; i != pTau.nWideTracks(); ++i) { + ptSum += tauTrk->pt(); + sumWeightedDR += deltaR * (tauTrk->pt()); + sumWeightedDR2 += deltaR * deltaR * (tauTrk->pt()); - double deltaR = Tau1P3PKineUtils::deltaR(leadTrkEta, leadTrkPhi, pTau.wideTrack(i)->eta(), pTau.wideTrack(i)->phi() ); - - ptSum += pTau.wideTrack(i)->pt(); - sumWeightedDR += deltaR * (pTau.wideTrack(i)->pt()); - sumWeightedDR2 += deltaR * deltaR * (pTau.wideTrack(i)->pt()); } double trkWidth2 = sumWeightedDR2 / ptSum - sumWeightedDR * sumWeightedDR / ptSum / ptSum; @@ -166,8 +157,9 @@ StatusCode TauCommonCalcVars::execute(xAOD::TauJet& pTau) { //FF: use now the 4-vector of the tau intermediate axis //P4EEtaPhiM P4CaloSeed(1., pDetails->seedCalo_eta(), pDetails->seedCalo_phi(), 0.); - - if ((pTau.nWideTracks() + pTau.nTracks()) > 0) { + + // if ((pTau.nWideTracks() + pTau.nTracks()) > 0) { + if (tauTracks.size()> 0) { double ptSum = 0; double innerPtSum = 0; @@ -175,26 +167,20 @@ StatusCode TauCommonCalcVars::execute(xAOD::TauJet& pTau) { double innerSumWeightedDR = 0; double sumWeightedDR2 = 0; - for (unsigned int i = 0; i != pTau.nTracks(); ++i) { - double deltaR = Tau1P3PKineUtils::deltaR(pTau.eta(), pTau.phi(), pTau.track(i)->eta(), pTau.track(i)->phi() ); + for (const xAOD::TauTrack* tauTrk : tauTracks){ + + double deltaR = Tau1P3PKineUtils::deltaR(pTau.eta(), pTau.phi(), tauTrk->eta(), tauTrk->phi() ); - ptSum += pTau.track(i)->pt(); - sumWeightedDR += deltaR * (pTau.track(i)->pt()); - sumWeightedDR2 += deltaR * deltaR * (pTau.track(i)->pt()); + ptSum += tauTrk->pt(); + sumWeightedDR += deltaR * (tauTrk->pt()); + sumWeightedDR2 += deltaR * deltaR * (tauTrk->pt()); //add calculation of innerTrkAvgDist - innerPtSum += pTau.track(i)->pt(); - innerSumWeightedDR += deltaR * (pTau.track(i)->pt()); - } - - for (unsigned int i = 0; i != pTau.nWideTracks(); ++i) { - - double deltaR = Tau1P3PKineUtils::deltaR(pTau.eta(), pTau.phi(), pTau.wideTrack(i)->eta(), pTau.wideTrack(i)->phi() ); - - ptSum += pTau.wideTrack(i)->pt(); - sumWeightedDR += deltaR * (pTau.wideTrack(i)->pt()); - sumWeightedDR2 += deltaR * deltaR * (pTau.wideTrack(i)->pt()); + if(tauTrk->flag(xAOD::TauJetParameters::TauTrackFlag::classifiedCharged)){ + innerPtSum += tauTrk->pt(); + innerSumWeightedDR += deltaR * (tauTrk->pt()); + } } if (ptSum > 0.0001) { diff --git a/Reconstruction/tauRecTools/Root/TauProcessorTool.cxx b/Reconstruction/tauRecTools/Root/TauProcessorTool.cxx index b304805894a6a03fe090637b2556a4a132a19346..d7928222168fea42fa63bf5dd324fc37a4da0894 100644 --- a/Reconstruction/tauRecTools/Root/TauProcessorTool.cxx +++ b/Reconstruction/tauRecTools/Root/TauProcessorTool.cxx @@ -5,6 +5,8 @@ #include "tauRecTools/TauProcessorTool.h" #include "xAODTau/TauJetContainer.h" +#include "xAODTau/TauTrackContainer.h" +#include "xAODTau/TauTrackAuxContainer.h" #include "xAODTracking/VertexContainer.h" #include "xAODTracking/VertexAuxContainer.h" #include "xAODPFlow/PFOContainer.h" @@ -41,6 +43,7 @@ TauProcessorTool::TauProcessorTool(const std::string& type) : declareProperty("deepCopyHadronicPFOContainer", m_deep_copy_hadronicPFOContainer=true); declareProperty("deepCopyNeutralPFOContainer", m_deep_copy_neutralPFOContainer=true); declareProperty("deepCopySecVtxContainer", m_deep_copy_SecVtxContainer=false); + declareProperty("deepCopyTauTrackContainer", m_deep_copy_TauTrackContainer=true); } //________________________________________ @@ -140,6 +143,12 @@ StatusCode TauProcessorTool::execute(){ xAOD::Vertex* v(0); ATH_CHECK(deepCopy(pSecVtxContainer, pSecVtxAuxContainer, v, "TauSecondaryVertices")); } + if(m_deep_copy_TauTrackContainer){ + xAOD::TauTrackContainer* pTrackContainer(0); + xAOD::TauTrackAuxContainer* pTauTrackAuxContainer(0); + xAOD::TauTrack* v(0); + ATH_CHECK(deepCopy(pTrackContainer, pTauTrackAuxContainer, v, "TauTracks")); + } if(m_deep_copy_neutralPFOContainer){ xAOD::PFOContainer* neutralPFOContainer(0); diff --git a/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx b/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx index 89e740aac7f393124bab192540fc767072cb5192..cbfce175b1b54578ef9bf537b9b66b3a1af50aae 100644 --- a/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx +++ b/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx @@ -110,7 +110,7 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) { if (!taujetseed) ATH_MSG_DEBUG("Taujet->jet() pointer is NULL: calo cluster variables will be set to -1111"); else ATH_MSG_DEBUG("problem in calculating calo cluster variables -> will be set to -1111"); - pTau.setDetail(xAOD::TauJetParameters::numCells , static_cast<int>(0) ); + if(!m_inAODmode) pTau.setDetail(xAOD::TauJetParameters::numCells , static_cast<int>(0) ); pTau.setDetail(xAOD::TauJetParameters::numTopoClusters , static_cast<int>(DEFAULT) ); pTau.setDetail(xAOD::TauJetParameters::numEffTopoClusters , static_cast<float>(DEFAULT) ); pTau.setDetail(xAOD::TauJetParameters::topoInvMass, static_cast<float>(DEFAULT) ); @@ -312,7 +312,7 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) { } // set new approximate energy flow variables for tau ID - pTau.setDetail(xAOD::TauJetParameters::ptRatioEflowApprox, static_cast<float>(approxSubstructure4Vec.Pt()/ pTau.detail<float>(xAOD::TauJetParameters::LC_TES_precalib)) ); + pTau.setDetail(xAOD::TauJetParameters::ptRatioEflowApprox, static_cast<float>(approxSubstructure4Vec.Pt()/ pTau.ptDetectorAxis()) ); pTau.setDetail(xAOD::TauJetParameters::mEflowApprox, static_cast<float>(approxSubstructure4Vec.M()) ); diff --git a/Reconstruction/tauRecTools/Root/TauTrackClassifier.cxx b/Reconstruction/tauRecTools/Root/TauTrackClassifier.cxx new file mode 100644 index 0000000000000000000000000000000000000000..d00e9b3629f0405e117562af1dc08ccf7cbbb7d5 --- /dev/null +++ b/Reconstruction/tauRecTools/Root/TauTrackClassifier.cxx @@ -0,0 +1,270 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +// ASG include(s) +#include "PathResolver/PathResolver.h" + +// xAOD include(s) +#include "xAODTracking/TrackParticle.h" +#include "xAODTau/TauTrackContainer.h" +#include "xAODTau/TauxAODHelpers.h" + +// local include(s) +#include "tauRecTools/TauTrackClassifier.h" + +using namespace tauRecTools; + +//============================================================================== +// class TauTrackClassifier +//============================================================================== + +//______________________________________________________________________________ +TauTrackClassifier::TauTrackClassifier(const std::string& sName) + : TauRecToolBase(sName) +{ + declareProperty("Classifiers", m_vClassifier ); + declareProperty("TauTrackContainerName", m_tauTrackConName="TauTracks"); +} + +//______________________________________________________________________________ +TauTrackClassifier::~TauTrackClassifier() +{ +} + +//______________________________________________________________________________ +StatusCode TauTrackClassifier::initialize() +{ + ATH_MSG_DEBUG("intialize classifiers"); + + for (auto cClassifier : m_vClassifier){ + ATH_MSG_INFO("TauTrackClassifier tool : " << cClassifier ); + ATH_CHECK(cClassifier.retrieve()); + } + return StatusCode::SUCCESS; +} + +//______________________________________________________________________________ +StatusCode TauTrackClassifier::execute(xAOD::TauJet& xTau) +{ + xAOD::TauTrackContainer* tauTrackCon = 0; + ATH_CHECK(evtStore()->retrieve(tauTrackCon, m_tauTrackConName)); + std::vector<xAOD::TauTrack*> vTracks = xAOD::TauHelpers::allTauTracksNonConst(&xTau, tauTrackCon); + for (xAOD::TauTrack* xTrack : vTracks) + { + // reset all track flags and set status to unclassified + xTrack->setFlag(xAOD::TauJetParameters::classifiedCharged, false); + xTrack->setFlag(xAOD::TauJetParameters::classifiedConversion, false); + xTrack->setFlag(xAOD::TauJetParameters::classifiedIsolation, false); + xTrack->setFlag(xAOD::TauJetParameters::classifiedFake, false); + xTrack->setFlag(xAOD::TauJetParameters::unclassified, true); + + // execute the bdt track classifier + for (auto cClassifier : m_vClassifier) + CHECK(cClassifier->classifyTrack(*xTrack, xTau)); + } + std::vector< ElementLink< xAOD::TauTrackContainer > > &tauTrackLinks(xTau.tauTrackLinksNonConst()); + std::sort(tauTrackLinks.begin(), tauTrackLinks.end(), sortTracks); + float charge=0.0; + for( const xAOD::TauTrack* trk : xTau.tracks(xAOD::TauJetParameters::classifiedCharged) ){ + charge += trk->track()->charge(); + } + xTau.setCharge(charge); + return StatusCode::SUCCESS; +} + +//============================================================================== +// class TrackMVABDT +//============================================================================== + +//______________________________________________________________________________ +TrackMVABDT::TrackMVABDT(const std::string& sName) + : TauRecToolBase(sName) + , m_sInputWeightsPath("") + , m_fThreshold(0.) + , m_iSignalType(xAOD::TauJetParameters::classifiedCharged) + , m_iBackgroundType(xAOD::TauJetParameters::classifiedFake) + , m_iExpectedFlag(xAOD::TauJetParameters::unclassified) + , m_rReader(0) +{ + declareProperty( "InputWeightsPath", m_sInputWeightsPath ); + declareProperty( "Threshold", m_fThreshold ); + declareProperty( "BackgroundType" , m_iBackgroundType ); + declareProperty( "SignalType", m_iSignalType ); + declareProperty( "ExpectedFlag", m_iExpectedFlag ); + m_mParsedVarsBDT.clear(); +} + +//______________________________________________________________________________ +TrackMVABDT::~TrackMVABDT() +{ + delete m_rReader; +} + +//______________________________________________________________________________ +StatusCode TrackMVABDT::initialize() +{ + m_rReader = new TMVA::Reader( "!Silent" ); + m_mAvailableVars= + { + {"TracksAuxDyn.tauPt",tauPt} + , {"TracksAuxDyn.tauEta",tauEta} + , {"TracksAuxDyn.trackEta",trackEta} + , {"TracksAuxDyn.z0sinThetaTJVA",z0sinThetaTJVA} + , {"TracksAuxDyn.rConv",rConv} + , {"TracksAuxDyn.rConvII",rConvII} + , {"TracksAuxDyn.DRJetSeedAxis",DRJetSeedAxis} + , {"TracksAux.d0",d0} + , {"TracksAux.qOverP",qOverP} + , {"TracksAux.theta",theta} + , {"TracksAux.numberOfInnermostPixelLayerHits", numberOfInnermostPixelLayerHits} + , {"TracksAux.numberOfPixelHits",numberOfPixelHits} + , {"TracksAux.numberOfPixelDeadSensors",numberOfPixelDeadSensors} + , {"TracksAux.numberOfPixelSharedHits",numberOfPixelSharedHits} + , {"TracksAux.numberOfSCTHits",numberOfSCTHits} + , {"TracksAux.numberOfSCTDeadSensors",numberOfSCTDeadSensors} + , {"TracksAux.numberOfSCTSharedHits",numberOfSCTSharedHits} + , {"TracksAux.numberOfTRTHighThresholdHits",numberOfTRTHighThresholdHits} + , {"TracksAux.numberOfTRTHits",numberOfTRTHits} + , {"TracksAux.numberOfPixelHits+TracksAux.numberOfPixelDeadSensors", nPixHits} + , {"TracksAux.numberOfPixelHits+TracksAux.numberOfPixelDeadSensors+TracksAux.numberOfSCTHits+TracksAux.numberOfSCTDeadSensors", nSiHits} + }; + + for (auto fVar : m_mAvailableVars) + fVar.second = 0.; + + ATH_CHECK(addWeightsFile()); + + return StatusCode::SUCCESS; +} + +//______________________________________________________________________________ +StatusCode TrackMVABDT::classifyTrack(xAOD::TauTrack& xTrack, const xAOD::TauJet& xTau) +{ + setVars(xTrack, xTau); + + //why? + if (xTrack.flag((xAOD::TauJetParameters::TauTrackFlag) m_iExpectedFlag)==false) + return StatusCode::SUCCESS; + + double dValue = m_rReader->EvaluateMVA( m_sInputWeightsPath ); + + xTrack.setFlag((xAOD::TauJetParameters::TauTrackFlag) m_iExpectedFlag, false); + if (m_fThreshold < dValue) + xTrack.setFlag((xAOD::TauJetParameters::TauTrackFlag) m_iSignalType, true); + else + xTrack.setFlag((xAOD::TauJetParameters::TauTrackFlag) m_iBackgroundType, true); + + xTrack.addBdtScore(dValue); + + return StatusCode::SUCCESS; +} + +//______________________________________________________________________________ +StatusCode TrackMVABDT::addWeightsFile() +{ + // better to replace that part at some point + m_sInputWeightsPath = PathResolverFindCalibFile(m_sInputWeightsPath); + ATH_MSG_DEBUG("InputWeightsPath: " << m_sInputWeightsPath); + + if (m_mParsedVarsBDT.empty()) + { + ATH_CHECK(parseVariableContent()); + + for (auto i : m_mParsedVarsBDT) + ATH_MSG_DEBUG(i.first<<" "<<i.second); + + for (size_t i = 0; i < m_mParsedVarsBDT.size(); i++){ + std::string sVarName = m_mParsedVarsBDT[i]; + m_rReader->AddVariable( sVarName, &(m_mAvailableVars[sVarName])); + } + } + + m_rReader->BookMVA( m_sInputWeightsPath, m_sInputWeightsPath ); + return StatusCode::SUCCESS; +} + +//______________________________________________________________________________ +StatusCode TrackMVABDT::parseVariableContent() +{ + // example <Variable VarIndex="0" Expression="TracksAuxDyn.tauPt" Label="TracksAuxDyn.tauPt" Title="tauPt" Unit="" Internal="TracksAuxDyn.tauPt" Type="F" Min="1.50000762e+04" Max="5.32858625e+05"/> + std::string sLine; + std::ifstream sFileStream (m_sInputWeightsPath); + if (sFileStream.is_open()) + { + while ( getline (sFileStream,sLine) ) + { + size_t iPosVarindex = sLine.find("VarIndex="); + + if ( iPosVarindex == std::string::npos ) + continue; + + iPosVarindex += 10; + + size_t iPosVarindexEnd = sLine.find("\"",iPosVarindex); + size_t iPosExpression = sLine.find("Expression=")+12; + size_t iPosExpressionEnd = sLine.find("\"",iPosExpression); + + int iVarIndex = std::stoi(sLine.substr(iPosVarindex, iPosVarindexEnd-iPosVarindex)); + std::string sExpression = sLine.substr(iPosExpression, iPosExpressionEnd-iPosExpression); + + m_mParsedVarsBDT.insert(std::pair<int, std::string >(iVarIndex,sExpression)); + } + sFileStream.close(); + return StatusCode::SUCCESS; + } + + ATH_MSG_ERROR("Unable to open file "<<m_sInputWeightsPath); + return StatusCode::FAILURE; +} + +//______________________________________________________________________________ +void TrackMVABDT::setVars(const xAOD::TauTrack& xTrack, const xAOD::TauJet& xTau) +{ + const xAOD::TrackParticle* xTrackParticle = xTrack.track(); + uint8_t iTracksNumberOfInnermostPixelLayerHits = 0; xTrackParticle->summaryValue(iTracksNumberOfInnermostPixelLayerHits, xAOD::numberOfInnermostPixelLayerHits); + uint8_t iTracksNPixelHits = 0; xTrackParticle->summaryValue(iTracksNPixelHits, xAOD::numberOfPixelHits); + uint8_t iTracksNPixelSharedHits = 0; xTrackParticle->summaryValue(iTracksNPixelSharedHits, xAOD::numberOfPixelSharedHits); + uint8_t iTracksNPixelDeadSensors = 0; xTrackParticle->summaryValue(iTracksNPixelDeadSensors, xAOD::numberOfPixelDeadSensors); + uint8_t iTracksNSCTHits = 0; xTrackParticle->summaryValue(iTracksNSCTHits, xAOD::numberOfSCTHits); + uint8_t iTracksNSCTSharedHits = 0; xTrackParticle->summaryValue(iTracksNSCTSharedHits, xAOD::numberOfSCTSharedHits); + uint8_t iTracksNSCTDeadSensors = 0; xTrackParticle->summaryValue(iTracksNSCTDeadSensors, xAOD::numberOfSCTDeadSensors); + uint8_t iTracksNTRTHighThresholdHits = 0; xTrackParticle->summaryValue( iTracksNTRTHighThresholdHits, xAOD::numberOfTRTHighThresholdHits); + uint8_t iTracksNTRTHits = 0; xTrackParticle->summaryValue( iTracksNTRTHits, xAOD::numberOfTRTHits); + + float fTracksNumberOfInnermostPixelLayerHits = (float)iTracksNumberOfInnermostPixelLayerHits; + float fTracksNPixelHits = (float)iTracksNPixelHits; + float fTracksNPixelDeadSensors = (float)iTracksNPixelDeadSensors; + float fTracksNPixelSharedHits = (float)iTracksNPixelSharedHits; + float fTracksNSCTHits = (float)iTracksNSCTHits; + float fTracksNSCTDeadSensors = (float)iTracksNSCTDeadSensors; + float fTracksNSCTSharedHits = (float)iTracksNSCTSharedHits; + float fTracksNTRTHighThresholdHits = (float)iTracksNTRTHighThresholdHits; + float fTracksNTRTHits = (float)iTracksNTRTHits; + + float fTracksNPixHits = fTracksNPixelHits + fTracksNPixelDeadSensors; + float fTracksNSiHits = fTracksNPixelHits + fTracksNPixelDeadSensors + fTracksNSCTHits + fTracksNSCTDeadSensors; + + + m_mAvailableVars["TracksAuxDyn.tauPt"] = xTau.pt(); + m_mAvailableVars["TracksAuxDyn.tauEta"] = xTau.eta(); + m_mAvailableVars["TracksAuxDyn.z0sinThetaTJVA"] = xTrack.z0sinThetaTJVA(xTau); + m_mAvailableVars["TracksAuxDyn.rConv"] = xTrack.rConv(xTau); + m_mAvailableVars["TracksAuxDyn.rConvII"] = xTrack.rConvII(xTau); + m_mAvailableVars["TracksAuxDyn.DRJetSeedAxis"] = xTrack.dRJetSeedAxis(xTau); + m_mAvailableVars["TracksAuxDyn.trackEta"] = xTrackParticle->d0(); + m_mAvailableVars["TracksAux.d0"] = xTrackParticle->d0(); + m_mAvailableVars["TracksAux.qOverP"] = xTrackParticle->qOverP(); + m_mAvailableVars["TracksAux.theta"] = xTrackParticle->theta(); + m_mAvailableVars["TracksAux.numberOfInnermostPixelLayerHits"] = fTracksNumberOfInnermostPixelLayerHits; + m_mAvailableVars["TracksAux.numberOfPixelHits"] = fTracksNPixelHits; + m_mAvailableVars["TracksAux.numberOfPixelDeadSensors"] = fTracksNPixelDeadSensors; + m_mAvailableVars["TracksAux.numberOfPixelSharedHits"] = fTracksNPixelSharedHits; + m_mAvailableVars["TracksAux.numberOfSCTHits"] = fTracksNSCTHits; + m_mAvailableVars["TracksAux.numberOfSCTDeadSensors"] = fTracksNSCTDeadSensors; + m_mAvailableVars["TracksAux.numberOfSCTSharedHits"] = fTracksNSCTSharedHits; + m_mAvailableVars["TracksAux.numberOfTRTHighThresholdHits"] = fTracksNTRTHighThresholdHits; + m_mAvailableVars["TracksAux.numberOfTRTHits"] = fTracksNTRTHits; + m_mAvailableVars["TracksAux.numberOfPixelHits+TracksAux.numberOfPixelDeadSensors"] = fTracksNPixHits; + m_mAvailableVars["TracksAux.numberOfPixelHits+TracksAux.numberOfPixelDeadSensors+TracksAux.numberOfSCTHits+TracksAux.numberOfSCTDeadSensors"] = fTracksNSiHits; +} diff --git a/Reconstruction/tauRecTools/Root/TauTrackFilter.cxx b/Reconstruction/tauRecTools/Root/TauTrackFilter.cxx index df2cf026b42b7e9dfcc76d3a5114d0665e232791..d89b61b185c6c2c14af934f0d932221792519d0b 100644 --- a/Reconstruction/tauRecTools/Root/TauTrackFilter.cxx +++ b/Reconstruction/tauRecTools/Root/TauTrackFilter.cxx @@ -23,18 +23,19 @@ #include "tauRecTools/TauTrackFilterUtils.h" #include "xAODTracking/TrackParticleContainer.h" +#include "xAODTau/TauxAODHelpers.h" #include "TLorentzVector.h" void TrackFilterAlg(TLorentzVector tau, - std::vector<TLorentzVector>* inputtracks20, - std::vector<int>* inputtracks20charge, - std::vector<TLorentzVector>* inputtracks40, - std::vector<int>* inputtracks40charge, - std::vector<TLorentzVector>* outputtracksgood, - std::vector<int>* outputtracksgoodcharge, - std::vector<TLorentzVector>* outputtracksbad, - std::vector<int>* outputtracksbadcharge, + std::vector<TLorentzVector>& inputtracks20, + std::vector<int>& inputtracks20charge, + std::vector<TLorentzVector>& inputtracks40, + std::vector<int>& inputtracks40charge, + std::vector<TLorentzVector>& outputtracksgood, + std::vector<int>& outputtracksgoodcharge, + std::vector<TLorentzVector>& outputtracksbad, + std::vector<int>& outputtracksbadcharge, int& nProng, int& flag); @@ -49,6 +50,7 @@ TauTrackFilter::TauTrackFilter(const std::string &name ) : { declareProperty("ConfigPath", m_configPath); declareProperty("TrackContainerName", m_trackContainerName = "InDetTrackParticles"); + declareProperty("TauTrackContainerName", m_tauTrackConName = "TauTracks"); } //----------------------------------------------------------------------------- @@ -104,14 +106,14 @@ StatusCode TauTrackFilter::execute(xAOD::TauJet& pTau) { pTau.phi(), pTau.e()/1000); //GeV - std::vector<TLorentzVector>* inputtracks20 = new std::vector<TLorentzVector>; - std::vector<TLorentzVector>* inputtracks40 = new std::vector<TLorentzVector>; - std::vector<int>* inputtracks20charge = new std::vector<int>; - std::vector<int>* inputtracks40charge = new std::vector<int>; - std::vector<TLorentzVector>* outputtracksgood = new std::vector<TLorentzVector>; - std::vector<TLorentzVector>* outputtracksbad = new std::vector<TLorentzVector>; - std::vector<int>* outputtracksgoodcharge = new std::vector<int>; - std::vector<int>* outputtracksbadcharge = new std::vector<int>; + std::vector<TLorentzVector> inputtracks20; + std::vector<TLorentzVector> inputtracks40; + std::vector<int> inputtracks20charge; + std::vector<int> inputtracks40charge; + std::vector<TLorentzVector> outputtracksgood; + std::vector<TLorentzVector> outputtracksbad; + std::vector<int> outputtracksgoodcharge; + std::vector<int> outputtracksbadcharge; int nProng = 0; int flag = 0; @@ -121,7 +123,7 @@ StatusCode TauTrackFilter::execute(xAOD::TauJet& pTau) { m_TrkPass.clear(); for(unsigned int j=0; j<pTau.nTracks(); j++ ) { - const xAOD::TrackParticle *TauJetTrack = pTau.track(j); + const xAOD::TrackParticle *TauJetTrack = pTau.track(j)->track(); TLorentzVector inputTrack; inputTrack.SetPtEtaPhiE(TauJetTrack->pt()/1000, //GeV TauJetTrack->eta(), @@ -129,13 +131,13 @@ StatusCode TauTrackFilter::execute(xAOD::TauJet& pTau) { TauJetTrack->e()/1000); //GeV float dR = tau.DeltaR(inputTrack); if (dR < 0.2) { - inputtracks20->push_back(inputTrack); - inputtracks20charge->push_back(TauJetTrack->charge()); + inputtracks20.push_back(inputTrack); + inputtracks20charge.push_back(TauJetTrack->charge()); inputtracks20index.push_back(j); } else if (dR < 0.4) { - inputtracks40->push_back(inputTrack); - inputtracks40charge->push_back(TauJetTrack->charge()); + inputtracks40.push_back(inputTrack); + inputtracks40charge.push_back(TauJetTrack->charge()); inputtracks40index.push_back(j); } @@ -157,14 +159,14 @@ StatusCode TauTrackFilter::execute(xAOD::TauJet& pTau) { flag); // Store results - for (unsigned int j=0; j<outputtracksgood->size(); j++ ) { - for (unsigned int k=0; k<inputtracks20->size(); k++ ) { - if (outputtracksgood->at(j) == inputtracks20->at(k)) { + for (unsigned int j=0; j<outputtracksgood.size(); j++ ) { + for (unsigned int k=0; k<inputtracks20.size(); k++ ) { + if (outputtracksgood.at(j) == inputtracks20.at(k)) { m_TrkPass.at(inputtracks20index.at(k)) = true; } } - for (unsigned int k=0; k<inputtracks40->size(); k++ ) { - if (outputtracksgood->at(j) == inputtracks40->at(k)) { + for (unsigned int k=0; k<inputtracks40.size(); k++ ) { + if (outputtracksgood.at(j) == inputtracks40.at(k)) { m_TrkPass.at(inputtracks40index.at(k)) = true; } } @@ -172,19 +174,12 @@ StatusCode TauTrackFilter::execute(xAOD::TauJet& pTau) { m_nProng = nProng; m_flag = flag; - // Cleanup - delete inputtracks20; - delete inputtracks20charge; - delete inputtracks40; - delete inputtracks40charge; - delete outputtracksgood; - delete outputtracksbad; - delete outputtracksgoodcharge; - delete outputtracksbadcharge; - // Set values in EDM + xAOD::TauTrackContainer* tauTracks = 0; + ATH_CHECK(evtStore()->retrieve(tauTracks, m_tauTrackConName)); for (unsigned int numTrack=0; numTrack<m_TrkPass.size(); numTrack++) { - pTau.setTrackFlag(pTau.track(numTrack), xAOD::TauJetParameters::failTrackFilter, !m_TrkPass.at(numTrack)); + xAOD::TauTrack* tauTrk = xAOD::TauHelpers::tauTrackNonConst(&pTau, tauTracks, numTrack); //pTau.trackNonConst(numTrack); + tauTrk->setFlag(xAOD::TauJetParameters::failTrackFilter, !m_TrkPass.at(numTrack)); } pTau.setTrackFilterProngs(m_nProng); pTau.setTrackFilterQuality(m_flag); @@ -197,24 +192,24 @@ StatusCode TauTrackFilter::execute(xAOD::TauJet& pTau) { // Main algorithm //----------------------------------------------------------------------------- void TrackFilterAlg(TLorentzVector tau, - std::vector<TLorentzVector>* inputtracks20, - std::vector<int>* inputtracks20charge, - std::vector<TLorentzVector>* inputtracks40, - std::vector<int>* inputtracks40charge, - std::vector<TLorentzVector>* outputtracksgood, - std::vector<int>* outputtracksgoodcharge, - std::vector<TLorentzVector>* outputtracksbad, - std::vector<int>* outputtracksbadcharge, + std::vector<TLorentzVector>& inputtracks20, + std::vector<int>& inputtracks20charge, + std::vector<TLorentzVector>& inputtracks40, + std::vector<int>& inputtracks40charge, + std::vector<TLorentzVector>& outputtracksgood, + std::vector<int>& outputtracksgoodcharge, + std::vector<TLorentzVector>& outputtracksbad, + std::vector<int>& outputtracksbadcharge, int& nProng, int& flag) { std::vector<TauTrackFilterUtils::TrackInfo> unsorted,tracks,SScombination; TauTrackFilterUtils::TrackInfo track; - unsigned int tracknum = inputtracks20->size(); - unsigned int widetracknum = inputtracks40->size(); + unsigned int tracknum = inputtracks20.size(); + unsigned int widetracknum = inputtracks40.size(); for(unsigned int i=0;i<tracknum;i++){ - track.p4 = (*inputtracks20)[i]; - track.charge = (*inputtracks20charge)[i]; + track.p4 = (inputtracks20)[i]; + track.charge = (inputtracks20charge)[i]; unsorted.push_back(track); } while (unsorted.size() > 0){ @@ -235,8 +230,8 @@ void TrackFilterAlg(TLorentzVector tau, bool test3prong = true, test2prong = true, test1prong = true; for(unsigned int i=0;i<widetracknum;i++){ - outputtracksbad->push_back((*inputtracks40)[i]); - outputtracksbadcharge->push_back((*inputtracks40charge)[i]); + outputtracksbad.push_back((inputtracks40)[i]); + outputtracksbadcharge.push_back((inputtracks40charge)[i]); } if(tracknum > 4){ //Anything with more than 4 tracks is a fake. flag = 0; @@ -283,8 +278,8 @@ void TrackFilterAlg(TLorentzVector tau, if(goodcombo){ //A Combination is found which passes 3prong hypothesis for(unsigned int i=0;i<combination.size();i++){ if (isSS) SScombination.push_back(combination[i]); - outputtracksgood->push_back(combination[i].p4); - outputtracksgoodcharge->push_back(combination[i].charge); + outputtracksgood.push_back(combination[i].p4); + outputtracksgoodcharge.push_back(combination[i].charge); } if(isSS) flag = 0; else flag = 1; @@ -293,8 +288,8 @@ void TrackFilterAlg(TLorentzVector tau, test2prong = false; if(!test1prong){ //Fill Bad Track in the Case of 4 trk taus if(tracknum == 4){ - outputtracksbad->push_back(tracks[3].p4); - outputtracksbadcharge->push_back(tracks[3].charge); + outputtracksbad.push_back(tracks[3].p4); + outputtracksbadcharge.push_back(tracks[3].charge); } } } @@ -305,14 +300,14 @@ void TrackFilterAlg(TLorentzVector tau, if(TauTrackFilterUtils::pass2prong(pair,tau)){ nProng = 2; for(unsigned int i=0;i<pair.size();i++){ //Fill Good Tracks - outputtracksgood->push_back(pair[i].p4); - outputtracksgoodcharge->push_back(pair[i].charge); + outputtracksgood.push_back(pair[i].p4); + outputtracksgoodcharge.push_back(pair[i].charge); } test1prong = false; if(tracknum == 3){ flag = 2; - outputtracksbad->push_back(tracks[2].p4); //Fill Bad Track in Case of 3 trk Taus - outputtracksbadcharge->push_back(tracks[2].charge); + outputtracksbad.push_back(tracks[2].p4); //Fill Bad Track in Case of 3 trk Taus + outputtracksbadcharge.push_back(tracks[2].charge); } else flag = 1; //Good 2 Prong if only 2 trks } @@ -338,13 +333,13 @@ void TrackFilterAlg(TLorentzVector tau, if(tracks[1].p4.Pt()/tracks[0].p4.Pt() < ratio10) goodcase = true; //Test 2trk taus most likely to actually be 1pngs } if((TauTrackFilterUtils::pass1prong(tracks[0].p4,tau))&&(goodcase)){ //A track is found which passes 1prong hypothesis - outputtracksgood->push_back(tracks[0].p4); - outputtracksgoodcharge->push_back(tracks[0].charge); + outputtracksgood.push_back(tracks[0].p4); + outputtracksgoodcharge.push_back(tracks[0].charge); nProng = 1; if (tracknum == 2){ flag = 2; - outputtracksbad->push_back(tracks[1].p4); //Fill Bad Track in Case of 3 trk Taus - outputtracksbadcharge->push_back(tracks[1].charge); + outputtracksbad.push_back(tracks[1].p4); //Fill Bad Track in Case of 3 trk Taus + outputtracksbadcharge.push_back(tracks[1].charge); } else flag = 1; } diff --git a/Reconstruction/tauRecTools/cmt/requirements b/Reconstruction/tauRecTools/cmt/requirements index 1ba42c91abf4feab380aba43a681da3cda4cb710..3c53df7bfca38477a69eb1dbbd1158d63cf68472 100644 --- a/Reconstruction/tauRecTools/cmt/requirements +++ b/Reconstruction/tauRecTools/cmt/requirements @@ -51,7 +51,7 @@ use VxVertex VxVertex-* Tracking/TrkEven #use tauEvent tauEvent-* Reconstruction use xAODCaloEvent xAODCaloEvent-* Event/xAOD use xAODPFlow xAODPFlow-* Event/xAOD -use TauAnalysisTools TauAnalysisTools-* PhysicsAnalysis/TauID +use InDetTrackSelectionTool InDetTrackSelectionTool-* InnerDetector/InDetRecTools diff --git a/Reconstruction/tauRecTools/share/TF2pileupForOfflineID_1p.root b/Reconstruction/tauRecTools/share/TF2pileupForOfflineID_1p.root deleted file mode 100644 index 9c291900a685fdad8f49e16df169aa7b6b1dabec..0000000000000000000000000000000000000000 Binary files a/Reconstruction/tauRecTools/share/TF2pileupForOfflineID_1p.root and /dev/null differ diff --git a/Reconstruction/tauRecTools/share/TF2pileupForOfflineID_3p.root b/Reconstruction/tauRecTools/share/TF2pileupForOfflineID_3p.root deleted file mode 100644 index 3be443eab64f1fe1c8aa05d264f9ff8e625f3007..0000000000000000000000000000000000000000 Binary files a/Reconstruction/tauRecTools/share/TF2pileupForOfflineID_3p.root and /dev/null differ diff --git a/Reconstruction/tauRecTools/src/JetSeedBuilder.cxx b/Reconstruction/tauRecTools/src/JetSeedBuilder.cxx index f67c1dc36d25b4e042842723d0d822042eeb5f1b..6d23975ca841b813563637fd1d2133407d8b6930 100644 --- a/Reconstruction/tauRecTools/src/JetSeedBuilder.cxx +++ b/Reconstruction/tauRecTools/src/JetSeedBuilder.cxx @@ -179,8 +179,8 @@ StatusCode JetSeedBuilder::execute(xAOD::TauJet& pTau) { // SL/SX trigger mode with negative jet_seed - do not set TauJet eta and phi in JetSeedBuilder ATH_MSG_DEBUG("TauJet eta/phi will be set in Level2 Trigger for negative energy jet"); - pTau.setDetail(xAOD::TauJetParameters::seedCalo_eta , static_cast<float>( pTau.eta() ) ); - pTau.setDetail(xAOD::TauJetParameters::seedCalo_phi , static_cast<float>( pTau.phi() ) ); + // pTau.setDetail(xAOD::TauJetParameters::seedCalo_eta , static_cast<float>( pTau.eta() ) ); + // pTau.setDetail(xAOD::TauJetParameters::seedCalo_phi , static_cast<float>( pTau.phi() ) ); pTau.setP4(pJetSeed->pt(),pTau.eta(),pTau.phi(),0.0); @@ -192,12 +192,12 @@ StatusCode JetSeedBuilder::execute(xAOD::TauJet& pTau) { // sigstateH.controlObject(pJetSeed); } - pTau.setDetail(xAOD::TauJetParameters::seedCalo_eta , static_cast<float>( pJetSeed->eta() ) ); - pTau.setDetail(xAOD::TauJetParameters::seedCalo_phi , static_cast<float>( pJetSeed->phi() ) ); + // pTau.setDetail(xAOD::TauJetParameters::seedCalo_eta , static_cast<float>( pJetSeed->eta() ) ); + // pTau.setDetail(xAOD::TauJetParameters::seedCalo_phi , static_cast<float>( pJetSeed->phi() ) ); if ( pJetSeed->pt() > 1e-7) pTau.setP4(static_cast<float>( pJetSeed->pt() ) ,static_cast<float>( pJetSeed->eta() ) ,static_cast<float>( pJetSeed->phi() ) ,0.0 ); else - pTau.setP4(static_cast<float>( pJetSeed->pt() ) ,static_cast<float>( pJetSeed->eta() ) ,static_cast<float>( pJetSeed->phi() ) , 0.0 ); + pTau.setP4(static_cast<float>( 1e-7 ) ,static_cast<float>( pJetSeed->eta() ) ,static_cast<float>( pJetSeed->phi() ) , 0.0 ); //store 4-vector of seed pTau.setP4( xAOD::TauJetParameters::JetSeed, pJetSeed->pt(), pJetSeed->eta(), pJetSeed->phi(), pJetSeed->m() ); diff --git a/Reconstruction/tauRecTools/src/PhotonConversionVertex.cxx b/Reconstruction/tauRecTools/src/PhotonConversionVertex.cxx index 5acee948f63eaccf861f696d015b749776e8195e..03b5efb82625ede4756768f9f0b44239338938db 100644 --- a/Reconstruction/tauRecTools/src/PhotonConversionVertex.cxx +++ b/Reconstruction/tauRecTools/src/PhotonConversionVertex.cxx @@ -124,10 +124,14 @@ bool PhotonConversionVertex::openContainer(T* &container, std::string containerN template <class T> bool PhotonConversionVertex::saveContainer(T* &container, std::string containerName) { - StatusCode sc = evtStore()->record(container, containerName); - if (!container || sc.isFailure()) - ATH_MSG_FATAL("Container (" << containerName << ") cannot be saved in StoreGate"); - return container; + if(!container){ + ATH_MSG_FATAL("Container (" << containerName << ") cannot be saved in StoreGate"); + return false; + } + StatusCode sc = evtStore()->record(container, containerName); + if (sc.isFailure()) + ATH_MSG_FATAL("Container (" << containerName << ") cannot be saved in StoreGate"); + return container; } template <class T> diff --git a/Reconstruction/tauRecTools/src/TauAxisSetter.cxx b/Reconstruction/tauRecTools/src/TauAxisSetter.cxx index 81b7015147c5b6b42ac744886075b8fd16312dc1..360327ca3298a3466b98468552a32cc38c952cb2 100644 --- a/Reconstruction/tauRecTools/src/TauAxisSetter.cxx +++ b/Reconstruction/tauRecTools/src/TauAxisSetter.cxx @@ -88,6 +88,10 @@ StatusCode TauAxisSetter::execute(xAOD::TauJet& pTau) if (BaryCenter.DeltaR(tempClusterVector) > m_clusterCone) continue; + ElementLink<xAOD::IParticleContainer> linkToCluster; + linkToCluster.toContainedElement( *(static_cast<const xAOD::IParticleContainer*> ((*cItr)->rawConstituent()->container())), (*cItr)->rawConstituent() ); + pTau.addClusterLink(linkToCluster); + nConstituents++; tauDetectorAxis += tempClusterVector; } @@ -111,10 +115,6 @@ StatusCode TauAxisSetter::execute(xAOD::TauJet& pTau) ATH_MSG_VERBOSE("jet axis:" << (*pTau.jetLink())->pt()<< " " << (*pTau.jetLink())->eta() << " " << (*pTau.jetLink())->phi() << " " << (*pTau.jetLink())->e() ); // save values for detector axis. - // FixMe: consider dropping these details variables as they are duplicated in the detector axis 4 vector - pTau.setDetail(xAOD::TauJetParameters::LC_TES_precalib , static_cast<float>( tauDetectorAxis.Pt() ) ); - pTau.setDetail(xAOD::TauJetParameters::seedCalo_eta, static_cast<float>( tauDetectorAxis.Eta() ) ); - pTau.setDetail(xAOD::TauJetParameters::seedCalo_phi, static_cast<float>( tauDetectorAxis.Phi() ) ); ATH_MSG_VERBOSE("detector axis:" << tauDetectorAxis.Pt()<< " " << tauDetectorAxis.Eta() << " " << tauDetectorAxis.Phi() << " " << tauDetectorAxis.E() ); // detectorAxis (set default) diff --git a/Reconstruction/tauRecTools/src/TauCalibrateEM.cxx b/Reconstruction/tauRecTools/src/TauCalibrateEM.cxx index 5061009ebde172b4cf63abf870bff5b8330ab2b0..e0ee73b556c0b3c204a12095e566b8aa662e0bac 100644 --- a/Reconstruction/tauRecTools/src/TauCalibrateEM.cxx +++ b/Reconstruction/tauRecTools/src/TauCalibrateEM.cxx @@ -110,7 +110,7 @@ StatusCode TauCalibrateEM::execute(xAOD::TauJet& pTau) { ATH_MSG_DEBUG("input variables: em_pt " << emscale_pt << " eta " << emscale_eta << " ntrack " << pTau.nTracks() << " emfrac " << emfrac); - double new_pt = evaluate_new_pt(emscale_pt / GeV, fabs(emscale_eta), pTau.nTracks(), emfrac); + /*double new_pt = */evaluate_new_pt(emscale_pt / GeV, fabs(emscale_eta), pTau.nTracks(), emfrac); // do NOT set TauJet energy, as this will be done in tauCalibrateLC @@ -119,7 +119,9 @@ StatusCode TauCalibrateEM::execute(xAOD::TauJet& pTau) { // instead fill place holder in TauCommonDetails //pDetails->setSeedCalo_etEMCalib(new_pt * GeV); - pTau.setDetail( xAOD::TauJetParameters::EM_TES_scale, static_cast<float>( new_pt * GeV ) ); + //r21 cleanup + ATH_MSG_WARNING("EM_TES_scale removed"); + //pTau.setDetail( xAOD::TauJetParameters::EM_TES_scale, static_cast<float>( new_pt * GeV ) ); return StatusCode::SUCCESS; } diff --git a/Reconstruction/tauRecTools/src/TauConversionFinder.cxx b/Reconstruction/tauRecTools/src/TauConversionFinder.cxx index a435d146412fe98a9307849d2a48535b48b6385f..449b45b79391a7765173e09c5cc65e54c2bd8c11 100644 --- a/Reconstruction/tauRecTools/src/TauConversionFinder.cxx +++ b/Reconstruction/tauRecTools/src/TauConversionFinder.cxx @@ -119,21 +119,24 @@ StatusCode TauConversionFinder::eventFinalize() { // Find conversion in normal tau tracks if (m_do_normal) { for (unsigned int j = 0; j < numTracks; ++j) { - const xAOD::TrackParticle *pTauTrack = pTau.track(j); - const Trk::Track* tau_trk_def = pTauTrack->track(); + //const xAOD::TrackParticle *pTauTrack = pTau.track(j); + xAOD::TauTrack *pTauTrack = pTau.trackNonConst(j); + const Trk::Track* tau_trk_def = pTauTrack->track()->track(); if (conv_trk == tau_trk_def) { if (conv_trk->trackSummary()->getPID(Trk::eProbabilityHT) > m_eProb_cut) { - if (!pTau.trackFlag(pTauTrack, xAOD::TauJetParameters::isConversion)) { + //if (!pTau.trackFlag(pTauTrack, xAOD::TauJetParameters::isConversion)) { + if (!pTauTrack->flag(xAOD::TauJetParameters::isConversionOld)) { ElementLink<xAOD::TrackParticleContainer> phoConvLink ; - phoConvLink.setElement(pTauTrack) ; + //phoConvLink.setElement(pTauTrack) ; + phoConvLink.setElement(pTauTrack->track()) ; phoConvLink.setStorableObject( *trackContainer ) ; - phoConvLink.index(); - pTau.addTrackLink( phoConvLink ) ; - pTau.setTrackFlag(pTauTrack, xAOD::TauJetParameters::isConversion, true); + //phoConvLink.index(); + pTauTrack->addTrackLink( phoConvLink ); + pTauTrack->setFlag(xAOD::TauJetParameters::isConversionOld, true); if (m_adjust_tau_charge) - pTau.setCharge(pTau.charge() - pTau.track(j)->charge()); + pTau.setCharge(pTau.charge() - pTau.track(j)->track()->charge()); m_numProng--; } diff --git a/Reconstruction/tauRecTools/src/TauConversionTagger.cxx b/Reconstruction/tauRecTools/src/TauConversionTagger.cxx index b44a2ecb89490a00cf0b118984a1d9fc33c76e08..1f1907c5e6cf7f9282aec600ecfb9f7dcd871236 100644 --- a/Reconstruction/tauRecTools/src/TauConversionTagger.cxx +++ b/Reconstruction/tauRecTools/src/TauConversionTagger.cxx @@ -90,7 +90,7 @@ StatusCode TauConversionTagger::execute(xAOD::TauJet& pTau) { for(unsigned int j=0; j<pTau.nTracks(); j++ ) { - const xAOD::TrackParticle *TauJetTrack = pTau.track(j); + const xAOD::TrackParticle *TauJetTrack = pTau.track(j)->track(); const Trk::Perigee* perigee = m_trackToVertexTool->perigeeAtVertex(*TauJetTrack, (*pTau.vertexLink())->position()); // Declare TrackSummary info @@ -169,8 +169,11 @@ StatusCode TauConversionTagger::execute(xAOD::TauJet& pTau) { } ATH_MSG_VERBOSE("Is tau track a conversion? : " << m_TrkIsConv); - if (m_TrkIsConv && !pTau.trackFlag(TauJetTrack, xAOD::TauJetParameters::isConversion)) - pTau.setTrackFlag(TauJetTrack, xAOD::TauJetParameters::isConversion, true); + // if (m_TrkIsConv && !pTau.trackFlag(TauJetTrack, xAOD::TauJetParameters::isConversion)) + // pTau.setTrackFlag(TauJetTrack, xAOD::TauJetParameters::isConversion, true); + xAOD::TauTrack* tauTrack = pTau.trackNonConst(j); + if(m_TrkIsConv && !tauTrack->flag(xAOD::TauJetParameters::isConversionOld)) + tauTrack->setFlag( xAOD::TauJetParameters::isConversionOld, true); } return StatusCode::SUCCESS; diff --git a/Reconstruction/tauRecTools/src/TauElectronVetoVariables.cxx b/Reconstruction/tauRecTools/src/TauElectronVetoVariables.cxx index 60edf799ae435a175d711a3ae105a253951eb559..e5d949ded6a6813a41f9d65c3cfde24de26e3df0 100644 --- a/Reconstruction/tauRecTools/src/TauElectronVetoVariables.cxx +++ b/Reconstruction/tauRecTools/src/TauElectronVetoVariables.cxx @@ -166,7 +166,7 @@ StatusCode TauElectronVetoVariables::execute(xAOD::TauJet& pTau) // get the extrapolation into the calo const Trk::CaloExtension* caloExtension = 0; - if( !m_caloExtensionTool->caloExtension(*pTau.track(0),caloExtension) || caloExtension->caloLayerIntersections().empty() ){ + if( !m_caloExtensionTool->caloExtension(*pTau.track(0)->track(),caloExtension) || caloExtension->caloLayerIntersections().empty() ){ ATH_MSG_WARNING("extrapolation of leading track to calo surfaces failed " ); return StatusCode::SUCCESS; } @@ -340,23 +340,25 @@ StatusCode TauElectronVetoVariables::execute(xAOD::TauJet& pTau) uint8_t TRTHTOutliers; uint8_t TRTHits; uint8_t TRTOutliers; + + const xAOD::TrackParticle* leadTrack = pTau.track(0)->track(); - if ( !pTau.track(0)->summaryValue( TRTHits, xAOD::SummaryType::numberOfTRTHits ) ) + if ( !leadTrack->summaryValue( TRTHits, xAOD::SummaryType::numberOfTRTHits ) ) { ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate"); return StatusCode::SUCCESS; } - if ( !pTau.track(0)->summaryValue( TRTHTHits, xAOD::SummaryType::numberOfTRTHighThresholdHits ) ) + if ( !leadTrack->summaryValue( TRTHTHits, xAOD::SummaryType::numberOfTRTHighThresholdHits ) ) { ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate"); return StatusCode::SUCCESS; } - if ( !pTau.track(0)->summaryValue( TRTOutliers, xAOD::SummaryType::numberOfTRTOutliers ) ) + if ( !leadTrack->summaryValue( TRTOutliers, xAOD::SummaryType::numberOfTRTOutliers ) ) { ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate"); return StatusCode::SUCCESS; } - if ( !pTau.track(0)->summaryValue( TRTHTOutliers, xAOD::SummaryType::numberOfTRTHighThresholdOutliers ) ) + if ( !leadTrack->summaryValue( TRTHTOutliers, xAOD::SummaryType::numberOfTRTHighThresholdOutliers ) ) { ATH_MSG_DEBUG("retrieval of track summary value failed. Not filling electron veto variables for this one prong candidate"); return StatusCode::SUCCESS; diff --git a/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.cxx b/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.cxx index 0a9a9de73bcf14ebd6e56d57b4c7dd0d54ec609e..02bc803bf64341d92a26e0efc3b1ea0f43edfb84 100644 --- a/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.cxx +++ b/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.cxx @@ -37,6 +37,7 @@ TauPi0ClusterScaler::TauPi0ClusterScaler( const string& name ) : declareProperty("ParticleCaloExtensionTool", m_caloExtensionTool); declareProperty("ChargedPFOContainerName", m_chargedPFOContainerName); declareProperty("runOnAOD", m_AODmode=false); + declareProperty("storeCaloSamplings", m_storeCaloSamplings=true); } //------------------------------------------------------------------------- @@ -53,10 +54,6 @@ StatusCode TauPi0ClusterScaler::initialize() // retrieve tools ATH_MSG_DEBUG( "Retrieving tools" ); CHECK( m_caloExtensionTool.retrieve() ); - // Create vector with default values - for (int layer = 0 ; layer != CaloCell_ID::FCAL0; ++layer) { - m_defaultValues.push_back(-10.); - } return StatusCode::SUCCESS; } @@ -77,7 +74,19 @@ StatusCode TauPi0ClusterScaler::eventInitialize() { CHECK( evtStore()->retrieve(m_chargedPFOContainer, m_chargedPFOContainerName) ); CHECK( evtStore()->retrieve( m_chargedPFOAuxStore, m_chargedPFOContainerName + "Aux." ) ); } - return StatusCode::SUCCESS; + + //Check if TauTracks have sampling decorations + const xAOD::TauTrackContainer* tauTracks = 0; + ATH_CHECK( evtStore()->retrieve(tauTracks, "TauTracks") ); + for( const xAOD::TauTrack* trk : *tauTracks ){ + if( trk->isAvailable<float>("CaloSamplingEtaEM") ) { + m_caloSamplingsStored = true; + break; + } + m_caloSamplingsStored = false; + } + + return StatusCode::SUCCESS; } @@ -93,23 +102,11 @@ StatusCode TauPi0ClusterScaler::execute(xAOD::TauJet& pTau) // Clear vector of cell-based charged PFO Links. Required when rerunning on xAOD level. pTau.clearProtoChargedPFOLinks(); - //--------------------------------------------------------------------- - // only run on 1-5 prong taus - //--------------------------------------------------------------------- - if (pTau.nTracks() == 0 || pTau.nTracks() >5 ) { - return StatusCode::SUCCESS; - } - ATH_MSG_DEBUG("ClusterScaler: new tau. \tpt = " << pTau.pt() << "\teta = " << pTau.eta() << "\tphi = " << pTau.phi() << "\tnprongs = " << pTau.nTracks()); - + //incase tau is rejected, still fill vector of samplings + //so that TauTracks are consistently decorated //--------------------------------------------------------------------- // get tau tracks //--------------------------------------------------------------------- - vector<const xAOD::TrackParticle*> tracks; - for(unsigned iTrack = 0; iTrack<pTau.nTracks();++iTrack){ - const xAOD::TrackParticle* track = pTau.track(iTrack); - tracks.push_back(track); - } - //--------------------------------------------------------------------- // prepare extrapolation of tracks to calo layers //--------------------------------------------------------------------- @@ -118,16 +115,54 @@ StatusCode TauPi0ClusterScaler::execute(xAOD::TauJet& pTau) // which is called once per event (and not once per tau) m_tracksEtaAtSampling.clear(); m_tracksPhiAtSampling.clear(); - m_extrapolatedSamplings.clear(); - // Fill with default values - for(int layer = 0 ; layer != CaloCell_ID::FCAL0; ++layer) { - m_extrapolatedSamplings.push_back(false); - } - for(unsigned iTrack = 0; iTrack<tracks.size();++iTrack){ - m_tracksEtaAtSampling.push_back( m_defaultValues ); - m_tracksPhiAtSampling.push_back( m_defaultValues ); + vector<const xAOD::TauTrack*> tracks; + for(xAOD::TauTrack* track : pTau.allTracks()){ + + float extrap_eta_EM, extrap_phi_EM, extrap_eta_Had, extrap_phi_Had; + + if(m_caloSamplingsStored==false) { + //no decorations, so do calo extrapolations + int sampling_EM, sampling_Had; + if(fabs(track->eta())<1.45) sampling_EM = CaloSampling::EMB2;//2 + else sampling_EM = CaloSampling::EME2;//6 + if(fabs(track->eta())<1.0) sampling_Had = CaloSampling::TileBar1; + else if(fabs(track->eta())<1.5) sampling_Had = CaloSampling::TileExt1; + else sampling_Had = CaloSampling::HEC1; + + getExtrapolatedPositions(track, sampling_EM, extrap_eta_EM, extrap_phi_EM); + getExtrapolatedPositions(track, sampling_Had, extrap_eta_Had, extrap_phi_Had); + if(m_storeCaloSamplings) { + //store the values in tracks + track->setDetail(xAOD::TauJetParameters::CaloSamplingEtaEM, extrap_eta_EM); + track->setDetail(xAOD::TauJetParameters::CaloSamplingPhiEM, extrap_phi_EM); + track->setDetail(xAOD::TauJetParameters::CaloSamplingEtaHad, extrap_eta_Had); + track->setDetail(xAOD::TauJetParameters::CaloSamplingPhiHad, extrap_phi_Had); + } + } + else { + //no need to perform calo extrapolation + track->detail(xAOD::TauJetParameters::CaloSamplingEtaEM, extrap_eta_EM ); + track->detail(xAOD::TauJetParameters::CaloSamplingPhiEM, extrap_phi_EM ); + track->detail(xAOD::TauJetParameters::CaloSamplingEtaHad, extrap_eta_Had); + track->detail(xAOD::TauJetParameters::CaloSamplingPhiHad, extrap_phi_Had); + } + if(track->flag(xAOD::TauJetParameters::classifiedCharged)) { + //now fill the extrapolated values in the v<v<float> > + tracks.push_back(track); + m_tracksEtaAtSampling.push_back({extrap_eta_EM, extrap_eta_Had}); + m_tracksPhiAtSampling.push_back({extrap_phi_EM, extrap_phi_Had}); + } } + + //--------------------------------------------------------------------- + // only run on 1-5 prong taus + //--------------------------------------------------------------------- + if (pTau.nTracks() == 0 || pTau.nTracks() >5 ) { + return StatusCode::SUCCESS; + } + ATH_MSG_DEBUG("ClusterScaler: new tau. \tpt = " << pTau.pt() << "\teta = " << pTau.eta() << "\tphi = " << pTau.phi() << "\tnprongs = " << pTau.nTracks()); + //--------------------------------------------------------------------- // get energy in HCal associated to the different tracks //--------------------------------------------------------------------- @@ -138,10 +173,10 @@ StatusCode TauPi0ClusterScaler::execute(xAOD::TauJet& pTau) // Create charged PFOs //--------------------------------------------------------------------- for(unsigned iTrack = 0; iTrack<tracks.size();++iTrack){ - const xAOD::TrackParticle* track = tracks.at(iTrack); + const xAOD::TrackParticle* track = tracks.at(iTrack)->track(); xAOD::PFO* chargedPFO = new xAOD::PFO(); m_chargedPFOContainer->push_back(chargedPFO); - ElementLink<xAOD::TrackParticleContainer> myTrackLink = pTau.trackLinks().at(iTrack); + ElementLink<xAOD::TrackParticleContainer> myTrackLink = pTau.trackNonConst(iTrack)->trackLinks()[0]; if(!chargedPFO->setTrackLink(myTrackLink)) ATH_MSG_WARNING("Could not add Track to PFO"); chargedPFO->setCharge(track->charge()); chargedPFO->setP4(track->pt(),track->eta(),track->phi(),track->m()); @@ -169,19 +204,12 @@ StatusCode TauPi0ClusterScaler::execute(xAOD::TauJet& pTau) unsigned nNeutPFO = pTau.nProtoNeutralPFOs(); for(unsigned int iNeutPFO=0; iNeutPFO<nNeutPFO; iNeutPFO++, thisCluster++) { const xAOD::PFO* curNeutPFO_const = pTau.protoNeutralPFO( iNeutPFO ); - int maxESample = 2; - if (fabs(curNeutPFO_const->eta()) > 1.45) maxESample = 6; - // check if tracks have been extrapolated to this sampling. Do so if this is not the case - if(m_extrapolatedSamplings.at(maxESample)==false){ - this->getExtrapolatedPositions(tracks,maxESample); - m_extrapolatedSamplings.at(maxESample)=true; - } for(unsigned iTrack = 0; iTrack<tracks.size();++iTrack){ if(EestInEcal.at(iTrack)<0.001) continue; // No need to subtract TLorentzVector extTrack; - extTrack.SetPtEtaPhiE(tracks.at(iTrack)->pt(), m_tracksEtaAtSampling.at(iTrack).at(maxESample), m_tracksPhiAtSampling.at(iTrack).at(maxESample), tracks.at(iTrack)->e()); + extTrack.SetPtEtaPhiE(tracks.at(iTrack)->pt(), m_tracksEtaAtSampling.at(iTrack).at(0), m_tracksPhiAtSampling.at(iTrack).at(0), tracks.at(iTrack)->e()); // get eta/phi distance of cell to track double deltaEta = extTrack.Eta()-curNeutPFO_const->eta(); double deltaPhi = TVector2::Phi_mpi_pi( extTrack.Phi() - curNeutPFO_const->phi());; @@ -252,35 +280,33 @@ StatusCode TauPi0ClusterScaler::execute(xAOD::TauJet& pTau) } void TauPi0ClusterScaler::getExtrapolatedPositions( - vector<const xAOD::TrackParticle*> tracks, - int sampling) + const xAOD::TauTrack* track, + int sampling, float& extrap_eta, float& extrap_phi) { - for (unsigned iTrack = 0 ; iTrack < tracks.size(); ++iTrack ) { - // get the extrapolation into the calo - ATH_MSG_DEBUG( "Try extrapolation of track with pt = " << tracks.at(iTrack)->pt() << ", eta " << tracks.at(iTrack)->eta() << ", phi" << tracks.at(iTrack)->phi() - << " to layer " << sampling); - const Trk::CaloExtension* caloExtension = 0; - if (!m_caloExtensionTool->caloExtension(*tracks.at(iTrack),caloExtension) - || caloExtension->caloLayerIntersections().size() < (unsigned int)(sampling+1)) return; - - // store if track extrapolation successful, only use entry layer - const Trk::TrackParameters* param_at_calo = caloExtension->caloLayerIntersections().at(sampling); - if (param_at_calo) { - ATH_MSG_DEBUG( "Extrapolated track with eta=" << tracks.at(iTrack)->eta() - << " phi="<<tracks.at(iTrack)->phi() - << " to eta=" << param_at_calo->position().eta() - << " phi="<<param_at_calo->position().phi() - ); - m_tracksEtaAtSampling.at(iTrack).at(sampling)=param_at_calo->position().eta(); - m_tracksPhiAtSampling.at(iTrack).at(sampling)=param_at_calo->position().phi(); - } - else ATH_MSG_DEBUG("Could not extrapolate track with pt = " << tracks.at(iTrack)->pt() << ", eta " << tracks.at(iTrack)->eta() << ", phi" << tracks.at(iTrack)->phi() - << " to layer " << sampling); - } + extrap_eta=-10; + extrap_phi=-10; + ATH_MSG_DEBUG( "Try extrapolation of track with pt = " << track->pt() << ", eta " << track->eta() << ", phi" << track->phi() + << " to layer " << sampling); + const Trk::CaloExtension* caloExtension = 0; + if (!m_caloExtensionTool->caloExtension(*track->track(),caloExtension) + || caloExtension->caloLayerIntersections().size() < (unsigned int)(sampling+1)) return; + // store if track extrapolation successful, only use entry layer + const Trk::TrackParameters* param_at_calo = caloExtension->caloLayerIntersections().at(sampling); + if (param_at_calo) { + ATH_MSG_DEBUG( "Extrapolated track with eta=" << track->eta() + << " phi="<<track->phi() + << " to eta=" << param_at_calo->position().eta() + << " phi="<<param_at_calo->position().phi() + ); + extrap_eta=param_at_calo->position().eta(); + extrap_phi=param_at_calo->position().phi(); + } + else ATH_MSG_DEBUG("Could not extrapolate track with pt = " << track->pt() << ", eta " << track->eta() << ", phi" << track->phi() + << " to layer " << sampling); } vector<double> TauPi0ClusterScaler::getEstEcalEnergy( - vector<const xAOD::TrackParticle*> tracks, + vector<const xAOD::TauTrack*> tracks, const xAOD::TauJet& pTau, vector<vector<ElementLink<xAOD::IParticleContainer> > >& hadPFOLinks) { @@ -310,28 +336,16 @@ vector<double> TauPi0ClusterScaler::getEstEcalEnergy( "\t deltaEtaToTau = " << deltaEtaToTau << "\t deltaPhiToTau = " << deltaPhiToTau << "\t deltaRToTau_squared = " << deltaRToTau_squared ); */ - // Decide which sampling to extrapolate to. Choose Hcal samplings that usually contain most energy (|eta| dependent) - int sample = 13; // |eta| <= 1.0 - if ( fabs(curHadPFO->eta())>1.5 ) sample = 9; // 1.5 < |eta| - else if( fabs(curHadPFO->eta())>1.0 ) sample = 19; // 1.0 < |eta| <= 1.5 - - // check if tracks have been extrapolated to this sampling. Do so if this is not the case - if(m_extrapolatedSamplings.at(sample)==false){ - this->getExtrapolatedPositions(tracks,sample); - ATH_MSG_DEBUG("Extrapolate to layer " << sample << "\teta = " - << m_tracksEtaAtSampling.at(0).at(sample) << "\t phi = " << m_tracksPhiAtSampling.at(0).at(sample) ); - m_extrapolatedSamplings.at(sample)=true; - } // Assign PFO to track int closestTrack = -1; double dRToClosestTrack_squared = 0.16; // XXX can be tuned later for(unsigned iTrack = 0; iTrack<tracks.size();++iTrack){ - const xAOD::TrackParticle* track = tracks.at(iTrack); + const xAOD::TrackParticle* track = tracks.at(iTrack)->track(); // set extrapolated track direction TLorentzVector extTrack; - extTrack.SetPtEtaPhiE(track->pt(), m_tracksEtaAtSampling.at(iTrack).at(sample), m_tracksPhiAtSampling.at(iTrack).at(sample), track->e()); + extTrack.SetPtEtaPhiE(track->pt(), m_tracksEtaAtSampling.at(iTrack).at(1), m_tracksPhiAtSampling.at(iTrack).at(1), track->e()); // get eta/phi distance of cell to track double deltaEta = extTrack.Eta()-curHadPFO->eta(); diff --git a/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.h b/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.h index 2cc3ed5960158c2e0bb9ccaa761258c2896b466a..ac68d0225a3fedf3bc711052a96e01571f8ba26f 100644 --- a/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.h +++ b/Reconstruction/tauRecTools/src/TauPi0ClusterScaler.h @@ -42,14 +42,13 @@ private: ToolHandle<Trk::IParticleCaloExtensionTool> m_caloExtensionTool; /** @brief extrapolated position of tracks and vector of bools to keep track for which samplings this has already been done */ - std::vector<std::vector<double> > m_tracksEtaAtSampling; - std::vector<std::vector<double> > m_tracksPhiAtSampling; - std::vector<bool> m_extrapolatedSamplings; - std::vector<double> m_defaultValues; + //inner vector is size 2: EM, Had layers + std::vector<std::vector<float> > m_tracksEtaAtSampling; + std::vector<std::vector<float> > m_tracksPhiAtSampling; /** @brief get extrapolated track position at each layer */ - void getExtrapolatedPositions( std::vector<const xAOD::TrackParticle*>, - int sampling); + void getExtrapolatedPositions( const xAOD::TauTrack*, + int sampling, float& eta, float& phi); /** @brief new charged PFO container and name */ xAOD::PFOContainer* m_chargedPFOContainer; @@ -58,10 +57,12 @@ private: /** @brief run on AOD */ bool m_AODmode; + bool m_storeCaloSamplings;//configurable to store whether you want calo samplings on TauTracks + bool m_caloSamplingsStored=false;//check whether TauJets have sampling decorations or not /** @brief estimate energy deposited in Ecal by each charged pion */ std::vector<double> getEstEcalEnergy( - std::vector<const xAOD::TrackParticle*> tracks, + std::vector<const xAOD::TauTrack*> tracks, const xAOD::TauJet& pTau, std::vector<std::vector<ElementLink<xAOD::IParticleContainer> > >& hadPFOLinks); diff --git a/Reconstruction/tauRecTools/src/TauPi0CreateROI.cxx b/Reconstruction/tauRecTools/src/TauPi0CreateROI.cxx index 763b2d9c802e901c85f8fb2e305c7a3dd6fe61ae..0e73063fdbc593bbe9a2337503d1b2750eff02e9 100644 --- a/Reconstruction/tauRecTools/src/TauPi0CreateROI.cxx +++ b/Reconstruction/tauRecTools/src/TauPi0CreateROI.cxx @@ -166,7 +166,7 @@ void TauPi0CreateROI::storeCell(const CaloCell* cell){ if(isNewCell){ CaloCell* copyCell = cell->clone(); - m_pPi0CellContainer->push_back(const_cast<CaloCell*> (copyCell)); + m_pPi0CellContainer->push_back(copyCell); m_addedCellsMap[cellHash] = copyCell; } } diff --git a/Reconstruction/tauRecTools/src/TauPi0Selector.cxx b/Reconstruction/tauRecTools/src/TauPi0Selector.cxx index 67fa5d2e261364a27b56e7a84e8b043430f7cbe7..a0c9609bbc6d10b70479c79f6b277cac20fd7fe5 100644 --- a/Reconstruction/tauRecTools/src/TauPi0Selector.cxx +++ b/Reconstruction/tauRecTools/src/TauPi0Selector.cxx @@ -51,7 +51,7 @@ StatusCode TauPi0Selector::execute(xAOD::TauJet& pTau) // Clear vector of cell-based pi0 PFO Links. Required when rerunning on xAOD level. pTau.clearProtoPi0PFOLinks(); // Set proto decay mode to "not set". Will be overwritten for taus with 1-5 tracks - pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_NotSet); + pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_NotSet); //--------------------------------------------------------------------- // only run on 1-5 prong taus @@ -61,7 +61,7 @@ StatusCode TauPi0Selector::execute(xAOD::TauJet& pTau) } // Set proto decay mode to "other". Will be overwritten for taus with 1 or 3 tracks - pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_Other); + pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_Other); //--------------------------------------------------------------------- // retrieve neutral PFOs from tau. Apply selection and create links to @@ -115,13 +115,13 @@ StatusCode TauPi0Selector::execute(xAOD::TauJet& pTau) // pTau.setMPanTauCellBasedProto( p4.M()); if(pTau.nTracks()==1){ - if(nRecoPi0s==0) pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_1p0n); - else if(nRecoPi0s==1) pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_1p1n); - else pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_1pXn); + if(nRecoPi0s==0) pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_1p0n); + else if(nRecoPi0s==1) pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_1p1n); + else pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_1pXn); } if(pTau.nTracks()==3){ - if(nRecoPi0s==0) pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_3p0n); - else pTau.setPanTauDetail(xAOD::TauJetParameters::pantau_CellBasedInput_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_3pXn); + if(nRecoPi0s==0) pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_3p0n); + else pTau.setPanTauDetail(xAOD::TauJetParameters::PanTau_DecayModeProto, xAOD::TauJetParameters::DecayMode::Mode_3pXn); } diff --git a/Reconstruction/tauRecTools/src/TauShotVariableHelpers.cxx b/Reconstruction/tauRecTools/src/TauShotVariableHelpers.cxx index 3645dd3b5db95f8884e4c39f18458878b2e9db5b..87eeb517b97725d3388848c77765b4cb6a028dad 100644 --- a/Reconstruction/tauRecTools/src/TauShotVariableHelpers.cxx +++ b/Reconstruction/tauRecTools/src/TauShotVariableHelpers.cxx @@ -323,10 +323,10 @@ namespace TauShotVariableHelpers { float pt_largerWindow = 0.; float pt_smallerWindow = 0.; for(int iCell = 0; iCell != nCells_eta; ++iCell ){ - if(fabs(iCell-seedIndex)>largerWindow/2) continue; + if(abs(iCell-seedIndex)>largerWindow/2) continue; if(shotCells.at(0).at(iCell)!=NULL) pt_largerWindow+=shotCells.at(0).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(0).at(iCell)); if(shotCells.at(1).at(iCell)!=NULL) pt_largerWindow+=shotCells.at(1).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(1).at(iCell)); - if(fabs(iCell-seedIndex)>smallerWindow/2) continue; + if(abs(iCell-seedIndex)>smallerWindow/2) continue; if(shotCells.at(0).at(iCell)!=NULL) pt_smallerWindow+=shotCells.at(0).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(0).at(iCell)); if(shotCells.at(1).at(iCell)!=NULL) pt_smallerWindow+=shotCells.at(1).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(1).at(iCell)); } @@ -344,10 +344,10 @@ namespace TauShotVariableHelpers { float pt_largerWindow = 0.; float pt_smallerWindow = 0.; for(int iCell = 0; iCell != nCells_eta; ++iCell ){ - if(fabs(iCell-seedIndex)>largerWindow/2) continue; + if(abs(iCell-seedIndex)>largerWindow/2) continue; if(shotCells.at(0).at(iCell)!=NULL) pt_largerWindow+=shotCells.at(0).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(0).at(iCell)); if(shotCells.at(1).at(iCell)!=NULL) pt_largerWindow+=shotCells.at(1).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(1).at(iCell)); - if(fabs(iCell-seedIndex)>smallerWindow/2) continue; + if(abs(iCell-seedIndex)>smallerWindow/2) continue; if(shotCells.at(0).at(iCell)!=NULL) pt_smallerWindow+=shotCells.at(0).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(0).at(iCell)); if(shotCells.at(1).at(iCell)!=NULL) pt_smallerWindow+=shotCells.at(1).at(iCell)->pt()*m_caloWeightTool->wtCell(shotCells.at(1).at(iCell)); } diff --git a/Reconstruction/tauRecTools/src/TauTrackFinder.cxx b/Reconstruction/tauRecTools/src/TauTrackFinder.cxx index 9ab9fa9621f3d231ff9238c7081e533d0ffe7750..d3bd289915d6953793621234cfa0a2c295e08375 100644 --- a/Reconstruction/tauRecTools/src/TauTrackFinder.cxx +++ b/Reconstruction/tauRecTools/src/TauTrackFinder.cxx @@ -8,6 +8,7 @@ #include "xAODTau/TauJet.h" #include "xAODTau/TauJetContainer.h" +#include "xAODTau/TauTrackContainer.h" #include "TauTrackFinder.h" #include "tauRecTools/KineUtils.h" @@ -19,6 +20,7 @@ TauTrackFinder::TauTrackFinder(const std::string& name ) : m_caloExtensionTool("Trk::ParticleCaloExtensionTool/ParticleCaloExtensionTool"), m_trackSelectorTool_tau(""), m_trackToVertexTool("Reco::TrackToVertex"), + m_trackSelectorTool_tau_xAOD(""), m_z0maxDelta(1000), m_applyZ0cut(false), m_storeInOtherTrks(true), @@ -28,12 +30,15 @@ TauTrackFinder::TauTrackFinder(const std::string& name ) : declareProperty("MaxJetDrTau", m_maxJetDr_tau = 0.2); declareProperty("MaxJetDrWide", m_maxJetDr_wide = 0.4); declareProperty("TrackSelectorToolTau", m_trackSelectorTool_tau); + declareProperty("TrackSelectorToolTauxAOD", m_trackSelectorTool_tau_xAOD); declareProperty("TrackParticleContainer", m_inputTrackParticleContainerName = "InDetTrackParticles"); + declareProperty("TrackParticleContainer", m_inputTauTrackContainerName = "TauTracks"); declareProperty("ParticleCaloExtensionTool", m_caloExtensionTool ); declareProperty("TrackToVertexTool",m_trackToVertexTool); declareProperty("maxDeltaZ0wrtLeadTrk", m_z0maxDelta); declareProperty("removeTracksOutsideZ0wrtLeadTrk", m_applyZ0cut); declareProperty("StoreRemovedCoreWideTracksInOtherTracks", m_storeInOtherTrks = true); + declareProperty("removeDuplicateCoreTracks", m_removeDuplicateCoreTracks = true); declareProperty("BypassSelector", m_bypassSelector = false); declareProperty("BypassExtrapolator", m_bypassExtrapolator = false); } @@ -46,6 +51,7 @@ StatusCode TauTrackFinder::initialize() { // Get the TrackSelectorTool if (!retrieveTool(m_trackSelectorTool_tau)) return StatusCode::FAILURE; + //if (!retrieveTool(m_trackSelectorTool_tau_xAOD)) return StatusCode::FAILURE; // Get the TJVA if (!retrieveTool(m_trackToVertexTool)) return StatusCode::FAILURE; @@ -77,10 +83,14 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { StatusCode sc; // get the track particle container from StoreGate const xAOD::TrackParticleContainer* trackParticleCont = 0; + xAOD::TauTrackContainer* tauTrackCon = 0; //for tau trigger bool inTrigger = tauEventData()->inTrigger(); - if (inTrigger) sc = tauEventData()->getObject( "TrackContainer", trackParticleCont ); + if (inTrigger) { + ATH_CHECK(tauEventData()->getObject( "TrackContainer", trackParticleCont )); + ATH_CHECK(tauEventData()->getObject( "TauTrackContainer", tauTrackCon )); + } if( !inTrigger || !trackParticleCont || sc.isFailure() ) { // try standard @@ -90,6 +100,14 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { } } + if( !inTrigger || !tauTrackCon || sc.isFailure() ) { + // try standard + if (!openContainer(tauTrackCon, m_inputTauTrackContainerName)) { + if (!inTrigger) return StatusCode::FAILURE; // in offline we don't reconstruct tau candidates without having a track container + else return StatusCode::SUCCESS; // we don't want stop trigger if there is no track container + } + } + std::vector<const xAOD::TrackParticle*> tauTracks; std::vector<const xAOD::TrackParticle*> wideTracks; std::vector<const xAOD::TrackParticle*> otherTracks; @@ -106,12 +124,6 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { this->removeOffsideTracksWrtLeadTrk(tauTracks, wideTracks, otherTracks, pVertex, m_z0maxDelta); } - //clear tracks first (needed for "rerun mode" if called on AODs again) - pTau.clearTrackLinks(); - pTau.clearWideTrackLinks(); - pTau.clearOtherTrackLinks(); - - bool alreadyUsed = false; //check for tracks used in multiple taus xAOD::TauJetContainer* pContainer = tauEventData()->xAODTauContainer; if(pContainer==0){ @@ -122,34 +134,25 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { return StatusCode::FAILURE; } - - for (std::vector<const xAOD::TrackParticle*>::iterator track_it = tauTracks.begin(); track_it != tauTracks.end() ;) - { - alreadyUsed = false; - - //loop over all up-to-now reconstructed tau candidates - xAOD::TauJetContainer::const_iterator tau_it = pContainer->begin(); - xAOD::TauJetContainer::const_iterator tau_end = pContainer->end(); - for( ; tau_it != tau_end; tau_it++ ) - { - if( (*tau_it) == &pTau ) continue; - //loop over core tracks - for (unsigned int j = 0; j < (*tau_it)->nTracks(); ++j) - { - if ((*track_it) == (*tau_it)->track(j)) - { - ATH_MSG_WARNING("Found a track that is identical with a track already associated to another tau. Will not add this track to more than one tau candidate"); - alreadyUsed = true; - } - } - } - //if this track has already been used by another tau, don't associate it to this new one - if (alreadyUsed) track_it = tauTracks.erase(track_it); - else ++track_it; + if(m_removeDuplicateCoreTracks){ + bool alreadyUsed = false; + for (std::vector<const xAOD::TrackParticle*>::iterator track_it = tauTracks.begin(); track_it != tauTracks.end() ;) + { + alreadyUsed = false; + + //loop over all up-to-now core tracks + for( const xAOD::TauTrack* tau_trk : (*tauTrackCon) ) { + if(! tau_trk->flagWithMask( (1<<xAOD::TauJetParameters::TauTrackFlag::coreTrack) | (1<<xAOD::TauJetParameters::TauTrackFlag::passTrkSelector))) continue; //originally it was coreTrack&passTrkSelector + if( (*track_it) == tau_trk->track()) alreadyUsed = true; + } + //if this track has already been used by another tau, don't associate it to this new one + if(alreadyUsed) ATH_MSG_WARNING( "Found Already Used track new, now removing: " << *track_it ); + if (alreadyUsed) track_it = tauTracks.erase(track_it); + else ++track_it; + } } - // associated track to tau candidate and calculate charge float charge = 0; for (unsigned int i = 0; i < tauTracks.size(); ++i) { @@ -160,9 +163,24 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { << " phi " << trackParticle->phi() ); charge += trackParticle->charge(); + + xAOD::TauTrack* track = new xAOD::TauTrack(); + tauTrackCon->push_back(track); ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle; linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle); - pTau.addTrackLink(linkToTrackParticle); + track->addTrackLink(linkToTrackParticle); + track->setP4(trackParticle->pt(), trackParticle->eta(), trackParticle->phi(), trackParticle->m()); + //track->setCharge(track->charge()); + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::coreTrack, true); + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelector, true); + // in case TrackClassifier is not run, still get sensible results + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::classifiedCharged, true); + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::unclassified, true); + //track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelectionTight, m_trackSelectorTool_tau_xAOD->accept(trackParticle)); + ElementLink<xAOD::TauTrackContainer> linkToTauTrack; + linkToTauTrack.toContainedElement(*tauTrackCon, track); + pTau.addTauTrackLink(linkToTauTrack); + ATH_MSG_VERBOSE(name() << " added core track nr: " << i << " eta " << pTau.track(i)->eta() << " phi " << pTau.track(i)->phi() @@ -182,13 +200,24 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { << " eta " << trackParticle->eta() << " phi " << trackParticle->phi() ); + + xAOD::TauTrack* track = new xAOD::TauTrack(); + tauTrackCon->push_back(track); ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle; linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle); - pTau.addWideTrackLink(linkToTrackParticle); - ATH_MSG_VERBOSE(name() << " added wide track nr: " << i - << " eta " << pTau.wideTrack(i)->eta() - << " phi " << pTau.wideTrack(i)->phi() - ); + track->addTrackLink(linkToTrackParticle); + track->setP4(trackParticle->pt(), trackParticle->eta(), trackParticle->phi(), trackParticle->m()); + //track->setCharge(track->charge()); + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::wideTrack, true); + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelector, true); + // in case TrackClassifier is not run, still get sensible results + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::classifiedIsolation, true); // for sake of trigger, reset in TauTrackClassifier + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::unclassified, true); + //track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelectionTight, m_trackSelectorTool_tau_xAOD->accept(trackParticle)); + ElementLink<xAOD::TauTrackContainer> linkToTauTrack; + linkToTauTrack.toContainedElement(*tauTrackCon, track); + pTau.addTauTrackLink(linkToTauTrack); + } /// was @@ -201,13 +230,24 @@ StatusCode TauTrackFinder::execute(xAOD::TauJet& pTau) { << " eta " << trackParticle->eta() << " phi " << trackParticle->phi() ); - ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle; - linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle); - pTau.addOtherTrackLink(linkToTrackParticle); - ATH_MSG_VERBOSE(name() << " added other track nr: " << i - << " eta " << pTau.otherTrack(i)->eta() - << " phi " << pTau.otherTrack(i)->phi() - ); + + // bool accepted=m_trackSelectorTool_tau_xAOD->accept(trackParticle); + // if(accepted){ + xAOD::TauTrack* track = new xAOD::TauTrack(); + tauTrackCon->push_back(track); + ElementLink<xAOD::TrackParticleContainer> linkToTrackParticle; + linkToTrackParticle.toContainedElement(*trackParticleCont, trackParticle); + track->addTrackLink(linkToTrackParticle); + track->setP4(trackParticle->pt(), trackParticle->eta(), trackParticle->phi(), trackParticle->m()); + float dR = track->p4().DeltaR(pTau.p4()); + if(dR<=0.2) track->setFlag(xAOD::TauJetParameters::TauTrackFlag::coreTrack, true); + else track->setFlag(xAOD::TauJetParameters::TauTrackFlag::wideTrack, true); + track->setFlag(xAOD::TauJetParameters::TauTrackFlag::unclassified, true); + // track->setFlag(xAOD::TauJetParameters::TauTrackFlag::passTrkSelectionTight, accepted); + ElementLink<xAOD::TauTrackContainer> linkToTauTrack; + linkToTauTrack.toContainedElement(*tauTrackCon, track); + pTau.addTauTrackLink(linkToTauTrack); + // } } @@ -287,10 +327,11 @@ StatusCode TauTrackFinder::extrapolateToCaloSurface(xAOD::TauJet& pTau) { Trk::TrackParametersIdHelper parsIdHelper; - for (unsigned int itr = 0; itr < 10 && itr < pTau.nTracks(); ++itr) { - - const xAOD::TrackParticle *orgTrack = pTau.track(itr); - + // for (unsigned int itr = 0; itr < 10 && itr < pTau.nAllTracks(); ++itr) { + + for( xAOD::TauTrack* tauTrack : pTau.allTracks() ) { + const xAOD::TrackParticle *orgTrack = tauTrack->track(); + if( !orgTrack ) continue; // get the extrapolation into the calo @@ -306,8 +347,8 @@ StatusCode TauTrackFinder::extrapolateToCaloSurface(xAOD::TauJet& pTau) { CaloSampling::CaloSample sample = parsIdHelper.caloSample((*cur)->cIdentifier()); if( sample == CaloSampling::EME1 || sample == CaloSampling::EMB1 ){ - pTau.setTrackEtaStrip( itr, (*cur)->position().eta() ); - pTau.setTrackPhiStrip( itr, (*cur)->position().phi() ); + tauTrack->setEtaStrip((*cur)->position().eta()); + tauTrack->setPhiStrip((*cur)->position().phi()); break; } } diff --git a/Reconstruction/tauRecTools/src/TauTrackFinder.h b/Reconstruction/tauRecTools/src/TauTrackFinder.h index e6ad1425b930618c48a9a95dcd3db9e372015347..2206aec17acaafa44bec0e38f0a99e1d0d195031 100644 --- a/Reconstruction/tauRecTools/src/TauTrackFinder.h +++ b/Reconstruction/tauRecTools/src/TauTrackFinder.h @@ -15,6 +15,9 @@ #include "VxVertex/RecVertex.h" +// xAOD Tracking Tool +#include "InDetTrackSelectionTool/IInDetTrackSelectionTool.h" + namespace Trk { class ITrackSelectorTool; class IParticleCaloExtensionTool; @@ -110,6 +113,7 @@ private: std::string m_inputTauJetContainerName; std::string m_inputTrackParticleContainerName; std::string m_inputPrimaryVertexContainerName; + std::string m_inputTauTrackContainerName; //------------------------------------------------------------- //! tools @@ -117,6 +121,8 @@ private: ToolHandle< Trk::IParticleCaloExtensionTool > m_caloExtensionTool; ToolHandle<Trk::ITrackSelectorTool> m_trackSelectorTool_tau; ToolHandle<Reco::ITrackToVertex> m_trackToVertexTool; + ToolHandle<InDet::IInDetTrackSelectionTool> m_trackSelectorTool_tau_xAOD; + //------------------------------------------------------------- //! Input parameters for algorithm @@ -130,6 +136,7 @@ private: float m_z0maxDelta; bool m_applyZ0cut; bool m_storeInOtherTrks; + bool m_removeDuplicateCoreTracks; std::vector<float> m_vDeltaZ0coreTrks; std::vector<float> m_vDeltaZ0wideTrks; diff --git a/Reconstruction/tauRecTools/src/TauVertexFinder.cxx b/Reconstruction/tauRecTools/src/TauVertexFinder.cxx index 09fbbee02024165b2927d27bdc40310aa98b047d..1367033dd7001e21dcae06f1d5e5d731a37e04f3 100644 --- a/Reconstruction/tauRecTools/src/TauVertexFinder.cxx +++ b/Reconstruction/tauRecTools/src/TauVertexFinder.cxx @@ -4,8 +4,8 @@ #include "TauVertexFinder.h" -//#include "VxVertex/RecVertex.h" -//#include "VxVertex/VxCandidate.h" +#include "VxVertex/RecVertex.h" +#include "VxVertex/VxCandidate.h" #include "xAODTracking/VertexContainer.h" #include "xAODTracking/Vertex.h" @@ -15,16 +15,18 @@ #include "xAODTau/TauJet.h" TauVertexFinder::TauVertexFinder(const std::string& name ) : -TauRecToolBase(name), -m_printMissingContainerINFO(true), -m_maxJVF(-100.), -m_assocTracksName(""), -m_trackVertexAssocName("") + TauRecToolBase(name), + m_printMissingContainerINFO(true), + m_maxJVF(-100.), + m_TrackSelectionToolForTJVA(""), + m_assocTracksName(""), + m_trackVertexAssocName("") { - declareProperty("UseTJVA", m_useTJVA=true); - declareProperty("PrimaryVertexContainer", m_inputPrimaryVertexContainerName = "PrimaryVertices"); - declareProperty("AssociatedTracks",m_assocTracksName); - declareProperty("TrackVertexAssociation",m_trackVertexAssocName); + declareProperty("UseTJVA", m_useTJVA=true); + declareProperty("PrimaryVertexContainer", m_inputPrimaryVertexContainerName = "PrimaryVertices"); + declareProperty("AssociatedTracks",m_assocTracksName); + declareProperty("TrackVertexAssociation",m_trackVertexAssocName); + declareProperty("InDetTrackSelectionToolForTJVA",m_TrackSelectionToolForTJVA); } TauVertexFinder::~TauVertexFinder() { @@ -32,23 +34,24 @@ TauVertexFinder::~TauVertexFinder() { // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * StatusCode TauVertexFinder::initialize() { - if (m_useTJVA) ATH_MSG_INFO("using TJVA to determine tau vertex"); - return StatusCode::SUCCESS; + if (m_useTJVA) ATH_MSG_INFO("using TJVA to determine tau vertex"); + ATH_CHECK(retrieveTool(m_TrackSelectionToolForTJVA)); + return StatusCode::SUCCESS; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * StatusCode TauVertexFinder::finalize() { - return StatusCode::SUCCESS; + return StatusCode::SUCCESS; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * StatusCode TauVertexFinder::eventInitialize() { - return StatusCode::SUCCESS; + return StatusCode::SUCCESS; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * StatusCode TauVertexFinder::eventFinalize() { - return StatusCode::SUCCESS; + return StatusCode::SUCCESS; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * @@ -64,107 +67,117 @@ StatusCode TauVertexFinder::execute(xAOD::TauJet& pTau) { sc = tauEventData()->getObject("VxPrimaryCandidate", vxContainer); if (sc.isFailure() || !vxContainer) { //not in trigger mode or no vxContainer was set by trigger - ATH_MSG_DEBUG("no VxPrimaryCandidateContainer for trigger -> try standard way"); - if (!openContainer(vxContainer, m_inputPrimaryVertexContainerName)) { - if (m_printMissingContainerINFO) { - ATH_MSG_INFO(m_inputPrimaryVertexContainerName << " container not found --> skip TauVertexFinder (no further info)"); - m_printMissingContainerINFO=false; - } - return StatusCode::SUCCESS; - } - - // find default PrimaryVertex (needed if TJVA is switched off or fails) - // see: https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/VertexReselectionOnAOD - // code adapted from - // https://svnweb.cern.ch/trac/atlasoff/browser/Tracking/TrkEvent/VxVertex/trunk/VxVertex/PrimaryVertexSelector.h - if (vxContainer->size()>0) { - // simple loop through and get the primary vertex - xAOD::VertexContainer::const_iterator vxIter = vxContainer->begin(); - xAOD::VertexContainer::const_iterator vxIterEnd = vxContainer->end(); - for ( size_t ivtx = 0; vxIter != vxIterEnd; ++vxIter, ++ivtx ){ - // the first and only primary vertex candidate is picked - if ( (*vxIter)->vertexType() == xAOD::VxType::PriVtx){ - primaryVertex = (*vxIter); - break; - } - } - } + ATH_MSG_DEBUG("no VxPrimaryCandidateContainer for trigger -> try standard way"); + if (!openContainer(vxContainer, m_inputPrimaryVertexContainerName)) { + if (m_printMissingContainerINFO) { + ATH_MSG_INFO(m_inputPrimaryVertexContainerName << " container not found --> skip TauVertexFinder (no further info)"); + m_printMissingContainerINFO=false; + } + return StatusCode::SUCCESS; } - else { // trigger mode - // find default PrimaryVertex (highest sum pt^2) - if (vxContainer->size()>0) primaryVertex = (*vxContainer)[0]; + + // find default PrimaryVertex (needed if TJVA is switched off or fails) + // see: https://twiki.cern.ch/twiki/bin/viewauth/AtlasProtected/VertexReselectionOnAOD + // code adapted from + // https://svnweb.cern.ch/trac/atlasoff/browser/Tracking/TrkEvent/VxVertex/trunk/VxVertex/PrimaryVertexSelector.h + if (vxContainer->size()>0) { + // simple loop through and get the primary vertex + xAOD::VertexContainer::const_iterator vxIter = vxContainer->begin(); + xAOD::VertexContainer::const_iterator vxIterEnd = vxContainer->end(); + for ( size_t ivtx = 0; vxIter != vxIterEnd; ++vxIter, ++ivtx ){ + // the first and only primary vertex candidate is picked + if ( (*vxIter)->vertexType() == xAOD::VxType::PriVtx){ + primaryVertex = (*vxIter); + break; + } + } } + } + else { // trigger mode + // find default PrimaryVertex (highest sum pt^2) + if (vxContainer->size()>0) primaryVertex = (*vxContainer)[0]; + } - ATH_MSG_VERBOSE("size of VxPrimaryContainer is: " << vxContainer->size() ); + ATH_MSG_VERBOSE("size of VxPrimaryContainer is: " << vxContainer->size() ); - // associate vertex to tau - if (primaryVertex) pTau.setVertex(vxContainer, primaryVertex); + // associate vertex to tau + if (primaryVertex) pTau.setVertex(vxContainer, primaryVertex); - //stop here if TJVA is disabled or vertex container is empty - if (!m_useTJVA || vxContainer->size()==0) return StatusCode::SUCCESS; - - // try to find new PV with TJVA - ATH_MSG_DEBUG("TJVA enabled -> try to find new PV for the tau candidate"); - - ElementLink<xAOD::VertexContainer> newPrimaryVertexLink = getPV_TJVA(pTau, *vxContainer ); - if (newPrimaryVertexLink.isValid()) { - // set new primary vertex - // will overwrite default one which was set above - pTau.setVertexLink(newPrimaryVertexLink); - // save highest JVF value - pTau.setDetail(xAOD::TauJetParameters::TauJetVtxFraction,static_cast<float>(m_maxJVF)); - ATH_MSG_DEBUG("TJVA vertex found and set"); - } - else { - ATH_MSG_DEBUG("couldn't find new PV for TJVA"); - } - - return StatusCode::SUCCESS; + //stop here if TJVA is disabled or vertex container is empty + if (!m_useTJVA || vxContainer->size()==0) return StatusCode::SUCCESS; + + // try to find new PV with TJVA + ATH_MSG_DEBUG("TJVA enabled -> try to find new PV for the tau candidate"); + + ElementLink<xAOD::VertexContainer> newPrimaryVertexLink = getPV_TJVA(pTau, *vxContainer ); + if (newPrimaryVertexLink.isValid()) { + // set new primary vertex + // will overwrite default one which was set above + pTau.setVertexLink(newPrimaryVertexLink); + // save highest JVF value + pTau.setDetail(xAOD::TauJetParameters::TauJetVtxFraction,static_cast<float>(m_maxJVF)); + ATH_MSG_DEBUG("TJVA vertex found and set"); + } + else { + ATH_MSG_DEBUG("couldn't find new PV for TJVA"); + } + + return StatusCode::SUCCESS; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ElementLink<xAOD::VertexContainer> TauVertexFinder::getPV_TJVA(const xAOD::TauJet& pTau, const xAOD::VertexContainer& vertices) { - const xAOD::Jet* pJetSeed = (*pTau.jetLink()); - - // the implementation follows closely the example given in modifyJet(...) in https://svnweb.cern.ch/trac/atlasoff/browser/Reconstruction/Jet/JetMomentTools/trunk/Root/JetVertexFractionTool.cxx#15 - - // Get the tracks associated to the jet - std::vector<const xAOD::TrackParticle*> assocTracks; - if (! pJetSeed->getAssociatedObjects(m_assocTracksName, assocTracks)) { - ATH_MSG_ERROR("Could not retrieve the AssociatedObjects named \""<< m_assocTracksName <<"\" from jet"); - return ElementLink<xAOD::VertexContainer>(); - } - - // Get the TVA object - const jet::TrackVertexAssociation* tva = NULL; - if (evtStore()->retrieve(tva,m_trackVertexAssocName).isFailure()) { - ATH_MSG_ERROR("Could not retrieve the TrackVertexAssociation from evtStore: " << m_trackVertexAssocName); - return ElementLink<xAOD::VertexContainer>(); - } - - // Calculate Jet Vertex Fraction - std::vector<float> jvf; - jvf.resize(vertices.size()); - for (size_t iVertex = 0; iVertex < vertices.size(); ++iVertex) { - jvf.at(iVertex) = getJetVertexFraction(vertices.at(iVertex),assocTracks,tva); - } + const xAOD::Jet* pJetSeed = (*pTau.jetLink()); + std::vector<const xAOD::TrackParticle*> tracksForTJVA; + const double dDeltaRMax(0.2); + + // the implementation follows closely the example given in modifyJet(...) in https://svnweb.cern.ch/trac/atlasoff/browser/Reconstruction/Jet/JetMomentTools/trunk/Root/JetVertexFractionTool.cxx#15 + + // Get the tracks associated to the jet + std::vector<const xAOD::TrackParticle*> assocTracks; + if (! pJetSeed->getAssociatedObjects(m_assocTracksName, assocTracks)) { + ATH_MSG_ERROR("Could not retrieve the AssociatedObjects named \""<< m_assocTracksName <<"\" from jet"); + return ElementLink<xAOD::VertexContainer>(); + } + + // Store tracks that meet TJVA track selection criteria and are between deltaR of 0.2 with the jet seed + // To be included in the TJVA calculation + // Maybe not as efficient as deleting unwanted tracks from assocTrack but quicker and safer for now. + for ( auto xTrack : assocTracks ){ + if ( (xTrack->p4().DeltaR(pJetSeed->p4())<dDeltaRMax) && m_TrackSelectionToolForTJVA->accept(*xTrack) ) + tracksForTJVA.push_back(xTrack); + } + + // Get the TVA object + const jet::TrackVertexAssociation* tva = NULL; + if (evtStore()->retrieve(tva,m_trackVertexAssocName).isFailure()) { + ATH_MSG_ERROR("Could not retrieve the TrackVertexAssociation from evtStore: " << m_trackVertexAssocName); + return ElementLink<xAOD::VertexContainer>(); + } + + // Calculate Jet Vertex Fraction + std::vector<float> jvf; + jvf.resize(vertices.size()); + for (size_t iVertex = 0; iVertex < vertices.size(); ++iVertex) { + jvf.at(iVertex) = getJetVertexFraction(vertices.at(iVertex),tracksForTJVA,tva); + } - // Get the highest JVF vertex and store maxJVF for later use - // Note: the official JetMomentTools/JetVertexFractionTool doesn't provide any possibility to access the JVF value, but just the vertex. - m_maxJVF=-100.; - size_t maxIndex = 0; - for (size_t iVertex = 0; iVertex < jvf.size(); ++iVertex) { - if (jvf.at(iVertex) > m_maxJVF) { - m_maxJVF = jvf.at(iVertex); - maxIndex = iVertex; - } + // Get the highest JVF vertex and store maxJVF for later use + // Note: the official JetMomentTools/JetVertexFractionTool doesn't provide any possibility to access the JVF value, but just the vertex. + m_maxJVF=-100.; + size_t maxIndex = 0; + for (size_t iVertex = 0; iVertex < jvf.size(); ++iVertex) { + if (jvf.at(iVertex) > m_maxJVF) { + m_maxJVF = jvf.at(iVertex); + maxIndex = iVertex; } + } - // Set the highest JVF vertex - ElementLink<xAOD::VertexContainer> vtxlink = ElementLink<xAOD::VertexContainer>(vertices,vertices.at(maxIndex)->index()); + // Set the highest JVF vertex + ElementLink<xAOD::VertexContainer> vtxlink = ElementLink<xAOD::VertexContainer>(vertices,vertices.at(maxIndex)->index()); - return vtxlink; + return vtxlink; } // reimplementation of JetVertexFractionTool::getJetVertexFraction(const xAOD::Vertex* vertex, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const @@ -172,40 +185,40 @@ ElementLink<xAOD::VertexContainer> TauVertexFinder::getPV_TJVA(const xAOD::TauJe // see https://svnweb.cern.ch/trac/atlasoff/browser/Reconstruction/Jet/JetMomentTools/trunk/Root/JetVertexFractionTool.cxx float TauVertexFinder::getJetVertexFraction(const xAOD::Vertex* vertex, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const { - float sumTrackPV = 0; - float sumTrackAll = 0; - for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) + float sumTrackPV = 0; + float sumTrackAll = 0; + for (size_t iTrack = 0; iTrack < tracks.size(); ++iTrack) { - const xAOD::TrackParticle* track = tracks.at(iTrack); - const xAOD::Vertex* ptvtx = tva->associatedVertex(track); - if (ptvtx != nullptr) { // C++11 feature - if (ptvtx->index() == vertex->index()) sumTrackPV += track->pt(); - } - sumTrackAll += track->pt(); + const xAOD::TrackParticle* track = tracks.at(iTrack); + const xAOD::Vertex* ptvtx = tva->associatedVertex(track); + if (ptvtx != nullptr) { // C++11 feature + if (ptvtx->index() == vertex->index()) sumTrackPV += track->pt(); + } + sumTrackAll += track->pt(); } - return sumTrackAll!=0 ? sumTrackPV/sumTrackAll : 0; + return sumTrackAll!=0 ? sumTrackPV/sumTrackAll : 0; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Helpers template <class T> bool TauVertexFinder::openContainer(T* &container, std::string containerName, bool printFATAL) { - StatusCode sc = evtStore()->retrieve(container, containerName); - if (sc.isFailure() || !container) { - if (printFATAL) ATH_MSG_FATAL("Container (" << containerName << ") not found in StoreGate"); - return 0; - } - return container; + StatusCode sc = evtStore()->retrieve(container, containerName); + if (sc.isFailure() || !container) { + if (printFATAL) ATH_MSG_FATAL("Container (" << containerName << ") not found in StoreGate"); + return 0; + } + return container; } template <class T> bool TauVertexFinder::retrieveTool(T & tool) { - if (tool.retrieve().isFailure()) { - ATH_MSG_FATAL("Failed to retrieve tool " << tool); - return false; - } else { - ATH_MSG_VERBOSE("Retrieved tool " << tool); - } - return true; + if (tool.retrieve().isFailure()) { + ATH_MSG_FATAL("Failed to retrieve tool " << tool); + return false; + } else { + ATH_MSG_VERBOSE("Retrieved tool " << tool); + } + return true; } diff --git a/Reconstruction/tauRecTools/src/TauVertexFinder.h b/Reconstruction/tauRecTools/src/TauVertexFinder.h index de30101ac0ea1b522a52a58993f13f7d5c819f75..8aca69d0bffe08df9a9376b938423920461c0bde 100644 --- a/Reconstruction/tauRecTools/src/TauVertexFinder.h +++ b/Reconstruction/tauRecTools/src/TauVertexFinder.h @@ -9,6 +9,7 @@ #include "GaudiKernel/ToolHandle.h" #include "xAODTracking/VertexContainer.h" #include "JetEDM/TrackVertexAssociation.h" +#include "InDetTrackSelectionTool/InDetTrackSelectionTool.h" ///////////////////////////////////////////////////////////////////////////// @@ -25,50 +26,51 @@ class TauVertexFinder : virtual public TauRecToolBase { public: - //------------------------------------------------------------- - //! Constructor and Destructor - //------------------------------------------------------------- - TauVertexFinder(const std::string& name); - ASG_TOOL_CLASS2(TauVertexFinder, TauRecToolBase, ITauToolBase); - ~TauVertexFinder(); + //------------------------------------------------------------- + //! Constructor and Destructor + //------------------------------------------------------------- + TauVertexFinder(const std::string& name); + ASG_TOOL_CLASS2(TauVertexFinder, TauRecToolBase, ITauToolBase); + ~TauVertexFinder(); - //------------------------------------------------------------- - //! Algorithm functions - //------------------------------------------------------------- - virtual StatusCode initialize(); - virtual StatusCode eventInitialize(); - virtual StatusCode execute(xAOD::TauJet& pTau); - virtual StatusCode eventFinalize(); - virtual StatusCode finalize(); + //------------------------------------------------------------- + //! Algorithm functions + //------------------------------------------------------------- + virtual StatusCode initialize(); + virtual StatusCode eventInitialize(); + virtual StatusCode execute(xAOD::TauJet& pTau); + virtual StatusCode eventFinalize(); + virtual StatusCode finalize(); - virtual void cleanup(xAOD::TauJet* ) { } - virtual void print() const { } + virtual void cleanup(xAOD::TauJet* ) { } + virtual void print() const { } - ElementLink<xAOD::VertexContainer> getPV_TJVA(const xAOD::TauJet& tauJet, const xAOD::VertexContainer& vertices); + ElementLink<xAOD::VertexContainer> getPV_TJVA(const xAOD::TauJet& tauJet, const xAOD::VertexContainer& vertices); private: - float getJetVertexFraction(const xAOD::Vertex* vertex, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const; + float getJetVertexFraction(const xAOD::Vertex* vertex, const std::vector<const xAOD::TrackParticle*>& tracks, const jet::TrackVertexAssociation* tva) const; - //------------------------------------------------------------- - //! Convenience functions to handle storegate objects - //------------------------------------------------------------- - template <class T> - bool openContainer(T* &container, std::string containerName, bool printFATAL=false); + //------------------------------------------------------------- + //! Convenience functions to handle storegate objects + //------------------------------------------------------------- + template <class T> + bool openContainer(T* &container, std::string containerName, bool printFATAL=false); - template <class T> - bool retrieveTool(T &tool); + template <class T> + bool retrieveTool(T &tool); private: - bool m_printMissingContainerINFO; - float m_maxJVF; + bool m_printMissingContainerINFO; + float m_maxJVF; + ToolHandle< InDet::IInDetTrackSelectionTool > m_TrackSelectionToolForTJVA; - //------------------------------------------------------------- - //! Configureables - //------------------------------------------------------------- - bool m_useTJVA; - std::string m_inputPrimaryVertexContainerName; - std::string m_assocTracksName; - std::string m_trackVertexAssocName; + //------------------------------------------------------------- + //! Configureables + //------------------------------------------------------------- + bool m_useTJVA; + std::string m_inputPrimaryVertexContainerName; + std::string m_assocTracksName; + std::string m_trackVertexAssocName; }; #endif // not TAUREC_TAUVERTEXFINDER_H diff --git a/Reconstruction/tauRecTools/src/TauVertexVariables.cxx b/Reconstruction/tauRecTools/src/TauVertexVariables.cxx index 5960a223fee02e17c4abb38a0ee4ef46234618e0..07ae82fd651c4c11d1be6607750eb7e32a2e2103 100644 --- a/Reconstruction/tauRecTools/src/TauVertexVariables.cxx +++ b/Reconstruction/tauRecTools/src/TauVertexVariables.cxx @@ -131,7 +131,7 @@ StatusCode TauVertexVariables::execute(xAOD::TauJet& pTau) { if (tauEventData()->hasObject("Beamspot")) scBeam = tauEventData()->getObject("Beamspot", vxcand); if(scBeam){ - myIPandSigma = m_trackToVertexIPEstimator->estimate(pTau.track(0), vxcand); + myIPandSigma = m_trackToVertexIPEstimator->estimate(pTau.track(0)->track(), vxcand); } else { ATH_MSG_DEBUG("No Beamspot object in tau candidate"); } @@ -142,7 +142,7 @@ StatusCode TauVertexVariables::execute(xAOD::TauJet& pTau) { //check if vertex has a valid type (skip if vertex has type NoVtx) if (vxcand->vertexType() > 0){ - myIPandSigma = m_trackToVertexIPEstimator->estimate(pTau.track(0), vxcand); + myIPandSigma = m_trackToVertexIPEstimator->estimate(pTau.track(0)->track(), vxcand); } } @@ -204,10 +204,10 @@ StatusCode TauVertexVariables::execute(xAOD::TauJet& pTau) { std::vector<const xAOD::TrackParticle*> xaodTracks; std::vector<const Trk::Track*> origTracks; for (unsigned i = 0; i < pTau.nTracks(); ++i) { - xaodTracks.push_back(pTau.track(i)); + xaodTracks.push_back(pTau.track(i)->track()); ATH_MSG_VERBOSE("xAOD::TrackParticle " <<i<<": "<< pTau.track(i)->pt() << " " << pTau.track(i)->eta() << " " << pTau.track(i)->phi()); if (pTau.track(i)->track()) { - origTracks.push_back(pTau.track(i)->track()); + origTracks.push_back(pTau.track(i)->track()->track()); // for debugging /* diff --git a/Reconstruction/tauRecTools/src/components/tauRecTools_entries.cxx b/Reconstruction/tauRecTools/src/components/tauRecTools_entries.cxx index 7f25d5e6416d9819dfd32d9b8a7a699ca6f7f397..a635f917a8385a0a24524b81898b76a0db95fec6 100644 --- a/Reconstruction/tauRecTools/src/components/tauRecTools_entries.cxx +++ b/Reconstruction/tauRecTools/src/components/tauRecTools_entries.cxx @@ -29,6 +29,7 @@ #include "tauRecTools/TauBuilderTool.h" #include "tauRecTools/MvaTESEvaluator.h" #include "tauRecTools/MvaTESVariableDecorator.h" +#include "tauRecTools/TauTrackClassifier.h" #include "tauRecTools/CombinedP4FromRecoTaus.h" @@ -65,4 +66,6 @@ DECLARE_TOOL_FACTORY( TauProcessorTool ) DECLARE_TOOL_FACTORY( TauBuilderTool ) DECLARE_TOOL_FACTORY( MvaTESVariableDecorator ) DECLARE_TOOL_FACTORY( MvaTESEvaluator ) +DECLARE_NAMESPACE_TOOL_FACTORY( tauRecTools, TauTrackClassifier ) +DECLARE_NAMESPACE_TOOL_FACTORY( tauRecTools, TrackMVABDT ) DECLARE_TOOL_FACTORY( CombinedP4FromRecoTaus ) diff --git a/Reconstruction/tauRecTools/tauRecTools/ITauToolBase.h b/Reconstruction/tauRecTools/tauRecTools/ITauToolBase.h index 2125058976c27e4a35f884dec99616e0d359ef05..953f17d1b1f25bf5a05b082d5e0ac422c61b1011 100644 --- a/Reconstruction/tauRecTools/tauRecTools/ITauToolBase.h +++ b/Reconstruction/tauRecTools/tauRecTools/ITauToolBase.h @@ -51,8 +51,8 @@ class ITauToolBase : virtual public asg::IAsgTool virtual void setTauEventData(TauEventData* data) = 0; - //make virtual - virtual StatusCode readConfig() {return StatusCode::SUCCESS; } + //make pure + virtual StatusCode readConfig() = 0; }; diff --git a/Reconstruction/tauRecTools/tauRecTools/MvaTESVariableDecorator.h b/Reconstruction/tauRecTools/tauRecTools/MvaTESVariableDecorator.h index 980d2bb90848d5b6fd00ab5d5902b92e41a8ec21..94eab92a82f3959ec668a6661d4f9599386428e7 100644 --- a/Reconstruction/tauRecTools/tauRecTools/MvaTESVariableDecorator.h +++ b/Reconstruction/tauRecTools/tauRecTools/MvaTESVariableDecorator.h @@ -30,6 +30,8 @@ class MvaTESVariableDecorator const xAOD::VertexContainer* m_xVertexContainer; //! int m_mu; //! int m_nVtx; //! + + }; diff --git a/Reconstruction/tauRecTools/tauRecTools/TauBuilderTool.h b/Reconstruction/tauRecTools/tauRecTools/TauBuilderTool.h index c8776c84da2aafdcbefc7d47903980906cbdb2b9..73229c053429f20a2ca2ba885c4f0b00f0eeb5ef 100644 --- a/Reconstruction/tauRecTools/tauRecTools/TauBuilderTool.h +++ b/Reconstruction/tauRecTools/tauRecTools/TauBuilderTool.h @@ -32,6 +32,8 @@ class TauBuilderTool : public asg::AsgTool, virtual public ITauToolExecBase { private: std :: string m_tauContainerName; //!< tau output container std :: string m_tauAuxContainerName; //!< tau output aux store container + std::string m_tauTrackContainerName; + std::string m_tauTrackAuxContainerName; std::string m_seedContainerName; //!< seed input container double m_maxEta; //!< only build taus with eta_seed < m_maxeta double m_minPt; //!< only build taus with pt_seed > m_minpt diff --git a/Reconstruction/tauRecTools/tauRecTools/TauCalibrateLC.h b/Reconstruction/tauRecTools/tauRecTools/TauCalibrateLC.h index 5c76e9702effb0d0dde951068a61faee0cfd9d67..0599d63f1fb7b875f3fc5d4b20a59bb55fe1e469 100644 --- a/Reconstruction/tauRecTools/tauRecTools/TauCalibrateLC.h +++ b/Reconstruction/tauRecTools/tauRecTools/TauCalibrateLC.h @@ -25,8 +25,7 @@ public: ASG_TOOL_CLASS2( TauCalibrateLC, TauRecToolBase, ITauToolBase ) - TauCalibrateLC(const std::string& type); - TauCalibrateLC(); + TauCalibrateLC(const std::string& name="TauCalibrateLC"); ~TauCalibrateLC(); virtual StatusCode initialize(); diff --git a/Reconstruction/tauRecTools/tauRecTools/TauProcessorTool.h b/Reconstruction/tauRecTools/tauRecTools/TauProcessorTool.h index 5b4ec84d8e474ca3f1776c27d40d532ca2f315bf..adfe60f6b6a54fed30e00c0cf014dc19ae1a914b 100644 --- a/Reconstruction/tauRecTools/tauRecTools/TauProcessorTool.h +++ b/Reconstruction/tauRecTools/tauRecTools/TauProcessorTool.h @@ -41,6 +41,7 @@ class TauProcessorTool : public asg::AsgTool, virtual public ITauToolExecBase { bool m_deep_copy_hadronicPFOContainer; bool m_deep_copy_neutralPFOContainer; bool m_deep_copy_SecVtxContainer; + bool m_deep_copy_TauTrackContainer; TauEventData m_data; ToolHandleArray<ITauToolBase> m_tools; diff --git a/Reconstruction/tauRecTools/tauRecTools/TauTrackClassifier.h b/Reconstruction/tauRecTools/tauRecTools/TauTrackClassifier.h new file mode 100644 index 0000000000000000000000000000000000000000..e87aaaea93e30c10a39446c5e4e42b751d3780be --- /dev/null +++ b/Reconstruction/tauRecTools/tauRecTools/TauTrackClassifier.h @@ -0,0 +1,159 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TAUREC_TAUTRACKCLASSIFIER_H +#define TAUREC_TAUTRACKCLASSIFIER_H + +// ASG include(s) +#include "AsgTools/AsgTool.h" +#include "AsgTools/ToolHandleArray.h" + +// xAOD include(s) +#include "xAODTau/TauJet.h" +#include "xAODTau/TauTrack.h" + +// local include(s) +#include "tauRecTools/TauRecToolBase.h" + +// ROOT include(s) +#include "TMVA/Reader.h" + +/** + * @brief Implementation of a TrackClassifier based on an MVA + * + * @author Dirk Duschinger + * + */ + +namespace tauRecTools +{ + +class TrackMVABDT; + +//______________________________________________________________________________ +class TauTrackClassifier + : virtual public TauRecToolBase +{ +public: + + ASG_TOOL_CLASS2( TauTrackClassifier, TauRecToolBase, ITauToolBase ) + + TauTrackClassifier(const std::string& sName="TauTrackClassifier"); + ~TauTrackClassifier(); + + // retrieve all track classifier sub tools + virtual StatusCode initialize(); + // pass all tracks in the tau cone to all track classifier sub tools + virtual StatusCode execute(xAOD::TauJet& pTau); + +private: + ToolHandleArray<TrackMVABDT> m_vClassifier; + std::string m_tauTrackConName; + +}; // class TauTrackClassifier + +//______________________________________________________________________________ +class TrackMVABDT + : public TauRecToolBase +{ + /// Create a proper constructor for Athena + ASG_TOOL_CLASS2( TrackMVABDT, + TauRecToolBase, + ITauToolBase) + + public: + + TrackMVABDT(const std::string& sName); + ~TrackMVABDT(); + + // configure the TMVA reader object and build a general map to store variables + // for possible TMVA inputs. Only Variables defined in the xml weights file + // are passed to the TMVA reader + StatusCode initialize(); + + // executes TMVA reader to get the BDT score, makes the decision and resets + // classification flags + StatusCode classifyTrack(xAOD::TauTrack& xTrack, const xAOD::TauJet& xTau); + // set BDT input variables in the corresponding map entries + void setVars(const xAOD::TauTrack& xTrack, const xAOD::TauJet& xTau); + + // load the xml weights file and configure the TMVA reader with the correct + // variable addresses + StatusCode addWeightsFile(); + // parse the weights file for the line showing the input variable used by that + // particular BDT names and store them + StatusCode parseVariableContent(); + +private: + // configurable variables + std::string m_sInputWeightsPath; + float m_fThreshold; + int m_iSignalType; + int m_iBackgroundType; + int m_iExpectedFlag; + +private: + TMVA::Reader* m_rReader; //! + + std::map<int, std::string> m_mParsedVarsBDT; //! + std::map<std::string, float> m_mAvailableVars; //! + + // possible bdt input variables + float tauPt; + float tauEta; + float trackEta; + float z0sinThetaTJVA; + float rConv; + float rConvII; + float DRJetSeedAxis; + float d0; + float qOverP; + float theta; + float numberOfInnermostPixelLayerHits; + float numberOfPixelHits; + float numberOfPixelDeadSensors; + float numberOfPixelSharedHits; + float numberOfSCTHits; + float numberOfSCTDeadSensors; + float numberOfSCTSharedHits; + float numberOfTRTHighThresholdHits; + float numberOfTRTHits; + float nPixHits; + float nSiHits; + +}; // class TrackMVABDT + + + xAOD::TauTrack::TrackFlagType isolateClassifiedBits(xAOD::TauTrack::TrackFlagType flag){ + static int flagsize=sizeof(flag)*8; + flag=flag<<(flagsize-xAOD::TauJetParameters::classifiedFake-1); + flag=flag>>(flagsize-xAOD::TauJetParameters::classifiedCharged+1); + return flag; + } + +//bool sortTracks(xAOD::TauTrack* xTrack1, xAOD::TauTrack* xTrack2) + bool sortTracks(const ElementLink<xAOD::TauTrackContainer> &l1, const ElementLink<xAOD::TauTrackContainer> &l2) +{ + //should we be safe and ask if the links are available? + const xAOD::TauTrack* xTrack1 = *l1; + const xAOD::TauTrack* xTrack2 = *l2; + + //return classified charged, then isolation (wide tracks), else by pt + xAOD::TauTrack::TrackFlagType f1 = isolateClassifiedBits(xTrack1->flagSet()); + xAOD::TauTrack::TrackFlagType f2 = isolateClassifiedBits(xTrack2->flagSet()); + + if(f1==f2) + return xTrack1->pt()>xTrack2->pt(); + return f1<f2; + + //this somehow causes a crash + /* static uint16_t flag1 = xTrack1->flagSet() >> (xAOD::TauJetParameters::classifiedCharged - 1); */ + /* static uint16_t flag2 = xTrack2->flagSet() >> (uint16_t(xAOD::TauJetParameters::classifiedCharged) - 1); */ + /* return (flag1<flag2) || // sort by type, true tracks first */ + /* ((flag1==flag2) && (xTrack1->pt()>xTrack2->pt())); // sort by pt if equal types */ +} + +} // namespace tauRecTools + +#endif // not TAUTRACKCLASSIFIER diff --git a/Reconstruction/tauRecTools/tauRecTools/TauTrackFilter.h b/Reconstruction/tauRecTools/tauRecTools/TauTrackFilter.h index 907309a8e8f70095a5d72eb4abe95af5db5dd63b..a7d39dda7ee7aadf1ad524c2e765bf07f21db33f 100644 --- a/Reconstruction/tauRecTools/tauRecTools/TauTrackFilter.h +++ b/Reconstruction/tauRecTools/tauRecTools/TauTrackFilter.h @@ -40,6 +40,7 @@ public: private: std::string m_configPath; std::string m_trackContainerName; + std::string m_tauTrackConName; std::vector<bool> m_TrkPass; int m_nProng; int m_flag;