Adding tracktwoMVABDT chains to master

This is a MR to add new flavour of tau chains tracktwoMVABDT where the track based Roi updater makes use of a BDT selected tau track instead of just taking the leading one. This MR also superseed the one from @adsalvad and takes most of the code from !35872 (closed)

Trigger/TrigHypothesis/TrigTauHypo/src/TrigTauTrackRoiUpdaterMT.cxx

@@ -19,6 +19,8 @@
 #include "TrkTrackSummary/TrackSummary.h"
 #include "xAODTrigger/TrigPassBits.h"
 
+#include "PathResolver/PathResolver.h"
+#include "tauRecTools/HelperFunctions.h"
 
 TrigTauTrackRoiUpdaterMT::TrigTauTrackRoiUpdaterMT(const std::string & name, ISvcLocator* pSvcLocator) :
   AthReentrantAlgorithm(name, pSvcLocator)
@@ -41,13 +43,22 @@ StatusCode TrigTauTrackRoiUpdaterMT::initialize()
   ATH_MSG_DEBUG( " nSiHoles            		" 	<< m_nSiHoles );
   ATH_MSG_DEBUG( " UpdateEta            		" 	<< m_updateEta );
   ATH_MSG_DEBUG( " UpdatePhi            		" 	<< m_updatePhi );
+  ATH_MSG_DEBUG(  " useBDT               		" 	<< m_useBDT );
 
+  if(m_useBDT) {
+    ATH_MSG_DEBUG( " BDTweights     	        	" << m_BDTweights );
+
+    std::string inputWeightsPath = PathResolverFindCalibFile(m_BDTweights);
+    ATH_MSG_INFO(" InputWeightsPath: " << inputWeightsPath);
+    m_reader = std::make_unique<tauRecTools::BDTHelper>();
+    ATH_CHECK(m_reader->initialize(inputWeightsPath));
+  }
 
   ATH_MSG_DEBUG( "Initialising HandleKeys" );
   CHECK( m_roIInputKey.initialize()        );
-  CHECK( m_tracksKey.initialize() );
-  CHECK( m_roIOutputKey.initialize()  );
-
+  CHECK( m_tracksKey.initialize()          );
+  CHECK( m_roIOutputKey.initialize()       );
+  CHECK( m_tauJetKey.initialize(!m_tauJetKey.key().empty()));
   return StatusCode::SUCCESS;
 
 }
