SiSPSeededTrackFinder.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
 
#include "SiSPSeededTrackFinder/SiSPSeededTrackFinder.h"
 
#include "RoiDescriptor/RoiDescriptor.h"
#include "SiSPSeededTrackFinderData/SiSpacePointsSeedMakerEventData.h"
#include "SiSPSeededTrackFinderData/SiTrackMakerEventData_xk.h"
#include "TrkPatternParameters/PatternTrackParameters.h"
#include "TrkRIO_OnTrack/RIO_OnTrack.h"
#include "TrkTrackSummary/TrackSummary.h"
 
#include "InDetIdentifier/PixelID.h"
#include "InDetIdentifier/SCT_ID.h"
#include "Identifier/Identifier.h"
 
#include <set>
#include <fstream>
 
namespace {
double trackQuality(const Trk::Track* Tr) {
 
  double quality = 0.;
  double baseScorePerHit = 17.;
  for (const Trk::TrackStateOnSurface* m : *(Tr->trackStateOnSurfaces())) {
    if (not m->type(Trk::TrackStateOnSurface::Measurement))
      continue;
    const Trk::FitQualityOnSurface fq = m->fitQualityOnSurface();
    if (!fq)
      continue;
 
    double x2 = fq.chiSquared();
    double hitQualityScore;
    if (fq.numberDoF() == 2)
      hitQualityScore = (1.2 * (baseScorePerHit - x2 * .5));  // pix
    else
      hitQualityScore = (baseScorePerHit - x2);  // sct
    if (hitQualityScore < 0.)
      hitQualityScore =
          0.;  // do not allow a bad hit to decrement the overall score
    quality += hitQualityScore;
  }
  if (Tr->info().trackProperties(Trk::TrackInfo::BremFit))
    quality *= 0.7;
 
  return quality;
}
}  // namespace
 
 
InDet::SiSPSeededTrackFinder::SiSPSeededTrackFinder
(const std::string& name,ISvcLocator* pSvcLocator) : AthReentrantAlgorithm(name, pSvcLocator)
{
}
 
 
StatusCode InDet::SiSPSeededTrackFinder::initialize() 
{
  ATH_CHECK(m_evtKey.initialize());
  ATH_CHECK(m_mbtsKey.initialize(m_useMBTS));
  ATH_CHECK(m_SpacePointsPixelKey.initialize(SG::AllowEmpty));
  ATH_CHECK(m_SpacePointsSCTKey.initialize(SG::AllowEmpty));
  ATH_CHECK(m_outputTracksKey.initialize());
 
  ATH_CHECK( m_prdToTrackMap.initialize( !m_prdToTrackMap.key().empty() ) );
 
  ATH_CHECK( m_caloClusterROIKey.initialize(m_useITkConvSeeded) );
 
  ATH_CHECK( m_seedsmaker.retrieve() );
  ATH_CHECK( m_zvertexmaker.retrieve( DisableTool{ not m_useZvertexTool } ));
 
  ATH_CHECK( m_trackmaker.retrieve());
 
  ATH_CHECK( m_regsel_strip.retrieve( DisableTool{ not m_useITkConvSeeded } ) );
 
  if (m_ITKGeometry and m_doFastTracking) ATH_CHECK(m_etaDependentCutsSvc.retrieve());
 
  ATH_CHECK( m_trackSummaryTool.retrieve( DisableTool{ m_trackSummaryTool.name().empty()} ));
 
  if (m_useNewStrategy and m_beamSpotKey.key().empty()) {
    m_useNewStrategy = false;
    m_useZBoundaryFinding = false;
  }
 
  ATH_CHECK(m_beamSpotKey.initialize( (m_useNewStrategy or m_useZBoundaryFinding or m_ITKGeometry or m_useITkConvSeeded) && not m_beamSpotKey.key().empty() ));
  if (m_useNewStrategy or m_useZBoundaryFinding or m_ITKGeometry or m_useITkConvSeeded) {
 
    if (not m_beamSpotKey.key().empty()) {
      ATH_CHECK( m_proptool.retrieve() );
 
      magneticFieldInit();
      
      if (m_histsize < 100) m_histsize = 100;
      m_zstep = static_cast<double>(m_histsize)/(2.*m_zcut);
    } else {
      m_proptool.disable();
      m_useNewStrategy = false;
      m_useZBoundaryFinding = false;
    }
  } else {
    m_proptool.disable();
  }
 
  if (msgLvl(MSG::DEBUG)) {
    dump(MSG::DEBUG, nullptr);
  }
  m_neventsTotal   = 0;
  m_neventsTotalV  = 0;
  m_problemsTotal  = 0;
  m_problemsTotalV = 0;
  return StatusCode::SUCCESS;
}
 
 
StatusCode InDet::SiSPSeededTrackFinder::execute(const EventContext& ctx) const
{ 
  if (m_ITKGeometry and m_doFastTracking) return itkFastTrackingStrategy(ctx);
  else if (m_useITkConvSeeded) return itkConvStrategy(ctx);
  else if (not m_useNewStrategy and not m_useZBoundaryFinding and not m_ITKGeometry) {
    return oldStrategy(ctx);
  }
  return newStrategy(ctx);
}
 
 
namespace InDet {
  class ExtendedSiTrackMakerEventData_xk final: public InDet::SiTrackMakerEventData_xk
  {
  public:
    explicit ExtendedSiTrackMakerEventData_xk(const SG::ReadHandleKey<Trk::PRDtoTrackMap> &key) { 
      if (!key.key().empty()) {
        m_prdToTrackMap = SG::ReadHandle<Trk::PRDtoTrackMap>(key);
        setPRDtoTrackMap(m_prdToTrackMap.cptr());
      }
    }
  private:
    SG::ReadHandle<Trk::PRDtoTrackMap> m_prdToTrackMap;
  };
}
 
 
StatusCode InDet::SiSPSeededTrackFinder::oldStrategy(const EventContext& ctx) const
{
  SG::WriteHandle<TrackCollection> outputTracks{m_outputTracksKey, ctx};
  ATH_CHECK(outputTracks.record(std::make_unique<TrackCollection>()));
  if (not isGoodEvent(ctx)) {
    return StatusCode::SUCCESS;
  }
 
  SiSpacePointsSeedMakerEventData seedEventData;
  bool ZVE = false;
  if (m_useZvertexTool) {
    std::list<Trk::Vertex> vertices = m_zvertexmaker->newEvent(ctx, seedEventData);
    if (not vertices.empty()) ZVE = true;
    m_seedsmaker->find3Sp(ctx, seedEventData, vertices);
  } else {
    m_seedsmaker->newEvent(ctx, seedEventData, -1);
    std::list<Trk::Vertex> vertexList;
    m_seedsmaker->find3Sp(ctx, seedEventData, vertexList);
  }
 
  const bool PIX = true;
  const bool SCT = true;
  InDet::ExtendedSiTrackMakerEventData_xk trackEventData(m_prdToTrackMap);
  m_trackmaker->newEvent(ctx, trackEventData, PIX, SCT);
 
  bool ERR = false;
  Counter_t counter{};
  const InDet::SiSpacePointsSeed* seed = nullptr;
  std::multimap<double, Trk::Track*> qualitySortedTrackCandidates;
  // Loop through all seed and reconsrtucted tracks collection preparation
  //
  while ((seed = m_seedsmaker->next(ctx, seedEventData))) {
    ++counter[kNSeeds];
    std::list<Trk::Track*> trackList = m_trackmaker->getTracks(ctx, trackEventData, seed->spacePoints());
    for (Trk::Track* t: trackList) {
      qualitySortedTrackCandidates.insert(std::make_pair(-trackQuality(t), t));
    }
    if (not ZVE and (counter[kNSeeds] >= m_maxNumberSeeds)) {
      ERR = true;
      ++m_problemsTotal;
      break;
    }
  }
  m_trackmaker->endEvent(trackEventData);
 
  // Remove shared tracks with worse quality
  //
  filterSharedTracks(qualitySortedTrackCandidates);
 
  // Save good tracks in track collection
  //
  for (const std::pair<const double, Trk::Track*> & qualityAndTrack: qualitySortedTrackCandidates) {
    ++counter[kNTracks];
    if (m_trackSummaryTool.isEnabled()) {
       m_trackSummaryTool->computeAndReplaceTrackSummary(*(qualityAndTrack.second),
                                                         false /* DO NOT suppress hole search*/);
    }
    outputTracks->push_back(qualityAndTrack.second);
  }
 
  m_counterTotal[kNSeeds] += counter[kNSeeds];
 
  if (ZVE) ++m_neventsTotalV;
  else     ++m_neventsTotal;
 
  if (ERR) {
    outputTracks->clear();
  } else {
    m_counterTotal[kNTracks] += counter[kNTracks];
  }
 
  // Print common event information
  //
  if (msgLvl(MSG::DEBUG)) {
    dump(MSG::DEBUG, &counter);
  }
 
  return StatusCode::SUCCESS;
}
 
