TRT_SeededTrackFinder.cxx

Go to the documentation of this file.

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

//   Implementation file for class InDet::TRT_SeededTrackFinder
// (c) ATLAS Detector software
// Version 1.0 02/20/2006 Thomas Koffas
// modified by Markus ELsing
 
#include "GaudiKernel/MsgStream.h"
#include "TRT_SeededTrackFinder/TRT_SeededTrackFinder.h"
#include "TrkTrack/TrackInfo.h"
#include "TrkPseudoMeasurementOnTrack/PseudoMeasurementOnTrack.h"

#include "InDetRIO_OnTrack/SiClusterOnTrack.h"
 
#include "SiSPSeededTrackFinderData/SiCombinatorialTrackFinderData_xk.h"
 
#include "TrkEventPrimitives/FitQualityOnSurface.h"
#include "TrkEventPrimitives/FitQuality.h"
 
#include "RoiDescriptor/RoiDescriptor.h"
 
#include <memory>
#include <cmath>
using namespace std;

// Constructor
 
InDet::TRT_SeededTrackFinder::TRT_SeededTrackFinder
(const std::string& name, ISvcLocator* pSvcLocator)
  : AthReentrantAlgorithm(name, pSvcLocator)
{ }

// Initialisation
 
StatusCode InDet::TRT_SeededTrackFinder::initialize()
{
 
  //Get the TRT seeded track maker tool
  //
  ATH_CHECK(m_trackmaker.retrieve());
 
  //Get the refitting tool
  //
  ATH_CHECK( m_fitterTool.retrieve(      DisableTool{ !m_doRefit } ));
  ATH_CHECK( m_extrapolator.retrieve(    DisableTool{ !m_SiExtensionCuts } ));
  ATH_CHECK( m_beamSpotKey.initialize(   m_SiExtensionCuts) );
 
  // optional PRD to track association map
  ATH_CHECK( m_prdToTrackMap.initialize( !m_prdToTrackMap.key().empty() ) );
 
  ATH_CHECK( m_trackSummaryTool.retrieve( DisableTool{ m_trackSummaryTool.name().empty() } ) );
 
  // Get tool for track extension to TRT
  //
  ATH_CHECK( m_trtExtension.retrieve(    DisableTool{ !m_doExtension} ));
 
 
  ATH_CHECK(  m_SegmentsKey.initialize()) ;  
  ATH_CHECK( m_outTracksKey.initialize());
 
  ATH_CHECK( m_caloClusterROIKey.initialize(m_caloSeededRoI) );
 
  if(m_caloSeededRoI){
    ATH_CHECK( m_regionSelector.retrieve());
  } else {
    m_regionSelector.disable();
  }
 
  // Get output print level
  //
  if(msgLvl(MSG::DEBUG)) {
     dumptools(msg(MSG::DEBUG));
  }
   //Global counters. See the include file for definitions
  m_totalStat = Stat_t();
  return StatusCode::SUCCESS;
}
 
namespace InDet {
  class MyPRDtoTrackMap  : public Trk::PRDtoTrackMap {
  public:
    unsigned int size() const {
      return m_prepRawDataTrackMap.size();
    }
  };
 
  class ExtendedSiCombinatorialTrackFinderData_xk : public InDet::SiCombinatorialTrackFinderData_xk
  {
  public:
    explicit ExtendedSiCombinatorialTrackFinderData_xk(const SG::ReadHandleKey<Trk::PRDtoTrackMap> &key) {
      if (!key.key().empty()) {
        m_prdToTrackMap = SG::ReadHandle<Trk::PRDtoTrackMap>(key);
        setPRDtoTrackMap(m_prdToTrackMap.cptr());
      }
    }
  private:
    void dummy() {}
    SG::ReadHandle<Trk::PRDtoTrackMap> m_prdToTrackMap;
  };
}

