TrigFastTrackFinder.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/

// TrigFastTrackFinder.cxx
// -------------------------------
// ATLAS Collaboration
//
// package created 16/04/2013 by Dmitry Emeliyanov (see ChangeLog for more details)
//
 
 
 
#include "TrigFastTrackFinder.h"
 
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "InDetPrepRawData/SCT_Cluster.h"
#include "InDetPrepRawData/PixelCluster.h"
#include "InDetRIO_OnTrack/SiClusterOnTrack.h"
#include "xAODTrigger/TrigCompositeAuxContainer.h"
#include "TrigInDetPattRecoEvent/TrigInDetTriplet.h"
#include "TrigInDetToolInterfaces/TrigL2HitResidual.h"
#include "TrigInDetEvent/TrigSiSpacePointBase.h"
#include "TrigInDetPattRecoEvent/TrigInDetTracklet.h"
#include "TrigInDetPattRecoTools/TrigTrackSeedGenerator.h"
 
//
#include "InDetIdentifier/SCT_ID.h"
#include "InDetIdentifier/PixelID.h"
#include "TrkTrack/Track.h"
 
#include "TrkParameters/TrackParameters.h"
#include "TrkTrack/Track.h"
 
#include "CxxUtils/phihelper.h"
 
#include "AthenaMonitoringKernel/Monitored.h"
 
#include "TrigInDetPattRecoTools/TrigInDetUtils.h"
 
#include "PathResolver/PathResolver.h"
 
//for UTT
#include "InDetRIO_OnTrack/PixelClusterOnTrack.h"
#include "InDetRIO_OnTrack/SCT_ClusterOnTrack.h"
 
#include <cmath>
#include <iostream>
#include <algorithm>
#include <memory>
 
TrigFastTrackFinder::TrigFastTrackFinder(const std::string& name, ISvcLocator* pSvcLocator) :
 
  AthReentrantAlgorithm(name, pSvcLocator),
  m_trackMaker("InDet::SiTrackMaker_xk/InDetTrigSiTrackMaker"),
  m_trigInDetTrackFitter("TrigInDetTrackFitter"),
  m_trigZFinder("TrigZFinder/TrigZFinder", this ),
  m_trackSummaryTool("Trk::ITrackSummaryTool/ITrackSummaryTool"),
  m_disTrkFitter("Trk::GlobalChi2Fitter/InDetTrackFitter"),
  m_useBeamSpot(true),
  m_doZFinder(false),
  m_doZFinderOnly(false),
  m_storeZFinderVertices(false),
  m_nfreeCut(5),
  m_countTotalRoI(0),
  m_countRoIwithEnoughHits(0),
  m_countRoIwithTracks(0),
  m_pixelId(0),
  m_sctId(0),
  m_idHelper(0),
  m_particleHypothesis(Trk::pion),
  m_useNewLayerNumberScheme(false),
  m_useGPU(false),
  m_LRTmode(false),
  m_LRTD0Min(0.0),
  m_LRTHardMinPt(0.0),
  m_trigseedML_LUT(""),
  m_dodEdxTrk(false),
  m_doDisappearingTrk(false),
  m_ITkMode(false),
  m_standaloneMode(false)
{

  declareProperty( "MinHits",               m_minHits = 5,"Minimum number of hits needed to perform tracking"  );
 
  //** Zfinder mode
  declareProperty( "doZFinder",               m_doZFinder = true,"Use fast ZFinder to find z of primary vertices");
  declareProperty( "doZFinderOnly",           m_doZFinderOnly = false,"stop processing after ZFinder - no tracking performed");
  declareProperty( "VertexSeededMode",        m_vertexSeededMode = false); //** NOT USED Obsolete? ATR-24242
  declareProperty( "doFastZVertexSeeding",    m_doFastZVseeding = true,"Use ZFinder vertex information to filter seeds");
  declareProperty( "zVertexResolution",       m_tcs.m_zvError = 10.0," Half-width (mm) in z of z region used to filter seeds when doFastZVertexSeeding enabled" );
  declareProperty( "zVertexResolutionEndcap", m_tcs.m_zvErrorEndcap = -1," Half-width (mm) in z of region used to filter seeds when doFastZVertexSeeding enabled, for endcap pixels; set to zVertexResolution value later if left negative" );
  declareProperty( "StoreZFinderVertices",    m_storeZFinderVertices = false ); //** NOT USED - to be implemented ATR-24242
  

  declareProperty("useNewLayerNumberScheme", m_useNewLayerNumberScheme = false,"Use LayerNumberTool for layer numbers");

  declareProperty("Doublet_FilterRZ",          m_tcs.m_doubletFilterRZ = true,"Enable check that doublet is consistent with the RoI in the RZ plane");
  declareProperty("DoubletDR_Max",             m_tcs.m_doublet_dR_Max = 270.0,"Maximum Radial distance between spacepoints forming Doublets");
  declareProperty("SeedRadBinWidth",           m_tcs.m_seedRadBinWidth = 2.0);

  declareProperty("Triplet_D0Max",            m_tcs.m_tripletD0Max      = 4.0,"Maximum d0 for triplet");
  declareProperty("Triplet_D0_PPS_Max",       m_tcs.m_tripletD0_PPS_Max = 1.7,"Maximin d0 for PPS doublets");
  declareProperty("Triplet_nMaxPhiSlice",     m_tcs.m_nMaxPhiSlice = 53,"Number of phi-slices used for seeding");
  declareProperty("Triplet_MaxBufferLength",     m_tcs.m_maxTripletBufferLength = 3,"Maximum number of triplets sharing a common middle spacepoint");
  declareProperty("TripletDoPSS",            m_tcs.m_tripletDoPSS = false,"Allow PSS Triplet seeds");
  declareProperty("TripletDoPPS",            m_tcs.m_tripletDoPPS = true,"Allow PPS triplet seeds");
  declareProperty("TripletDoConfirm",        m_tcs.m_tripletDoConfirm = false,"Enable triplet confirmation");
  declareProperty("DoubletDR_Max_Confirm",     m_tcs.m_doublet_dR_Max_Confirm = 150.0,"doublet max DR when TripletDoConfirm enabled");
  declareProperty("TripletMaxCurvatureDiff",    m_tcs.m_curv_delta = 0.001,"Maximum curvature difference allowed in seed confirmation");//for the triplet confirmation
  declareProperty("Triplet_DtCut",            m_tcs.m_tripletDtCut      = 10.0);//i.e. 10*sigma_MS
  declareProperty("pTmin",                    m_pTmin = 1000.0,"Triplet pT threshold is pTmin*Triplet_MinPtFrac" );
  declareProperty("Triplet_MinPtFrac",        m_tripletMinPtFrac = 0.3,"Triplet pT threshold is pTmin*Triplet_MinPtFrac");
  declareProperty("doSeedRedundancyCheck",    m_checkSeedRedundancy = false,"skip Triplets already used in a track");

  declareProperty("UseTrigSeedML",              m_tcs.m_useTrigSeedML = 0,"set ML-based seed selection mode (0 disables)" );
  declareProperty("TrigSeedML_LUT",             m_trigseedML_LUT = "trigseed_ml_pixel_barrel_kde.lut","LUT used by ML-based seed selection");
  declareProperty("maxEC_Pixel_cluster_length", m_tcs.m_maxEC_len = 1.5,"Maximum Endcap Pixel cluster length for ML-based seed selection" );
 
 
  //* Clone removal (removal of tracks sharing too many hits */
  
  declareProperty( "FreeClustersCut"   ,m_nfreeCut,"Minimum number of unshared clusters");
 
  //** Cuts applied to final tracks after Fit
  declareProperty("TrackInitialD0Max",            m_initialD0Max      = 10.0,"Maximum d0 of track");
  declareProperty("TrackZ0Max",                   m_Z0Max      = 300.0,"Maximum z0 of track");
 
 
  /* Monitoring */
  declareProperty( "doResMon",       m_doResMonitoring = true,"enable unbiased residual monitoring");
  declareProperty( "UseBeamSpot",           m_useBeamSpot = true,"Monitor d0 with respect to beamspot");
 
  //* Collection Names */
  declareProperty("TracksName",
                  m_outputTracksKey = std::string("TrigFastTrackFinder_Tracks"),
                  "TrackCollection name");
 
  declareProperty("inputTracksName",
                  m_inputTracksKey = std::string(""),
                  "TrackCollection name");
 
  declareProperty("RoIs", m_roiCollectionKey = std::string("OutputRoIs"), "RoIs to read in");
 
  //* Tools */
  declareProperty( "initialTrackMaker", m_trackMaker);
  declareProperty( "trigInDetTrackFitter",   m_trigInDetTrackFitter );
  declareProperty( "trigZFinder",   m_trigZFinder );
  declareProperty("TrackSummaryTool", m_trackSummaryTool);
  
  // Accleration
  declareProperty("useGPU", m_useGPU = false,"Use GPU acceleration");
 
  // Large Radius Tracking
  declareProperty("LRT_Mode", m_LRTmode,"Enable Large Radius Tracking mode" );
  declareProperty("LRT_D0Min", m_LRTD0Min=0.0,"Minimum d0 for tracks to be saved in LRT Mode" );
  declareProperty("LRT_HardMinPt", m_LRTHardMinPt=0.0,"Minimum pT for tracks to be saved in LRT Mode");
 
  // UTT
  declareProperty("dodEdxTrk",         m_dodEdxTrk         = false);
  declareProperty("doDisappearingTrk", m_doDisappearingTrk = false);
  declareProperty("DisTrackFitter",    m_disTrkFitter );
 
  // Phase II
  declareProperty("ITkMode",           m_ITkMode = false);
  declareProperty("StandaloneMode",    m_standaloneMode = false);
  
}
 
//--------------------------------------------------------------------------
 
TrigFastTrackFinder::~TrigFastTrackFinder() {}
 
//-----------------------------------------------------------------------
 
StatusCode TrigFastTrackFinder::initialize() {
 
  ATH_CHECK(m_roiCollectionKey.initialize(!m_standaloneMode) );
  ATH_CHECK(m_outputTracksKey.initialize() );
 
  // optional input tracks collection if present the clusters on previously found tracks are not used to form seeds
  if (m_LRTmode) {
    ATH_CHECK(m_inputTracksKey.initialize( !m_inputTracksKey.key().empty() ) );
  }
  // optional PRD to track association map
  ATH_CHECK(m_prdToTrackMap.initialize( !m_prdToTrackMap.key().empty() ));
 
 
  ATH_CHECK(m_beamSpotKey.initialize());
 
  ATH_MSG_DEBUG(" TrigFastTrackFinder : MinHits set to " << m_minHits);
 
  ATH_CHECK(m_numberingTool.retrieve());
 
  if(!m_ITkMode) {
    ATH_CHECK(m_spacePointTool.retrieve());
  }
  else {
    m_spacePointTool.disable();
  }
  
  ATH_CHECK(m_trackMaker.retrieve());
 
  ATH_CHECK(m_trigInDetTrackFitter.retrieve());
  
  if (m_doZFinder) {
    ATH_CHECK(m_trigZFinder.retrieve());
    // If m_tcs.m_zvErrorEndcap has negative default value, it was not set by user,
    // so set it to the same value as m_tcs.m_zvError
    if (m_tcs.m_zvErrorEndcap < 0) m_tcs.m_zvErrorEndcap = m_tcs.m_zvError;
  } else {
    m_trigZFinder.disable();
  }
  
  ATH_CHECK(m_trackSummaryTool.retrieve());
 
  //Get ID helper
  ATH_CHECK(detStore()->retrieve(m_idHelper, "AtlasID"));
 
  ATH_CHECK(detStore()->retrieve(m_pixelId, "PixelID"));
 
  ATH_CHECK(detStore()->retrieve(m_sctId, "SCT_ID"));
 
  // monitoring
  if ( !m_monTool.empty() ) {
     ATH_CHECK(m_monTool.retrieve() );
  }
  else {
     ATH_MSG_INFO("Monitoring tool is empty");
  }
 
  ATH_MSG_DEBUG(" doResMon " << m_doResMonitoring);
 
  if(m_useGPU) {//for GPU acceleration
    ATH_CHECK(m_accelSvc.retrieve());
    ATH_CHECK(m_accelTool.retrieve());
  }
 
  ATH_MSG_INFO("Use GPU acceleration : "<<std::boolalpha<<m_useGPU);
 
  if (m_LRTmode) {
    ATH_MSG_INFO(" FTF configures in Large Radius Tracking Mode");
    // set TrigTrackSeedGenerator to LRTmode
    m_tcs.m_LRTmode=m_LRTmode;
 
  }
 
  if(m_tcs.m_useTrigSeedML > 0) {
    //LUT params
    int lut_w = 30;
    int lut_h = 45;
    float lut_range[4] = {0.0,3.0,0.0,9.0};
    TrigSeedML_LUT L(1,lut_w,lut_h,lut_range);
    //read data from LUT file
    std::string lut_fileName = PathResolver::find_file(m_trigseedML_LUT, "DATAPATH");
    if (lut_fileName.empty()) {
      ATH_MSG_ERROR("Cannot find TrigSeedML LUT file " << lut_fileName);
      return StatusCode::FAILURE;
    }
    else {
      ATH_MSG_INFO(lut_fileName);
      std::ifstream ifs(lut_fileName.c_str());
      int row, col0, col1;
      while(!ifs.eof()) {
        ifs >> row >> col0 >> col1;
        if(ifs.eof()) break;
        for(int c=col0;c<=col1;c++) L.setBin(row, c);
      }
      ifs.close();
      ATH_MSG_INFO("TrigSeedML LUT initialized from file " << m_trigseedML_LUT);
      m_tcs.m_vLUT.push_back(L);
    }
  }
  if (m_ITkMode) {
    ATH_CHECK(m_seedingTool.retrieve());
  } else {
    m_seedingTool.disable();
  }
 
  // UTT tools
  if( m_doDisappearingTrk ) {
     ATH_CHECK(m_extrapolator.retrieve());
     ATH_MSG_DEBUG("Retrieved tool " << m_extrapolator);
 
     ATH_CHECK(m_disTrkFitter.retrieve());
     ATH_MSG_DEBUG("Retrieved tool " << m_disTrkFitter);
  } else {
    m_extrapolator.disable();
    m_disTrkFitter.disable();
  }
 
  // UTT read/write handles
  ATH_CHECK( m_dEdxTrkKey.initialize(m_dodEdxTrk) );
  ATH_CHECK( m_dEdxHitKey.initialize(m_dodEdxTrk) );
  ATH_CHECK( m_disTrkCandKey.initialize(m_doDisappearingTrk) );
 
  //
  ATH_MSG_DEBUG("FTF : " << name()  );
  ATH_MSG_DEBUG("   m_tcs.m_doubletFilterRZ    : " <<  m_tcs.m_doubletFilterRZ     );
  ATH_MSG_DEBUG("   m_tcs.m_doublet_dR_Max     : " <<  m_tcs.m_doublet_dR_Max      );
  ATH_MSG_DEBUG("   m_tcs.m_doublet_dR_Max_Confirm  : " << m_tcs.m_doublet_dR_Max_Confirm   );
  ATH_MSG_DEBUG("   m_tcs.m_seedRadBinWidth    : " <<  m_tcs.m_seedRadBinWidth     );
  ATH_MSG_DEBUG("   m_tcs.m_tripletD0Max       : " <<  m_tcs.m_tripletD0Max        );
  ATH_MSG_DEBUG("   m_tcs.m_tripletD0_PPS_Max  : " <<  m_tcs.m_tripletD0_PPS_Max   );
  ATH_MSG_DEBUG("   m_tcs.m_nMaxPhiSlice       : " <<  m_tcs.m_nMaxPhiSlice        );
  ATH_MSG_DEBUG("   m_tcs.m_maxTripletBufferLength  : " << m_tcs.m_maxTripletBufferLength   );
  ATH_MSG_DEBUG("   m_tcs.m_tripletDoPSS       : " <<  m_tcs.m_tripletDoPSS        );
  ATH_MSG_DEBUG("   m_tcs.m_tripletDoPPS       : " <<  m_tcs.m_tripletDoPPS        );
  ATH_MSG_DEBUG("   m_tcs.m_tripletDoConfirm   : " <<  m_tcs.m_tripletDoConfirm    );
  ATH_MSG_DEBUG("   m_tcs.m_curv_delta         : " <<  m_tcs.m_curv_delta          );
  ATH_MSG_DEBUG("   m_tcs.m_tripletDtCut       : " <<  m_tcs.m_tripletDtCut        );
  ATH_MSG_DEBUG("   m_tcs.m_useTrigSeedML      : " <<  m_tcs.m_useTrigSeedML       );
  ATH_MSG_DEBUG("   m_trigseedML_LUT           : " <<  m_trigseedML_LUT     );
  ATH_MSG_DEBUG("   m_tcs.m_maxEC_len          : " <<  m_tcs.m_maxEC_len    );
  ATH_MSG_DEBUG("   m_vertexSeededMode         : " <<  m_vertexSeededMode   );
  ATH_MSG_DEBUG("   m_doZFinder                : " <<  m_doZFinder          );
  ATH_MSG_DEBUG("   m_doZFinderOnly            : " <<  m_doZFinderOnly      );
  ATH_MSG_DEBUG("   m_doFastZVseeding          : " <<  m_doFastZVseeding    );
  ATH_MSG_DEBUG("   m_tcs.m_zvError            : " <<  m_tcs.m_zvError      );
  ATH_MSG_DEBUG("   m_tcs.m_zvErrorEndcap      : " <<  m_tcs.m_zvErrorEndcap    );
  ATH_MSG_DEBUG("   m_storeZFinderVertices     : " <<  m_storeZFinderVertices   );
  ATH_MSG_DEBUG("   m_tripletMinPtFrac         : " <<  m_tripletMinPtFrac   );
  ATH_MSG_DEBUG("   m_pTmin                    : " <<  m_pTmin              );
  ATH_MSG_DEBUG("   m_initialD0Max             : " <<  m_initialD0Max       );
  ATH_MSG_DEBUG("   m_Z0Max                    : " <<  m_Z0Max              );
  ATH_MSG_DEBUG("   m_checkSeedRedundancy      : " <<  m_checkSeedRedundancy   );
  ATH_MSG_DEBUG("   m_minHits                  : " <<  m_minHits            );
  ATH_MSG_DEBUG("   "  );
  ATH_MSG_DEBUG("   m_useBeamSpot              : " <<  m_useBeamSpot       );
  ATH_MSG_DEBUG("   m_nfreeCut                 : " <<  m_nfreeCut          );
  ATH_MSG_DEBUG("   m_spacePointTool           : " <<  m_spacePointTool    );
  ATH_MSG_DEBUG("   m_numberingTool            : " <<  m_numberingTool     );
  ATH_MSG_DEBUG("   m_trackMaker               : " <<  m_trackMaker        );
  ATH_MSG_DEBUG("   m_trigInDetTrackFitter     : " <<  m_trigInDetTrackFitter   );
  ATH_MSG_DEBUG("   m_trigZFinder              : " <<  m_trigZFinder       );
  ATH_MSG_DEBUG("   m_trackSummaryTool         : " <<  m_trackSummaryTool  );
  ATH_MSG_DEBUG("   m_doResMonitoring          : " <<  m_doResMonitoring   );
  ATH_MSG_DEBUG("   m_doCloneRemoval           : " <<  m_doCloneRemoval    );
  ATH_MSG_DEBUG("   m_useNewLayerNumberScheme  : " <<  m_useNewLayerNumberScheme   );
  ATH_MSG_DEBUG("   m_useGPU                   : " <<  m_useGPU            );
  ATH_MSG_DEBUG("   m_LRTmode                  : " <<  m_LRTmode           );
  ATH_MSG_DEBUG("   m_dodEdxTrk                : " <<  m_dodEdxTrk         );
  ATH_MSG_DEBUG("   m_ITkMode                  : " <<  m_ITkMode         );
  ATH_MSG_DEBUG("   m_useTracklets             : " <<  m_useTracklets         );
 
  ATH_MSG_DEBUG(" Initialized successfully");
 
 
  return StatusCode::SUCCESS;
}
 
 
//-------------------------------------------------------------------------
 