// Execute with new strategy
 
StatusCode InDet::SiSPSeededTrackFinder::newStrategy(const EventContext& ctx) const
{
  SG::WriteHandle<TrackCollection> outputTracks{m_outputTracksKey, ctx};
  ATH_CHECK(outputTracks.record(std::make_unique<TrackCollection>()));
  if (not isGoodEvent(ctx)) {
    return StatusCode::SUCCESS;
  }
 
  SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle{m_beamSpotKey, ctx};
  Trk::PerigeeSurface beamPosPerigee(beamSpotHandle->beamPos());
 
  SiSpacePointsSeedMakerEventData seedEventData;
  
  m_seedsmaker->newEvent(ctx, seedEventData, 0);
  std::list<Trk::Vertex> vertexList;
  m_seedsmaker->find3Sp(ctx, seedEventData, vertexList);
 
  const bool PIX = true ;
  const bool SCT = true ;
  InDet::ExtendedSiTrackMakerEventData_xk trackEventData(m_prdToTrackMap);
  m_trackmaker->newEvent(ctx, trackEventData, PIX, SCT);
 
  std::vector<int> numberHistogram(m_histsize, 0);
  std::vector<double> zWeightedHistogram(m_histsize, 0.);
  std::vector<double> ptWeightedHistogram(m_histsize, 0.);
 
  bool ERR = false;
  Counter_t counter{};
  const InDet::SiSpacePointsSeed* seed = nullptr;
 
  std::multimap<double, Trk::Track*> qualitySortedTrackCandidates;
 
  bool doWriteNtuple = m_seedsmaker->getWriteNtupleBoolProperty();
  long EvNumber = 0.;            //Event number variable to be used for the validation ntuple 
 
  if (doWriteNtuple) {
    SG::ReadHandle<xAOD::EventInfo> eventInfo(m_evtKey,ctx);
    if(!eventInfo.isValid()) {EvNumber = -1.0;} else {EvNumber = eventInfo->eventNumber();}
  }
 
  while ((seed = m_seedsmaker->next(ctx, seedEventData))) {
 
    ++counter[kNSeeds];
    bool firstTrack{true};
 
      std::list<Trk::Track*> trackList = m_trackmaker->getTracks(ctx, trackEventData, seed->spacePoints());
      for (Trk::Track* t: trackList) {
 
        qualitySortedTrackCandidates.insert(std::make_pair(-trackQuality(t), t));
 
        if (firstTrack && m_doDumpGBTSTrainingDataLRT) {
          collectGBTSTrainingData(t);
        }
 
        if (firstTrack and not m_ITKGeometry) {
          fillZHistogram(t, beamPosPerigee, numberHistogram, zWeightedHistogram, ptWeightedHistogram);
        }
        firstTrack = false;
      }  
      if(doWriteNtuple) { m_seedsmaker->writeNtuple(seed, !trackList.empty() ? trackList.front() : nullptr, ISiSpacePointsSeedMaker::StripSeed, EvNumber) ; } 
        
 
    if (counter[kNSeeds] >= m_maxNumberSeeds) {
      ERR = true;
      ++m_problemsTotal;
      break;
    }
  }
 
  if(not m_SpacePointsPixelKey.empty()) {
    m_seedsmaker->newEvent(ctx, seedEventData, 1);
  } else {
    ATH_MSG_WARNING("SpacePointsPixelKey is empty. Skipping the second seeding pass that uses pixel seeds.");
  }
 
  std::pair<double,double> zBoundaries;
  if (not m_ITKGeometry) {
    findZvertex(vertexList, zBoundaries, numberHistogram, zWeightedHistogram, ptWeightedHistogram);
    m_seedsmaker->find3Sp(ctx, seedEventData, vertexList, &(zBoundaries.first));
  } else {
    m_seedsmaker->find3Sp(ctx, seedEventData, vertexList);
  }
 
  while ((seed = m_seedsmaker->next(ctx, seedEventData))) {
 
    ++counter[kNSeeds];
 
    std::list<Trk::Track*> trackList = m_trackmaker->getTracks(ctx, trackEventData, seed->spacePoints());
 
    for (Trk::Track* t: trackList) {
      qualitySortedTrackCandidates.insert(std::make_pair(-trackQuality(t), t));
    }
 
    if(doWriteNtuple) { m_seedsmaker->writeNtuple(seed, !trackList.empty() ? trackList.front() : nullptr, ISiSpacePointsSeedMaker::PixelSeed, EvNumber); }
 
    if (counter[kNSeeds] >= m_maxNumberSeeds) {
      ERR = true;
      ++m_problemsTotal;
      break;
    }
  }
 
  m_trackmaker->endEvent(trackEventData);
 
  filterSharedTracks(qualitySortedTrackCandidates);
 
  for (const std::pair<const double, Trk::Track*> & qualityAndTrack: qualitySortedTrackCandidates) {
    ++counter[kNTracks];
    if (m_trackSummaryTool.isEnabled()) {
       m_trackSummaryTool->computeAndReplaceTrackSummary(*qualityAndTrack.second,
                                                         false /* DO NOT suppress hole search*/);
       InDet::PatternHoleSearchOutcome theOutcome; 
       if (m_writeHolesFromPattern && trackEventData.combinatorialData().findPatternHoleSearchOutcome(qualityAndTrack.second,theOutcome)){
         qualityAndTrack.second->trackSummary()->update(Trk::numberOfPixelHoles, theOutcome.nPixelHoles); 
         qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTHoles, theOutcome.nSCTHoles); 
         qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTDoubleHoles, theOutcome.nSCTDoubleHoles); 
         qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTDeadSensors, theOutcome.nSCTDeads); 
         qualityAndTrack.second->trackSummary()->update(Trk::numberOfPixelDeadSensors, theOutcome.nPixelDeads); 
       }
    }
    outputTracks->push_back(qualityAndTrack.second);
    if (m_doDumpGBTSTrainingData) {
      collectGBTSTrainingData(qualityAndTrack.second);
    }
  }
 
  m_counterTotal[kNSeeds] += counter[kNSeeds] ;
 
  ++m_neventsTotal;
 
  if (ERR) {
    outputTracks->clear();
  } else {
    m_counterTotal[kNTracks] += counter[kNTracks];
  }
 
  // Print common event information
  //
  if (msgLvl(MSG::DEBUG)) {
    dump(MSG::DEBUG, &counter);
  }
  return StatusCode::SUCCESS;
}
 
 
 