// Execute
StatusCode 
InDet::TRT_SeededTrackFinder::execute(const EventContext& ctx) const{
  //Counters. See the include file for definitions
  Stat_t ev_stat;
  // counter
  int nTrtSegCur = 0;
  SG::ReadHandle<Trk::SegmentCollection>  segments(m_SegmentsKey,ctx);
  if(!segments.isValid()){
    ATH_MSG_FATAL ("No segment with name " << segments.name() << " found in StoreGate!");
    return StatusCode::FAILURE;
  } else {
    ATH_MSG_DEBUG ("Found segments collection " << segments.name() << " in StoreGate!");
  }
  // number of segments + statistics
  ev_stat.m_counter[Stat_t::kNTrtSeg] = int(segments->size());
  ATH_MSG_DEBUG ("TRT track container size " << ev_stat.m_counter[Stat_t::kNTrtSeg]);
  if(ev_stat.m_counter[Stat_t::kNTrtSeg]>m_MaxSegNum) {
    ATH_MSG_DEBUG ("TRT track container size huge; will process event partially if number of max segments reached !!!");
  }
  // Event dependent data of SiCombinatorialTrackFinder_xk
  InDet::ExtendedSiCombinatorialTrackFinderData_xk combinatorialData(m_prdToTrackMap);
  std::unique_ptr<InDet::ITRT_SeededTrackFinder::IEventData> event_data_p;
  if(m_caloSeededRoI ) {
    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();
      for (const ROIPhiRZ &the_roi : *calo_rois) {
          double eta = the_roi.eta();
          double phi = the_roi.phi();
          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, beamZ, roiZMin,roiZMax);
          roiComp->push_back(roi);
      }
 
    }
    std::vector<IdentifierHash> listOfSCTIds;
    std::vector<IdentifierHash> listOfPixIds;
    
    m_regionSelector->lookup(ctx)->HashIDList( *roiComp, listOfSCTIds );
 
    event_data_p = m_trackmaker->newRegion(ctx, combinatorialData, listOfPixIds, listOfSCTIds);
  } else {
    event_data_p = m_trackmaker->newEvent(ctx, combinatorialData);
  }
 
  std::unique_ptr<InDet::ITRT_TrackExtensionTool::IEventData> ext_event_data_p( m_trtExtension->newEvent(ctx) );
  std::unique_ptr<TrackCollection> outTracks = std::make_unique<TrackCollection>();
  std::vector<Trk::Track*> tempTracks;                           //Temporary track collection
  tempTracks.reserve(128);
  // loop over event
  ATH_MSG_DEBUG ("Begin looping over all TRT segments in the event");
  Trk::SegmentCollection::const_iterator iseg    = segments->begin();
  Trk::SegmentCollection::const_iterator isegEnd = segments->end();
  InDet::TRT_DetElementLink_xk::TRT_DetElemUsedMap map;
  for(; iseg != isegEnd; ++ iseg) {
    // Get the track segment
    const Trk::TrackSegment *trackTRT = dynamic_cast<const Trk::TrackSegment*>(*iseg);
    if(!trackTRT){
      ATH_MSG_ERROR ("No pointer to segment !");
      continue;
    } else {
      // the segment finder is applying a TRT(eta) cut and a pt preselection, so we don't do that here
      //Ask for at least 10 TRT hits in order to process
      if(trackTRT->numberOfMeasurementBases() < m_minTRTonSegment) {
        ATH_MSG_DEBUG ("TRT segment fails nTRT hit cut, reject.");
        // statistics
        ev_stat.m_counter[Stat_t::kNTrtFailSel]++;
      } else {
        // do we continue to process ?
        nTrtSegCur++;
        if(nTrtSegCur>=m_MaxSegNum) {
          ATH_MSG_DEBUG ("====> Reached maximal number of segments in event, stop !!!");
          // statistics
          ev_stat.m_counter[Stat_t::kNTrtLimit]++;
          break;
        }
        // Get the number of the TRT track segment ROTs
        ATH_MSG_DEBUG ("=> New segment to process, number Of TRT ROTs : " << (trackTRT->numberOfMeasurementBases()));
        // statistics
        ev_stat.m_counter[Stat_t::Stat_t::kNTrtSegGood]++;
        // ok, call track maker and get list of possible track candidates
        std::list<Trk::Track*> trackSi = m_trackmaker->getTrack(ctx, *event_data_p, *trackTRT); //Get the possible Si extensions
        if (trackSi.empty()) {
          ATH_MSG_DEBUG ("No Si track candidates associated to the TRT track ");
          // statistics
          ev_stat.m_counter[Stat_t::kNTrtNoSiExt]++;
          // obsolete backup of TRT only
          if(m_saveTRT && trackTRT->numberOfMeasurementBases() > m_minTRTonly){
            Trk::Track* trtSeg = nullptr;trtSeg = segToTrack(ctx, *trackTRT);
            if(!trtSeg) {
              ATH_MSG_DEBUG ("Failed to make a track out of the TRT segment!");
              continue;
            }
            // statistics
            ev_stat.m_counter[Stat_t::kNBckTrk]++; ev_stat.m_counter[Stat_t::Stat_t::kNBckTrkTrt]++;
            // add track to output list
            outTracks->push_back(trtSeg);
          }
          continue;
        } else {
          // Found useful extensions
          ATH_MSG_DEBUG ("Found " << (trackSi.size()) << " Si tracks associated to the TRT track ");
          // Merge the resolved Si extensions with the original TRT track segment
          std::list<Trk::Track*>::const_iterator itt    = trackSi.begin();
          std::list<Trk::Track*>::const_iterator ittEnd = trackSi.end();
          for (; itt != ittEnd ; ++itt){
            tempTracks.push_back(*itt);
            // get list of TSOS
            const Trk::TrackStates* temptsos = (*itt)->trackStateOnSurfaces();
            if (!temptsos) {
              ATH_MSG_DEBUG ("Silicon extension empty ???");
              continue;
            }
            // Add the track to the list of tracks in the event
            ATH_MSG_DEBUG ("Silicon extension found has length of : " << temptsos->size());
            // do we do a preselection ?
            if (m_SiExtensionCuts) {
              // get parameters without errors
              auto input = (*itt)->trackParameters()->front()->uniqueClone();
              // cuts on parameters
              if (std::abs(input->pT()) < m_minPt) {
                ATH_MSG_DEBUG ("Track pt < "<<m_minPt<<", reject it");
                // statistics
                ev_stat.m_counter[Stat_t::kNExtCut]++;
                continue;
              }
              if (std::abs(input->eta()) > m_maxEta) {
                ATH_MSG_DEBUG ("Track eta > "<<m_maxEta<<", reject it");
                // statistics
                ev_stat.m_counter[Stat_t::kNExtCut]++;
                continue;
              }
              // --- beam spot position
              Amg::Vector3D beamSpotPosition(0,0,0);
              if (m_SiExtensionCuts){
                       SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey, ctx };
                       beamSpotPosition = beamSpotHandle->beamVtx().position();
              }
              // --- create surface
              Trk::PerigeeSurface perigeeSurface(beamSpotPosition);
 
              // uses perigee on track or extrapolates, no material in any case, we cut on impacts
              std::unique_ptr<const Trk::TrackParameters> parm =
                m_extrapolator->extrapolateDirectly(
                  ctx, *input, perigeeSurface);
              std::unique_ptr<const Trk::Perigee> extrapolatedPerigee = nullptr;
              if (parm && parm->associatedSurface().type() == Trk::SurfaceType::Perigee) {
                extrapolatedPerigee.reset(static_cast<const Trk::Perigee*>(parm.release()));
              }
              if (!extrapolatedPerigee) {
                ATH_MSG_WARNING("Extrapolation of perigee failed, this should never happen" );
                // statistics
                ev_stat.m_counter[Stat_t::kNExtCut]++;
                continue;
              }
 
              ATH_MSG_VERBOSE ("extrapolated perigee: "<<*extrapolatedPerigee);
              if (std::abs(extrapolatedPerigee->parameters()[Trk::d0]) > m_maxRPhiImp) {
                ATH_MSG_DEBUG ("Track Rphi impact > "<<m_maxRPhiImp<<", reject it");
                // statistics
                ev_stat.m_counter[Stat_t::kNExtCut]++;
                continue;
              }
              if (std::abs(extrapolatedPerigee->parameters()[Trk::z0]) > m_maxZImp) {
                ATH_MSG_DEBUG ("Track Z impact > "<<m_maxZImp<<", reject it");
                // statistics
                ev_stat.m_counter[Stat_t::kNExtCut]++;
                continue;
              }
            }
 
            // do re run a Track extension into TRT ?
            Trk::Track* globalTrackNew = nullptr;
            // do we have 4 and extension is enabled ?
            if(int(temptsos->size())>=4 && m_doExtension){
              // Add the track to the list of tracks in the event
              ATH_MSG_DEBUG ("Try to improve TRT calling extension tool.");
              // statistics
              ev_stat.m_counter[Stat_t::Stat_t::kNTrtExtCalls]++;
              // call extension tool
              std::vector<const Trk::MeasurementBase*>& tn =
                  m_trtExtension->extendTrack(ctx, *(*itt), *ext_event_data_p, map);
              if (tn.empty()) {
                // Fallback if extension failed
                ATH_MSG_DEBUG ("No new segment found, use input segment as fallback.");
                // statistics
                ev_stat.m_counter[Stat_t::Stat_t::kNTrtExtFail]++;
                // merge Si with input track segments
                globalTrackNew = mergeSegments(**itt,*trackTRT);
              } else if (!m_rejectShortExten ||
                         tn.size() >= trackTRT->numberOfMeasurementBases()) {
                // Use the extension to instead of the segment
                ATH_MSG_DEBUG ("Successful extension, number of TRT hits : " << tn.size() << " was : " << (trackTRT->numberOfMeasurementBases()));
                // merge the extension with the Si track
                globalTrackNew = mergeExtension(**itt,tn);
                // Add the track to the list of tracks in the event
                ATH_MSG_DEBUG ("Merged extension with Si segment");
                // statistics
                ev_stat.m_counter[Stat_t::kNTrtExt]++;
                // clean up
                std::vector<const Trk::MeasurementBase*>::const_iterator iv, ive=tn.end();
                for(iv=tn.begin(); iv!=ive; ++iv) delete (*iv);
              } else {
                // Extension is shorter, let's fall back onto the original
                ATH_MSG_DEBUG ("Extension too short, number of TRT hits : " << tn.size() << " was : " << (trackTRT->numberOfMeasurementBases()) << ". Use Segement !");
                // merge segments
                globalTrackNew = mergeSegments(**itt,*trackTRT);
                // Add the track to the list of tracks in the event
                ATH_MSG_DEBUG ("Merged TRT segment with Si segment");
                // statistics
                ev_stat.m_counter[Stat_t::Stat_t::kNTrtExtBad]++;
                // clean up
                std::vector<const Trk::MeasurementBase*>::const_iterator iv, ive=tn.end();
                for(iv=tn.begin(); iv!=ive; ++iv) delete (*iv);
              }
            } else {
              // no extension tool, jsut add the two
              ATH_MSG_DEBUG ("Do not try to extend Si track, merging it with input TRT.");
              // merge segments
              globalTrackNew = mergeSegments(**itt,*trackTRT);
            }
            // do we have an track candidate ?
            if(!globalTrackNew){
              ATH_MSG_DEBUG ("Failed to merge TRT+Si track segment !");
              if(m_saveTRT && trackTRT->numberOfMeasurementBases() > m_minTRTonly) {
                Trk::Track* trtSeg = nullptr;trtSeg = segToTrack(ctx, *trackTRT);
                if(!trtSeg){
                  ATH_MSG_DEBUG ("Failed to make a track out of the  TRT segment!");
                  continue;
                }
                ATH_MSG_DEBUG ("Add TRT only to output list");
                // statistis
                ev_stat.m_counter[Stat_t::kNBckTrk]++; ev_stat.m_counter[Stat_t::Stat_t::kNBckTrkTrt]++;
                // add it to output list
                if (m_trackSummaryTool.isEnabled()) {
                  m_trackSummaryTool->computeAndReplaceTrackSummary(*trtSeg,
                                      false /* DO NOT suppress hole search*/);
                }
                outTracks->push_back(trtSeg);
              }
            } else {
              ATH_MSG_DEBUG ("Save merged TRT+Si track segment!");
              // statistics
              ev_stat.m_counter[Stat_t::kNBckTrk]++; ev_stat.m_counter[Stat_t::Stat_t::kNBckTrkSi]++;
              // add it to output list
              if (m_trackSummaryTool.isEnabled()) {
                m_trackSummaryTool->computeAndReplaceTrackSummary(*globalTrackNew,
                                    false /* DO NOT suppress hole search*/);
              }
              outTracks->push_back(globalTrackNew);
            }
          }
        }
      }
    }
  }
  // further debugging of results
  if(m_doStat){
    Analyze(outTracks.get());
  }
  if (SG::WriteHandle<TrackCollection>(m_outTracksKey,ctx).record(std::move(outTracks)).isFailure()){
     ATH_MSG_ERROR("Failed to record " << m_outTracksKey.key());
     return StatusCode::FAILURE;
  }
  // Update the total counters
  {
     std::lock_guard<std::mutex> lock(m_statMutex);
     m_totalStat += ev_stat;
  }
 
  for (auto *p : tempTracks){
    delete p;
  }
  m_trackmaker->endEvent(*event_data_p);
 
  //Print common event information
  if(msgLvl(MSG::DEBUG)){
     dumpevent(msg(MSG::DEBUG), ev_stat);
  }
 
  return StatusCode::SUCCESS;
}