StatusCode TrigFastTrackFinder::start()
{
  //getting magic numbers from the layer numbering tool
 
  m_tcs.m_maxBarrelPix    = m_numberingTool->offsetBarrelSCT();
  m_tcs.m_minEndcapPix    = m_numberingTool->offsetEndcapPixels();
  m_tcs.m_maxEndcapPix    = m_numberingTool->offsetEndcapSCT();
  m_tcs.m_maxSiliconLayer = m_numberingTool->maxSiliconLayerNum();
  m_tcs.m_layerGeometry.clear();
 
  if(m_useNewLayerNumberScheme) {
    const std::vector<TrigInDetSiLayer>* pVL = m_numberingTool->layerGeometry();
    std::copy(pVL->begin(),pVL->end(),std::back_inserter(m_tcs.m_layerGeometry));
  }
 
  m_tcs.m_tripletPtMin = m_tripletMinPtFrac*m_pTmin;
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode TrigFastTrackFinder::execute(const EventContext& ctx) const {
 
  SG::WriteHandle<TrackCollection> outputTracks(m_outputTracksKey, ctx);
 
  outputTracks = std::make_unique<TrackCollection>();
 
  const TrackCollection* inputTracks = nullptr;
  if (m_LRTmode) {
    if (!m_inputTracksKey.key().empty()){
      SG::ReadHandle<TrackCollection> inputTrackHandle(m_inputTracksKey, ctx);
      ATH_CHECK(inputTrackHandle.isValid());
      inputTracks = inputTrackHandle.ptr();
    }
  }
  InDet::ExtendedSiTrackMakerEventData_xk trackEventData(m_prdToTrackMap, ctx);
 
  //RoI preparation/update
 
  if (m_standaloneMode) {
    //the default fullscan TrigRoiDescriptor settings for beamspot width (z-range) are incorrect
    const TrigRoiDescriptor internalRoI = TrigRoiDescriptor(0, -4.5, 4.5, 0, -M_PI, M_PI, 0, -150.0, 150.0);
 
    ATH_CHECK(findTracks(trackEventData, internalRoI, inputTracks, *outputTracks, ctx));
 
  } else {
    SG::ReadHandle<TrigRoiDescriptorCollection> roiCollection(m_roiCollectionKey, ctx);
    
    ATH_CHECK(roiCollection.isValid());
    
    if ( roiCollection->size()>1 ) ATH_MSG_WARNING( "More than one Roi in the collection: " << m_roiCollectionKey << ", this is not supported - use a composite Roi: Using the first Roi ONLY" );
    
    if ( roiCollection->size()==0) {
      ATH_MSG_ERROR("No Roi found for " << m_roiCollectionKey.key() );
      return StatusCode::FAILURE;
    }
 
    TrigRoiDescriptor internalRoI = **roiCollection->begin();
 
    ATH_CHECK(findTracks(trackEventData, internalRoI, inputTracks, *outputTracks, ctx));
  }


  m_countTotalRoI++;
 
  return StatusCode::SUCCESS;
}
 
 
 
StatusCode TrigFastTrackFinder::findTracks(InDet::SiTrackMakerEventData_xk &trackEventData,
                                           const TrigRoiDescriptor& roi,
                                           const TrackCollection* inputTracks,
                                           TrackCollection& outputTracks,
                                           const EventContext& ctx) const {
  ATH_MSG_DEBUG( "Input RoI " << roi );
 
  auto mnt_roi_nTracks = Monitored::Scalar<int>("roi_nTracks", 0);
  std::vector<int> vec_seedSize;
  auto mnt_seedSize = Monitored::Collection("trk_seedSize", vec_seedSize);
  auto mnt_roi_nSPs    = Monitored::Scalar<int>("roi_nSPs",    0);
  auto mnt_roi_nSPsPIX = Monitored::Scalar<int>("roi_nSPsPIX", 0);
  auto mnt_roi_nSPsSCT = Monitored::Scalar<int>("roi_nSPsSCT", 0);
  auto monSP = Monitored::Group(m_monTool, mnt_roi_nSPsPIX, mnt_roi_nSPsSCT, mnt_seedSize);
 
  auto mnt_timer_Total                 = Monitored::Timer<std::chrono::milliseconds>("TIME_Total");
  auto mnt_timer_SpacePointConversion  = Monitored::Timer<std::chrono::milliseconds>("TIME_SpacePointConversion");
  auto mnt_timer_PatternReco           = Monitored::Timer<std::chrono::milliseconds>("TIME_PattReco");
  auto mnt_timer_TripletMaking         = Monitored::Timer<std::chrono::milliseconds>("TIME_Triplets");
  auto mnt_timer_CombTracking          = Monitored::Timer<std::chrono::milliseconds>("TIME_CmbTrack");
  auto mnt_timer_TrackFitter           = Monitored::Timer<std::chrono::milliseconds>("TIME_TrackFitter");
  auto mnt_timer_dEdxTrk               = Monitored::Timer<std::chrono::milliseconds>("TIME_dEdxTrk");
  auto mnt_timer_disTrkZVertex         = Monitored::Timer<std::chrono::milliseconds>("TIME_disTrkZVertex");
  auto mnt_timer_disTrk                = Monitored::Timer<std::chrono::milliseconds>("TIME_disappearingTrack");
  auto monTime = Monitored::Group(m_monTool, mnt_roi_nTracks, mnt_roi_nSPs, mnt_timer_Total, mnt_timer_SpacePointConversion,
                  mnt_timer_PatternReco, mnt_timer_TripletMaking, mnt_timer_CombTracking, mnt_timer_TrackFitter,
                  mnt_timer_dEdxTrk, mnt_timer_disTrkZVertex, mnt_timer_disTrk);
 
  auto mnt_roi_lastStageExecuted = Monitored::Scalar<int>("roi_lastStageExecuted", 0);
  auto monDataError              = Monitored::Group(m_monTool, mnt_roi_lastStageExecuted);
 
  mnt_timer_Total.start();
  mnt_timer_SpacePointConversion.start();
 
 
  mnt_roi_lastStageExecuted = 1;
 
  std::vector<TrigSiSpacePointBase> convertedSpacePoints;
  
  convertedSpacePoints.reserve(5000);
 
  std::map<Identifier, std::vector<long int> > siClusterMap;
 
  if (!m_ITkMode) {
 
    if (m_LRTmode) {
      // In LRT mode read the input track collection and enter the clusters on track into the cluster map so these are not used for seeding
      if (!m_inputTracksKey.key().empty()) {
    ATH_MSG_DEBUG("LRT Mode: Got input track collection with "<<inputTracks->size()<< "tracks");
    long int trackIndex=0;
    for (auto t:*inputTracks) {
      updateClusterMap(trackIndex++, t, siClusterMap);
    }
      }
      ATH_CHECK(m_spacePointTool->getSpacePoints(roi, convertedSpacePoints, mnt_roi_nSPsPIX, mnt_roi_nSPsSCT, ctx, &siClusterMap));
    }
    else {
      ATH_CHECK(m_spacePointTool->getSpacePoints(roi, convertedSpacePoints, mnt_roi_nSPsPIX, mnt_roi_nSPsSCT, ctx));
    }
    
  }
 
  mnt_timer_SpacePointConversion.stop();
  mnt_roi_nSPs    = mnt_roi_nSPsPIX + mnt_roi_nSPsSCT;
 
  if (!m_ITkMode) {
    
    if( mnt_roi_nSPs >= m_minHits ) {
      ATH_MSG_DEBUG("REGTEST / Found " << mnt_roi_nSPs << " space points.");
      ATH_MSG_DEBUG("REGTEST / Found " << mnt_roi_nSPsPIX << " Pixel space points.");
      ATH_MSG_DEBUG("REGTEST / Found " << mnt_roi_nSPsSCT << " SCT space points.");
      ATH_MSG_DEBUG("REGTEST / converted space points size = " << convertedSpacePoints.size());
      m_countRoIwithEnoughHits++;
    }
    else {
      ATH_MSG_DEBUG("No tracks found - too few hits in ROI to run " << mnt_roi_nSPs);
      ATH_CHECK( createEmptyUTTEDMs(ctx) );
      return StatusCode::SUCCESS;
    }
  }
  
  mnt_roi_lastStageExecuted = 2;
 
  std::unique_ptr<TrigRoiDescriptor> tmpRoi = std::make_unique<TrigRoiDescriptor>(roi);
 
  auto vertices = std::make_unique<TrigVertexCollection>();
  std::vector<float> vZv;
 
  if (m_doZFinder) {
    auto mnt_timer_ZFinder = Monitored::Timer<std::chrono::milliseconds>("TIME_ZFinder");
    auto monTimeZFinder    = Monitored::Group(m_monTool, mnt_timer_ZFinder);
    mnt_timer_ZFinder.start();

    tmpRoi = std::make_unique<TrigRoiDescriptor>(true);
    tmpRoi->setComposite(true);
 
    vertices = std::make_unique<TrigVertexCollection>(*m_trigZFinder->findZ( convertedSpacePoints, roi));
 
    ATH_MSG_DEBUG("vertices->size(): " << vertices->size());
 
 
    if ( m_doFastZVseeding ) {
      vZv.reserve(vertices->size());
      for (const auto vertex : *vertices) {
        ATH_MSG_DEBUG("REGTEST / ZFinder vertex: " << *vertex);
        float z      = vertex->z();
        float zMinus = z - 7.0;
        float zPlus  = z + 7.0;
        TrigRoiDescriptor* newRoi =  new TrigRoiDescriptor(roi.eta(), roi.etaMinus(), roi.etaPlus(),
            roi.phi(), roi.phiMinus(), roi.phiPlus(), z, zMinus, zPlus);
        tmpRoi->push_back(newRoi);
        vZv.push_back(z);
      }
 
      ATH_MSG_DEBUG("REGTEST / tmpRoi: " << *tmpRoi);
    }
 
    mnt_timer_ZFinder.stop();
 
    if (  m_doZFinderOnly ) {
      ATH_CHECK( createEmptyUTTEDMs(ctx) );
      return StatusCode::SUCCESS;
    }
  }
 
 
  mnt_roi_lastStageExecuted = 3;
 
  mnt_timer_PatternReco.start();
 
  mnt_timer_TripletMaking.start();
 
  std::vector<TrigInDetTriplet> triplets;
 
  std::vector<TrigInDetTracklet> new_tracklets;
 
  if(!m_useGPU) {
 
    if (m_ITkMode) {
 
      TrigInDetTrackSeedingResult seed_result = m_seedingTool->findSeeds(roi, new_tracklets, ctx);
      
      mnt_roi_nSPsPIX = seed_result.m_nPixelSPs;
      mnt_roi_nSPsSCT = seed_result.m_nStripSPs;
      mnt_roi_nSPs    = mnt_roi_nSPsPIX + mnt_roi_nSPsSCT;
 
      if( mnt_roi_nSPs >= m_minHits ) {
    ATH_MSG_DEBUG("REGTEST / Found " << mnt_roi_nSPs << " space points.");
    ATH_MSG_DEBUG("REGTEST / Found " << mnt_roi_nSPsPIX << " Pixel space points.");
    ATH_MSG_DEBUG("REGTEST / Found " << mnt_roi_nSPsSCT << " SCT space points.");
    m_countRoIwithEnoughHits++;
      }
      else {
    ATH_MSG_DEBUG("No tracks found - too few hits in ROI to run " << mnt_roi_nSPs);
    ATH_CHECK( createEmptyUTTEDMs(ctx) );
    return StatusCode::SUCCESS;
      }
      
    } else {
      TRIG_TRACK_SEED_GENERATOR seedGen(m_tcs);
      
      seedGen.loadSpacePoints(convertedSpacePoints);
      
      if (m_doZFinder && m_doFastZVseeding) {
    seedGen.createSeeds(tmpRoi.get(), vZv);
      }
      else {
    seedGen.createSeeds(tmpRoi.get());
      }
      
      seedGen.getSeeds(triplets);
    }
  }
  else {
    //GPU offloading begins ...
    ATH_CHECK(m_accelSvc->isReady()); // or remake it into afterfork incident
    makeSeedsOnGPU(m_tcs, tmpRoi.get(), convertedSpacePoints, triplets);
 
    //GPU offloading ends ...
  }
 
  unsigned int nTrackSeeds = m_useTracklets ? new_tracklets.size() : triplets.size();
  
  ATH_MSG_DEBUG("number of triplets: " << nTrackSeeds);
  
  mnt_timer_TripletMaking.stop();
  mnt_roi_lastStageExecuted = 4;
 
  mnt_timer_CombTracking.start();
 
  // 8. Combinatorial tracking
 
  std::vector<std::tuple<bool, double,Trk::Track*>> qualityTracks; //bool used for later filtering
  qualityTracks.reserve(nTrackSeeds);
 
  auto mnt_roi_nSeeds  = Monitored::Scalar<int>("roi_nSeeds",  0);
  auto monTrk_seed = Monitored::Group(m_monTool, mnt_roi_nSeeds);
 
  long int trackIndex=0;
 
  bool PIX = true;
  bool SCT = true;
 
  m_trackMaker->newTrigEvent(ctx, trackEventData, PIX, SCT);
  if( m_doDisappearingTrk ) trackEventData.combinatorialData().setFlagToReturnFailedTrack(true);
 
  std::vector<Trk::Track*> disFailTrks;
  std::vector<Trk::Track*> disCombTrks;
  std::vector<std::tuple<bool, double, Trk::Track*>> qualityDisFailTrks;
  std::vector<std::tuple<bool, double, Trk::Track*>> qualityDisCombTrks;
  int disTrk_n_disCombTrks=0;
  int disTrk_n_disCombTrks_cleaning=0;
  int disTrk_n_disFailTrks=0;
  int disTrk_n_disFailTrks_cleaning=0;
 
  
  
  for(unsigned int seedIdx=0;seedIdx!=nTrackSeeds;seedIdx++) {
 
    std::vector<const Trk::SpacePoint*> spVec;
  
    if( m_useTracklets && (!new_tracklets.empty())) { //create an n-SP seed
      spVec = new_tracklets[seedIdx].seed();
    }
    else {
    
      const TrigInDetTriplet &seed = triplets[seedIdx];
      const Trk::SpacePoint* osp1 = seed.s1().offlineSpacePoint();
      const Trk::SpacePoint* osp2 = seed.s2().offlineSpacePoint();
      const Trk::SpacePoint* osp3 = seed.s3().offlineSpacePoint();
 
      spVec = {osp1, osp2, osp3};//create a 3-SP seed
    }
 
    vec_seedSize.push_back(spVec.size());//monitoring seed length for GBTS seeding
 
    if(m_checkSeedRedundancy) {
      //check if clusters do not belong to any track
      std::vector<Identifier> clusterIds;
      extractClusterIds(spVec.at(0), clusterIds);
      extractClusterIds(spVec.at(1), clusterIds);
      extractClusterIds(spVec.at(2), clusterIds);
      if(usedByAnyTrack(clusterIds, siClusterMap)) {
    continue;
      }
    }
    ++mnt_roi_nSeeds;
 
    std::list<Trk::Track*> tracks;
    std::list<Trk::Track*> tracksFail;
    std::list<Trk::Track*> tracksAll = m_trackMaker->getTracks(ctx, trackEventData, spVec);
    auto resultCode = trackEventData.combinatorialData().resultCode();
    if( ! m_doDisappearingTrk || (m_doDisappearingTrk && resultCode==InDet::SiCombinatorialTrackFinderData_xk::ResultCode::Success) ) {
       tracks = tracksAll;
    }
    else {
       tracksFail = tracksAll;
    }
 
    if( m_doDisappearingTrk ) {
       ATH_MSG_VERBOSE("size of tracks=" << tracks.size() << ", tracksFail=" << tracksFail.size() << ": resultCode=" << resultCode);
       const TrigInDetTriplet &seed = triplets[seedIdx];
       for(std::list<Trk::Track*>::const_iterator t=tracks.begin(); t!=tracks.end(); ++t) {
      if( ! (*t) ) continue;
      m_trackSummaryTool->updateTrack(ctx, **t);
      disTrk_n_disCombTrks++;
      if( (*t)->perigeeParameters()!=0 && isCleaningPassDisTrack(seed, (*t), false) ) {
         ATH_MSG_VERBOSE("... combTrk, cleaningPass");
         disTrk_n_disCombTrks_cleaning++;
         disCombTrks.push_back((*t));
         qualityDisCombTrks.emplace_back(std::make_tuple(true, -disTrackQuality((*t)), (*t)));
      }
       }
       for(std::list<Trk::Track*>::const_iterator t=tracksFail.begin(); t!=tracksFail.end(); ++t) {
      if( ! (*t) ) continue;
      m_trackSummaryTool->updateTrack(ctx, **t);
      disTrk_n_disFailTrks++;
      if( (*t)->perigeeParameters()!=0 && isCleaningPassDisTrack(seed, (*t), true) ) {
         ATH_MSG_VERBOSE("... failTrk, cleaningPass");
         disTrk_n_disFailTrks_cleaning++;
         disFailTrks.push_back((*t));
         qualityDisFailTrks.emplace_back(std::make_tuple(true, -disTrackQuality((*t)), (*t)));
      }
      else {
         delete(*t); // delete failed trk but not disFailTrk candidate
      }
       }
    }
 
    for(std::list<Trk::Track*>::const_iterator t=tracks.begin(); t!=tracks.end(); ++t) {
      if((*t)) {
        float d0 = (*t)->perigeeParameters()==0 ? 10000.0 : (*t)->perigeeParameters()->parameters()[Trk::d0];
        if (std::abs(d0) > m_initialD0Max) {
          ATH_MSG_DEBUG("REGTEST / Reject track with d0 = " << d0 << " > " << m_initialD0Max);
          qualityTracks.push_back(std::make_tuple(false,0,(*t)));//Flag track as bad, but keep in vector for later deletion
          continue;
        }
        if(m_checkSeedRedundancy) {
          //update clusterMap
          updateClusterMap(trackIndex++, (*t), siClusterMap);
        }
        if(m_doCloneRemoval) {
          qualityTracks.push_back(std::make_tuple(true, -trackQuality((*t)), (*t)));
        }
        else {
          qualityTracks.push_back(std::make_tuple(true, 0, (*t)));
        }
      }
    }
    ATH_MSG_VERBOSE("Found "<<tracks.size()<<" tracks using triplet");
  }
 
  if( m_doDisappearingTrk ) {
     ATH_MSG_DEBUG("===> nr of disFailTrks=" << disTrk_n_disFailTrks << " -> cleaning pass=" << disTrk_n_disFailTrks_cleaning);
     ATH_MSG_DEBUG("===> nr of disCombTrks=" << disTrk_n_disCombTrks << " -> cleaning pass=" << disTrk_n_disCombTrks_cleaning);
  }
 
  m_trackMaker->endEvent(trackEventData);
 
  //clone removal
  if(m_doCloneRemoval) {
    filterSharedTracks(qualityTracks);
  }
 
  // filter shared hits
  if( m_doDisappearingTrk ) {
    filterSharedDisTracks(qualityDisFailTrks);
    filterSharedDisTracks(qualityDisCombTrks);
  }
 
  TrackCollection initialTracks;
  initialTracks.reserve(qualityTracks.size());
 
  TrackCollection extraDisCombTracks;
  // if( m_doDisappearingTrk ) extraDisCombTracks.reserve(qualityTracks.size());
 
  unsigned int idx=0;
  std::vector<unsigned int> indexDisCombTrk;
  for(const auto& q : qualityTracks) {
    bool needed_for_disCombTrk = false;
    if( m_doDisappearingTrk ) {
       Trk::Track* trk_q = std::get<2>(q);
       for(const auto& qdis : qualityDisCombTrks ) {
      if( std::get<2>(qdis) == trk_q ) {
         needed_for_disCombTrk = std::get<0>(qdis);
         break;
      }
       }
       if( needed_for_disCombTrk) ATH_MSG_VERBOSE("idx=" << idx << " ===> neded for disCombTrk");
    }
    if (std::get<0>(q)==true) {
      initialTracks.push_back(std::get<2>(q));
      if( m_doDisappearingTrk && needed_for_disCombTrk ) indexDisCombTrk.push_back(idx);
    }
    else {
      if( ! m_doDisappearingTrk ) {
     delete std::get<2>(q);
      }
      else {
     if( needed_for_disCombTrk ) {
        ATH_MSG_VERBOSE("... adding to extraDisCombTracks");
        extraDisCombTracks.push_back(std::get<2>(q));
     }
     else {
        delete std::get<2>(q);
     }
      }
    }
    idx++;
  }
  qualityTracks.clear();
 
  ATH_MSG_DEBUG("After clone removal "<<initialTracks.size()<<" tracks left");
 
 
  mnt_timer_CombTracking.stop();
  mnt_timer_PatternReco.stop();
 
  mnt_roi_lastStageExecuted = 5;
 
  mnt_timer_TrackFitter.start();
 
  if( ! m_dodEdxTrk ) {
 
    if(m_doTrackRefit) {
      m_trigInDetTrackFitter->fit(initialTracks, outputTracks, ctx, m_particleHypothesis);
    }
    else {
      outputTracks = std::move(initialTracks);
    }
  }
  else {
     TrackCollection outputTrackswTP;
     
     m_trigInDetTrackFitter->fit(initialTracks, outputTracks, outputTrackswTP, ctx, m_particleHypothesis, true); // add TP to TSoS for dEdx
     
     // large dEdx finding
     mnt_timer_dEdxTrk.start();
     for(auto t=outputTrackswTP.begin(); t!=outputTrackswTP.end();t++) { m_trackSummaryTool->updateTrack(ctx, **t); }
     ATH_CHECK( finddEdxTrk(ctx,outputTrackswTP) );
     
  }
 
  if( m_dodEdxTrk ) mnt_timer_dEdxTrk.stop(); // to include timing to destroy TrackCollection object
 
  if( outputTracks.empty() ) {
     ATH_MSG_DEBUG("REGTEST / No tracks fitted");
  }
  
 
  bool do_recoverDisCombTrk = true;
  if( m_doDisappearingTrk && (initialTracks.size()!=outputTracks.size()) ) {
     ATH_MSG_DEBUG("part of initialTracks fails in fitting. do not try to recover DisCombTracks");
     do_recoverDisCombTrk = false;
  }
 
  TrackCollection fittedExtraDisCombTracks;
  fittedExtraDisCombTracks.reserve(extraDisCombTracks.size());
  TrackCollection fittedDisCombTrks(SG::VIEW_ELEMENTS);
  if( m_doDisappearingTrk) {
     ATH_MSG_VERBOSE("nr of extraDisCombTracks=" << extraDisCombTracks.size());
     if( extraDisCombTracks.size() > 0 ) {
    ATH_MSG_VERBOSE("fitting extraDisCombTracks ...");
    m_trigInDetTrackFitter->fit(extraDisCombTracks, fittedExtraDisCombTracks, ctx, m_particleHypothesis);
    for (auto fittedTrack = fittedExtraDisCombTracks.begin(); fittedTrack!=fittedExtraDisCombTracks.end(); ++fittedTrack) {
       (*fittedTrack)->info().setPatternRecognitionInfo(Trk::TrackInfo::FastTrackFinderSeed);
       m_trackSummaryTool->updateTrack(ctx, **fittedTrack);
       fittedDisCombTrks.push_back(*fittedTrack);
    }
     }
  }
 
  //check track parameters
 
  for ( auto fittedTrack = outputTracks.begin(); fittedTrack!=outputTracks.end(); ) {
    if ((*fittedTrack)->perigeeParameters()) {
      float d0 = (*fittedTrack)->perigeeParameters()->parameters()[Trk::d0];
      float z0 = (*fittedTrack)->perigeeParameters()->parameters()[Trk::z0];
      if (std::abs(d0) > m_initialD0Max || std::abs(z0) > m_Z0Max) {
        if(m_LRTmode){
          ATH_MSG_DEBUG("REGTEST / Reject track after fit with d0 = " << d0 << " z0= "  << z0
              << " larger than limits (" << m_initialD0Max << ", " << m_Z0Max << ")");
        }else{
          ATH_MSG_WARNING("REGTEST / Reject track after fit with d0 = " << d0 << " z0= "  << z0
              << " larger than limits (" << m_initialD0Max << ", " << m_Z0Max << ")");
        }
        ATH_MSG_DEBUG(**fittedTrack);
        fittedTrack = outputTracks.erase(fittedTrack);
        continue;
      }
 
      if(m_LRTmode){
        //reject tracks which have a d0 below a cut but only when an input track collection (from ftf) is also present
        if(m_LRTD0Min>0.0){
          if(std::abs(d0) < m_LRTD0Min && !m_inputTracksKey.key().empty()){
            ATH_MSG_DEBUG("REGTEST / Reject track after fit for min d0 (" << d0 << " < " << m_LRTD0Min <<")");
            fittedTrack = outputTracks.erase(fittedTrack);
            continue;
          }
        }
 
        //calculate pt
        float trkPt = 0.0;
        if(m_LRTHardMinPt > 0.0){
          //avoid a floating poitn error
          if(std::abs((*fittedTrack)->perigeeParameters()->parameters()[Trk::qOverP]) >= 1e-9){
            trkPt = std::sin((*fittedTrack)->perigeeParameters()->parameters()[Trk::theta])/std::abs((*fittedTrack)->perigeeParameters()->parameters()[Trk::qOverP]);
 
            if(trkPt < m_LRTHardMinPt){
              ATH_MSG_DEBUG("REGTEST / Reject track after fit for min pt (" << trkPt << " < " << m_LRTHardMinPt <<")");
              fittedTrack = outputTracks.erase(fittedTrack);
              continue;
            }
          }
        }
      }
    }
    ++fittedTrack;
  }
 
  mnt_timer_TrackFitter.stop();
 
  //make track summary
 
  size_t counter(1);
  idx = 0;
 
  for ( auto fittedTrack  = outputTracks.begin();fittedTrack!=outputTracks.end();++fittedTrack) {
 
    (*fittedTrack)->info().setPatternRecognitionInfo(Trk::TrackInfo::FastTrackFinderSeed);
    ATH_MSG_VERBOSE("Updating fitted track: " << counter);
    ATH_MSG_VERBOSE(**fittedTrack);
    m_trackSummaryTool->updateTrack(ctx, **fittedTrack);
    ATH_MSG_VERBOSE("Updated track: " << counter);
    ATH_MSG_VERBOSE(**fittedTrack);
 
    if( m_doDisappearingTrk && do_recoverDisCombTrk ) {
       if( std::find(indexDisCombTrk.begin(),indexDisCombTrk.end(),idx)!=indexDisCombTrk.end() ) {
      ATH_MSG_VERBOSE("fittedTrack idx=" << idx << ": recovers also for DisCombTrack");
      fittedDisCombTrks.push_back(*fittedTrack);
       }
    }
 
    ++counter;
    idx++;
  }
  
  if( outputTracks.empty() ) {
    ATH_MSG_DEBUG("REGTEST / No tracks reconstructed");
  }
  mnt_roi_lastStageExecuted = 6;
 
  mnt_roi_nTracks = outputTracks.size();
 
  // z-vertex for UTT
  std::vector<double> disTrk_v_xVtx;
  std::vector<double> disTrk_v_yVtx;
  std::vector<double> disTrk_v_zVtx;
  if( m_doDisappearingTrk ) {
     mnt_timer_disTrkZVertex.start();
     recoVertexForDisTrack(ctx, outputTracks, disTrk_v_xVtx, disTrk_v_yVtx, disTrk_v_zVtx);
     mnt_timer_disTrkZVertex.stop();
  }
 
  // disappearing track reco
  if( m_doDisappearingTrk ) {
     mnt_timer_disTrk.start();
     ATH_CHECK( findDisTracks(ctx,outputTracks,qualityDisFailTrks,qualityDisCombTrks,fittedDisCombTrks,disTrk_v_xVtx,disTrk_v_yVtx,disTrk_v_zVtx) );
     mnt_timer_disTrk.stop();
  }
  //monitor Z-vertexing
 
  //monitor number of tracks
  ATH_MSG_DEBUG("REGTEST / Found " << outputTracks.size() << " tracks");
  if( !outputTracks.empty() )
    m_countRoIwithTracks++;

  fillMon(outputTracks, *vertices, roi, ctx);
 
  mnt_roi_lastStageExecuted = 7;
 
  mnt_timer_Total.stop();
 
  return StatusCode::SUCCESS;
}
 
double TrigFastTrackFinder::trackQuality(const Trk::Track* Tr) const {
 
  Trk::TrackStates::const_iterator
    m  = Tr->trackStateOnSurfaces()->begin(),
       me = Tr->trackStateOnSurfaces()->end  ();
 
  double quality = 0. ;
  const double W       = 17.;
 
  for(; m!=me; ++m) {
    const Trk::FitQualityOnSurface fq =  (*m)->fitQualityOnSurface();
    if(!fq) continue;
 
    double x2 = fq.chiSquared();
    double q;
    if(fq.numberDoF() == 2) q = (1.2*(W-x2*.5));
    else                     q =      (W-x2    );
    if(q < 0.) q = 0.;
    quality+=q;
  }
  return quality;
}
 
void TrigFastTrackFinder::filterSharedTracks(std::vector<std::tuple<bool, double,Trk::Track*>>& QT) const {
 
  std::set<const Trk::PrepRawData*> clusters;
 
  const Trk::PrepRawData* prd[100];
 
  std::sort(QT.begin(), QT.end(),
      [](const std::tuple<bool, double, Trk::Track*>& lhs, const std::tuple<bool, double, Trk::Track*>& rhs) {
      return std::get<1>(lhs) < std::get<1>(rhs); } );
 
  for (auto& q : QT) {
    DataVector<const Trk::MeasurementBase>::const_iterator
      m  = std::get<2>(q)->measurementsOnTrack()->begin(),
         me = std::get<2>(q)->measurementsOnTrack()->end  ();
 
    int nf = 0, nc = 0;
    for(; m!=me; ++m ) {
 
      const Trk::PrepRawData* pr = ((const Trk::RIO_OnTrack*)(*m))->prepRawData();
      if(pr) {
        ++nc;
        if(clusters.find(pr)==clusters.end()) {prd[nf++]=pr; if(nf==100) break;}
      }
    }
    if((nf >= m_nfreeCut) || (nf == nc) ) {
      for(int n=0; n!=nf; ++n) clusters.insert(prd[n]);
    }
    else  {
      std::get<0>(q) = false;
    }
  }
}
 
//---------------------------------------------------------------------------
 
StatusCode TrigFastTrackFinder::finalize()
{
 
  ATH_MSG_INFO("=========================================================");
  ATH_MSG_INFO("TrigFastTrackFinder::finalize() - TrigFastTrackFinder Statistics: ");
  ATH_MSG_INFO("RoI processed: " <<  m_countTotalRoI);
  ATH_MSG_INFO("RoI with enough SPs : " <<  m_countRoIwithEnoughHits);
  ATH_MSG_INFO("RoI with Track(s)   : " << m_countRoIwithTracks);
  ATH_MSG_INFO("=========================================================");
 
  return StatusCode::SUCCESS;
}
 
void TrigFastTrackFinder::updateClusterMap(long int trackIdx, const Trk::Track* pTrack, std::map<Identifier, std::vector<long int> >& clusterMap) const {
  //loop over clusters
 
  for(auto tMOT = pTrack->measurementsOnTrack()->begin(); tMOT != pTrack->measurementsOnTrack()->end(); ++tMOT) {
 
    const InDet::SiClusterOnTrack* siCLOT = dynamic_cast<const InDet::SiClusterOnTrack*>(*tMOT);
    if (siCLOT==nullptr) continue;
    const InDet::SiCluster* siCL = dynamic_cast<const InDet::SiCluster*>(siCLOT->prepRawData());
    if (siCL==nullptr) continue;
 
    if(m_ITkMode) {
      //skip non-pixel clusters because all seeds are PPP in ITK mode
      const InDet::PixelCluster* pixCL = dynamic_cast<const InDet::PixelCluster*>(siCL);
      if(pixCL==nullptr) continue;
    }
 
    Identifier id = siCL->identify();
    clusterMap[id].push_back(trackIdx);
  }
}
 
void TrigFastTrackFinder::extractClusterIds(const Trk::SpacePoint* pSP, std::vector<Identifier>& vIds) const {
  const InDet::SiCluster* pCL = dynamic_cast<const InDet::SiCluster*>(pSP->clusterList().first);
  if(pCL!=nullptr) vIds.push_back(pCL->identify());
  //check second cluster : SCT uv clusters only !
  pCL = dynamic_cast<const InDet::SiCluster*>(pSP->clusterList().second);
  if(pCL!=nullptr) vIds.push_back(pCL->identify());
}
 
bool TrigFastTrackFinder::usedByAnyTrack(const std::vector<Identifier>& vIds, std::map<Identifier, std::vector<long int> >& clusterMap) const {
 
  std::vector<long int> xSection;
  //initializing
  std::map<Identifier, std::vector<long int> >::iterator itm0 = clusterMap.find(*vIds.begin());
  if(itm0 == clusterMap.end()) return false;
  xSection.reserve((*itm0).second.size());
  std::copy((*itm0).second.begin(), (*itm0).second.end(), std::back_inserter(xSection));
  std::vector<Identifier>::const_iterator it = vIds.begin();++it;
  for(;it!=vIds.end();++it) {
    std::map<Identifier, std::vector<long int> >::iterator itm1 = clusterMap.find(*it);
    if(itm1 == clusterMap.end()) return false;
    std::vector<long int> tmp;
    std::set_intersection(xSection.begin(), xSection.end(), (*itm1).second.begin(),(*itm1).second.end(), std::back_inserter(tmp));
    if(tmp.empty()) return false;
    //update xSection
    xSection.clear();
    xSection.reserve(tmp.size());
    std::copy(tmp.begin(), tmp.end(), std::back_inserter(xSection));
  }
  return !xSection.empty();
}
 
void TrigFastTrackFinder::fillMon(const TrackCollection& tracks, const TrigVertexCollection& vertices,
                                  const TrigRoiDescriptor& roi, const EventContext& ctx) const {
  float shift_x = 0;
  float shift_y = 0;
  if(m_useBeamSpot) {
      SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
      FTF::getBeamSpotShift(shift_x, shift_y, **beamSpotHandle);
  }
  auto mnt_roi_eta      = Monitored::Scalar<float>("roi_eta",      0.0);
  auto mnt_roi_phi      = Monitored::Scalar<float>("roi_phi",      0.0);
  auto mnt_roi_etaWidth = Monitored::Scalar<float>("roi_etaWidth", 0.0);
  auto mnt_roi_phiWidth = Monitored::Scalar<float>("roi_phiWidth", 0.0);
  auto mnt_roi_z        = Monitored::Scalar<float>("roi_z",        0.0);
  auto mnt_roi_zWidth   = Monitored::Scalar<float>("roi_zWidth",  0.0);
  auto monRoI           = Monitored::Group(m_monTool, mnt_roi_eta, mnt_roi_phi, mnt_roi_etaWidth, mnt_roi_phiWidth, mnt_roi_z, mnt_roi_zWidth);
 
  if (roi.composite()){
    for(unsigned int i=0; i<roi.size(); i++) {
      const IRoiDescriptor *subroi = roi.at(i);
      if (subroi){
    mnt_roi_eta      = subroi->eta();
    mnt_roi_phi      = subroi->phi();
    mnt_roi_etaWidth = subroi->etaPlus() - subroi->etaMinus();
    mnt_roi_phiWidth = CxxUtils::wrapToPi(subroi->phiPlus() - subroi->phiMinus());
    mnt_roi_z        = subroi->zed();
    mnt_roi_zWidth   = subroi->zedPlus() - subroi->zedMinus();
    monRoI.fill();
      }
    }
  } 
  else {
    mnt_roi_eta      = roi.eta();
    mnt_roi_phi      = roi.phi();
    mnt_roi_etaWidth = roi.etaPlus() - roi.etaMinus();
    mnt_roi_phiWidth = CxxUtils::wrapToPi(roi.phiPlus() - roi.phiMinus());
    mnt_roi_z        = roi.zed();
    mnt_roi_zWidth   = roi.zedPlus() - roi.zedMinus();
  }
 
  std::vector<float> mnt_trk_pt;
  std::vector<float> mnt_trk_a0;
  std::vector<float> mnt_trk_z0;
  std::vector<float> mnt_trk_phi0;
  std::vector<float> mnt_trk_eta;
  std::vector<float> mnt_trk_chi2dof;
  std::vector<float> mnt_trk_nSiHits;
  std::vector<float> mnt_trk_nPIXHits;
  std::vector<float> mnt_trk_nSCTHits;
  std::vector<float> mnt_trk_a0beam;
  std::vector<float> mnt_trk_z0beam;
  std::vector<float> mnt_trk_dPhi0;
  std::vector<float> mnt_trk_dEta;
 
  auto mon_pt       = Monitored::Collection("trk_pt",       mnt_trk_pt);
  auto mon_a0       = Monitored::Collection("trk_a0",       mnt_trk_a0);
  auto mon_z0       = Monitored::Collection("trk_z0",       mnt_trk_z0);
  auto mon_phi0     = Monitored::Collection("trk_phi0",     mnt_trk_phi0);
  auto mon_eta      = Monitored::Collection("trk_eta",      mnt_trk_eta);
  auto mon_chi2dof  = Monitored::Collection("trk_chi2dof",  mnt_trk_chi2dof);
  auto mon_nSiHits  = Monitored::Collection("trk_nSiHits",  mnt_trk_nSiHits);
  auto mon_nPIXHits = Monitored::Collection("trk_nPIXHits", mnt_trk_nPIXHits);
  auto mon_nSCTHits = Monitored::Collection("trk_nSCTHits", mnt_trk_nSCTHits);
  auto mon_a0beam   = Monitored::Collection("trk_a0beam",   mnt_trk_a0beam);
  auto mon_z0beam   = Monitored::Collection("trk_z0beam",   mnt_trk_z0beam);
  auto mon_dPhi0    = Monitored::Collection("trk_dPhi0",    mnt_trk_dPhi0);
  auto mon_dEta     = Monitored::Collection("trk_dEta",     mnt_trk_dEta);
  auto monTrk       = Monitored::Group(m_monTool, mon_pt, mon_a0, mon_z0, mon_phi0, mon_eta, mon_chi2dof,
                       mon_nSiHits, mon_nPIXHits, mon_nSCTHits, mon_a0beam, mon_z0beam, mon_dPhi0, mon_dEta);
 
  std::vector<float> mnt_roi_zVertices;
  auto mon_roi_nZvertices = Monitored::Scalar<int>("roi_nZvertices", 0);
  auto mon_roi_zVertices  = Monitored::Collection("roi_zVertices", mnt_roi_zVertices);
  auto monVtx = Monitored::Group(m_monTool, mon_roi_nZvertices, mon_roi_zVertices);
  mon_roi_nZvertices = vertices.size();
  for (const auto vertex : vertices) {
    mnt_roi_zVertices.push_back(vertex->z());
  }
 
  for (auto track : tracks) {
    const Trk::TrackParameters* trackPars = track->perigeeParameters();
    if(trackPars==nullptr) {
      continue;
    }
 
    if(trackPars->covariance()==nullptr) {
      continue;
    }
 
    float a0 = trackPars->parameters()[Trk::d0];
    float z0 = trackPars->parameters()[Trk::z0];
    float phi0 = trackPars->parameters()[Trk::phi0];
    float theta = trackPars->parameters()[Trk::theta];
    float eta = -log(tan(0.5*theta));
    mnt_trk_a0.push_back(a0);
    mnt_trk_z0.push_back(z0);
    mnt_trk_phi0.push_back(phi0);
    mnt_trk_a0beam.push_back(a0+shift_x*std::sin(phi0)-shift_y*std::cos(phi0));
    mnt_trk_z0beam.push_back(z0+(shift_x*std::cos(phi0)+shift_y*std::sin(phi0))/std::tan(theta));
    mnt_trk_eta.push_back(eta);
    for(unsigned int i=0; i<roi.size(); i++) {
       mnt_trk_dPhi0.push_back(CxxUtils::wrapToPi(phi0 - (roi.at(i))->phi()));
       mnt_trk_dEta.push_back(eta - (roi.at(i))->eta());
    }
 
    float qOverP = trackPars->parameters()[Trk::qOverP];
    if (qOverP==0) {
      ATH_MSG_DEBUG("REGTEST / q/p == 0, adjusting to 1e-12");
      qOverP = 1e-12;
    }
    float pT=sin(theta)/qOverP;
 
    const Trk::FitQuality* fq = track->fitQuality();
    float chi2 = 1e8;
    if (fq) {
      ATH_MSG_VERBOSE("Fitted chi2: " << fq->chiSquared());
      ATH_MSG_VERBOSE("Fitted ndof: " << fq->numberDoF());
      if(fq->numberDoF()!=0) {
        chi2 = fq->chiSquared()/fq->numberDoF();
      }
    }
    mnt_trk_pt.push_back(pT);
    mnt_trk_chi2dof.push_back(chi2);
 
    int nPix=0, nSct=0;
 
    for(auto tSOS = track->trackStateOnSurfaces()->begin();
        tSOS!=track->trackStateOnSurfaces()->end(); ++tSOS) {
      if ((*tSOS)->type(Trk::TrackStateOnSurface::Perigee) == false) {
        const Trk::FitQualityOnSurface fq =  (*tSOS)->fitQualityOnSurface();
        if(!fq) continue;
        int nd = fq.numberDoF();
        if(nd==2) nPix++;
        if(nd==1) nSct++;
      }
    }
    mnt_trk_nPIXHits.push_back(nPix);
    mnt_trk_nSCTHits.push_back(nSct/2);
    mnt_trk_nSiHits.push_back(nPix + nSct/2);
 
    ATH_MSG_DEBUG("REGTEST / track npix/nsct/phi0/pt/eta/d0/z0/chi2: " <<
        nPix   << " / "  <<
        nSct/2 << " / "  <<
        phi0   << " / "  <<
        pT     << " / " <<
        eta    << " / " <<
        a0     << " / " <<
        z0     << " / " <<
        chi2);
    // tighter selection for unbiased residuals
    bool goodTrack = std::fabs(pT)>1000. && (nPix + nSct/2) > 3 && nSct > 0;
    if (goodTrack && m_doResMonitoring) {
      runResidualMonitoring(*track, ctx);
    }
  }
}
 
void TrigFastTrackFinder::runResidualMonitoring(const Trk::Track& track, const EventContext& ctx) const {
 
  std::vector<float> mnt_layer_IBL;
  std::vector<float> mnt_layer_PixB;
  std::vector<float> mnt_layer_PixE;
  std::vector<float> mnt_layer_SCTB;
  std::vector<float> mnt_layer_SCTE;
  std::vector<float> mnt_hit_IBLPhiResidual;
  std::vector<float> mnt_hit_IBLEtaResidual;
  std::vector<float> mnt_hit_IBLPhiPull;
  std::vector<float> mnt_hit_IBLEtaPull;
  std::vector<float> mnt_hit_PIXBarrelPhiResidual;
  std::vector<float> mnt_hit_PIXBarrelEtaResidual;
  std::vector<float> mnt_hit_PIXBarrelPhiPull;
  std::vector<float> mnt_hit_PIXBarrelEtaPull;
  std::vector<float> mnt_hit_SCTBarrelResidual;
  std::vector<float> mnt_hit_SCTBarrelPull;
  std::vector<float> mnt_hit_PIXEndcapPhiResidual;
  std::vector<float> mnt_hit_PIXEndcapEtaResidual;
  std::vector<float> mnt_hit_PIXEndcapPhiPull;
  std::vector<float> mnt_hit_PIXEndcapEtaPull;
  std::vector<float> mnt_hit_SCTEndcapResidual;
  std::vector<float> mnt_hit_SCTEndcapPull;
  std::vector<float> mnt_hit_PIXBarrelL1PhiResidual;
  std::vector<float> mnt_hit_PIXBarrelL1EtaResidual;
  std::vector<float> mnt_hit_PIXBarrelL2PhiResidual;
  std::vector<float> mnt_hit_PIXBarrelL2EtaResidual;
  std::vector<float> mnt_hit_PIXBarrelL3PhiResidual;
  std::vector<float> mnt_hit_PIXBarrelL3EtaResidual;
  std::vector<float> mnt_hit_PIXEndcapL1PhiResidual;
  std::vector<float> mnt_hit_PIXEndcapL1EtaResidual;
  std::vector<float> mnt_hit_PIXEndcapL2PhiResidual;
  std::vector<float> mnt_hit_PIXEndcapL2EtaResidual;
  std::vector<float> mnt_hit_PIXEndcapL3PhiResidual;
  std::vector<float> mnt_hit_PIXEndcapL3EtaResidual;
  std::vector<float> mnt_hit_SCTBarrelL1PhiResidual;
  std::vector<float> mnt_hit_SCTBarrelL2PhiResidual;
  std::vector<float> mnt_hit_SCTBarrelL3PhiResidual;
  std::vector<float> mnt_hit_SCTBarrelL4PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL1PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL2PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL3PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL4PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL5PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL6PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL7PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL8PhiResidual;
  std::vector<float> mnt_hit_SCTEndcapL9PhiResidual;
  auto mon_layer_IBL  = Monitored::Collection("layer_IBL", mnt_layer_IBL);
  auto mon_layer_PixB = Monitored::Collection("layer_PixB",mnt_layer_PixB);
  auto mon_layer_PixE = Monitored::Collection("layer_PixE",mnt_layer_PixE);
  auto mon_layer_SCTB = Monitored::Collection("layer_SCTB",mnt_layer_SCTB);
  auto mon_layer_SCTE = Monitored::Collection("layer_SCTE",mnt_layer_SCTE);
  auto mon_hit_IBLPhiResidual = Monitored::Collection("hit_IBLPhiResidual",mnt_hit_IBLPhiResidual);
  auto mon_hit_IBLEtaResidual = Monitored::Collection("hit_IBLEtaResidual",mnt_hit_IBLEtaResidual);
  auto mon_hit_IBLPhiPull = Monitored::Collection("hit_IBLPhiPull",mnt_hit_IBLPhiPull);
  auto mon_hit_IBLEtaPull = Monitored::Collection("hit_IBLEtaPull",mnt_hit_IBLEtaPull);
  auto mon_hit_PIXBarrelPhiResidual = Monitored::Collection("hit_PIXBarrelPhiResidual",mnt_hit_PIXBarrelPhiResidual);
  auto mon_hit_PIXBarrelEtaResidual = Monitored::Collection("hit_PIXBarrelEtaResidual",mnt_hit_PIXBarrelEtaResidual);
  auto mon_hit_PIXBarrelPhiPull = Monitored::Collection("hit_PIXBarrelPhiPull",mnt_hit_PIXBarrelPhiPull);
  auto mon_hit_PIXBarrelEtaPull = Monitored::Collection("hit_PIXBarrelEtaPull",mnt_hit_PIXBarrelEtaPull);
  auto mon_hit_SCTBarrelResidual = Monitored::Collection("hit_SCTBarrelResidual",mnt_hit_SCTBarrelResidual);
  auto mon_hit_SCTBarrelPull = Monitored::Collection("hit_SCTBarrelPull",mnt_hit_SCTBarrelPull);
  auto mon_hit_PIXEndcapPhiResidual = Monitored::Collection("hit_PIXEndcapPhiResidual",mnt_hit_PIXEndcapPhiResidual);
  auto mon_hit_PIXEndcapEtaResidual = Monitored::Collection("hit_PIXEndcapEtaResidual",mnt_hit_PIXEndcapEtaResidual);
  auto mon_hit_PIXEndcapPhiPull = Monitored::Collection("hit_PIXEndcapPhiPull",mnt_hit_PIXEndcapPhiPull);
  auto mon_hit_PIXEndcapEtaPull = Monitored::Collection("hit_PIXEndcapEtaPull",mnt_hit_PIXEndcapEtaPull);
  auto mon_hit_SCTEndcapResidual = Monitored::Collection("hit_SCTEndcapResidual",mnt_hit_SCTEndcapResidual);
  auto mon_hit_SCTEndcapPull = Monitored::Collection("hit_SCTEndcapPull",mnt_hit_SCTEndcapPull);
  auto mon_hit_PIXBarrelL1PhiResidual = Monitored::Collection("hit_PIXBarrelL1PhiResidual",mnt_hit_PIXBarrelL1PhiResidual);
  auto mon_hit_PIXBarrelL1EtaResidual = Monitored::Collection("hit_PIXBarrelL1EtaResidual",mnt_hit_PIXBarrelL1EtaResidual);
  auto mon_hit_PIXBarrelL2PhiResidual = Monitored::Collection("hit_PIXBarrelL2PhiResidual",mnt_hit_PIXBarrelL2PhiResidual);
  auto mon_hit_PIXBarrelL2EtaResidual = Monitored::Collection("hit_PIXBarrelL2EtaResidual",mnt_hit_PIXBarrelL2EtaResidual);
  auto mon_hit_PIXBarrelL3PhiResidual = Monitored::Collection("hit_PIXBarrelL3PhiResidual",mnt_hit_PIXBarrelL3PhiResidual);
  auto mon_hit_PIXBarrelL3EtaResidual = Monitored::Collection("hit_PIXBarrelL3EtaResidual",mnt_hit_PIXBarrelL3EtaResidual);
  auto mon_hit_PIXEndcapL1PhiResidual = Monitored::Collection("hit_PIXEndcapL1PhiResidual",mnt_hit_PIXEndcapL1PhiResidual);
  auto mon_hit_PIXEndcapL1EtaResidual = Monitored::Collection("hit_PIXEndcapL1EtaResidual",mnt_hit_PIXEndcapL1EtaResidual);
  auto mon_hit_PIXEndcapL2PhiResidual = Monitored::Collection("hit_PIXEndcapL2PhiResidual",mnt_hit_PIXEndcapL2PhiResidual);
  auto mon_hit_PIXEndcapL2EtaResidual = Monitored::Collection("hit_PIXEndcapL2EtaResidual",mnt_hit_PIXEndcapL2EtaResidual);
  auto mon_hit_PIXEndcapL3PhiResidual = Monitored::Collection("hit_PIXEndcapL3PhiResidual",mnt_hit_PIXEndcapL3PhiResidual);
  auto mon_hit_PIXEndcapL3EtaResidual = Monitored::Collection("hit_PIXEndcapL3EtaResidual",mnt_hit_PIXEndcapL3EtaResidual);
  auto mon_hit_SCTBarrelL1PhiResidual = Monitored::Collection("hit_SCTBarrelL1PhiResidual",mnt_hit_SCTBarrelL1PhiResidual);
  auto mon_hit_SCTBarrelL2PhiResidual = Monitored::Collection("hit_SCTBarrelL2PhiResidual",mnt_hit_SCTBarrelL2PhiResidual);
  auto mon_hit_SCTBarrelL3PhiResidual = Monitored::Collection("hit_SCTBarrelL3PhiResidual",mnt_hit_SCTBarrelL3PhiResidual);
  auto mon_hit_SCTBarrelL4PhiResidual = Monitored::Collection("hit_SCTBarrelL4PhiResidual",mnt_hit_SCTBarrelL4PhiResidual);
  auto mon_hit_SCTEndcapL1PhiResidual = Monitored::Collection("hit_SCTEndcapL1PhiResidual",mnt_hit_SCTEndcapL1PhiResidual);
  auto mon_hit_SCTEndcapL2PhiResidual = Monitored::Collection("hit_SCTEndcapL2PhiResidual",mnt_hit_SCTEndcapL2PhiResidual);
  auto mon_hit_SCTEndcapL3PhiResidual = Monitored::Collection("hit_SCTEndcapL3PhiResidual",mnt_hit_SCTEndcapL3PhiResidual);
  auto mon_hit_SCTEndcapL4PhiResidual = Monitored::Collection("hit_SCTEndcapL4PhiResidual",mnt_hit_SCTEndcapL4PhiResidual);
  auto mon_hit_SCTEndcapL5PhiResidual = Monitored::Collection("hit_SCTEndcapL5PhiResidual",mnt_hit_SCTEndcapL5PhiResidual);
  auto mon_hit_SCTEndcapL6PhiResidual = Monitored::Collection("hit_SCTEndcapL6PhiResidual",mnt_hit_SCTEndcapL6PhiResidual);
  auto mon_hit_SCTEndcapL7PhiResidual = Monitored::Collection("hit_SCTEndcapL7PhiResidual",mnt_hit_SCTEndcapL7PhiResidual);
  auto mon_hit_SCTEndcapL8PhiResidual = Monitored::Collection("hit_SCTEndcapL8PhiResidual",mnt_hit_SCTEndcapL8PhiResidual);
  auto mon_hit_SCTEndcapL9PhiResidual = Monitored::Collection("hit_SCTEndcapL9PhiResidual",mnt_hit_SCTEndcapL9PhiResidual);
 
  auto monRes = Monitored::Group(m_monTool, mon_layer_IBL, mon_layer_PixB, mon_layer_PixE, mon_layer_SCTB, mon_layer_SCTE, mon_hit_IBLPhiResidual, mon_hit_IBLEtaResidual, mon_hit_IBLPhiPull, mon_hit_IBLEtaPull, mon_hit_PIXBarrelPhiResidual, mon_hit_PIXBarrelEtaResidual, mon_hit_PIXBarrelPhiPull, mon_hit_PIXBarrelEtaPull, mon_hit_SCTBarrelResidual, mon_hit_SCTBarrelPull, mon_hit_PIXEndcapPhiResidual, mon_hit_PIXEndcapEtaResidual, mon_hit_PIXEndcapPhiPull, mon_hit_PIXEndcapEtaPull, mon_hit_SCTEndcapResidual, mon_hit_SCTEndcapPull, mon_hit_PIXBarrelL1PhiResidual, mon_hit_PIXBarrelL1EtaResidual, mon_hit_PIXBarrelL2PhiResidual, mon_hit_PIXBarrelL2EtaResidual, mon_hit_PIXBarrelL3PhiResidual, mon_hit_PIXBarrelL3EtaResidual, mon_hit_PIXEndcapL1PhiResidual, mon_hit_PIXEndcapL1EtaResidual, mon_hit_PIXEndcapL2PhiResidual, mon_hit_PIXEndcapL2EtaResidual, mon_hit_PIXEndcapL3PhiResidual, mon_hit_PIXEndcapL3EtaResidual, mon_hit_SCTBarrelL1PhiResidual, mon_hit_SCTBarrelL2PhiResidual, mon_hit_SCTBarrelL3PhiResidual, mon_hit_SCTBarrelL4PhiResidual, mon_hit_SCTEndcapL1PhiResidual, mon_hit_SCTEndcapL2PhiResidual, mon_hit_SCTEndcapL3PhiResidual, mon_hit_SCTEndcapL4PhiResidual, mon_hit_SCTEndcapL5PhiResidual, mon_hit_SCTEndcapL6PhiResidual, mon_hit_SCTEndcapL7PhiResidual, mon_hit_SCTEndcapL8PhiResidual, mon_hit_SCTEndcapL9PhiResidual);
 
  std::vector<TrigL2HitResidual> vResid;
  vResid.clear();
  StatusCode scRes = m_trigInDetTrackFitter->getUnbiasedResiduals(track,vResid, ctx);
  if(!scRes.isSuccess()) return;
  for(std::vector<TrigL2HitResidual>::iterator it=vResid.begin();it!=vResid.end();++it) {
    Identifier id = it->identify();
    int pixlayer= (m_pixelId->layer_disk(id) );
    int sctlayer= (m_sctId->layer_disk(id) );
 
    switch(it->regionId()) {
      case Region::PixBarrel :
        mnt_layer_PixB.push_back(pixlayer);
        mnt_hit_PIXBarrelPhiResidual.push_back(it->phiResidual());
        mnt_hit_PIXBarrelPhiPull.push_back(it->phiPull());
        mnt_hit_PIXBarrelEtaResidual.push_back(it->etaResidual());
        mnt_hit_PIXBarrelEtaPull.push_back(it->etaPull());
        if (pixlayer == 1) {
          mnt_hit_PIXBarrelL1PhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXBarrelL1EtaResidual.push_back(it->etaResidual());
        }
        if (pixlayer == 2) {
          mnt_hit_PIXBarrelL2PhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXBarrelL2EtaResidual.push_back(it->etaResidual());
        }
        if (pixlayer == 3) {
          mnt_hit_PIXBarrelL3PhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXBarrelL3EtaResidual.push_back(it->etaResidual());
        }
        break;
      case Region::PixEndcap :
        ATH_MSG_DEBUG("Pixel Endcap "  );
        mnt_layer_PixE.push_back(pixlayer);
        mnt_hit_PIXEndcapPhiResidual.push_back(it->phiResidual());
        mnt_hit_PIXEndcapPhiPull.push_back(it->phiPull());
        mnt_hit_PIXEndcapEtaResidual.push_back(it->etaResidual());
        mnt_hit_PIXEndcapEtaPull.push_back(it->etaPull());
        if (pixlayer == 0) {
          mnt_hit_PIXEndcapL1PhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXEndcapL1EtaResidual.push_back(it->etaResidual());
        }
        if (pixlayer == 1) {
          mnt_hit_PIXEndcapL2PhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXEndcapL2EtaResidual.push_back(it->etaResidual());
        }
        if (pixlayer == 2) {
          mnt_hit_PIXEndcapL3PhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXEndcapL3EtaResidual.push_back(it->etaResidual());
        }
        break;
      case Region::SctBarrel :
        mnt_layer_SCTB.push_back(sctlayer);
        mnt_hit_SCTBarrelResidual.push_back(it->phiResidual());
        mnt_hit_SCTBarrelPull.push_back(it->phiPull());
        if (sctlayer == 0) {
          mnt_hit_SCTBarrelL1PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 1) {
          mnt_hit_SCTBarrelL2PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 2) {
          mnt_hit_SCTBarrelL3PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 3) {
          mnt_hit_SCTBarrelL4PhiResidual.push_back(it->phiResidual());
        }
        break;
      case Region::SctEndcap :
        ATH_MSG_DEBUG("SCT Endcap"  );
        mnt_layer_SCTE.push_back(sctlayer);
        mnt_hit_SCTEndcapResidual.push_back(it->phiResidual());
        mnt_hit_SCTEndcapPull.push_back(it->phiPull());
        if (sctlayer == 0) {
          mnt_hit_SCTEndcapL1PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 1) {
          mnt_hit_SCTEndcapL2PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 2) {
          mnt_hit_SCTEndcapL3PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 3) {
          mnt_hit_SCTEndcapL4PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 4) {
          mnt_hit_SCTEndcapL5PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 5) {
          mnt_hit_SCTEndcapL6PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 6) {
          mnt_hit_SCTEndcapL7PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 7) {
          mnt_hit_SCTEndcapL8PhiResidual.push_back(it->phiResidual());
        }
        if (sctlayer == 8) {
          mnt_hit_SCTEndcapL9PhiResidual.push_back(it->phiResidual());
        }
        break;
      case Region::IBL :
        mnt_layer_IBL.push_back(pixlayer);
        if (m_tcs.m_maxSiliconLayer==32) {
          mnt_hit_IBLPhiResidual.push_back(it->phiResidual());
          mnt_hit_IBLPhiPull.push_back(it->phiPull());
          mnt_hit_IBLEtaResidual.push_back(it->etaResidual());
          mnt_hit_IBLEtaPull.push_back(it->etaPull());
        }
        else {//No IBL, fill pixel histograms instead
          ATH_MSG_DEBUG("IBL wrong region"  );
          mnt_hit_PIXBarrelPhiResidual.push_back(it->phiResidual());
          mnt_hit_PIXBarrelPhiPull.push_back(it->phiPull());
          mnt_hit_PIXBarrelEtaResidual.push_back(it->etaResidual());
          mnt_hit_PIXBarrelEtaPull.push_back(it->etaPull());
        }
        break;
      case Region::Undefined :
        ATH_MSG_DEBUG("Undefined ID region");
        break;
    }
  }
}
 
void TrigFastTrackFinder::makeSeedsOnGPU(const TrigCombinatorialSettings& tcs, const IRoiDescriptor* roi, const std
::vector<TrigSiSpacePointBase>& vsp, std::vector<TrigInDetTriplet>& output) const {
 
  output.clear();
 
  TrigAccel::DATA_EXPORT_BUFFER* dataBuffer = new TrigAccel::DATA_EXPORT_BUFFER(5000);//i.e. 5KB
 
  size_t actualSize = m_accelTool->exportSeedMakingJob(tcs, roi, vsp, *dataBuffer);
 
  ATH_MSG_DEBUG("SeedMakingJob is ready, data size for transfer = " <<actualSize);
 
  std::shared_ptr<TrigAccel::OffloadBuffer> pBuff = std::make_shared<TrigAccel::OffloadBuffer>(dataBuffer);
 
  TrigAccel::Work* pJob = m_accelSvc->createWork(TrigAccel::InDetJobControlCode::MAKE_SEEDS, pBuff);
 
  if(pJob) {
    ATH_MSG_DEBUG("Work item created for task "<<TrigAccel::InDetJobControlCode::MAKE_SEEDS);
 
    pJob->run();
 
    std::shared_ptr<TrigAccel::OffloadBuffer> pOB = pJob->getOutput();
 
    int nTriplets = m_accelTool->extractTripletsFromOutput(pOB,vsp, output);
 
    ATH_MSG_DEBUG("Found "<<nTriplets<<" triplets on GPU");
  }
 
  delete pJob;
  delete dataBuffer;
}
 
StatusCode TrigFastTrackFinder::createEmptyUTTEDMs(const EventContext& ctx) const
{
   if( m_dodEdxTrk ) {
      SG::WriteHandle<xAOD::TrigCompositeContainer> dEdxTrkHandle(m_dEdxTrkKey, ctx);
      ATH_CHECK( dEdxTrkHandle.record(std::make_unique<xAOD::TrigCompositeContainer>(), std::make_unique<xAOD::TrigCompositeAuxContainer>()) );
 
      SG::WriteHandle<xAOD::TrigCompositeContainer> dEdxHitHandle(m_dEdxHitKey, ctx);
      ATH_CHECK( dEdxHitHandle.record(std::make_unique<xAOD::TrigCompositeContainer>(), std::make_unique<xAOD::TrigCompositeAuxContainer>()) );
   }
   if( m_doDisappearingTrk ) {
      SG::WriteHandle<xAOD::TrigCompositeContainer> disTrkCandHandle(m_disTrkCandKey, ctx);
      ATH_CHECK( disTrkCandHandle.record(std::make_unique<xAOD::TrigCompositeContainer>(), std::make_unique<xAOD::TrigCompositeAuxContainer>()) );
   }
   return StatusCode::SUCCESS;
}
 
 
int TrigFastTrackFinder::getSPLayer(int layer, float eta) const
{
   float abseta = std::fabs(eta);
 
   // Pixel barrel or SCT barrel
   if( 0<=layer && layer <=7 ) {
      ATH_MSG_VERBOSE("layer=" << layer << ", eta=" << abseta);
      return layer;
   }
 
   int base = 0;
 
   //
   const float PixBR6limit = 1.29612;
   const float PixBR5limit = 1.45204;
   const float PixBR4limit = 1.64909;
   const float PixBR3limit = 1.90036;
   const float PixBR2limit = 2.2146;
 
   // Pixel Endcap #1
   base = 8;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Pix EC1, eta=" << abseta);
      if( abseta > PixBR2limit ) return 2;
      return 3;
   }
 
   // Pixel Endcap #2
   base = 9;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Pix EC2, eta=" << abseta);
      if( abseta > PixBR2limit ) return 2;
      return 3;
   }
 
   // Pixel Endcap #3
   base = 10;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Pix EC3, eta=" << abseta);
      return 3;
   }
 
   // SCT Endcap #1
   base = 11;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC1, eta=" << abseta);
      if( abseta < PixBR6limit )      return 7;
      else if( abseta < PixBR5limit ) return 6;
      return 5;
   }
 
   // SCT Endcap #2
   base = 12;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC2, eta=" << abseta);
      if( abseta < PixBR5limit )      return 7;
      else if( abseta < PixBR4limit ) return 6;
      return 4;
   }
 
   // SCT Endcap #3
   base = 13;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC3, eta=" << abseta);
      if( abseta < PixBR4limit ) return 7;
      return 5;
   }
 
   // SCT Endcap #4
   base = 14;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC4, eta=" << abseta);
      if( abseta < PixBR4limit ) return 6;
      else if( abseta < PixBR3limit ) return 6;
      return 4;
   }
 
   // SCT Endcap #5
   base = 15;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC5, eta=" << abseta);
      if( abseta < PixBR3limit ) return 7;
      return 5;
   }
 
   // SCT Endcap #6
   base = 16;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC6, eta=" << abseta);
      if( abseta < PixBR3limit ) return 6;
      return 4;
   }
 
   // SCT Endcap #7
   base = 17;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC7, eta=" << abseta);
      if( abseta < PixBR3limit ) return 7;
      return 5;
   }
 
   // SCT Endcap #8
   base = 18;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC8, eta=" << abseta);
      if( abseta < PixBR3limit ) return 7;
      return 6;
   }
 
   // SCT Endcap #9
   base = 19;
   if( layer==base || layer==(base+12) ) {
      ATH_MSG_VERBOSE("Sct EC9, eta=" << abseta);
      return 7;
   }
 
   return 0;
}
 
 
StatusCode TrigFastTrackFinder::finddEdxTrk(const EventContext& ctx, const TrackCollection& outputTracks) const
{
   SG::WriteHandle<xAOD::TrigCompositeContainer> dEdxTrkHandle(m_dEdxTrkKey, ctx);
   ATH_CHECK( dEdxTrkHandle.record(std::make_unique<xAOD::TrigCompositeContainer>(), std::make_unique<xAOD::TrigCompositeAuxContainer>()) );
   auto dEdxTrkContainer = dEdxTrkHandle.ptr();
   dEdxTrkContainer->reserve(outputTracks.size());
 
   SG::WriteHandle<xAOD::TrigCompositeContainer> dEdxHitHandle(m_dEdxHitKey, ctx);
   ATH_CHECK( dEdxHitHandle.record(std::make_unique<xAOD::TrigCompositeContainer>(), std::make_unique<xAOD::TrigCompositeAuxContainer>()) );
   auto dEdxHitContainer = dEdxHitHandle.ptr();
 
   std::vector<float> mnt_dedx;
   std::vector<int>   mnt_dedx_nusedhits;
   auto mon_dedx           = Monitored::Collection("trk_dedx",           mnt_dedx);
   auto mon_dedx_nusedhits = Monitored::Collection("trk_dedx_nusedhits", mnt_dedx_nusedhits);
   auto mondEdx            = Monitored::Group(m_monTool, mon_dedx, mon_dedx_nusedhits);
 
   int i_track=0;
 
   ATH_MSG_VERBOSE("========== in finddEdxTrk ==========");
 
 
   static constexpr float TRKCUT_PTGEV_LOOSE  =  3.0;
   static constexpr float TRKCUT_PTGEV_TIGHT  = 10.0;
   static constexpr float TRKCUT_DEDX_LOOSE   =  1.25;
   static constexpr float TRKCUT_DEDX_TIGHT   =  1.55;
 
   for (const auto track: outputTracks) {
 
      float shift_x = 0; float shift_y = 0;
      if(m_useBeamSpot) {
        SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
        FTF::getBeamSpotShift(shift_x, shift_y, **beamSpotHandle);
      }
 
      ATH_MSG_VERBOSE("+++++++ i_track: " << i_track << " +++++++");
      i_track++;
 
      trackInfo theTrackInfo;
      bool igt = FTF::isGoodTrackUTT(track, theTrackInfo, shift_x, shift_y, TRKCUT_PTGEV_LOOSE);
      if (not igt) {continue;}
 
      ATH_MSG_VERBOSE("calculate dEdx -->");
      int pixelhits=0; int n_usedhits=0;
      std::vector<float> v_pixhit_dedx; std::vector<float> v_pixhit_tot; std::vector<float> v_pixhit_trkchi2; std::vector<float> v_pixhit_trkndof;
      std::vector<int> v_pixhit_iblovfl; std::vector<int> v_pixhit_loc; std::vector<int> v_pixhit_layer;
      float dedx = dEdx(track,pixelhits,n_usedhits,v_pixhit_dedx,v_pixhit_tot,v_pixhit_trkchi2,v_pixhit_trkndof,
            v_pixhit_iblovfl,v_pixhit_loc,v_pixhit_layer);
      ATH_MSG_VERBOSE("--> dedx = " << dedx);
 
      mnt_dedx.push_back(dedx);
      mnt_dedx_nusedhits.push_back(n_usedhits);
 
      bool hpt = (theTrackInfo.ptGeV >= TRKCUT_PTGEV_TIGHT && dedx >= TRKCUT_DEDX_LOOSE);
      bool lpt = (theTrackInfo.ptGeV >= TRKCUT_PTGEV_LOOSE && dedx >= TRKCUT_DEDX_TIGHT);
      if( ! hpt && ! lpt ) continue;
 
      xAOD::TrigComposite *dEdxTrk = new xAOD::TrigComposite();
      dEdxTrkContainer->push_back(dEdxTrk);
      dEdxTrk->setDetail<int>  ("dEdxTrk_id",     i_track);
      dEdxTrk->setDetail<float>("dEdxTrk_pt",     theTrackInfo.ptGeV*Gaudi::Units::GeV);
      dEdxTrk->setDetail<float>("dEdxTrk_eta",    theTrackInfo.eta);
      dEdxTrk->setDetail<float>("dEdxTrk_phi",    theTrackInfo.phi0);
      dEdxTrk->setDetail<float>("dEdxTrk_a0beam", theTrackInfo.a0beam);
      dEdxTrk->setDetail<float>("dEdxTrk_dedx",   dedx);
      dEdxTrk->setDetail<int>  ("dEdxTrk_dedx_n_usedhits",  n_usedhits);
      dEdxTrk->setDetail<int>  ("dEdxTrk_n_hits_innermost", theTrackInfo.n_hits_innermost);
      dEdxTrk->setDetail<int>  ("dEdxTrk_n_hits_inner",     theTrackInfo.n_hits_inner);
      dEdxTrk->setDetail<int>  ("dEdxTrk_n_hits_pix",       theTrackInfo.n_hits_pix);
      dEdxTrk->setDetail<int>  ("dEdxTrk_n_hits_sct",       theTrackInfo.n_hits_sct);
 
      for(unsigned int i=0; i<v_pixhit_dedx.size(); i++) {
     xAOD::TrigComposite *dEdxHit = new xAOD::TrigComposite();
     dEdxHitContainer->push_back(dEdxHit);
     dEdxHit->setDetail<int>  ("dEdxHit_trkid",   i_track);
     dEdxHit->setDetail<float>("dEdxHit_dedx",    v_pixhit_dedx[i]);
     dEdxHit->setDetail<float>("dEdxHit_tot",     v_pixhit_tot[i]);
     dEdxHit->setDetail<float>("dEdxHit_trkchi2", v_pixhit_trkchi2[i]);
     dEdxHit->setDetail<float>("dEdxHit_trkndof", v_pixhit_trkndof[i]);
     dEdxHit->setDetail<int>  ("dEdxHit_iblovfl", v_pixhit_iblovfl[i]);
     dEdxHit->setDetail<int>  ("dEdxHit_loc",     v_pixhit_loc[i]);
     dEdxHit->setDetail<int>  ("dEdxHit_layer",   v_pixhit_layer[i]);
      }
   }
   return StatusCode::SUCCESS;
}
 
 
 