// ITk fast tracking strategy
 
StatusCode InDet::SiSPSeededTrackFinder::itkFastTrackingStrategy(const EventContext& ctx) const
{
  SG::WriteHandle<TrackCollection> outputTracks{m_outputTracksKey, ctx};
  ATH_CHECK(outputTracks.record(std::make_unique<TrackCollection>()));
 
  const bool PIX = true ;
  const bool STRIP = true ;
  InDet::ExtendedSiTrackMakerEventData_xk trackEventData(m_prdToTrackMap);
  m_trackmaker->newTrigEvent(ctx, trackEventData, PIX, STRIP);
 
  SiSpacePointsSeedMakerEventData seedEventData;
 
  m_seedsmaker->newEvent(ctx, seedEventData, 0);
  std::list<Trk::Vertex> vertexList;
  m_seedsmaker->find3Sp(ctx, seedEventData, vertexList);
 
  bool ERR = false;
  Counter_t counter{};
  const InDet::SiSpacePointsSeed* seed = nullptr;
 
  std::multimap<double, Trk::Track*> qualitySortedTrackCandidates;
 
  bool doWriteNtuple = m_seedsmaker->getWriteNtupleBoolProperty();
  long EvNumber = 0.;            //Event number variable to be used for the validation ntuple
 
  if (doWriteNtuple) {
    SG::ReadHandle<xAOD::EventInfo> eventInfo(m_evtKey,ctx);
    if(!eventInfo.isValid()) {EvNumber = -1.0;} else {EvNumber = eventInfo->eventNumber();}
  }
 
  while ((seed = m_seedsmaker->next(ctx, seedEventData))) {
 
    ++counter[kNSeeds];
 
    std::list<Trk::Track*> trackList = m_trackmaker->getTracks(ctx, trackEventData, seed->spacePoints());
    for (Trk::Track* t: trackList) {
 
      qualitySortedTrackCandidates.insert(std::make_pair(-trackQuality(t), t));
 
    }
    if(doWriteNtuple) { m_seedsmaker->writeNtuple(seed, !trackList.empty() ? trackList.front() : nullptr, ISiSpacePointsSeedMaker::PixelSeed, EvNumber) ; }
 
    if (counter[kNSeeds] >= m_maxNumberSeeds) {
      ERR = true;
      ++m_problemsTotal;
      break;
    }
  }
 
  m_trackmaker->endEvent(trackEventData);
 
  filterSharedTracksFast(qualitySortedTrackCandidates);
 
  for (const std::pair<const double, Trk::Track*> & qualityAndTrack: qualitySortedTrackCandidates) {
    ++counter[kNTracks];
 
    if (m_trackSummaryTool.isEnabled()) {
      m_trackSummaryTool->computeAndReplaceTrackSummary(*qualityAndTrack.second);
      InDet::PatternHoleSearchOutcome theOutcome;
      if (m_writeHolesFromPattern && trackEventData.combinatorialData().findPatternHoleSearchOutcome(qualityAndTrack.second,theOutcome)){
        qualityAndTrack.second->trackSummary()->update(Trk::numberOfPixelHoles, theOutcome.nPixelHoles);
        qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTHoles, theOutcome.nSCTHoles);
        qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTDoubleHoles, theOutcome.nSCTDoubleHoles);
        qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTDeadSensors, theOutcome.nSCTDeads);
        qualityAndTrack.second->trackSummary()->update(Trk::numberOfPixelDeadSensors, theOutcome.nPixelDeads);
      }
    }
 
    outputTracks->push_back(qualityAndTrack.second);
  }
 
  m_counterTotal[kNSeeds] += counter[kNSeeds] ;
 
  ++m_neventsTotal;
 
  if (ERR) {
    outputTracks->clear();
  } else {
    m_counterTotal[kNTracks] += counter[kNTracks];
  }
 
  // Print common event information
  //
  if (msgLvl(MSG::DEBUG)) {
    dump(MSG::DEBUG, &counter);
  }
  return StatusCode::SUCCESS;
}
 
 
// Conversion Strategy for ITk
 