// Finalize
 
StatusCode InDet::TRT_SeededTrackFinder::finalize()
{
  if(msgLvl(MSG::INFO)){
     msg(MSG::INFO) << "\n";
     dumpevent(msg(MSG::INFO), m_totalStat);
     msg(MSG::INFO) << endmsg;
  }
  return StatusCode::SUCCESS;
}
 

// Dumps conditions information into the MsgStream
 
MsgStream& InDet::TRT_SeededTrackFinder::dumptools( MsgStream& out ) const
{
  int n = 65-m_trackmaker.type().size();
  std::string s1; for(int i=0; i<n; ++i) s1.append(" "); s1.append("|");
  n     = 65-m_fitterTool.type().size();
  std::string s2; for(int i=0; i<n; ++i) s2.append(" "); s2.append("|");
  n     = 65-m_trtExtension.type().size();
  std::string s3; for(int i=0; i<n; ++i) s3.append(" "); s3.append("|");
  n     = 65-m_SegmentsKey.key().size();
  std::string s4; for(int i=0; i<n; ++i) s4.append(" "); s4.append("|");
  n     = 65-m_outTracksKey.key().size();
  std::string s5; for(int i=0; i<n; ++i) s5.append(" "); s5.append("|");
 
  out<<"|----------------------------------------------------------------------"
     <<"-------------------|"
       <<std::endl;
  out<<"| Tool for TRT seeded track finding | "<<m_trackmaker.type()    <<s1<<std::endl;
  out<<"| Tool for final track refitting    | "<<m_fitterTool.type()    <<s2<<std::endl;
  out<<"| Tool for TRT trac extension       | "<<m_trtExtension.type()  <<s3<<std::endl;
  out<<"| Location of input tracks          | "<<m_SegmentsKey.key()      <<s4<<std::endl;
  out<<"| Location of output tracks         | "<<m_outTracksKey.key()     <<s5<<std::endl;
  out<<"|----------------------------------------------------------------------"
     <<"-------------------|"
      <<std::endl;
  return out;
}