float TrigFastTrackFinder::dEdx(const Trk::Track* track, int& pixelhits, int& n_usedhits,
                std::vector<float>& v_pixhit_dedx,    std::vector<float>& v_pixhit_tot,
                std::vector<float>& v_pixhit_trkchi2, std::vector<float>& v_pixhit_trkndof,
                std::vector<int>&   v_pixhit_iblovfl, std::vector<int>&   v_pixhit_loc,     std::vector<int>& v_pixhit_layer) const
{
   const float Pixel_sensorthickness=.025;      // 250 microns Pixel Planars
   const float IBL_3D_sensorthickness=.023;     // 230 microns IBL 3D
   const float IBL_PLANAR_sensorthickness=.020; // 200 microns IBL Planars
 
   const float energyPair = 3.68e-6; // Energy in MeV to create an electron-hole pair in silicon
   const float sidensity = 2.329;    // silicon density in g cm^-3
 
   float conversion_factor=energyPair/sidensity;
 
   // Loop over pixel hits on track, and calculate dE/dx:
 
   pixelhits  = 0;
   n_usedhits = 0;
 
   v_pixhit_dedx.clear();
   v_pixhit_tot.clear();
   v_pixhit_trkchi2.clear();
   v_pixhit_trkndof.clear();
   v_pixhit_iblovfl.clear();
   v_pixhit_loc.clear();
   v_pixhit_layer.clear();
 
   const int PIXLOC_IBL_PL = 0;
   const int PIXLOC_IBL_3D = 1;
   const int PIXLOC_PIX_LY = 2;
   const int PIXLOC_PIX_EC = 3;
   const int PIXLOC_IBL_UNKNOWN = 4;
   const int PIXLOC_PIX_UNKNOWN = 5;
 
   std::multimap<float,int> dEdxMap;
   float dEdxValue = 0;
 
   // Check for track states:
   const Trk::TrackStates* recoTrackStates = track->trackStateOnSurfaces();
   if (recoTrackStates) {
 
      Trk::TrackStates::const_iterator tsosIter    = recoTrackStates->begin();
      Trk::TrackStates::const_iterator tsosIterEnd = recoTrackStates->end();
 
      int i_tsos=0;
 
      // Loop over track states on surfaces (i.e. generalized hits):
      for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
 
     ATH_MSG_VERBOSE("-------- TSoS: " << i_tsos++  << " --------");
 
     const Trk::MeasurementBase *measurement = (*tsosIter)->measurementOnTrack();
     if ( measurement == nullptr ) {
        ATH_MSG_VERBOSE("no measurement on this TSoS, skip it");
        continue;
     }
     const Trk::TrackParameters* tp = (*tsosIter)->trackParameters();
     if( tp == nullptr ) {
        ATH_MSG_VERBOSE("no trackParameters() to this TSoS, skip it");
        continue;
     }
     const InDet::PixelClusterOnTrack *pixclus = dynamic_cast<const InDet::PixelClusterOnTrack*>(measurement);
     if ( pixclus == nullptr ) {
        ATH_MSG_VERBOSE("this TSoS is not Pixel, skip it");
        continue;
     }
     const InDet::PixelCluster* prd = pixclus->prepRawData();
     if( prd == nullptr ) {
        ATH_MSG_VERBOSE("no PrepRawData(), skip it");
        continue;
     }
 
     float dotProd  = tp->momentum().dot(tp->associatedSurface().normal());
     float cosalpha = std::abs(dotProd/tp->momentum().mag());
     ATH_MSG_VERBOSE("dotProd / cosalpha = " << dotProd  << " / " << cosalpha);
     if (std::abs(cosalpha)<.16) continue;
 
     const std::vector<int>& v_tots = prd->totList();
     float charge = prd->totalCharge();
     float tot    = prd->totalToT();
     ATH_MSG_VERBOSE("charge / ToT = " << charge << " / " << tot);
     charge *= cosalpha; // path length correction
 
     double locx = pixclus->localParameters()[Trk::locX];
     double locy = pixclus->localParameters()[Trk::locY];
     int    bec        = m_pixelId->barrel_ec(pixclus->identify());
     int    layer      = m_pixelId->layer_disk(pixclus->identify());
     int    eta_module = m_pixelId->eta_module(pixclus->identify()); //check eta module to select thickness
 
     float chi2  = 0.;
     float ndof  = 0.;
     const Trk::FitQualityOnSurface fq = (*tsosIter)->fitQualityOnSurface();
     if(fq) {
        chi2  = fq.chiSquared();
        ndof  = fq.doubleNumberDoF();
     }
 
     int iblOverflow=0; // keep track if this is IBL overflow
     float thickness=0;
     int loc=-1;
 
     if ( (bec==0) and (layer==0) ){ // IBL
        const float overflowIBLToT = 16; // m_overflowIBLToT = m_offlineCalibSvc->getIBLToToverflow();
        for (int pixToT : v_tots) {
           if (pixToT >= overflowIBLToT) {
          iblOverflow = 1; // overflow pixel hit -- flag cluster
          break; //no need to check other hits of this cluster
           }
        }
        if(iblOverflow==1) ATH_MSG_VERBOSE("IBL overflow");
 
        if(((eta_module>=-10 && eta_module<=-7)||(eta_module>=6 && eta_module<=9)) && (std::abs(locy)<10. && (locx>-8.33 && locx <8.3)) ){ // IBL 3D
           thickness = IBL_3D_sensorthickness;
           loc = PIXLOC_IBL_3D;
        }
        else if((eta_module>=-6 && eta_module<=5) && (std::abs(locy)<20. &&( locx >-8.33 && locx <8.3 )) ){ // IBL planer
           thickness = IBL_PLANAR_sensorthickness;
           loc = PIXLOC_IBL_PL;
        }
        else {
           ATH_MSG_VERBOSE("unknown IBL module");
           loc = PIXLOC_IBL_UNKNOWN;
        }
     }
     else if(bec==0 && std::abs(locy)<30. &&  (( locx > -8.20 && locx < -0.60 ) || ( locx > 0.50 && locx < 8.10 ) ) ){ //PIXEL layer
        thickness = Pixel_sensorthickness;
        loc = PIXLOC_PIX_LY;
     }
     else if(std::abs(bec) == 2 && std::abs(locy)<30. && ( ( locx > -8.15 && locx < -0.55 ) || ( locx > 0.55 && locx < 8.15 ) ) ) { //PIXEL endcap
        thickness = Pixel_sensorthickness;
        loc = PIXLOC_PIX_EC;
     }
     else {
        ATH_MSG_VERBOSE("unknown Pixel module");
        loc = PIXLOC_IBL_UNKNOWN;
     }
 
     dEdxValue = 0;
     if( loc != PIXLOC_IBL_UNKNOWN && loc != PIXLOC_PIX_UNKNOWN ) {
        dEdxValue = charge*conversion_factor/thickness;
        dEdxMap.insert(std::pair<float,int>(dEdxValue, iblOverflow));
        pixelhits++;
     }
     ATH_MSG_VERBOSE("dEdx=" << dEdxValue);
     v_pixhit_dedx.push_back(dEdxValue); v_pixhit_tot.push_back(tot);
     v_pixhit_trkchi2.push_back(chi2); v_pixhit_trkndof.push_back(ndof);
     v_pixhit_iblovfl.push_back(iblOverflow); v_pixhit_loc.push_back(loc); v_pixhit_layer.push_back(layer);
      }
   }
 
   // Now calculate dEdx, multimap is already sorted in ascending order
   // float averageCharge=-1;
 
   float averagedEdx=0.;
   int IBLOverflow=0;
 
   int i_map=0;
 
   for (std::pair<float,int> itdEdx : dEdxMap) {
      ATH_MSG_VERBOSE("++++++++ i_map: " << i_map++  << " ++++++++");
      if(itdEdx.second==0){
     ATH_MSG_VERBOSE("usedhits, dEdx=" << itdEdx.first);
         averagedEdx += itdEdx.first;
         n_usedhits++;
      }
      if(itdEdx.second > 0){
     ATH_MSG_VERBOSE("IBLOverflow");
         IBLOverflow++;
      }
      // break, skipping last or the two last elements depending on total measurements
      if (((int)pixelhits >= 5) and ((int)n_usedhits >= (int)pixelhits-2)) {
     ATH_MSG_VERBOSE("break, skipping last or two last elements");
     break;
      }
 
      // break, IBL Overflow case pixelhits==3 and 4
      if((int)IBLOverflow>0 and ((int)pixelhits==3) and (int)n_usedhits==1) {
     ATH_MSG_VERBOSE("break, IBL overflow case, pixel hits=3");
     break;
      }
      if((int)IBLOverflow>0 and ((int)pixelhits==4) and (int)n_usedhits==2) {
     ATH_MSG_VERBOSE("break, IBL overflow case, pixel hits=4");
     break;
      }
 
      if (((int)pixelhits > 1) and ((int)n_usedhits >=(int)pixelhits-1)) {
     ATH_MSG_VERBOSE("break, skipping last??");
     break;
      }
 
      if((int)IBLOverflow>0 and (int)pixelhits==1){ // only IBL in overflow
     ATH_MSG_VERBOSE("break, only IBL in overflow");
         averagedEdx=itdEdx.first;
         break;
      }
   }
 
   if (n_usedhits > 0 or (n_usedhits==0 and(int)IBLOverflow>0 and (int)pixelhits==1)) {
      if(n_usedhits > 0) averagedEdx = averagedEdx / n_usedhits;
      //if(n_usedhits == 0 and (int)IBLOverflow > 0 and (int)pixelhits == 1) averagedEdx = averagedEdx;  //no-op
      ATH_MSG_DEBUG("=====> averaged dEdx = " << averagedEdx << " =====>");;
      ATH_MSG_DEBUG("   +++ Used hits: " << n_usedhits << ", IBL overflows: " << IBLOverflow );;
      ATH_MSG_DEBUG("   +++ Original number of measurements = " << pixelhits << " (map size = " << dEdxMap.size() << ") ");
      return averagedEdx;
   }
 
   // -- false return
   ATH_MSG_DEBUG("dEdx not calculated, return 0");
   return 0.;
}
 