StatusCode InDet::SiSPSeededTrackFinder::itkConvStrategy(const EventContext& ctx) const
{
  ATH_MSG_DEBUG("Executing " << name() << "::itkConvStrategy");
 
  SG::WriteHandle<TrackCollection> outputTracks{m_outputTracksKey, ctx};
  ATH_CHECK(outputTracks.record(std::make_unique<TrackCollection>()));
  if (not isGoodEvent(ctx)) {
    return StatusCode::SUCCESS;
  }
 
  SiSpacePointsSeedMakerEventData seedEventData;
 
  SG::ReadHandle<ROIPhiRZContainer> calo_rois(m_caloClusterROIKey, ctx);
  std::unique_ptr<RoiDescriptor> roiComp = std::make_unique<RoiDescriptor>(true);
 
  if(calo_rois.isValid()) {
    RoiDescriptor * roi =nullptr;
    SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle{m_beamSpotKey, ctx};
    double beamZ = beamSpotHandle->beamVtx().position().z();
    roiComp->clear();
    roiComp->setComposite();
 
    const ROIPhiRZContainer &calo_rois_ref=*calo_rois;
    for (const ROIPhiRZ &calo_roi : calo_rois_ref) {
        double phi = calo_roi.phi();
        if (std::abs(phi)>=M_PI && phi!=-M_PI) continue; // skip duplicates < -pi and >pi
        double eta = calo_roi.eta();
        double z = beamZ;
        double roiPhiMin = phi - m_deltaPhi;
        double roiPhiMax = phi + m_deltaPhi;
        double roiEtaMin = eta - m_deltaEta;
        double roiEtaMax = eta + m_deltaEta;
        double roiZMin = beamZ - m_deltaZ;
        double roiZMax = beamZ + m_deltaZ;
        roi = new RoiDescriptor( eta, roiEtaMin, roiEtaMax,phi, roiPhiMin ,roiPhiMax,z,roiZMin,roiZMax);
        roiComp->push_back(roi);
    }
  }
  else {
    ATH_MSG_ERROR("Calo RoI is not valid: " << m_caloClusterROIKey.key());
    return StatusCode::FAILURE;
  }
 
  std::vector<IdentifierHash> listOfStripIds;
  std::vector<IdentifierHash> listOfPixIds;
 
  m_regsel_strip->lookup(ctx)->HashIDList( *roiComp, listOfStripIds );
  
  m_seedsmaker->newRegion(ctx, seedEventData, listOfPixIds, listOfStripIds);
  std::list<Trk::Vertex> vertexList;
  m_seedsmaker->find3Sp(ctx, seedEventData, vertexList);
 
  const bool PIX = true ;
  const bool STRIP = true ;
  InDet::ExtendedSiTrackMakerEventData_xk trackEventData(m_prdToTrackMap);
  m_trackmaker->newEvent(ctx, trackEventData, PIX, STRIP);
 
  bool ERR = false;
  Counter_t counter{};
  const InDet::SiSpacePointsSeed* seed = nullptr;
 
  std::multimap<double, Trk::Track*> qualitySortedTrackCandidates;
 
  bool doWriteNtuple = m_seedsmaker->getWriteNtupleBoolProperty();
  long EvNumber = 0.;            //Event number variable to be used for the validation ntuple
 
  if (doWriteNtuple) {
    SG::ReadHandle<xAOD::EventInfo> eventInfo(m_evtKey,ctx);
    if(!eventInfo.isValid()) {EvNumber = -1.0;} else {EvNumber = eventInfo->eventNumber();}
  }
 
  while ((seed = m_seedsmaker->next(ctx, seedEventData))) {
 
    ++counter[kNSeeds];
 
    std::list<Trk::Track*> trackList = m_trackmaker->getTracks(ctx, trackEventData, seed->spacePoints());
    for (Trk::Track* t: trackList) {
      qualitySortedTrackCandidates.insert(std::make_pair(-trackQuality(t), t));
    }
 
    if(doWriteNtuple) { m_seedsmaker->writeNtuple(seed, !trackList.empty() ? trackList.front() : nullptr, ISiSpacePointsSeedMaker::StripSeed, EvNumber) ; }
 
    if (counter[kNSeeds] >= m_maxNumberSeeds) {
      ERR = true;
      ++m_problemsTotal;
      break;
    }
  }
 
  m_trackmaker->endEvent(trackEventData);
 
  filterSharedTracks(qualitySortedTrackCandidates);
 
  for (const std::pair<const double, Trk::Track*> & qualityAndTrack: qualitySortedTrackCandidates) {
    ++counter[kNTracks];
    if (m_trackSummaryTool.isEnabled()) {
      m_trackSummaryTool->computeAndReplaceTrackSummary(*qualityAndTrack.second,
                            false /* DO NOT suppress hole search*/);
      InDet::PatternHoleSearchOutcome theOutcome;
      if (m_writeHolesFromPattern && trackEventData.combinatorialData().findPatternHoleSearchOutcome(qualityAndTrack.second,theOutcome)){
    qualityAndTrack.second->trackSummary()->update(Trk::numberOfPixelHoles, theOutcome.nPixelHoles);
    qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTHoles, theOutcome.nSCTHoles);
    qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTDoubleHoles, theOutcome.nSCTDoubleHoles);
    qualityAndTrack.second->trackSummary()->update(Trk::numberOfSCTDeadSensors, theOutcome.nSCTDeads);
    qualityAndTrack.second->trackSummary()->update(Trk::numberOfPixelDeadSensors, theOutcome.nPixelDeads);
      }
    }
    outputTracks->push_back(qualityAndTrack.second);
  }
 
  m_counterTotal[kNSeeds] += counter[kNSeeds] ;
 
  ++m_neventsTotal;
 
  if (ERR) {
    outputTracks->clear();
  } else {
    m_counterTotal[kNTracks] += counter[kNTracks];
  }
 
  // Print common event information
  //
  if (msgLvl(MSG::DEBUG)) {
    dump(MSG::DEBUG, &counter);
  }
  return StatusCode::SUCCESS;
}
 
 
 