// Dumps event information into the MsgStream
 
MsgStream& InDet::TRT_SeededTrackFinder::dumpevent( MsgStream& out, const InDet::TRT_SeededTrackFinder::Stat_t &stat ) const
{
  out<<"|-------------------------------------------------------------------"<<std::endl;
  out<<"|  Investigated :"<<std::endl
     <<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNTrtSeg]       <<" TRT Segments on input"<<std::endl
     <<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNTrtFailSel]   <<" TRT Segments fail selection on input"<<std::endl
     <<"| "<<std::setw(7)<<stat.m_counter[Stat_t::Stat_t::kNTrtSegGood]   <<" TRT Segments after selection"<<std::endl;
  if (stat.m_counter[Stat_t::kNTrtLimit]>0) {
    out<<"|-------------------------------------------------------------------"<<std::endl;
    out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNTrtLimit]   <<" TRT segments lost because of processing limit"<<std::endl;
  }
  out<<"|-------------------------------------------------------------------"<<std::endl;
  out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNTrtNoSiExt]   <<" TRT segments without Si extension"<<std::endl;
  if (m_SiExtensionCuts) {
    out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNExtCut]     <<" number of Si extension failing cuts"<<std::endl;
  }
  if (stat.m_counter[Stat_t::Stat_t::kNBckTrkTrt]>0) {
    out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::Stat_t::kNBckTrkTrt]  <<" number ot TRT only tracks created"<<std::endl;
  }
  if (m_doExtension) {
    out<<"|-------------------------------------------------------------------"<<std::endl;
    out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::Stat_t::kNTrtExtCalls]<<" number of times TRT extension is called"<<std::endl
       <<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNTrtExt]     <<" number of good TRT extension"<<std::endl;
    if (m_rejectShortExten) {
      out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::Stat_t::kNTrtExtBad]<<" number of bad TRT extension"<<std::endl;
    }
    out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::Stat_t::kNTrtExtFail]<<" number of failed TRT extension and fallback"<<std::endl;
  }
  out<<"|-------------------------------------------------------------------"<<std::endl;
  out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::Stat_t::kNBckTrkSi]    <<" TRT+Si tracks created of output"<<std::endl;
  if (stat.m_counter[Stat_t::kNBckTrkSi] != stat.m_counter[Stat_t::kNBckTrk]) {
    out<<"| "<<std::setw(7)<<stat.m_counter[Stat_t::kNBckTrk]    <<" total number of tracks on output"<<std::endl;
  }
  out<<"|-------------------------------------------------------------------";
  return out;
}
 