bool TrigFastTrackFinder::isCleaningPassDisTrack(const TrigInDetTriplet& seed, Trk::Track* trk, bool isFail) const
{
   const float  PT_CUT                            = 3.0;
   //
   const double FAIL_CHI2_OV_NDOF_CUT             = 20.0;
   //
   const int    COMB_N_HITS_IBL_OR_BL_CUT         = 1;
   const int    COMB_N_HITS_PIX_BR_CUT            = 3;
   const double COMB_CHI2_OV_NDOF_PIX_BR_CUT      = 3.0;
   const int    COMB_N_GOOD_HITS_SCT_BR_CUT       = 2;
   const int    COMB_N_GOOD_DOUBLEHITS_SCT_BR_CUT = 0;
 
   // sanity check
   if( trk==nullptr )                      return false;
   if( trk->perigeeParameters()==nullptr ) return false;
   if( trk->fitQuality()==nullptr )        return false;
 
   // if allpix/barrel
   if( !seed.s1().isPixel() || !seed.s2().isPixel() || !seed.s3().isPixel() ) return false;
   float s1_z = seed.s1().z();
   float s2_z = seed.s2().z();
   float s3_z = seed.s3().z();
   const float PIXEL_BARREL_Z = 410.0;
   if( std::abs(s1_z) > PIXEL_BARREL_Z || std::abs(s2_z) > PIXEL_BARREL_Z || std::abs(s3_z) > PIXEL_BARREL_Z ) return false;
   ATH_MSG_VERBOSE("FTF::isCleaningPassDisTrack> ... barrel cut passed");
 
   // pt cut
   double theta  = trk->perigeeParameters()->parameters()[Trk::theta];
   double qOverP = std::abs(trk->perigeeParameters()->parameters()[Trk::qOverP]);
   if ( qOverP < 1e-12 ) qOverP = 1e-12;
   double pt = sin(theta)/qOverP;
   if( pt/1000.0 < PT_CUT ) return false;
   ATH_MSG_VERBOSE("FTF::isCleaningPassDisTrack> ... pT cut passed");
 
   // fail/comb dependent cuts
   if( isFail ) {
      double chi2 = trk->fitQuality()->chiSquared();
      double ndof = trk->fitQuality()->doubleNumberDoF();
      if( std::abs(ndof) < 1e-12 ) return false;
      if( chi2/ndof > FAIL_CHI2_OV_NDOF_CUT ) return false;
      ATH_MSG_VERBOSE("FTF::isCleaningPassDisTrack> ... (failTrk) Chi2 cut passed");
   }
   else {
 
      const auto & barrelInfo = getTrkBarrelLayerInfo(trk);
 
      int n_hits_iblbl = barrelInfo[0].nHits + barrelInfo[1].nHits;
      if( n_hits_iblbl < COMB_N_HITS_IBL_OR_BL_CUT ) return false;
 
      // PIX cuts
      int  n_hits_pixbr = 0;
      double chi2_pixbr = 0;
      int    ndof_pixbr = 0;
      for(unsigned int ily=0; ily<=3; ily++) {
          n_hits_pixbr += barrelInfo[ily].nHits;
          chi2_pixbr   += barrelInfo[ily].chiSq;
          ndof_pixbr   += barrelInfo[ily].nDof;
      }
      if( n_hits_pixbr < COMB_N_HITS_PIX_BR_CUT ) return false;
      if( ndof_pixbr < 1 )  return false;
      double chi2_ov_ndof_pixbr = chi2_pixbr / ndof_pixbr;
      if( chi2_ov_ndof_pixbr > COMB_CHI2_OV_NDOF_PIX_BR_CUT ) return false;
      ATH_MSG_VERBOSE("FTF::isCleaningPassDisTrack> ... (combTrk) Pix cut passed");
 
      // SCT cuts
      int n_hits_sctbr_good = 0;
      int n_doublehits_sctbr_good = 0;
      for(unsigned int ily=4; ily<=7; ily++) {
          n_hits_sctbr_good += barrelInfo[ily].nGood;
          if( barrelInfo[ily].nGood >= 2 ) n_doublehits_sctbr_good++;
      }
      if( n_hits_sctbr_good > COMB_N_GOOD_HITS_SCT_BR_CUT ) return false;
      if( n_doublehits_sctbr_good > COMB_N_GOOD_DOUBLEHITS_SCT_BR_CUT ) return false;
      ATH_MSG_VERBOSE("FTF::isCleaningPassDisTrack> ... (combTrk) SCT cut passed");
   }
 
   // cut passed
   return true;
}
 