@@ -100,38 +111,67 @@ StatusCode TrigTauTrackRoiUpdaterMT::execute(const EventContext& ctx) const
   double trkPtMax = 0;
   
   if(foundTracks!=0){
-    
+
+    //Find the track with the highest BDT score
+    if(m_useBDT) {
+      
+      // retrieve TauJet from TrigTauRecCaloOnlyMVASequence
+      SG::ReadHandle< xAOD::TauJetContainer > tauJetHandle = SG::makeHandle( m_tauJetKey,ctx );
+      const xAOD::TauJetContainer *foundTaus = tauJetHandle.get();
+
+      if (foundTaus->size()==0) {
+	      ATH_MSG_ERROR( "No TauJetContainer feature was found." );
+      }
+
+      const Trk::Track *leadTrackBDT = nullptr;
+      for(const Trk::Track* trk_it : *foundTracks)
+        if(trk_it->perigeeParameters()->pT() > trkPtMax) {
+          trkPtMax = trk_it->perigeeParameters()->pT();
+          leadTrackBDT = trk_it;
+        }
+      double trkBDTMax = 0.;
+
+      for(const Trk::Track* trk_it : *foundTracks) {
+        double BDTscore = getBDTscore( foundTaus->at(0), trk_it, leadTrackBDT );
+
+        if(BDTscore > trkBDTMax) {
+          trkBDTMax = BDTscore;
+          leadTrack = trk_it;
+        }
+      }
+    }
     // Find leading track
-    for (const Trk::Track* track : *foundTracks){
-      trackPer = track->perigeeParameters();
-      summary = track->trackSummary();
-      if(summary==nullptr) ATH_MSG_ERROR(" track summary not available in RoI updater " << trkPtMax << " " << name() );
-      if(trackPer && summary){
- 	     float trackPt = trackPer->pT();
-	     if ( trackPt > trkPtMax ) {
-	       int nPix  = summary->get(Trk::numberOfPixelHits);
-	       if(nPix<0) nPix=0;
-	       if(nPix < m_nHitPix) {
-            if(msgLvl(MSG::DEBUG))  msg() << MSG::DEBUG <<" Track rejected because of nHitPix " << nPix << " < " << m_nHitPix
-            << " " << name() << ": in execute()" << endmsg;
-	         continue;
-	       }
-	       int nPixHole = summary->get(Trk::numberOfPixelHoles);
-	       if (nPixHole < 0) nPixHole = 0;
-	       int nSCTHole = summary->get(Trk::numberOfSCTHoles);
-	       if (nSCTHole < 0) nSCTHole = 0;
-	       if((nPixHole + nSCTHole) > m_nSiHoles) {
-            if(msgLvl(MSG::DEBUG) )  msg() << MSG::DEBUG <<" Track rejected because of nSiHoles " << nPixHole+nSCTHole << " > " << m_nSiHoles 
-            << " " << name() << ": in execute()" << endmsg;
-	         continue;
-	       }
-	       leadTrack = track;
-	       trkPtMax = trackPt;
-	     }
+    else{
+      for (const Trk::Track* track : *foundTracks){
+        trackPer = track->perigeeParameters();
+        summary = track->trackSummary();
+        if(summary==nullptr) ATH_MSG_ERROR(" track summary not available in RoI updater " << trkPtMax << " " << name() );
+        if(trackPer && summary){
+        float trackPt = trackPer->pT();
+        if ( trackPt > trkPtMax ) {
+          int nPix  = summary->get(Trk::numberOfPixelHits);
+          if(nPix<0) nPix=0;
+          if(nPix < m_nHitPix) {
+              if(msgLvl(MSG::DEBUG))  msg() << MSG::DEBUG <<" Track rejected because of nHitPix " << nPix << " < " << m_nHitPix
+              << " " << name() << ": in execute()" << endmsg;
+            continue;
+          }
+          int nPixHole = summary->get(Trk::numberOfPixelHoles);
+          if (nPixHole < 0) nPixHole = 0;
+          int nSCTHole = summary->get(Trk::numberOfSCTHoles);
+          if (nSCTHole < 0) nSCTHole = 0;
+          if((nPixHole + nSCTHole) > m_nSiHoles) {
+              if(msgLvl(MSG::DEBUG) )  msg() << MSG::DEBUG <<" Track rejected because of nSiHoles " << nPixHole+nSCTHole << " > " << m_nSiHoles 
+              << " " << name() << ": in execute()" << endmsg;
+            continue;
+          }
+          leadTrack = track;
+          trkPtMax = trackPt;
+        }
+        }
       }
     }
-
-
+    
     if(leadTrack) {
       ATH_MSG_DEBUG(" leading track pT " << trkPtMax );
       leadTrkZ0  = leadTrack->perigeeParameters()->parameters()[Trk::z0];
@@ -191,3 +231,67 @@ StatusCode TrigTauTrackRoiUpdaterMT::execute(const EventContext& ctx) const
   return StatusCode::SUCCESS;
 
 }
+
+double TrigTauTrackRoiUpdaterMT::getBDTscore(const xAOD::TauJet* tau, const Trk::Track* track, const Trk::Track* leadtrack ) const
+{
+
+  std::map<TString, float*> BDTvars;
+  BDTInputVariables vars;
+    
+  BDTvars = {
+    {"log(Coretrack_pt)", &vars.logtrk_pt}, 
+    {"fabs(Coretrack_z0)", &vars.abstrck_z0},
+    {"fabs(Coretrack_d0)", &vars.abstrk_d0},
+    {"Coretrack_nPiHits", &vars.trk_nPiHits},
+    {"Coretrack_nSiHoles", &vars.trk_nSiHoles},
+    {"log(Coretrack_ratioptCalo)", &vars.logtrk_ratiopt},
+    {"Coretrack_dR", &vars.trk_dR},
+    {"Coretrack_dRleadtrk", &vars.trk_dRtoleadtrk},
+    {"Coretrack_CaloHadpt", &vars.CaloHad_pt},
+    {"Coretrack_CaloEMpt", &vars.CaloEM_pt}
+  };
+
+  const Trk::Perigee* trkPerigee = track->perigeeParameters();
+  const Trk::TrackSummary* trkSummary = track->trackSummary();
+
+  int nPixHit = trkSummary->get(Trk::numberOfPixelHits);
+  int nPixDead = trkSummary->get(Trk::numberOfPixelDeadSensors);
+
+  int nPixHole = trkSummary->get(Trk::numberOfPixelHoles);
+  int nSCTHole = trkSummary->get(Trk::numberOfSCTHoles);
+
+  float ratio_pt = tau->ptTrigCaloOnly() ? trkPerigee->pT()/tau->ptTrigCaloOnly() : 0.;
+
+  float dEta = tau->eta() - trkPerigee->eta();
+  float dPhi = CxxUtils::wrapToPi(tau->phi() - trkPerigee->parameters()[Trk::phi0]);
+  float dR = sqrt(dEta*dEta + dPhi*dPhi);
+
+  const Trk::Perigee* leadtrkPerigee = leadtrack->perigeeParameters();
+  float dEta_leadtrk = trkPerigee->eta() - leadtrkPerigee->eta();
+  float dPhi_leadtrk = CxxUtils::wrapToPi(trkPerigee->parameters()[Trk::phi0] - leadtrkPerigee->parameters()[Trk::phi0]);
+  float dR_leadtrk = sqrt(dEta_leadtrk*dEta_leadtrk + dPhi_leadtrk*dPhi_leadtrk);
+
+  float tau_emscale_ptEM = 0;
+  float tau_emscale_ptHad = 0;
+  if ( !tau->detail( xAOD::TauJetParameters::etEMAtEMScale, tau_emscale_ptEM ) ) {
+    ATH_MSG_WARNING("Retrieval of tau etEMAtEMScale detail failed.");
+  }
+  if ( !tau->detail( xAOD::TauJetParameters::etHadAtEMScale, tau_emscale_ptHad ) ) {
+    ATH_MSG_WARNING("Retrieval of tau etHadAtEMScale detail failed.");
+  }
+
+  vars.logtrk_pt = log( trkPerigee->pT() );
+  vars.abstrck_z0 = std::abs( trkPerigee->parameters()[Trk::z0] );
+  vars.abstrk_d0 = std::abs( trkPerigee->parameters()[Trk::d0] );
+  vars.trk_nPiHits = nPixHit + nPixDead;
+  vars.trk_nSiHoles = nPixHole + nSCTHole;
+  vars.logtrk_ratiopt = log( ratio_pt );
+  vars.trk_dR = dR;
+  vars.trk_dRtoleadtrk = dR_leadtrk;
+  vars.CaloHad_pt = tau_emscale_ptHad;
+  vars.CaloEM_pt = tau_emscale_ptEM;
+
+  double BDTval = m_reader->getResponse(BDTvars);
+
+  return BDTval;
+}