// Finalize
 
StatusCode InDet::SiSPSeededTrackFinder::finalize()
{
  
    dump(MSG::INFO, &m_counterTotal);
 
  if (m_doDumpGBTSTrainingData || m_doDumpGBTSTrainingDataLRT) {
    dumpGBTSTrainingData();
  }
  return StatusCode::SUCCESS;
}
 
// Dumps relevant information into the MsgStream
 
MsgStream& InDet::SiSPSeededTrackFinder::dump(MSG::Level assign_level, const InDet::SiSPSeededTrackFinder::Counter_t* counter) const
{
  msg(assign_level) <<std::endl;
  MsgStream& out_msg=msg();
  if (counter) dumpevent(out_msg ,*counter);
  else dumptools(out_msg);
  out_msg << endmsg;
  return out_msg;
}
 
// Dumps conditions information into the MsgStream
 
MsgStream& InDet::SiSPSeededTrackFinder::dumptools(MsgStream& out) const
{
  int n = 65-m_zvertexmaker.type().size();
  std::string s1; for (int i=0; i<n; ++i) s1.append(" "); s1.append("|");
  n     = 65-m_seedsmaker.type().size();
  std::string s2; for (int i=0; i<n; ++i) s2.append(" "); s2.append("|");
  n     = 65-m_trackmaker.type().size();
  std::string s3; for (int i=0; i<n; ++i) s3.append(" "); s3.append("|");
  n     = 65-m_outputTracksKey.key().size();
  std::string s4; for (int i=0; i<n; ++i) s4.append(" "); s4.append("|");
 
  std::string s5; 
  if (m_useZvertexTool) s5= "Yes"; else s5 = "No";
  n     = 65-s5.size(); for (int i=0; i<n; ++i) s5.append(" "); s5.append("|");
  
  out<<"|----------------------------------------------------------------"
     <<"----------------------------------------------------|"
     <<std::endl;
  out<<"| Use primary vertices z-coordinates      finding?| "<<s5
     <<std::endl;
  if (m_useZvertexTool) {
    out<<"| Tool for primary vertices z-coordinates finding | "<<m_zvertexmaker.type()<<s1
       <<std::endl;
  }
  out<<"| Tool for space points seeds             finding | "<<m_seedsmaker.type()<<s2
     <<std::endl;
  out<<"| Tool for space points seeded track      finding | "<<m_trackmaker.type()<<s3
     <<std::endl;
  out<<"| Location of output tracks                       | "<<m_outputTracksKey.key()<<s4
     <<std::endl;
  out<<"|----------------------------------------------------------------"
     <<"----------------------------------------------------|"
     <<std::endl;
  return out;
}
 
// Dumps event information into the ostream
 
MsgStream& InDet::SiSPSeededTrackFinder::dumpevent(MsgStream& out, const InDet::SiSPSeededTrackFinder::Counter_t& counter) const
{
  out<<"|-------------------------------------------------------------------";
  out<<"---------------------------------|"
     <<std::endl;
  out<<"|  Investigated "
     <<std::setw(9)<<counter[kNSeeds]<<" space points seeds and found ";
  out<<std::setw(9)<<counter[kNTracks];
  if      (m_ITKGeometry        ) out<<" tracks using new strategy for ITK   |"<<std::endl; 
  else if (m_useNewStrategy     ) out<<" tracks using new strategy ("<<std::setw(2)<<  m_nvertex <<")      |"<<std::endl;
  else if (m_useZBoundaryFinding) out<<" tracks using old strategy with Zb   |"<<std::endl;
  else                            out<<" tracks using old strategy           |"<<std::endl;
 
  out<<"|-------------------------------------------------------------------";
  out<<"---------------------------------|"
     <<std::endl;
  if (m_problemsTotal or m_problemsTotalV) {
    out<<"|  Events       "
       <<std::setw(7)<<m_neventsTotal   <<"            without Z-vertz   "
       <<std::setw(7)<<m_neventsTotalV  <<"            with    Z-vertex         |"
       <<std::endl;
    out<<"|  Problems     "
       <<std::setw(7)<<m_problemsTotal  <<"            without Z-vertz   "
       <<std::setw(7)<<m_problemsTotalV <<"            with    Z-vertex         |"
       <<std::endl;
    out<<"|-------------------------------------------------------------------";
    out<<"-----------------------------|"
       <<std::endl;
  }
  return out;
}
 
// Test is it good event for reconstruction (mainly for HI events)
 