std::array<TrigFastTrackFinder::OneLayerInfo_t, TrigFastTrackFinder::N_BARREL_LAYERS>
TrigFastTrackFinder::getTrkBarrelLayerInfo(Trk::Track* t) const
{
   static constexpr double CHI2_GOOD_CUT = 3.0;
   //nHits, chiSq, nDof, nGood in array
   std::array<TrigFastTrackFinder::OneLayerInfo_t, TrigFastTrackFinder::N_BARREL_LAYERS> result{};
   if (not t) return result;
 
   const Trk::TrackStates* recoTrackStates = t->trackStateOnSurfaces();
   if (recoTrackStates) {
      Trk::TrackStates::const_iterator tsosIter    = recoTrackStates->begin();
      Trk::TrackStates::const_iterator tsosIterEnd = recoTrackStates->end();
      for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
          const Trk::FitQualityOnSurface fq = (*tsosIter)->fitQualityOnSurface();
          double x2 = 0;
          int  ndof = 0;
          if(fq) {
            x2 = fq.chiSquared();
            ndof = fq.numberDoF();
          }
          bool chi2_good = (x2 <= CHI2_GOOD_CUT);
          const Trk::MeasurementBase *measurement = (*tsosIter)->measurementOnTrack();
          const InDet::PixelClusterOnTrack *pixclus = dynamic_cast<const InDet::PixelClusterOnTrack*>(measurement);
          if (pixclus!=nullptr) {
            int bec   = m_pixelId->barrel_ec(pixclus->identify());
            int layer = m_pixelId->layer_disk(pixclus->identify());
            if ( bec==0 ) {
              if( layer < 0 || 3 < layer ) {
                  ATH_MSG_WARNING("Pixel layer out of range, layer=" << layer);
                  continue;
              }
              result[layer].nHits++;
              result[layer].chiSq += x2;
              result[layer].nDof += ndof;
              if(chi2_good) result[layer].nGood++;
            }
          }
          const InDet::SCT_ClusterOnTrack *sctclus = dynamic_cast<const InDet::SCT_ClusterOnTrack*>(measurement);
          if (sctclus!=nullptr) {
          int bec   = m_sctId->barrel_ec(sctclus->identify());
          int layer = m_sctId->layer_disk(sctclus->identify());
          if ( bec==0 ) {
            if( layer < 0 || 3 < layer ) {
                ATH_MSG_WARNING("SCT layer out of range, layer=" << layer);
                continue;
            }
          layer += 4;
          result[layer].nHits++;
          result[layer].chiSq += x2;
          result[layer].nDof += ndof;
          if(chi2_good) result[layer].nGood++;
            }
          }
      }//end loop on TSoS
   }
   return result;
}
 