// Merge a Si extension and a TRT segment.Refit at the end
 
Trk::Track* InDet::TRT_SeededTrackFinder::mergeSegments(const Trk::Track& tT, const Trk::TrackSegment& tS) const {
  // TSOS from the track
  const Trk::TrackStates* stsos = tT.trackStateOnSurfaces();
  // fitQuality from track
  auto fq = tT.fitQuality()->uniqueClone();
  // output datavector of TSOS
  auto ntsos = std::make_unique<Trk::TrackStates>();
  int siHits = 0;
  // copy track Si states into track
  Trk::TrackStates::const_iterator p_stsos;
  for (p_stsos=stsos->begin(); p_stsos != stsos->end(); ++p_stsos) {
    ntsos->push_back( (*p_stsos)->clone() );
    if ((*p_stsos)->type(Trk::TrackStateOnSurface::Measurement)) siHits++;
  }
  // loop over segment
  for (int it = 0; it < int(tS.numberOfMeasurementBases()); it++) {
    //test if it is a pseudo measurement
    if ( dynamic_cast<const Trk::PseudoMeasurementOnTrack*>(tS.measurement(it)) ) {
      if (siHits < 4) {
        ATH_MSG_DEBUG ("Too few Si hits.Will keep pseudomeasurement...");
        const Trk::TrackStateOnSurface* seg_tsos = new Trk::TrackStateOnSurface(tS.measurement(it)->uniqueClone(), nullptr);
        ntsos->push_back(seg_tsos);
      }
    } else {
      std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
      typePattern.set(Trk::TrackStateOnSurface::Measurement);
      const Trk::TrackStateOnSurface* seg_tsos = new Trk::TrackStateOnSurface(tS.measurement(it)->uniqueClone(), nullptr, nullptr, typePattern);
      ntsos->push_back(seg_tsos);
    }
  }

  Trk::TrackInfo info;
  info.setPatternRecognitionInfo(Trk::TrackInfo::TRTSeededTrackFinder);
  std::unique_ptr<Trk::Track> newTrack(std::make_unique<Trk::Track>(info, std::move(ntsos), std::move(fq)));
 
  //Careful refitting at the end
  if (m_doRefit) {
    newTrack=m_fitterTool->fit(Gaudi::Hive::currentContext(),*newTrack, false, Trk::pion);
    if (!newTrack) {
      ATH_MSG_DEBUG ("Refit of TRT+Si track segment failed!");
      return nullptr;
    }
    const Trk::Perigee* perTrack=newTrack->perigeeParameters();
    if (perTrack) {
      const AmgSymMatrix(5)* CM = perTrack->covariance();
      if (!CM || std::sqrt((*CM)(1,1)) == 0. || std::sqrt((*CM)(3,3)) == 0.) {
        return nullptr;
      }
    }
  }
  return newTrack.release();
}