bool InDet::SiSPSeededTrackFinder::isGoodEvent(const EventContext& ctx) const {
 
  if ( not m_alwaysProtectAgainstBusyEvent ) { // if not enabled the protection is only applied to specific events
    // Test MBTS information from calorimeter
    //
    if (not m_useMBTS) return true;
 
    SG::ReadHandle<xAOD::EventInfo> eventInfo{m_evtKey, ctx};
    if (not eventInfo->isEventFlagBitSet(xAOD::EventInfo::Background, xAOD::EventInfo::MBTSTimeDiffHalo) ) {
      return true;
    }
  }
 
  // Test total number pixels space points
  //
  unsigned int nsp = 0;
  if (not m_SpacePointsPixelKey.empty()) {
    SG::ReadHandle<SpacePointContainer> spacePointsPixel{m_SpacePointsPixelKey, ctx};
    if (spacePointsPixel.isValid()) {
      for (const SpacePointCollection* spc: *spacePointsPixel) {
        nsp += spc->size();
      }
      if (static_cast<int>(nsp) > m_maxPIXsp) {
        ATH_MSG_WARNING("Found more than "<<m_maxPIXsp<<" pixels space points in background event. Skip track finding");
        return false;
      }
    }
  }
  // Test total number sct space points
  //
  nsp = 0;
  if (not m_SpacePointsSCTKey.empty()) {
    SG::ReadHandle<SpacePointContainer> spacePointsSCT{m_SpacePointsSCTKey, ctx};
    if (spacePointsSCT.isValid()) {
      for (const SpacePointCollection* spc: *spacePointsSCT) {
        nsp += spc->size();
      }
      if (static_cast<int>(nsp) > m_maxSCTsp) {
        ATH_MSG_WARNING("Found more than "<<m_maxSCTsp<<" sct space points in background event. Skip track finding");
        return false;
      }
    }
  }
 
  return true;
}
 
// Filer shared tracks
 
void InDet::SiSPSeededTrackFinder::filterSharedTracks(std::multimap<double, Trk::Track*>& qualitySortedTracks) const
{
  std::set<const Trk::PrepRawData*> clusters;
  
  std::vector<const Trk::PrepRawData*> freeClusters;
  freeClusters.reserve(15);    
  
  std::multimap<double, Trk::Track*>::iterator it_qualityAndTrack = qualitySortedTracks.begin();
 
  while (it_qualityAndTrack!=qualitySortedTracks.end()) {
    freeClusters.clear();
 
    std::set<const Trk::PrepRawData*>::iterator it_clustersEnd = clusters.end();
 
    int nClusters = 0; 
    for (const Trk::MeasurementBase* m: *((*it_qualityAndTrack).second->measurementsOnTrack())) {
 
      const Trk::PrepRawData* pr = (static_cast<const Trk::RIO_OnTrack*>(m))->prepRawData();
      if (pr) {
        ++nClusters;
        if (clusters.find(pr)==it_clustersEnd) {
          freeClusters.push_back(pr); 
        }
      }
    }
 
    int nFreeClusters = static_cast<int>(freeClusters.size()); 
    if (nFreeClusters >= m_nfreeCut || nFreeClusters==nClusters) {
      clusters.insert(freeClusters.begin(), freeClusters.end());
      ++it_qualityAndTrack;
    } else {
      delete (*it_qualityAndTrack).second;
      qualitySortedTracks.erase(it_qualityAndTrack++);
    }
  }
}
 
 
void InDet::SiSPSeededTrackFinder::filterSharedTracksFast(std::multimap<double, Trk::Track*>& qualitySortedTracks) const
{
  std::set<const Trk::PrepRawData*> clusters;
 
  std::vector<const Trk::PrepRawData*> freeClusters;
  freeClusters.reserve(15);
 
  std::multimap<double, Trk::Track*>::iterator it_qualityAndTrack = qualitySortedTracks.begin();
 
  while (it_qualityAndTrack!=qualitySortedTracks.end()) {
    freeClusters.clear();
 
    std::set<const Trk::PrepRawData*>::iterator it_clustersEnd = clusters.end();
 
    int nClusters = 0;
    int nPixels = 0;
    for (const Trk::TrackStateOnSurface* tsos: *((*it_qualityAndTrack).second->trackStateOnSurfaces())) {
 
      if(!tsos->type(Trk::TrackStateOnSurface::Measurement)) continue;
      const Trk::FitQualityOnSurface fq =  tsos->fitQualityOnSurface();
      if(!fq) continue;
      if(fq.numberDoF() == 2) ++nPixels;
 
      const Trk::MeasurementBase* mb = tsos->measurementOnTrack();
      const Trk::RIO_OnTrack*     ri = dynamic_cast<const Trk::RIO_OnTrack*>(mb);
      if(!ri) continue;
      const Trk::PrepRawData* pr = ri->prepRawData();
      if (pr) {
        ++nClusters;
        if (clusters.find(pr)==it_clustersEnd) {
          freeClusters.push_back(pr);
        }
      }
    }
 
    clusters.insert(freeClusters.begin(), freeClusters.end());
 
    int nFreeClusters = static_cast<int>(freeClusters.size());
    if( passEtaDepCuts( (*it_qualityAndTrack).second, nClusters, nFreeClusters, nPixels) ){
      ++it_qualityAndTrack;
    } else {
      delete (*it_qualityAndTrack).second;
      qualitySortedTracks.erase(it_qualityAndTrack++);
    }
  }
}
 
// Fill z coordinate histogram
 