double TrigFastTrackFinder::disTrackQuality(const Trk::Track* Tr) const
{
   Trk::TrackStates::const_iterator
      m  = Tr->trackStateOnSurfaces()->begin(),
      me = Tr->trackStateOnSurfaces()->end  ();
 
   double quality_pixel = 0. ;
   double quality_sct   = 0. ;
 
   const double W = 17.;
 
   for(; m!=me; ++m) {
      const Trk::FitQualityOnSurface fq =  (*m)->fitQualityOnSurface();
      if(!fq) continue;
 
      double x2 = fq.chiSquared();
      double q;
      if(fq.numberDoF() == 2) q = (1.2*(W-x2*.5));
      else                     q =      (W-x2    );
      if(q < 0.) q = 0.;
 
      const Trk::MeasurementBase *measurement = (*m)->measurementOnTrack();
      if (measurement) {
     const InDet::PixelClusterOnTrack *pixclus = dynamic_cast<const InDet::PixelClusterOnTrack*>(measurement);
     if( pixclus !=0 ) quality_pixel += q;
     else quality_sct += q;
      }
   }
 
   //
   double quality = quality_pixel;
   quality -= quality_sct;
   if( quality < 0. ) quality = 0.;
 
   return quality;
}
 
void TrigFastTrackFinder::recoVertexForDisTrack(const EventContext& ctx, TrackCollection& tracks, std::vector<double>& v_xvtx, std::vector<double>& v_yvtx, std::vector<double>& v_zvtx) const
{
   v_xvtx.clear();
   v_yvtx.clear();
   v_zvtx.clear();
 
   // beamspot and tilt
   SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
   Amg::Vector3D vertex = beamSpotHandle->beamPos();
   double xVTX = vertex.x();
   double yVTX = vertex.y();
   double tiltXZ = beamSpotHandle->beamTilt(0);
   double tiltYZ = beamSpotHandle->beamTilt(1);
 
   // zvertex
   const double CLUSTCUT_DIST_SIGMA = 5.0;
   const double CLUSTCUT_DIST       = 2.5;
   const double CLUSTCUT_SEED_PT    = 3.0;
 
   const int VTXCUT_N_TRACKS = 3;
   const int VTXCUT_ALGO     = 1; // 0: highest pT track, 1: sumPt
 
   std::vector<std::tuple<int,double,double,Trk::Track*>> QT;
   QT.reserve(tracks.size());
 
   for (auto t=tracks.begin(); t!=tracks.end(); ++t) {
      if( ! isGoodForDisTrackVertex(*t, ctx) ) continue;
      // consider for vertex fitting (idx, sort, weight, Trk)
      double theta = (*t)->perigeeParameters()->parameters()[Trk::theta];
      double qOverP = std::abs((*t)->perigeeParameters()->parameters()[Trk::qOverP]);
      if ( qOverP < 1e-12 ) qOverP = 1e-12;
      double pt = sin(theta)/qOverP;
      pt /= 1000.0;
      QT.emplace_back(std::make_tuple(-1,pt,pt,(*t)));
   }
 
   // sort
   std::sort(QT.begin(), QT.end(),
         [](const std::tuple<int,double,double,Trk::Track*>& lhs, const std::tuple<int,double,double,Trk::Track*>& rhs) {
        return std::get<1>(lhs) > std::get<1>(rhs); } );
 
   // clustering
   std::vector<int>    cluster_ntrk;
   std::vector<double> cluster_ptsum;
   std::vector<double> cluster_z;
   std::vector<double> cluster_wsum;
   std::vector<double> cluster_zerr;
   std::vector<double> cluster_w2sum;
 
   for(unsigned int i=0; i<QT.size(); ++i) {
      Trk::Track* t = std::get<3>(QT[i]);
      double z = t->perigeeParameters()->parameters()[Trk::z0];
      double zerr = sqrt((*(t->perigeeParameters()->covariance()))(Trk::z0,Trk::z0));
      double w  = std::get<2>(QT[i]);
      double pt = std::get<1>(QT[i]);
      if(i==0) {
     cluster_ntrk.push_back(1);             cluster_ptsum.push_back(pt);
     cluster_z.push_back(w*z);              cluster_wsum.push_back(w);
     cluster_zerr.push_back(w*zerr*w*zerr); cluster_w2sum.push_back(w*w);
     continue;
      }
      //
      const int IDX_INITIAL = 100;
      double dist_min = 100.0;
      int idx_min     = IDX_INITIAL;
      for(unsigned j=0; j<cluster_z.size(); ++j) {
     double dist = std::abs(z - cluster_z[j]/cluster_wsum[j]);
     if( dist < dist_min ) {
        dist_min = dist;
        idx_min  = j;
     }
      }
      int match_idx = IDX_INITIAL;
      if( idx_min != IDX_INITIAL ) {
     double c_zerr_min = std::sqrt(cluster_zerr[idx_min]/cluster_w2sum[idx_min]);
     double err = std::sqrt(zerr*zerr+c_zerr_min*c_zerr_min);
     if( std::abs(err) < 1e-12 ) err = 1e-12;
     double dist = dist_min / err;
     if( dist < CLUSTCUT_DIST_SIGMA && dist_min < CLUSTCUT_DIST ) { match_idx = idx_min; }
      }
      //
      if( match_idx == IDX_INITIAL ) {
     if( pt > CLUSTCUT_SEED_PT && dist_min > CLUSTCUT_DIST ) {
        cluster_ntrk.push_back(1);             cluster_ptsum.push_back(pt);
        cluster_z.push_back(w*z);              cluster_wsum.push_back(w);
        cluster_zerr.push_back(w*zerr*w*zerr); cluster_w2sum.push_back(w*w);
     }
     continue;
      }
      int    new_n     = cluster_ntrk[match_idx]  + 1;
      double new_ptsum = cluster_ptsum[match_idx] + pt;
      double new_z     = cluster_z[match_idx]     + w*z;
      double new_wsum  = cluster_wsum[match_idx]  + w;
      double new_zerr  = cluster_zerr[match_idx]  + w*zerr*w*zerr;
      double new_w2sum = cluster_w2sum[match_idx] + w*w;
      cluster_ntrk[match_idx]  = new_n;
      cluster_ptsum[match_idx] = new_ptsum;
      cluster_z[match_idx]     = new_z;
      cluster_wsum[match_idx]  = new_wsum;
      cluster_zerr[match_idx]  = new_zerr;
      cluster_w2sum[match_idx] = new_w2sum;
   }
   QT.clear();
 
   // determine zvtx (pt sort)
   std::vector<std::tuple<double,double,double,int>> zVtx;
   zVtx.reserve(tracks.size());
   for(unsigned int i=0; i<cluster_z.size(); i++) {
      if( cluster_ntrk[i] < VTXCUT_N_TRACKS ) continue;
      double z = cluster_z[i] / cluster_wsum[i];
      double zerr = std::sqrt(cluster_zerr[i] / cluster_w2sum[i]);
      zVtx.push_back(std::make_tuple(cluster_ptsum[i],z,zerr,cluster_ntrk[i]));
   }
   // ptsum sort
   if( VTXCUT_ALGO == 1 ) {
      std::sort(zVtx.begin(), zVtx.end(),
        [](const std::tuple<double,double,double,int>& lhs, const std::tuple<double,double,double,int>& rhs) {
           return std::get<0>(lhs) > std::get<0>(rhs); } );
   }
   ATH_MSG_VERBOSE("disTrkZVtertex> =====  looping zVtx size: " << zVtx.size());
   for(unsigned int i=0; i<zVtx.size(); i++) {
      double z    = std::get<1>(zVtx[i]);
      double zerr = std::get<2>(zVtx[i]);
      double pt   = std::get<0>(zVtx[i]);
      int    n    = std::get<3>(zVtx[i]);
      v_zvtx.push_back(z);
      v_xvtx.push_back(xVTX - tiltXZ*z); //correction for tilt
      v_yvtx.push_back(yVTX - tiltYZ*z); //correction for tilt
      ATH_MSG_VERBOSE("disTrkZVtertex> Vertex cand i=" << i << ":  z = " << z << " +- " << zerr << ", sum n / pt = " << n << " / " << pt);
   }
 
   // monitoring
   auto mnt_disTrk_nVtx  = Monitored::Scalar<int>  ("disTrk_nVtx",  0);
   auto mnt_disTrk_xVtx  = Monitored::Scalar<float>("disTrk_xVtx",  0);
   auto mnt_disTrk_yVtx  = Monitored::Scalar<float>("disTrk_yVtx",  0);
   auto mnt_disTrk_zVtx  = Monitored::Scalar<float>("disTrk_zVtx",  0);
   auto monDisTrk = Monitored::Group(m_monTool, mnt_disTrk_nVtx, mnt_disTrk_xVtx, mnt_disTrk_yVtx, mnt_disTrk_zVtx);
   mnt_disTrk_nVtx = v_zvtx.size();
   if(v_zvtx.size()>0) {
      mnt_disTrk_xVtx = v_xvtx[0];
      mnt_disTrk_yVtx = v_yvtx[0];
      mnt_disTrk_zVtx = v_zvtx[0];
   }
}
 