// Transform a TRT segment to track
 
Trk::Track* InDet::TRT_SeededTrackFinder::segToTrack(const EventContext& ctx, const Trk::TrackSegment& tS) const {
  ATH_MSG_DEBUG ("Transforming the TRT segment into a track...");
 
  //Get the track segment information and build the initial track parameters
  const Trk::StraightLineSurface* surf = dynamic_cast<const Trk::StraightLineSurface*>(&(tS.associatedSurface()));
  if (!surf) {
    throw std::logic_error("Unhandled surface.");
  }
  const AmgVector(5)& p = tS.localParameters();
  AmgSymMatrix(5) ep = AmgSymMatrix(5)(tS.localCovariance());
  auto ntsos = std::make_unique<Trk::TrackStates>();
  std::unique_ptr<Trk::TrackParameters> segPar =
    surf->createUniqueParameters<5, Trk::Charged>(
      p(0), p(1), p(2), p(3), p(4), std::move(ep));
  if (segPar) {
    ATH_MSG_DEBUG( "Initial TRT Segment Parameters for refitting " << (*segPar) );
  } else {
    ATH_MSG_DEBUG( "Could not get initial TRT segment parameters! " );
    return nullptr;
  }
  for (int it = 0; it < int(tS.numberOfMeasurementBases()); it++) {
    // on first measurement add parameters
    const Trk::TrackStateOnSurface* seg_tsos = nullptr;
    std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
    typePattern.set(Trk::TrackStateOnSurface::Measurement);
    if (it == 0){
      seg_tsos = new Trk::TrackStateOnSurface(tS.measurement(it)->uniqueClone(), std::move(segPar), nullptr, typePattern);
    } else {
      seg_tsos = new Trk::TrackStateOnSurface(tS.measurement(it)->uniqueClone(), nullptr, nullptr, typePattern);
    }
    ntsos->push_back(seg_tsos);
  }
  Trk::TrackInfo info;
  info.setPatternRecognitionInfo(Trk::TrackInfo::TRTSeededTrackFinder);
  std::unique_ptr<Trk::Track> newTrack = std::make_unique<Trk::Track>(info, std::move(ntsos), nullptr);
  // Careful refitting of the TRT stand alone track
  if (m_doRefit) {
    newTrack = m_fitterTool->fit(ctx,*newTrack, false, Trk::pion);
    if (!newTrack) {
      ATH_MSG_DEBUG ("Refit of TRT track segment failed!");
      return nullptr;
    }
    //Protect for tracks that have no really defined locz and theta parameters
    const Trk::Perigee* perTrack=newTrack->perigeeParameters();
    if (perTrack) {
      const AmgSymMatrix(5)* CM = perTrack->covariance();
      if (!CM || std::sqrt((*CM)(1,1)) == 0. || std::sqrt((*CM)(3,3)) == 0.) {
         return nullptr;
      }
    }
  }
  return newTrack.release();
}