void InDet::SiSPSeededTrackFinder::fillZHistogram(const Trk::Track* Tr,
                                                  const Trk::PerigeeSurface& beamPosPerigee,
                                                  std::vector<int>& numberHistogram,
                                                  std::vector<double>& zWeightedHistogram,
                                                  std::vector<double>& ptWeightedHistogram) const
{
  
  if (Tr->measurementsOnTrack()->size() < 10) return;
    
  const Trk::TrackParameters* paramsAtFirstSurface = Tr->trackStateOnSurfaces()->front()->trackParameters();
  Amg::Vector3D              position = paramsAtFirstSurface->position()          ;
  Amg::Vector3D              momentum = paramsAtFirstSurface->momentum()          ;
  
  constexpr double rSquare_max_forZHisto = 60.*60.; 
  if (position.x()*position.x()+position.y()*position.y() >= rSquare_max_forZHisto) return;
 
  double pT = sqrt(momentum.x()*momentum.x()+momentum.y()*momentum.y());
  if (pT < m_pTcut) return;
 
  Trk::PatternTrackParameters TP;
  if (not TP.production(paramsAtFirstSurface)) return;
      
  double step;
  if (not m_proptool->propagate(Gaudi::Hive::currentContext(),
                                TP, beamPosPerigee, TP, Trk::anyDirection, m_fieldprop, step, Trk::pion)) return;
          
  const AmgVector(5)& parsAtBeamSpot = TP.parameters();
  if (std::abs(parsAtBeamSpot[0]) > m_imcut) return;
  int z = static_cast<int>((parsAtBeamSpot[1]+m_zcut)*m_zstep);
  if (z >=0 and z < m_histsize) {
    ++numberHistogram[z];
    zWeightedHistogram[z] += parsAtBeamSpot[1];
    ptWeightedHistogram[z] += pT;
  }
  
}
 
// Find verteex  z coordinates
 
void  InDet::SiSPSeededTrackFinder::findZvertex(std::list<Trk::Vertex>& vertexZList,
                                                std::pair<double, double> & zBoundaries,
                                                const std::vector<int>& numberHistogram,
                                                const std::vector<double>& zWeightedHistogram,
                                                const std::vector<double>& ptWeightedHistogram) const
{
  zBoundaries = {1000., -1000};
 
  std::multimap<int   ,double> vertexZ_sortedByNtracks;
  std::multimap<double,double> vertexZ_sortedBySumPt;
 
  int lastBin = m_histsize-1;
  int minBinContentSum = 3;
 
  for (int binIndex=1; binIndex<lastBin; ++binIndex) { 
 
    int vertexNtracks = numberHistogram.at(binIndex-1)+numberHistogram.at(binIndex)+numberHistogram.at(binIndex+1);
 
    if (vertexNtracks>=minBinContentSum and (numberHistogram.at(binIndex) >= numberHistogram.at(binIndex-1) and numberHistogram.at(binIndex) >= numberHistogram.at(binIndex+1))) {
      double vertexZestimate = (zWeightedHistogram.at(binIndex-1)+zWeightedHistogram.at(binIndex)+zWeightedHistogram.at(binIndex+1))/static_cast<double>(vertexNtracks);
 
      if (vertexZestimate < zBoundaries.first) zBoundaries.first = vertexZestimate;
      if (vertexZestimate > zBoundaries.second) zBoundaries.second = vertexZestimate;
 
      if (m_useNewStrategy) { 
        double vertexSumPt = ptWeightedHistogram.at(binIndex-1)+ptWeightedHistogram.at(binIndex)+ptWeightedHistogram.at(binIndex+1);
        vertexZ_sortedByNtracks.insert(std::make_pair(-vertexNtracks, vertexZestimate));
        vertexZ_sortedBySumPt.insert(std::make_pair(-vertexSumPt, vertexZestimate));
      }
    }
  }
 
  if (m_useNewStrategy) {
 
    std::set<double> leadingVertices;
    int n = 0;
    std::multimap<double, double>::iterator vertex_pt_and_z = vertexZ_sortedBySumPt.begin();
    for (std::pair<int, double> nTrackAndZ: vertexZ_sortedByNtracks) {
      if (n++ >= m_nvertex) break; 
      leadingVertices.insert(nTrackAndZ.second);
      leadingVertices.insert((*vertex_pt_and_z++).second);
    }
 
    for (double v: leadingVertices) {
      vertexZList.emplace_back(Amg::Vector3D{0.,0.,v});
    }
  }
  if (zBoundaries.first > zBoundaries.second) {
    zBoundaries.first = -1000.;
    zBoundaries.second = +1000.;
  } else {
    zBoundaries.first -= 20.;
    zBoundaries.second += 20.;
  }
} 
 
 
// Callback function - get the magnetic field /
 
void InDet::SiSPSeededTrackFinder::magneticFieldInit() 
{
  // Build MagneticFieldProperties 
  //
  if(m_fieldmode == "NoField") {
    m_fieldprop = Trk::MagneticFieldProperties(Trk::NoField);
  } else {
    m_fieldprop = Trk::MagneticFieldProperties(Trk::FastField);
  }
}
 
 
// Check if track passes eta-dependent cuts for fast tracking
 
bool InDet::SiSPSeededTrackFinder::passEtaDepCuts(const Trk::Track* track,
                          int nClusters,
                          int nFreeClusters,
                          int nPixels) const
{
  Trk::TrackStates::const_iterator  m = track->trackStateOnSurfaces()->begin();
  const Trk::TrackParameters* par = (*m)->trackParameters();
  if(!par) return false;
 
  double eta = std::abs(par->eta());
  if(nClusters               < m_etaDependentCutsSvc->getMinSiHitsAtEta(eta)) return false;
  if(nFreeClusters           < m_etaDependentCutsSvc->getMinSiNotSharedAtEta(eta)) return false;
  if(nClusters-nFreeClusters > m_etaDependentCutsSvc->getMaxSharedAtEta(eta)) return false;
  if(nPixels                 < m_etaDependentCutsSvc->getMinPixelHitsAtEta(eta)) return false;
 
  if(par->pT() < m_etaDependentCutsSvc->getMinPtAtEta(eta)) return false;
  if(!(*m)->type(Trk::TrackStateOnSurface::Perigee)) return true ;
  if(std::abs(par->localPosition()[0]) > m_etaDependentCutsSvc->getMaxPrimaryImpactAtEta(eta)) return false;
  return true;
}
 