bool TrigFastTrackFinder::isGoodForDisTrackVertex(Trk::Track* t, const EventContext& ctx) const
{
   const double TRKCUT_CHI2_OV_NDOF        = 3.0;
   const double TRKCUT_PT                  = 1.0;
   const double TRKCUT_D0                  = 2.0;
   const int    TRKCUT_N_HITS_INNER        = 1;
   const int    TRKCUT_N_HITS_PIX          = 3;
   const int    TRKCUT_N_HITS              = 7;
 
   // sanity check
   if ( ! t->perigeeParameters() ) return false;
   if ( ! t->fitQuality() )        return false;
   if ( t->trackSummary()==0 ) {
      m_trackSummaryTool->updateTrack(ctx, *t);
      if ( t->trackSummary()==0 )  return false;
   }
 
   // chi2
   double chi2  = t->fitQuality()->chiSquared();
   double ndof  = t->fitQuality()->doubleNumberDoF();
   if( std::abs(ndof) < 1e-2 ) return false;
   double chi2_ov_ndof = chi2/ndof;
   if( chi2_ov_ndof > TRKCUT_CHI2_OV_NDOF )  return false;
 
   // pt
   double theta = t->perigeeParameters()->parameters()[Trk::theta];
   double qOverP = std::abs(t->perigeeParameters()->parameters()[Trk::qOverP]);
   if ( qOverP < 1e-12 ) qOverP = 1e-12;
   double pt    = std::sin(theta)/qOverP;
   pt /= 1000.0;
   if( pt < TRKCUT_PT ) return false;
 
   // d0
   double d0 = t->perigeeParameters()->parameters()[Trk::d0];
   if( std::abs(d0) > TRKCUT_D0 ) return false;
 
   // nr hits
   int n_hits_innermost = t->trackSummary()->get(Trk::SummaryType::numberOfInnermostPixelLayerHits);
   int n_hits_next_to_innermost = t->trackSummary()->get(Trk::SummaryType::numberOfNextToInnermostPixelLayerHits);
   int n_hits_inner = n_hits_innermost + n_hits_next_to_innermost;
   int n_hits_pix = t->trackSummary()->get(Trk::SummaryType::numberOfPixelHits);
   int n_hits_sct = t->trackSummary()->get(Trk::SummaryType::numberOfSCTHits);
   if( n_hits_inner  < TRKCUT_N_HITS_INNER )     return false;
   if( n_hits_pix    < TRKCUT_N_HITS_PIX )       return false;
   if( (n_hits_pix+n_hits_sct) < TRKCUT_N_HITS ) return false;
 
   // ok
   return true;
}
 
void TrigFastTrackFinder::filterSharedDisTracks(std::vector<std::tuple<bool, double,Trk::Track*>>& QT) const
{
   const int N_FREE_PIX_HITS_CUT = 2;
 
   std::set<const Trk::PrepRawData*> clusters;
 
   const Trk::PrepRawData* prd[100];
 
   std::sort(QT.begin(), QT.end(),
         [](const std::tuple<bool, double, Trk::Track*>& lhs, const std::tuple<bool, double, Trk::Track*>& rhs) {
        return std::get<1>(lhs) < std::get<1>(rhs); } );
 
   for (auto& q : QT) {
      DataVector<const Trk::MeasurementBase>::const_iterator
     m  = std::get<2>(q)->measurementsOnTrack()->begin(),
         me = std::get<2>(q)->measurementsOnTrack()->end  ();
 
      int nf = 0, nc = 0;
      for(; m!=me; ++m ) {
     const Trk::PrepRawData* pr = ((const Trk::RIO_OnTrack*)(*m))->prepRawData();
     if(pr) {
        const InDet::PixelClusterOnTrack *pixclus = dynamic_cast<const InDet::PixelClusterOnTrack*>((*m));
        if (pixclus) {
           ++nc;
           if(clusters.find(pr)==clusters.end()) {prd[nf++]=pr; if(nf==100) break;}
        }
     }
      }
      if((nf >= N_FREE_PIX_HITS_CUT) || (nf == nc) ) {
     for(int n=0; n!=nf; ++n) clusters.insert(prd[n]);
      }
      else  {
     std::get<0>(q) = false;
      }
   }
}
 
StatusCode TrigFastTrackFinder::findDisTracks(const EventContext& ctx,
                          TrackCollection& tracks,
                          std::vector<std::tuple<bool, double, Trk::Track*>>& qualityDisFailTrks,
                          std::vector<std::tuple<bool, double, Trk::Track*>>& qualityDisCombTrks,
                          TrackCollection& fittedDisCombTrks,
                          const std::vector<double>& v_xvtx,
                          const std::vector<double>& v_yvtx,
                          const std::vector<double>& v_zvtx) const
{
   SG::WriteHandle<xAOD::TrigCompositeContainer> disTrkCandHandle(m_disTrkCandKey, ctx);
   ATH_CHECK( disTrkCandHandle.record(std::make_unique<xAOD::TrigCompositeContainer>(), std::make_unique<xAOD::TrigCompositeAuxContainer>()) );
   auto disTrkCandContainer = disTrkCandHandle.ptr();
 
   // monitoring
   auto mnt_disFailTrk_n      = Monitored::Scalar<int>("disFailTrk_n",      0);
   auto mnt_disFailTrk_nclone = Monitored::Scalar<int>("disFailTrk_nclone", 0);
   auto mnt_disFailTrk_ncand  = Monitored::Scalar<int>("disFailTrk_ncand",  0);
   auto mnt_disCombTrk_n      = Monitored::Scalar<int>("disCombTrk_n",      0);
   auto mnt_disCombTrk_nclone = Monitored::Scalar<int>("disCombTrk_nclone", 0);
   auto mnt_disCombTrk_ncand  = Monitored::Scalar<int>("disCombTrk_ncand",  0);
   auto monDisTrk             = Monitored::Group(m_monTool, mnt_disFailTrk_n, mnt_disFailTrk_nclone, mnt_disFailTrk_ncand, mnt_disCombTrk_n, mnt_disCombTrk_nclone, mnt_disCombTrk_ncand);
 
   // select tracks to be used for isolation calculation
   std::vector<Trk::Track*> tracksForIso;
   for (auto t=tracks.begin(); t!=tracks.end(); ++t) {
      if( isGoodForDisTrackVertex(*t,ctx) ) tracksForIso.push_back(*t);
   }
 
   //
   const std::string prefix = "disTrkCand";
 
   // disFailTrk
   TrackCollection initialDisFailTrks;
   initialDisFailTrks.reserve(qualityDisFailTrks.size());
   std::vector<int> resultCodes;
   for(const auto& q : qualityDisFailTrks) {
      if (std::get<0>(q)==true) {
     initialDisFailTrks.emplace_back(std::get<2>(q));
      }
      else {
     delete std::get<2>(q);
      }
   }
   ATH_MSG_VERBOSE("===> nr of disFailTrk=" << qualityDisFailTrks.size() << " -> clone removal=" << initialDisFailTrks.size());
 
   TrackCollection fittedDisFailTrks;
   m_trigInDetTrackFitter->fit(initialDisFailTrks, fittedDisFailTrks, ctx, m_particleHypothesis);
   int n_disFailTrkCands = recoAndFillDisTrkCand(prefix, &fittedDisFailTrks, tracksForIso, disTrkCandContainer, v_xvtx, v_yvtx, v_zvtx, true, ctx);
   ATH_MSG_VERBOSE("disFailTrk: nr of cands = " << n_disFailTrkCands);
 
   mnt_disFailTrk_n      = qualityDisFailTrks.size();
   mnt_disFailTrk_nclone = initialDisFailTrks.size();
   mnt_disFailTrk_ncand  = n_disFailTrkCands;
 
   // disCombTrk
   ATH_MSG_VERBOSE("===> nr of disCombTrk=" << qualityDisCombTrks.size() << " -> clone removal=" << fittedDisCombTrks.size());
   int n_disCombTrkCands = recoAndFillDisTrkCand(prefix, &fittedDisCombTrks, tracksForIso, disTrkCandContainer, v_xvtx, v_yvtx, v_zvtx, false, ctx);
   ATH_MSG_VERBOSE("disCombTrk: nr of cands = " << n_disCombTrkCands);
 
   mnt_disCombTrk_n      = qualityDisCombTrks.size();
   mnt_disCombTrk_nclone = fittedDisCombTrks.size();
   mnt_disCombTrk_ncand  = n_disCombTrkCands;
 
   return StatusCode::SUCCESS;
}
 
bool TrigFastTrackFinder::isPreselPassDisTrackAfterRefit(Trk::Track* trk, Trk::Track* refitTrk, double refit_d0_wrtVtx, double refit_z0_wrtVtx) const
{
   const float  PRESEL_PT_GEV            =  5.0;
   const float  PRESEL_REFIT_PT_GEV_P3S1 = 10.0;
   const double PRESEL_D0_WRTVTX         =  5.0;
   const double PRESEL_Z0_WRTVTX         = 50.0;
 
   // sanity check
   if( trk == nullptr ) return false;
 
   DisTrkCategory cat = getDisTrkCategory(trk);
   if( cat==DisTrkCategory::Pix4l_Sct1p || cat==DisTrkCategory::Pix3l_Sct1p ) { if( refitTrk == nullptr ) return false; }
 
   // refit d0
   if( std::abs(refit_d0_wrtVtx) > PRESEL_D0_WRTVTX ) return false;
 
   // refit z0
   if( std::abs(refit_z0_wrtVtx) > PRESEL_Z0_WRTVTX ) return false;
 
   // pt (either trk or refit trk should have pt beyond cut)
   std::vector<float>       v_ptGeV;
   std::vector<Trk::Track*> v_trk;
   v_trk.push_back(trk);
   if( refitTrk != nullptr ) v_trk.push_back(refitTrk);
   for(auto t : v_trk) {
      float theta  = t->perigeeParameters()->parameters()[Trk::theta];
      float qOverP = std::abs(t->perigeeParameters()->parameters()[Trk::qOverP]);
      if ( qOverP < 1e-12 ) qOverP = 1e-12;
      float ptGeV = sin(theta)/qOverP/Gaudi::Units::GeV;
      v_ptGeV.push_back(ptGeV);
   }
   bool isLowPt = true;
   for(auto pt : v_ptGeV) {
      if( pt > PRESEL_PT_GEV ) { isLowPt = false; break; }
   }
   if( isLowPt ) return false;
 
   // refit pt cut for Pix3l_Sct1p which dominates rates
   if( cat==DisTrkCategory::Pix3l_Sct1p ) {
      float refitPt = v_ptGeV[1];
      if( refitPt < PRESEL_REFIT_PT_GEV_P3S1 ) return false;
   }
 
   // cut passed
   return true;
}
 
bool TrigFastTrackFinder::isPreselPassDisTrackBeforeRefit(Trk::Track* trk, double d0_wrtVtx, double z0_wrtVtx) const
{
   const double PRESEL_D0_WRTVTX =  5.0;
   const double PRESEL_Z0_WRTVTX = 50.0;
   const int    PRESEL_N_GOOD_BR_LAYERS_PIX    = 3;
   const double PRESEL_CHI2_OV_NDOF_PIX_BR_CUT = 5.0;
 
   // sanity check
   if( trk == nullptr ) return false;
   if( trk->perigeeParameters() == nullptr ) return false;
 
   // barrel hits
   const auto & barrelLayerInfo = getTrkBarrelLayerInfo(trk);
 
   // PIX cuts
   double chi2_pixbr = 0.0;
   int    ndof_pixbr = 0;
   int    n_good_brlayers_pix = 0;
   for(unsigned int ily=0; ily<=3; ily++) {
      if( barrelLayerInfo[ily].nGood >= 1 ) n_good_brlayers_pix++;
      chi2_pixbr += barrelLayerInfo[ily].chiSq;
      ndof_pixbr += barrelLayerInfo[ily].nDof;
   }
   if( n_good_brlayers_pix < PRESEL_N_GOOD_BR_LAYERS_PIX ) return false;
 
   if( ndof_pixbr < 1 )  return false;
   double chi2_ov_ndof_pixbr = chi2_pixbr / ndof_pixbr;
   if( chi2_ov_ndof_pixbr > PRESEL_CHI2_OV_NDOF_PIX_BR_CUT ) return false;
 
   // d0
   if( std::abs(d0_wrtVtx) > PRESEL_D0_WRTVTX ) return false;
 
   // z0
   if( std::abs(z0_wrtVtx) > PRESEL_Z0_WRTVTX ) return false;
 
   // cut passed
   return true;
}
 
std::unique_ptr<const Trk::TrackParameters> TrigFastTrackFinder::extrapolateDisTrackToBS(
   Trk::Track* t, const std::vector<double>& v_xvtx, const std::vector<double>& v_yvtx,  const std::vector<double>& v_zvtx, const EventContext& ctx) const
{
   float vtx_x  = 0;
   float vtx_y  = 0;
   float vtx_z  = 9999;
   float trk_z0 = t->perigeeParameters()->parameters()[Trk::z0];
   float z0_min = 9999;
   for(unsigned int i_vtx=0; i_vtx<v_zvtx.size(); i_vtx++) {
      float z = v_zvtx[i_vtx];
      if( std::abs(trk_z0-z) < z0_min ) {
     z0_min = std::abs(trk_z0-z);
     vtx_z  = z;
     vtx_x  = v_xvtx[i_vtx];
     vtx_y  = v_yvtx[i_vtx];
      }
   }
 
   Amg::Vector3D gp(vtx_x, vtx_y, vtx_z);
   Trk::PerigeeSurface persf(gp);
   std::unique_ptr<const Trk::TrackParameters> tmp =
      m_extrapolator->extrapolateDirectly(ctx, (*(t->perigeeParameters())), persf);
   if (tmp && tmp->associatedSurface().type() == Trk::SurfaceType::Perigee) {
      return tmp;
   }
   return nullptr;
}
 