// Merge a Si extension and a TRT track extension.Refit at the end
 
Trk::Track* InDet::TRT_SeededTrackFinder::
mergeExtension(const Trk::Track& tT, std::vector<const Trk::MeasurementBase*>& tS) const {
  // TSOS from the track
  const Trk::TrackStates* stsos = tT.trackStateOnSurfaces();
  // fitQuality from track
  auto fq = tT.fitQuality()->uniqueClone();
  // output datavector of TSOS
  auto ntsos = std::make_unique<Trk::TrackStates>();
  // copy track Si states into track
  Trk::TrackStates::const_iterator p_stsos;
  for (p_stsos = stsos->begin(); p_stsos != stsos->end(); ++p_stsos) {
    ntsos->push_back((*p_stsos)->clone());
  }
  // loop over TRT track extension
  for (auto & it : tS) {
    std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes> typePattern;
    typePattern.set(Trk::TrackStateOnSurface::Measurement);
    const Trk::TrackStateOnSurface* seg_tsos = new Trk::TrackStateOnSurface(it->uniqueClone(), nullptr, nullptr, typePattern);
    ntsos->push_back(seg_tsos);
  }
  Trk::TrackInfo info;
  info.setPatternRecognitionInfo(Trk::TrackInfo::TRTSeededTrackFinder);
  std::unique_ptr<Trk::Track> newTrack( std::make_unique<Trk::Track>(info, std::move(ntsos), std::move(fq)) );
  //Careful refitting at the end
  if (m_doRefit) {
    newTrack = (m_fitterTool->fit(Gaudi::Hive::currentContext(),*newTrack, false, Trk::pion) ) ;
    if (!newTrack) {
      ATH_MSG_DEBUG ("Refit of TRT+Si track segment failed!");
      return nullptr;
    }
    //Protect for tracks that have no really defined locz and theta parameters
    const Trk::Perigee* perTrack=newTrack->perigeeParameters();
    if (perTrack) {
      const AmgSymMatrix(5)* CM = perTrack->covariance();
      if (!CM || std::sqrt((*CM)(1,1)) == 0. || std::sqrt((*CM)(3,3)) == 0.) {
        return nullptr;
      }
    }
  }
 
  return newTrack.release();
}