void InDet::SiSPSeededTrackFinder::collectGBTSTrainingData(const Trk::Track* track) const {
 
  struct VLM_Data {
    int vol_id, lay_id, mod_id;
    float m_x, m_y, m_z;
  };
 
  const PixelID* IDp = 0;
  const SCT_ID* IDs = 0;
 
  if (detStore()->retrieve(IDp, "PixelID").isFailure()) {
    ATH_MSG_FATAL("Could not get Pixel ID helper");
  }
 
  if (detStore()->retrieve(IDs, "SCT_ID").isFailure()) {
    ATH_MSG_FATAL("Could not get SCT ID helper");
  }
 
  if (!IDs && !IDp) return;
 
  ++m_numGBTSTrainingData;
 
  std::vector<VLM_Data> vlm;
 
  for (const auto* s : *track->trackStateOnSurfaces()) {
    if (!s->type(Trk::TrackStateOnSurface::Measurement)) continue;
 
    const Trk::MeasurementBase* mb = s->measurementOnTrack();
    if (!mb) continue;
 
    const Trk::RIO_OnTrack* ri = dynamic_cast<const Trk::RIO_OnTrack*>(mb);
    if (!ri) continue;
 
    const Trk::PrepRawData* rd = ri->prepRawData();
    if (!rd) continue;
 
    const InDet::SiCluster* si = dynamic_cast<const InDet::SiCluster*>(rd);
    if (!si) continue;
 
    const Amg::Vector3D& pos = s->trackParameters()->position();
 
    if (dynamic_cast<const InDet::PixelCluster*>(si)) {  // Pixel
 
      Identifier id = si->identify();
 
      int bec = IDp->barrel_ec(id);
 
      int vol_id = 8;
 
      if (bec == -2) vol_id = 7;
      if (bec == 2) vol_id = 9;
 
      if (bec < -2 || bec > 2) continue;
 
      int lay_id = IDp->layer_disk(id);
      int eta_mod = IDp->eta_module(id);
      int phi_mod = IDp->phi_module(id);
 
      Identifier wafer_id = IDp->wafer_id(bec, lay_id, phi_mod, eta_mod);
 
      int mod_id = IDp->wafer_hash(wafer_id);
 
      int new_vol = 0, new_lay = 0;
 
      if (vol_id == 7 || vol_id == 9) {
        new_vol = 10 * vol_id + lay_id;
        new_lay = eta_mod;
      } else if (vol_id == 8) {
        new_lay = 0;
        new_vol = 10 * vol_id + lay_id;
      }
      if (vol_id != 0)
        vlm.emplace_back(new_vol, new_lay, mod_id, pos.x(), pos.y(), pos.z());
    }
 
    if (dynamic_cast<const InDet::SCT_Cluster*>(si)) {  // SCT
 
      Identifier id = si->identify();
 
      int bec = IDs->barrel_ec(id);
 
      int vol_id = 13;
 
      if (bec < 0) vol_id = 12;
      if (bec > 0) vol_id = 14;
 
      int lay_id = IDs->layer_disk(id);
      int eta_mod = IDs->eta_module(id);
      int phi_mod = IDs->phi_module(id);
      int side = IDs->side(id);
 
      Identifier wafer_id = IDs->wafer_id(bec, lay_id, phi_mod, eta_mod, side);
 
      int mod_id = IDs->wafer_hash(wafer_id);
 
      vlm.emplace_back(vol_id, lay_id, mod_id, pos.x(), pos.y(), pos.z());
    }
  }
 
  // remove single-strip cases where no spacepoint exists
 
  std::vector<VLM_Data> vlm2;
 
  for (std::size_t it1 = 0; it1 < vlm.size() - 1; it1++) {
    if (vlm.at(it1).vol_id > 14) {  // Pixels
      vlm2.push_back(vlm.at(it1));
      continue;
    }
 
    std::size_t it2 = it1 + 1;
 
    int src = vlm.at(it1).vol_id * 1000 + vlm.at(it1).lay_id;
    int dst = vlm.at(it2).vol_id * 1000 + vlm.at(it2).lay_id;
 
    if (src == dst) {  // a spacepoint can be formed
      vlm2.push_back(vlm.at(it1));
      vlm2.push_back(vlm.at(it2));
      it1 = it2;
      continue;
    }
  }
 
  // remove track segments which are too short
  if (m_removeShortSegments) {
 
    constexpr float minDist = 20.0;
 
    for (auto it = std::next(vlm2.begin()); it != vlm2.end(); ) {
      auto jt = std::prev(it);
      float dx = it->m_x - jt->m_x;
      float dy = it->m_y - jt->m_y;
      float dz = it->m_z - jt->m_z;
 
      float dist = std::sqrt(dx*dx + dy*dy + dz*dz);
 
      if (dist < minDist) it = vlm2.erase(it);
      else ++it;
    }
  }
 
  std::scoped_lock trainingDataLock(m_GBTSTrainingDataMutex);
 
  for (std::size_t it1 = 0; it1 < vlm2.size() - 1; ++it1) {
    std::size_t it2 = it1 + 1;
 
    int src = vlm2.at(it1).vol_id * 1000 + vlm2.at(it1).lay_id;
    int dst = vlm2.at(it2).vol_id * 1000 + vlm2.at(it2).lay_id;
 
    if (src != dst) {  // skip the same layer
      auto [im1, new1] = m_GBTSTrainingData.insert({src, {}});
      auto [im2, new2] = im1->second.insert({dst, 1ul});
      if (!new2) im2->second++;
    }
  }
}
 
void InDet::SiSPSeededTrackFinder::dumpGBTSTrainingData() const {
  std::ofstream tableFile(m_GBTSTrainingDataFileName);
  tableFile << "from,to,probability,flow\n";
 
  unsigned long nTotal = 0;
  for (const auto& [src, conns] : m_GBTSTrainingData) {
    unsigned long nTotalDst = 0;
    for (const auto& [dst, n] : conns) {
      nTotalDst += n;
    }
    nTotal += nTotalDst;
    for (const auto& [dst, n] : conns) {
      double prob = double(n) / double(nTotalDst);
      tableFile << src << ", " << dst << ", " << std::fixed << std::setprecision(6) << prob << ", " << prob << '\n';
    }
  }
  ATH_MSG_INFO("GBTS training data from " << m_numGBTSTrainingData << " tracks with " << nTotal << " pairs written to " << m_GBTSTrainingDataFileName.value());
}