TrigFastTrackFinder::DisTrkCategory TrigFastTrackFinder::getDisTrkCategory(Trk::Track* trk) const
{
   const auto & result = getTrkBarrelLayerInfo(trk);
 
   int n_good_brlayers_pix = 0;
   int n_hits_sct = 0;
   for(unsigned int ily=0; ily<8; ily++) {
      if( ily<=3 && result[ily].nGood >= 1 ) n_good_brlayers_pix++;
      if( 4<=ily ) {
          n_hits_sct += result[ily].nHits;
      }
   }
   if( trk->trackSummary()!=0 ) { n_hits_sct = trk->trackSummary()->get(Trk::SummaryType::numberOfSCTHits); }
 
   // category
   DisTrkCategory cat = DisTrkCategory::Other;
   if( n_good_brlayers_pix == 4 ) {
      if( n_hits_sct==0 ) { cat=DisTrkCategory::Pix4l_Sct0; }
      else                { cat=DisTrkCategory::Pix4l_Sct1p; }
   }
   else if( n_good_brlayers_pix == 3 ) {
      if( n_hits_sct==0 ) { cat=DisTrkCategory::Pix3l_Sct0; }
      else                { cat=DisTrkCategory::Pix3l_Sct1p; }
   }
   return cat;
}
 
void TrigFastTrackFinder::fillDisTrkCand(xAOD::TrigComposite* comp, const std::string& prefix, Trk::Track* trk, const std::unique_ptr<const Trk::TrackParameters>& vertexPerigee) const
{
   std::vector<Trk::Track*> vtmp;
   fillDisTrkCand(comp,prefix,trk,vertexPerigee,false,vtmp);
}
 
void TrigFastTrackFinder::fillDisTrkCand(xAOD::TrigComposite* comp, const std::string& prefix, Trk::Track* trk, const std::unique_ptr<const Trk::TrackParameters>& vertexPerigee,
                     bool fillIso, std::vector<Trk::Track*>& tracksForIso) const
{
   // category
   int category = (trk != nullptr) ? (int)getDisTrkCategory(trk) : -1;
   if( prefix.find("refit") == std::string::npos ) comp->setDetail<int16_t>(prefix+"_category",(int16_t)category);
 
   // track
   float theta=0; float eta=0; float pt=0; float d0=0; float z0=0; float phi=0; float chi2=0; float ndof=0;
   int n_hits_innermost=-1; int n_hits_next_to_innermost=-1; int n_hits_inner=-1; int n_hits_pix=-1; int n_hits_sct=-1;
   if( trk != nullptr ) {
      theta = trk->perigeeParameters()->parameters()[Trk::theta];
      eta   = -std::log(std::tan(0.5*theta));
      float qOverP = std::abs(trk->perigeeParameters()->parameters()[Trk::qOverP]);
      if ( qOverP < 1e-12 ) qOverP = 1e-12;
      pt    = sin(theta)/qOverP;
      d0    = trk->perigeeParameters()->parameters()[Trk::d0];
      z0    = trk->perigeeParameters()->parameters()[Trk::z0];
      phi   = trk->perigeeParameters()->parameters()[Trk::phi];
      chi2  = trk->fitQuality()->chiSquared();
      ndof  = trk->fitQuality()->doubleNumberDoF();
      if( trk->trackSummary()!=0 ) {
     n_hits_pix = trk->trackSummary()->get(Trk::SummaryType::numberOfPixelHits);
     n_hits_sct = trk->trackSummary()->get(Trk::SummaryType::numberOfSCTHits);
     n_hits_innermost = trk->trackSummary()->get(Trk::SummaryType::numberOfInnermostPixelLayerHits);
     n_hits_next_to_innermost = trk->trackSummary()->get(Trk::SummaryType::numberOfNextToInnermostPixelLayerHits);
     n_hits_inner = n_hits_innermost + n_hits_next_to_innermost;
      }
   }
   comp->setDetail<float>(prefix+"_pt",   pt);
   comp->setDetail<float>(prefix+"_eta",  eta);
   comp->setDetail<float>(prefix+"_phi",  phi);
   comp->setDetail<float>(prefix+"_d0",   d0);
   comp->setDetail<float>(prefix+"_z0",   z0);
   comp->setDetail<float>(prefix+"_chi2", chi2);
   comp->setDetail<float>(prefix+"_ndof", ndof);
   comp->setDetail<int16_t>(prefix+"_n_hits_innermost", (int16_t)n_hits_innermost);
   comp->setDetail<int16_t>(prefix+"_n_hits_inner",     (int16_t)n_hits_inner);
   comp->setDetail<int16_t>(prefix+"_n_hits_pix",       (int16_t)n_hits_pix);
   comp->setDetail<int16_t>(prefix+"_n_hits_sct",       (int16_t)n_hits_sct);
 
   // extrapolate
   float theta_wrtVtx=0; float eta_wrtVtx=0; float pt_wrtVtx=0; float d0_wrtVtx=0; float z0_wrtVtx=0; float phi_wrtVtx=0;
   if( vertexPerigee != nullptr ) {
      theta_wrtVtx = vertexPerigee->parameters()[Trk::theta];
      eta_wrtVtx   = -std::log(std::tan(0.5*theta_wrtVtx));
      float qOverP_wrtVtx = std::abs(vertexPerigee->parameters()[Trk::qOverP]);
      if ( qOverP_wrtVtx < 1e-12 ) qOverP_wrtVtx = 1e-12;
      pt_wrtVtx    = std::sin(theta_wrtVtx)/qOverP_wrtVtx;
      d0_wrtVtx    = vertexPerigee->parameters()[Trk::d0];
      z0_wrtVtx    = vertexPerigee->parameters()[Trk::z0];
      phi_wrtVtx   = vertexPerigee->parameters()[Trk::phi];
   }
   comp->setDetail<float>(prefix+"_pt_wrtVtx",  pt_wrtVtx);
   comp->setDetail<float>(prefix+"_eta_wrtVtx", eta_wrtVtx);
   comp->setDetail<float>(prefix+"_phi_wrtVtx", phi_wrtVtx);
   comp->setDetail<float>(prefix+"_d0_wrtVtx",  d0_wrtVtx);
   comp->setDetail<float>(prefix+"_z0_wrtVtx",  z0_wrtVtx);
 
   // barrel hits
   std::array<OneLayerInfo_t, N_BARREL_LAYERS> barrelInfo{};
   barrelInfo = getTrkBarrelLayerInfo(trk);
   comp->setDetail<float>(prefix+"_chi2sum_br_ibl",     barrelInfo[0].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_pix1",    barrelInfo[1].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_pix2",    barrelInfo[2].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_pix3",    barrelInfo[3].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_sct1",    barrelInfo[4].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_sct2",    barrelInfo[5].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_sct3",    barrelInfo[6].chiSq);
   comp->setDetail<float>(prefix+"_chi2sum_br_sct4",    barrelInfo[7].chiSq);
   comp->setDetail<float>(prefix+"_ndofsum_br_ibl",     barrelInfo[0].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_pix1",    barrelInfo[1].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_pix2",    barrelInfo[2].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_pix3",    barrelInfo[3].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_sct1",    barrelInfo[4].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_sct2",    barrelInfo[5].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_sct3",    barrelInfo[6].nDof);
   comp->setDetail<float>(prefix+"_ndofsum_br_sct4",    barrelInfo[7].nDof);
 
   // isolation
   if( fillIso ) {
      const float ISOL_CALC_Z0_DIFF_CUT          = 2.5;
      const float ISOL_CALC_DR_CUT_TO_AVOID_ZERO = 0.015;
      float iso1_dr01=0; float iso1_dr02=0;
      float iso2_dr01=0; float iso2_dr02=0;
      float iso3_dr01=0; float iso3_dr02=0;
      for(auto t=tracksForIso.begin(); t!=tracksForIso.end(); t++) {
     float z0_t   = (*t)->perigeeParameters()->parameters()[Trk::z0];
     if( std::abs(z0_t - z0) <= ISOL_CALC_Z0_DIFF_CUT ) {
        float theta_t = (*t)->perigeeParameters()->parameters()[Trk::theta];
        float qOverP_t= std::abs((*t)->perigeeParameters()->parameters()[Trk::qOverP]);
        if ( qOverP_t < 1e-12 ) qOverP_t = 1e-12;
        float pt_t    = std::sin(theta_t)/qOverP_t;
        float phi_t   = (*t)->perigeeParameters()->parameters()[Trk::phi];
        float eta_t   = -std::log(std::tan(theta_t/2.0));
        float deta    = eta_t - eta;
        float dphi    = std::abs(phi_t - phi);
        if( dphi > CLHEP::pi ) dphi = CLHEP::pi*2 - dphi;
        float dr   = std::sqrt(deta*deta + dphi*dphi);
        if( dr > ISOL_CALC_DR_CUT_TO_AVOID_ZERO && dr<0.1 && pt_t > 1.0*Gaudi::Units::GeV ) iso1_dr01 += pt_t;
        if( dr > ISOL_CALC_DR_CUT_TO_AVOID_ZERO && dr<0.2 && pt_t > 1.0*Gaudi::Units::GeV ) iso1_dr02 += pt_t;
        //
        if( dr > ISOL_CALC_DR_CUT_TO_AVOID_ZERO && dr<0.1 && pt_t > 2.0*Gaudi::Units::GeV ) iso2_dr01 += pt_t;
        if( dr > ISOL_CALC_DR_CUT_TO_AVOID_ZERO && dr<0.2 && pt_t > 2.0*Gaudi::Units::GeV ) iso2_dr02 += pt_t;
        //
        if( dr > ISOL_CALC_DR_CUT_TO_AVOID_ZERO && dr<0.1 && pt_t > 3.0*Gaudi::Units::GeV ) iso3_dr01 += pt_t;
        if( dr > ISOL_CALC_DR_CUT_TO_AVOID_ZERO && dr<0.2 && pt_t > 3.0*Gaudi::Units::GeV ) iso3_dr02 += pt_t;
     }
      }
      comp->setDetail<float>(prefix+"_iso1_dr01", iso1_dr01);
      comp->setDetail<float>(prefix+"_iso1_dr02", iso1_dr02);
      comp->setDetail<float>(prefix+"_iso2_dr01", iso2_dr01);
      comp->setDetail<float>(prefix+"_iso2_dr02", iso2_dr02);
      comp->setDetail<float>(prefix+"_iso3_dr01", iso3_dr01);
      comp->setDetail<float>(prefix+"_iso3_dr02", iso3_dr02);
   }
}
 
int TrigFastTrackFinder::recoAndFillDisTrkCand(const std::string& base_prefix,
                           TrackCollection* tracks, std::vector<Trk::Track*>& tracksForIso,
                           xAOD::TrigCompositeContainer* trigCompositeContainer,
                           const std::vector<double>& v_xvtx,
                           const std::vector<double>& v_yvtx,
                           const std::vector<double>& v_zvtx,
                           bool isFail,
                 const EventContext& ctx) const
{
   std::string prefix;
 
   int n_stored_tracks = 0;
 
   for (auto trk = tracks->begin(); trk!=tracks->end(); ++trk) {
 
      Trk::Track* ptrk = *trk;
 
      if( ptrk == nullptr ) continue;
      if( ptrk->perigeeParameters()==nullptr ) continue;
 
      // extrapolate to vertex
      std::unique_ptr<const Trk::TrackParameters> vertexPerigee = extrapolateDisTrackToBS(ptrk,v_xvtx,v_yvtx,v_zvtx, ctx);
      double d0 = ptrk->perigeeParameters()->parameters()[Trk::d0];
      double z0 = ptrk->perigeeParameters()->parameters()[Trk::z0];
      double d0_wrtVtx = 0;
      double z0_wrtVtx = 0;
      if( vertexPerigee != nullptr ) {
     d0_wrtVtx = vertexPerigee->parameters()[Trk::d0];
     z0_wrtVtx = vertexPerigee->parameters()[Trk::z0];
     ATH_MSG_VERBOSE("d0 : " << d0 << " -> extrapolate -> " << d0_wrtVtx);
     ATH_MSG_VERBOSE("z0 : " << z0 << " -> extrapolate -> " << z0_wrtVtx);
      }
 
      m_trackSummaryTool->updateTrack(ctx, *ptrk);
 
      // pre-selection before refit
      if( ! isPreselPassDisTrackBeforeRefit(ptrk,d0_wrtVtx,z0_wrtVtx) ) continue;
 
      // refit
      std::unique_ptr<Trk::Track> refit_trk = disTrk_refit(ptrk, ctx);
      if( refit_trk != nullptr ) m_trackSummaryTool->updateTrack(ctx, *refit_trk);
 
      // extrapolate refitted track to vertex
      double refit_d0 = 0;
      double refit_z0 = 0;
      double refit_d0_wrtVtx = 0;
      double refit_z0_wrtVtx = 0;
      std::unique_ptr<const Trk::TrackParameters> refitVertexPerigee = nullptr;
      if( refit_trk != nullptr ) {
        refitVertexPerigee = extrapolateDisTrackToBS(refit_trk.get(),v_xvtx,v_yvtx,v_zvtx,ctx);
        if( refitVertexPerigee == nullptr ) {
          ATH_MSG_VERBOSE("extrapote to BS fails for refit track");
        }
        else {
          refit_d0 = refit_trk.get()->perigeeParameters()->parameters()[Trk::d0];
          refit_z0 = refit_trk.get()->perigeeParameters()->parameters()[Trk::z0];
          refit_d0_wrtVtx = refitVertexPerigee->parameters()[Trk::d0];
          refit_z0_wrtVtx = refitVertexPerigee->parameters()[Trk::z0];
          ATH_MSG_VERBOSE("refit trk d0 : " << refit_d0 << " -> extrapolate -> " << refit_d0_wrtVtx);
          ATH_MSG_VERBOSE("refit trk z0 : " << refit_z0 << " -> extrapolate -> " << refit_z0_wrtVtx);
        }
      }
 
      // pre-selection after refit
      if( ! isPreselPassDisTrackAfterRefit(ptrk,refit_trk.get(),refit_d0_wrtVtx,refit_z0_wrtVtx) ) continue;
 
      // store it!
      n_stored_tracks++;
 
      xAOD::TrigComposite *comp = new xAOD::TrigComposite();
      comp->makePrivateStore();
      trigCompositeContainer->push_back(comp);
 
      //
      int is_fail = isFail ? 1 : 0;
      comp->setDetail<int16_t>(base_prefix+"_is_fail",(int16_t)is_fail);
 
      // store trk info
      prefix = base_prefix;
      fillDisTrkCand(comp,prefix,ptrk,vertexPerigee,true,tracksForIso);
 
      // store refit trk info
      prefix = base_prefix + "_refit";
      if( refit_trk != nullptr ) {
     fillDisTrkCand(comp,prefix,refit_trk.get(),refitVertexPerigee);
      }
      else {
     fillDisTrkCand(comp,prefix,nullptr,refitVertexPerigee);
      }
   }
 
   //
   ATH_MSG_VERBOSE("========> filling trigcomposite for " << prefix << " end");
   ATH_MSG_VERBOSE("nr of " << prefix << " tracks / stored = " << tracks->size() << " / " << n_stored_tracks);
 
   //
   return n_stored_tracks;
}
 
std::unique_ptr<Trk::Track> TrigFastTrackFinder::disTrk_refit(Trk::Track* t, const EventContext& ctx) const
{
   std::unique_ptr<Trk::Track> newtrack = nullptr;
 
   if( t == nullptr ) return newtrack;
   if( t->trackSummary() == nullptr ) m_trackSummaryTool->updateTrack(ctx, *t);
 
   ATH_MSG_VERBOSE("refitting - input track:");
   print_disTrk(t);
 
   const Trk::Perigee* origPerigee = t->perigeeParameters();
   if( origPerigee == nullptr ) return newtrack;
   ATH_MSG_VERBOSE("... origPerigee is there");
 
   // remove SCT hits
   std::vector<const Trk::MeasurementBase*> vec;
   int n_measurements       = 0;
   int n_measurements_refit = 0;
   const Trk::TrackStates* recoTrackStates = t->trackStateOnSurfaces();
   if (recoTrackStates) {
      Trk::TrackStates::const_iterator tsosIter    = recoTrackStates->begin();
      Trk::TrackStates::const_iterator tsosIterEnd = recoTrackStates->end();
      for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
     const Trk::MeasurementBase *measurement = (*tsosIter)->measurementOnTrack();
     if (measurement) {
        n_measurements++;
        const InDet::PixelClusterOnTrack *pixclus = dynamic_cast<const InDet::PixelClusterOnTrack*>(measurement);
        const InDet::SCT_ClusterOnTrack  *sctclus = dynamic_cast<const InDet::SCT_ClusterOnTrack*>(measurement);
        bool to_add = true;
        if ( !pixclus && sctclus ) to_add = false;
        if( to_add ) {
           vec.push_back(measurement);
           n_measurements_refit++;
        }
     }
      }
   }
   ATH_MSG_VERBOSE("... Nr of measurments / refit = " << n_measurements << " / " << n_measurements_refit);
 
   // perform refit
   newtrack = m_disTrkFitter->fit(ctx,vec, *origPerigee, false, m_particleHypothesis); // false to run outlier switch
   ATH_MSG_VERBOSE("... ---> refit track:");
   if( newtrack!=0 && newtrack.get() ) {
      print_disTrk(dynamic_cast<const Trk::Track*>(newtrack.get()));
   }
   else {
      ATH_MSG_VERBOSE("... refit failed");
   }
 
   //
   return newtrack;
}
 
void TrigFastTrackFinder::print_disTrk(const Trk::Track* t) const
{
   float chi2=0; float ndof=0; float d0=0; float z0=0; float phi=0; float theta=0; float pt=0;
   if( t!=nullptr ) {
      chi2  = t->fitQuality()->chiSquared();
      ndof  = t->fitQuality()->doubleNumberDoF();
      d0    = t->perigeeParameters()->parameters()[Trk::d0];
      z0    = t->perigeeParameters()->parameters()[Trk::z0];
      phi   = t->perigeeParameters()->parameters()[Trk::phi];
      theta = t->perigeeParameters()->parameters()[Trk::theta];
      float qOverP = std::abs(t->perigeeParameters()->parameters()[Trk::qOverP]);
      if ( qOverP < 1e-12 ) qOverP = 1e-12;
      pt = sin(theta)/qOverP;
      pt /= Gaudi::Units::GeV;
   }
   ATH_MSG_DEBUG("... pt / theta / phi / d0 / z0 = " << pt << " / " << theta << " / " << phi << " / " << d0 << " / " << z0);
   ATH_MSG_DEBUG("... chi2 / ndof = " << chi2 << " / " << ndof);
}