// Analysis of tracks
 
void 
InDet::TRT_SeededTrackFinder::Analyze(TrackCollection* tC) const{
  if(msgLvl(MSG::DEBUG)) {
    ATH_MSG_DEBUG( "Analyzing tracks..." );
    ATH_MSG_DEBUG( "Number of back tracks " << (tC->size()) );
    int nsct1{}, nsct2{}, nsct3{}, nsct4{}; //SCT layer counters
    int nsctTot1{}, nsctTot2{}, nsctTot3{}, nsctTot4{}; //SCT layer counters
    int npix1{}, npix2{}, npix3{}; //Pixel layer counters
    int npixTot1{}, npixTot2{}, npixTot3{}; //Pixel layer counters
    TrackCollection::const_iterator r    = tC->begin();
    TrackCollection::const_iterator re = tC->end();
    for (; r != re ; ++r){
      nsct1=nsct2=nsct3=nsct4=0; npix1=npix2=npix3=0; 
      const Trk::TrackStates* newtsos = (*r)->trackStateOnSurfaces();
      if(!newtsos) continue;
      Trk::TrackStates::const_iterator itp, itpe=newtsos->end();
      for(itp=newtsos->begin(); itp!=itpe; ++itp){
        const InDet::SiClusterOnTrack* clus = dynamic_cast<const InDet::SiClusterOnTrack*>((*itp)->measurementOnTrack());
        if(clus && ((*itp)->type(Trk::TrackStateOnSurface::Measurement))){  //Count the number of hits used in the track
          double rc = clus->globalPosition().perp();
          if((40.<=rc)&&(rc<80.)){npix1++;}    //1st pixel layer
          if((80.<=rc)&&(rc<100.)){npix2++;}   //2nd pixel layer
          if((100.<=rc)&&(rc<150.)){npix3++;}  //3rd pixel layer
          if((280.<=rc)&&(rc<340.)){nsct1++;}  //1st SCT layer
          if((340.<=rc)&&(rc<390.)){nsct2++;}  //2nd SCT layer
          if((390.<=rc)&&(rc<460.)){nsct3++;}  //3rd SCT layer
          if((460.<=rc)&&(rc<550.)){nsct4++;}  //4th SCT layer
        }
      }
      nsctTot1+=nsct1; nsctTot2+=nsct2; nsctTot3+=nsct3; nsctTot4+=nsct4;
      npixTot1+=npix1; npixTot2+=npix2; npixTot3+=npix3;
    }
    ATH_MSG_DEBUG("Total hits on 1st SCT: "<<nsctTot1<<" 2nd SCT: "<<nsctTot2<<" 3rd SCT: "<<nsctTot3<<" 4th SCT: "<<nsctTot4);
    ATH_MSG_DEBUG("Total hits on 1st Pixel: "<<npixTot1<<" 2nd Pixel: "<<npixTot2<<" 3rd Pixel: "<<npixTot3);
  }
}