SiSpacePointsSeedMaker_ATLxk.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
     
//   Implementation file for class SiSpacePointsSeedMaker_ATLxk
// AlgTool used for TRT_DriftCircleOnTrack object production
// Version 1.0 21/04/2004 I.Gavrilenko
 
#include "SiSpacePointsSeedTool_xk/SiSpacePointsSeedMaker_ATLxk.h"
 
#include "InDetPrepRawData/SiCluster.h"
#include "InDetReadoutGeometry/SiDetectorElement.h"
 
#include "TrkTrack/Track.h"
#include "TrkParameters/TrackParameters.h"
 
#include "StoreGate/ReadCondHandle.h"
#include "CxxUtils/checker_macros.h"
#include "CxxUtils/trapping_fp.h"
#include <iomanip>
#include <ostream>

 
InDet::SiSpacePointsSeedMaker_ATLxk::SiSpacePointsSeedMaker_ATLxk
(const std::string& t,const std::string& n,const IInterface* p)
  : base_class(t, n, p),
    m_thistSvc("THistSvc", n),
    m_outputTree(nullptr),
    m_treeName(""),
    m_treeFolder("/valNtuples/")
{
}

 
StatusCode InDet::SiSpacePointsSeedMaker_ATLxk::initialize()
{
  StatusCode sc = AlgTool::initialize();
 
  ATH_CHECK(m_spacepointsPixel.initialize(m_pixel));
  ATH_CHECK(m_spacepointsSCT.initialize(m_sct));
  ATH_CHECK(m_spacepointsOverlap.initialize(m_useOverlap));

  if (not m_beamSpotKey.empty()) {
    ATH_CHECK(m_beamSpotKey.initialize());
  }
 
  ATH_CHECK( m_fieldCondObjInputKey.initialize() );

  ATH_CHECK( m_prdToTrackMap.initialize( !m_prdToTrackMap.key().empty()));
  
  buildFrameWork();
 
  if (msgLvl(MSG::DEBUG)) { 
    EventData data;
    initializeEventData(data);
    data.nprint=0;
    dump(data, msg(MSG::DEBUG));
  }
 
  if (m_writeNtuple) {
 
    ATH_CHECK( m_thistSvc.retrieve() );
 
    m_treeName = (std::string("SeedTree_")+name());
    std::replace( m_treeName.begin(), m_treeName.end(), '.', '_' );
 
    m_outputTree = new TTree( m_treeName.c_str() , "SeedMakerValTool"); 
 
    m_outputTree->Branch("eventNumber",    &m_eventNumber); 
    m_outputTree->Branch("d0",             &m_d0);
    m_outputTree->Branch("z0",             &m_z0);
    m_outputTree->Branch("pt",             &m_pt);
    m_outputTree->Branch("eta",            &m_eta);
    m_outputTree->Branch("x1",             &m_x1);
    m_outputTree->Branch("x2",             &m_x2);
    m_outputTree->Branch("x3",             &m_x3);
    m_outputTree->Branch("y1",             &m_y1);
    m_outputTree->Branch("y2",             &m_y2);
    m_outputTree->Branch("y3",             &m_y3);
    m_outputTree->Branch("z1",             &m_z1);
    m_outputTree->Branch("z2",             &m_z2);
    m_outputTree->Branch("z3",             &m_z3);
    m_outputTree->Branch("r1",             &m_r1);
    m_outputTree->Branch("r2",             &m_r2);
    m_outputTree->Branch("r3",             &m_r3);
    m_outputTree->Branch("quality",        &m_quality);
    m_outputTree->Branch("seedType",       &m_type);
    m_outputTree->Branch("givesTrack",     &m_givesTrack);
    m_outputTree->Branch("dzdr_b",         &m_dzdr_b);
    m_outputTree->Branch("dzdr_t",         &m_dzdr_t);
    m_outputTree->Branch("track_pt",       &m_trackPt);
    m_outputTree->Branch("track_eta",      &m_trackEta);
 
    TString fullTreeName = m_treeFolder + m_treeName;
 
    ATH_CHECK(  m_thistSvc->regTree( fullTreeName.Data(), m_outputTree )  );
 
  }
 
  return sc;
}

 
StatusCode InDet::SiSpacePointsSeedMaker_ATLxk::finalize()
{
  return AlgTool::finalize();
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::newEvent(const EventContext& ctx, EventData& data, int iteration) const
{
  if (not data.initialized) initializeEventData(data);
 
  data.trigger = false;
  if (!m_pixel && !m_sct) return;

  data.iteration = iteration;
  if (iteration <=0) data.iteration = 0;
  erase(data);
  data.dzdrmin = m_dzdrmin0;
  data.dzdrmax = m_dzdrmax0;
  data.maxScore = m_maxScore;     

  if (data.iteration == 0) {
    if (not m_beamSpotKey.empty()) {
      buildBeamFrameWork(data);
    }

    double magField[3]{0,0,0};
    double globalPos[3] ={10.,10.,0.};
 
    MagField::AtlasFieldCache    fieldCache;
    SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
    const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
    if (fieldCondObj == nullptr) {
      ATH_MSG_ERROR("SiSpacePointsSeedMaker_ATLxk: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
      return;
    }
    fieldCondObj->getInitializedCache (fieldCache);
 
    if (fieldCache.solenoidOn()) {
      fieldCache.getFieldZR(globalPos,magField);
      data.K = 2.f/(300.f*float(magField[2]));
    } else {
      data.K = 2.f/(300.f* 5.f );
    } 
    data.ipt2K = m_ipt2/(data.K*data.K);
    data.ipt2C = m_ipt2*m_COF;
    data.COFK = m_COF*(data.K*data.K);

    data.i_spforseed = data.l_spforseed.begin();
    // Set the seed multiplicity strategy of the event data to the one configured
    // by the user for strip seeds
    data.maxSeedsPerSP = m_maxOneSizeSSS;
    data.keepAllConfirmedSeeds = m_alwaysKeepConfirmedStripSeeds;
    
  } 
  else {  
    data.r_first = 0;     
    // Set the seed multiplicity strategy of the event data to the one configured
    // by the user for pixel seeds
    data.maxSeedsPerSP = m_maxOneSizePPP;
    data.keepAllConfirmedSeeds = m_alwaysKeepConfirmedPixelSeeds;

    fillLists(data);
    return;
  }

 
  data.checketa = data.dzdrmin > 1.;

  float oneOverBinSizeR = 1.f/m_binSizeR;
  int   maxBinR  = m_nBinsR-1;
 
  for (int i=0; i<data.nr; ++i) {
    int n = data.r_index[i];
    data.r_map[n] = 0;
    data.r_Sorted[n].clear();
  }
  data.ns = data.nr = 0;

  SG::ReadHandle<Trk::PRDtoTrackMap>  prd_to_track_map;
  const Trk::PRDtoTrackMap *prd_to_track_map_cptr = nullptr;
  if (!m_prdToTrackMap.key().empty()) {
    prd_to_track_map=SG::ReadHandle<Trk::PRDtoTrackMap>(m_prdToTrackMap, ctx);
    if (!prd_to_track_map.isValid()) {
      ATH_MSG_ERROR("Failed to read PRD to track association map: " << m_prdToTrackMap.key());
    }
    prd_to_track_map_cptr = prd_to_track_map.cptr();
  }


  data.r_first = 0;

  if (m_pixel) {
 
    SG::ReadHandle<SpacePointContainer> spacepointsPixel{m_spacepointsPixel, ctx};
    if (spacepointsPixel.isValid()) {
      for (const SpacePointCollection* spc: *spacepointsPixel) {
        for (const Trk::SpacePoint* sp: *spc) {
          
          if ((prd_to_track_map_cptr &&  isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax) continue;

          const InDetDD::SiDetectorElement* de= 
            static_cast<const InDetDD::SiDetectorElement*>(sp->clusterList().first->detectorElement());
          if (!de || de->isDBM()) continue;
 
          InDet::SiSpacePointForSeed* sps = newSpacePoint(data, sp);
          if (!sps) continue;

          int radiusBin = static_cast<int>(sps->radius()*oneOverBinSizeR);
          if (radiusBin>maxBinR) radiusBin = maxBinR;

          data.r_Sorted[radiusBin].push_back(sps);
          ++data.r_map[radiusBin];
          if (data.r_map[radiusBin]==1) data.r_index[data.nr++] = radiusBin;
          if (radiusBin > data.r_first) data.r_first = radiusBin;
          ++data.ns;
        } 
      }
    } 
    ++data.r_first; 
  } 
  
  if (m_sct) {
 
    SG::ReadHandle<SpacePointContainer> spacepointsSCT{m_spacepointsSCT, ctx};
    if (spacepointsSCT.isValid()) {
 
      for (const SpacePointCollection* spc: *spacepointsSCT) {
        for (const Trk::SpacePoint* sp: *spc) {
          if ((prd_to_track_map_cptr &&  isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax) continue;

          InDet::SiSpacePointForSeed* sps = newSpacePoint(data, sp);
          if (!sps) continue;

          int radiusBin = static_cast<int>(sps->radius()*oneOverBinSizeR);
          if (radiusBin>maxBinR) radiusBin = maxBinR;
          data.r_Sorted[radiusBin].push_back(sps);
          ++data.r_map[radiusBin];
          if (data.r_map[radiusBin]==1) data.r_index[data.nr++] = radiusBin;
          ++data.ns;
        }
      }
    }

    if (m_useOverlap && !data.checketa) {
 
      SG::ReadHandle<SpacePointOverlapCollection> spacepointsOverlap{m_spacepointsOverlap, ctx};
      if (spacepointsOverlap.isValid()) {
  
        for (const Trk::SpacePoint* sp: *spacepointsOverlap) {
          if ((prd_to_track_map_cptr &&  isUsed(sp, *prd_to_track_map_cptr)) || sp->r() > m_r_rmax) continue;

          InDet::SiSpacePointForSeed* sps = newSpacePoint(data, sp);
          if (!sps) continue;

          int radiusBin = static_cast<int>(sps->radius()*oneOverBinSizeR);
          if (radiusBin>maxBinR) radiusBin = maxBinR;
          data.r_Sorted[radiusBin].push_back(sps);
          ++data.r_map[radiusBin];
          if (data.r_map[radiusBin]==1) data.r_index[data.nr++] = radiusBin;
          ++data.ns;
        }
      }
    }
  }
  
  if (iteration < 0) data.r_first = 0;

  fillLists(data);
}


void InDet::SiSpacePointsSeedMaker_ATLxk::newRegion
(const EventContext& ctx, EventData& data,
 const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT) const
{
  if (not data.initialized) initializeEventData(data);
 
  data.iteration  = 0;
  data.trigger = false;
  erase(data);
  if (!m_pixel && !m_sct) return;
 
  data.dzdrmin = m_dzdrmin0;
  data.dzdrmax = m_dzdrmax0;
  data.maxScore = m_maxScore;
 
  if (not m_beamSpotKey.empty()) {
    buildBeamFrameWork(data);
  }
 
  double magField[3]{0,0,0};
  double globalPos[3] ={10.,10.,0.};
 
  MagField::AtlasFieldCache    fieldCache;
  // Get field cache object
  SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
  const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
  
  if (fieldCondObj == nullptr) {
    ATH_MSG_ERROR("SiSpacePointsSeedMaker_ATLxk: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
    return;
  }
  fieldCondObj->getInitializedCache (fieldCache);
 
  // initialise pt cuts and conversion factors 
  if (fieldCache.solenoidOn()) {
    fieldCache.getFieldZR(globalPos,magField);
 
    data.K = 2.f/(300.f*float(magField[2]));
  } else {
    data.K = 2.f/(300.f* 5.f );
  }
 
  data.ipt2K = m_ipt2/(data.K*data.K);
  data.ipt2C = m_ipt2*m_COF;
  data.COFK = m_COF*(data.K*data.K);
 
  data.i_spforseed = data.l_spforseed.begin();
 
  float oneOverBinSizeR = 1.f/m_binSizeR;
  int   maxBinR  = m_nBinsR-1;
 
  data.r_first = 0;
  data.checketa = false;
 
  for (int i=0; i<data.nr; ++i) {
    int n = data.r_index[i];
    data.r_map[n] = 0;
    data.r_Sorted[n].clear();
  }
  data.ns = data.nr = 0;
 
  // Get pixels space points containers from store gate 
  //
  if (m_pixel && !vPixel.empty()) {
 
    SG::ReadHandle<SpacePointContainer> spacepointsPixel{m_spacepointsPixel, ctx};
    if ( spacepointsPixel.isValid() ) {
 
      // Loop through all trigger collections
      //
      for (const IdentifierHash& l: vPixel) {
        const auto *w = spacepointsPixel->indexFindPtr(l);
        if (w==nullptr) continue;
        for (const Trk::SpacePoint* sp: *w) {
          float r = sp->r();
          if (r > m_r_rmax) continue;
          InDet::SiSpacePointForSeed* sps = newSpacePoint(data, sp);
          int ir = static_cast<int>(sps->radius()*oneOverBinSizeR);
          if (ir>maxBinR) ir = maxBinR;
          data.r_Sorted[ir].push_back(sps);
          ++data.r_map[ir];
          if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
          ++data.ns;
        }
      }
    }
  }
 
  // Get sct space points containers from store gate 
  //
  if (m_sct && !vSCT.empty()) {
 
    SG::ReadHandle<SpacePointContainer> spacepointsSCT{m_spacepointsSCT, ctx};
    if (spacepointsSCT.isValid()) {
 
      // Loop through all trigger collections
      //
      for (const IdentifierHash& l: vSCT) {
        const auto *w = spacepointsSCT->indexFindPtr(l);
        if (w==nullptr) continue;
        for (const Trk::SpacePoint* sp: *w) {
          float r = sp->r();
          if (r > m_r_rmax) continue;
          InDet::SiSpacePointForSeed* sps = newSpacePoint(data, sp);
          int ir = static_cast<int>(sps->radius()*oneOverBinSizeR);
          if (ir>maxBinR) ir = maxBinR;
          data.r_Sorted[ir].push_back(sps);
          ++data.r_map[ir];
          if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
          ++data.ns;
        }
      }
    }
  }
  fillLists(data);
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::newRegion
(const EventContext& ctx, EventData& data,
 const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT,
 const IRoiDescriptor& IRD) const
{
  constexpr float twoPi = 2.*M_PI;
 
  if (not data.initialized) initializeEventData(data);
 
  newRegion(ctx, data, vPixel, vSCT);
  data.trigger = true;
 
  double dzdrmin = 1./std::tan(2.*std::atan(std::exp(-IRD.etaMinus())));
  double dzdrmax = 1./std::tan(2.*std::atan(std::exp(-IRD.etaPlus ())));
 
  data.zminB = IRD.zedMinus()-data.zbeam[0]; // min bottom Z
  data.zmaxB = IRD.zedPlus ()-data.zbeam[0]; // max bottom Z
  data.zminU = data.zminB+550.f*float(dzdrmin);
  data.zmaxU = data.zmaxB+550.f*float(dzdrmax);
  float fmax    = IRD.phiPlus ();
  float fmin    = IRD.phiMinus();
  if (fmin > fmax) fmin -= twoPi;
  data.ftrig  = (fmin+fmax)*.5f;
  data.ftrigW = (fmax-fmin)*.5f;
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::find2Sp(EventData& data, const std::list<Trk::Vertex>& lv) const
{
  if (not data.initialized) initializeEventData(data);
 
  data.zminU = m_zmin;
  data.zmaxU = m_zmax;
 
  int mode = 0;
  if (lv.begin()!=lv.end()) mode = 1;
  bool newv = newVertices(data, lv);
  
  if (newv || !data.state || data.nspoint!=2 || data.mode!=mode || data.nlist) {
    data.i_seede_Pro = data.l_seeds_Pro.begin();
    data.state = 1;
    data.nspoint = 2;
    data.nlist = 0;
    data.mode = mode;
    data.endlist = true;
    data.fvNmin = 0;
    data.fNmin = 0;
    data.zMin = 0;
    production2Sp(data);
  }
  data.i_seed_Pro = data.l_seeds_Pro.begin();
  
 
  if (msgLvl(MSG::DEBUG)) {
    data.nprint=1;
    dump(data, msg(MSG::DEBUG));
  }
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
{
  if (not data.initialized) initializeEventData(data);

  data.zminU = m_zmin;
  data.zmaxU = m_zmax;
  int mode = 2;
  if (lv.begin()!=lv.end()) mode = 3;
  bool newv = newVertices(data, lv);
  if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
    data.i_seede_Pro = data.l_seeds_Pro.begin();
    data.state = 1;
    data.nspoint = 3;
    data.nlist = 0;
    data.mode = mode;
    data.endlist = true;
    data.fvNmin = 0;
    data.fNmin = 0;
    data.zMin = 0;
    production3Sp(data); 
  }
  data.i_seed_Pro = data.l_seeds_Pro.begin();
 
  
  if (msgLvl(MSG::DEBUG)) {
    data.nprint=1;
    dump(data, msg(MSG::DEBUG));
  }
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv, const double* ZVertex) const
{
  if (not data.initialized) initializeEventData(data);

  data.zminU = ZVertex[0];
  if (data.zminU < m_zmin) data.zminU = m_zmin; 
  data.zmaxU = ZVertex[1];
  if (data.zmaxU > m_zmax) data.zmaxU = m_zmax; 

  int mode = 2;
  if (lv.begin()!=lv.end()) mode = 3;
  bool newv = newVertices(data, lv);
  if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
    data.i_seede_Pro = data.l_seeds_Pro.begin();
    data.state = 1;
    data.nspoint = 3;
    data.nlist = 0;
    data.mode = mode;
    data.endlist = true;
    data.fvNmin = 0;
    data.fNmin = 0;
    data.zMin = 0;
    production3Sp(data); 
  }
  data.i_seed_Pro = data.l_seeds_Pro.begin();
 
  if (msgLvl(MSG::DEBUG)) {  
    data.nprint=1;
    dump(data, msg(MSG::DEBUG));
  }
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::findVSp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
{
  if (not data.initialized) initializeEventData(data);
 
  data.zminU = m_zmin;
  data.zmaxU = m_zmax;
 
  int mode = 5;
  if (lv.begin()!=lv.end()) mode = 6;
  bool newv = newVertices(data, lv);
  
  if (newv || !data.state || data.nspoint!=4 || data.mode!=mode || data.nlist) {
    data.i_seede_Pro = data.l_seeds_Pro.begin();
    data.state = 1;
    data.nspoint = 4;
    data.nlist = 0;
    data.mode = mode;
    data.endlist = true;
    data.fvNmin = 0;
    data.fNmin = 0;
    data.zMin = 0;
    production3Sp(data);
  }
  data.i_seed_Pro = data.l_seeds_Pro.begin();
 
  if (msgLvl(MSG::DEBUG)) { 
    data.nprint=1;
    dump(data, msg(MSG::DEBUG));
  }
}

 
MsgStream& InDet::SiSpacePointsSeedMaker_ATLxk::dump(EventData& data, MsgStream& out) const
{
  if (not data.initialized) initializeEventData(data);
 
  if (data.nprint) return dumpEvent(data, out);
  return dumpConditions(data, out);
}

 
MsgStream& InDet::SiSpacePointsSeedMaker_ATLxk::dumpConditions(EventData &data, MsgStream& out) const
{
  int n = 42-m_spacepointsPixel.key().size();
  std::string s2;
  for (int i=0; i<n; ++i) s2.append(" ");
  s2.append("|");
  n     = 42-m_spacepointsSCT.key().size();
  std::string s3;
  for (int i=0; i<n; ++i) s3.append(" ");
  s3.append("|");
  n     = 42-m_spacepointsOverlap.key().size();
  std::string s4;
  for (int i=0; i<n; ++i) s4.append(" ");
  s4.append("|");
  n     = 42-m_beamSpotKey.key().size();
  std::string s5;
  for (int i=0; i<n; ++i) s5.append(" ");
  s5.append("|");
 
  out<<"|---------------------------------------------------------------------|"
     <<endmsg;
  out<<"| Pixel    space points   | "<<m_spacepointsPixel.key() <<s2
     <<endmsg;
  out<<"| SCT      space points   | "<<m_spacepointsSCT.key()<<s3
     <<endmsg;
  out<<"| Overlap  space points   | "<<m_spacepointsOverlap.key()<<s4
     <<endmsg;
  out<<"| BeamConditionsService   | "<<m_beamSpotKey.key()<<s5
     <<endmsg;
  out<<"| usePixel                | "
     <<std::setw(12)<<m_pixel 
     <<"                              |"<<endmsg;
  out<<"| useSCT                  | "
     <<std::setw(12)<<m_sct 
     <<"                              |"<<endmsg;
  out<<"| maxSize                 | "
     <<std::setw(12)<<m_maxsize 
     <<"                              |"<<endmsg;
  out<<"| maxSizeSP               | "
     <<std::setw(12)<<m_maxsizeSP
     <<"                              |"<<endmsg;
  out<<"| pTmin  (mev)            | "
     <<std::setw(12)<<std::setprecision(5)<<m_ptmin
     <<"                              |"<<endmsg;
  out<<"| max radius SP           | "
     <<std::setw(12)<<std::setprecision(5)<<m_r_rmax 
     <<"                              |"<<endmsg;
  out<<"| radius step             | "
     <<std::setw(12)<<std::setprecision(5)<<m_binSizeR
     <<"                              |"<<endmsg;
  out<<"| min Z-vertex position   | "
     <<std::setw(12)<<std::setprecision(5)<<m_zmin
     <<"                              |"<<endmsg;
  out<<"| max Z-vertex position   | "
     <<std::setw(12)<<std::setprecision(5)<<m_zmax
     <<"                              |"<<endmsg;
  out<<"| min radius first  SP(2) | "
     <<std::setw(12)<<std::setprecision(5)<<m_r1minv
     <<"                              |"<<endmsg;
  out<<"| min radius second SP(2) | "
     <<std::setw(12)<<std::setprecision(5)<<m_r2minv
     <<"                              |"<<endmsg;
  out<<"| max radius first  SP(2) | "
     <<std::setw(12)<<std::setprecision(5)<<m_r1maxv
     <<"                              |"<<endmsg;
  out<<"| max radius second SP(2) | "
     <<std::setw(12)<<std::setprecision(5)<<m_r2maxv
     <<"                              |"<<endmsg;
  out<<"| min space points dR     | "
     <<std::setw(12)<<std::setprecision(5)<<m_drmin
     <<"                              |"<<endmsg;
  out<<"| max space points dR     | "
     <<std::setw(12)<<std::setprecision(5)<<m_drmax
     <<"                              |"<<endmsg;
  out<<"| max dZ    impact        | "
     <<std::setw(12)<<std::setprecision(5)<<m_dzver 
     <<"                              |"<<endmsg;
  out<<"| max dZ/dR impact        | "
     <<std::setw(12)<<std::setprecision(5)<<m_dzdrver 
     <<"                              |"<<endmsg;
  out<<"| max       impact        | "
     <<std::setw(12)<<std::setprecision(5)<<m_maxdImpact
     <<"                              |"<<endmsg;
  out<<"| max       impact sss    | "
     <<std::setw(12)<<std::setprecision(5)<<m_maxdImpactSSS
     <<"                              |"<<endmsg;
  out<<"|---------------------------------------------------------------------|"
     <<endmsg;
  out<<"| Beam X center           | "
     <<std::setw(12)<<std::setprecision(5)<<data.xbeam[0]
     <<"                              |"<<endmsg;
  out<<"| Beam Y center           | "
     <<std::setw(12)<<std::setprecision(5)<<data.ybeam[0]
     <<"                              |"<<endmsg;
  out<<"| Beam Z center           | "
     <<std::setw(12)<<std::setprecision(5)<<data.zbeam[0]
     <<"                              |"<<endmsg;
  out<<"| Beam X-axis direction   | "
     <<std::setw(12)<<std::setprecision(5)<<data.xbeam[1]
     <<std::setw(12)<<std::setprecision(5)<<data.xbeam[2]
     <<std::setw(12)<<std::setprecision(5)<<data.xbeam[3]
     <<"      |"<<endmsg;
  out<<"| Beam Y-axis direction   | "
     <<std::setw(12)<<std::setprecision(5)<<data.ybeam[1]
     <<std::setw(12)<<std::setprecision(5)<<data.ybeam[2]
     <<std::setw(12)<<std::setprecision(5)<<data.ybeam[3]
     <<"      |"<<endmsg;
  out<<"| Beam Z-axis direction   | "
     <<std::setw(12)<<std::setprecision(5)<<data.zbeam[1]
     <<std::setw(12)<<std::setprecision(5)<<data.zbeam[2]
     <<std::setw(12)<<std::setprecision(5)<<data.zbeam[3]
     <<"      |"<<endmsg;
  out<<"|---------------------------------------------------------------------|"
     <<endmsg;
  return out;
}

 
MsgStream& InDet::SiSpacePointsSeedMaker_ATLxk::dumpEvent(EventData& data, MsgStream& out) 
{
  out<<"|---------------------------------------------------------------------|"
     <<endmsg;
  out<<"| ns                    | "
     <<std::setw(12)<<data.ns
     <<"                              |"<<endmsg;
  out<<"| nsaz                  | "
     <<std::setw(12)<<data.nsaz
     <<"                              |"<<endmsg;
  out<<"| nsazv                 | "
     <<std::setw(12)<<data.nsazv
     <<"                              |"<<endmsg;
  out<<"| seeds                   | "
     <<std::setw(12)<<data.l_seeds_Pro.size()
     <<"                              |"<<endmsg;
  out<<"|---------------------------------------------------------------------|"
     <<endmsg;
  return out;
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::findNext(EventData& data) const
{
  if (data.endlist) return;
 
  data.i_seede_Pro = data.l_seeds_Pro.begin();
 
  if      (data.mode==0 or data.mode==1) production2Sp(data);
  else if (data.mode==2 or data.mode==3) production3Sp(data);
  else if (data.mode==5 or data.mode==6) production3Sp(data);
 
  data.i_seed_Pro = data.l_seeds_Pro.begin();
  ++data.nlist;
}                       

 
bool InDet::SiSpacePointsSeedMaker_ATLxk::newVertices(EventData& data, const std::list<Trk::Vertex>& lV) const
{
  unsigned int s1 = data.l_vertex.size();
  unsigned int s2 = lV.size();

  data.isvertex = false;
  if (s1==0 && s2==0) return false;

  data.l_vertex.clear();
  if (s2 == 0) return false;

  data.isvertex = true;
  for (const Trk::Vertex& v: lV) {
    data.l_vertex.insert(static_cast<float>(v.position().z()));
  }

  data.zminU = (*data.l_vertex. begin())-20.f;
  if (data.zminU < m_zmin) data.zminU = m_zmin;
  data.zmaxU = (*data.l_vertex.rbegin())+20.f;
  if (data.zmaxU > m_zmax) data.zmaxU = m_zmax;

  return false;
}

// Initiate frame work for seed generator
 
void InDet::SiSpacePointsSeedMaker_ATLxk::buildFrameWork() 
{
  m_ptmin = std::max( std::abs(m_ptmin), float(100.*m_fieldScale));
  if (m_maxdImpactSSS < m_maxdImpact   ) m_maxdImpactSSS = m_maxdImpact;
  if (m_maxdImpactDecays < m_maxdImpactSSS) m_maxdImpactDecays = m_maxdImpactSSS;

  if (std::abs(m_etamin) < .1) m_etamin = -m_etamax;
  m_dzdrmax0  = 1.f/std::tan(2.f*std::atan(std::exp(-m_etamax)));
  m_dzdrmin0  = 1.f/std::tan(2.f*std::atan(std::exp(-m_etamin)));
  
  m_ipt       = 1.f/std::abs(.9f*m_ptmin);
  m_ipt2      = m_ipt*m_ipt;

  m_seedScoreThresholdPPPConfirmationSeed = m_seedScoreBonusPPP - 1.f; 
  m_seedScoreThresholdSSSConfirmationSeed = m_seedScoreBonusSSS - 1.f; 

  m_nBinsR = static_cast<int>((m_r_rmax+.1f)/m_binSizeR);


  constexpr float twoPi = 2.f*M_PI;

  const int   nPhiBinsMax = arraySizePhi;
  const float inverseSizePhiMax = static_cast<float>(nPhiBinsMax)/twoPi; 

  constexpr float inverseSizePhiMin = 100./60.;     

  
  if (m_optimisePhiBinning){
    constexpr float radiusPixelStart = 33.; 
    constexpr float radiusPixelEnd = 150.;  
    const float binSizePhi_PPP = m_pixel ? azimuthalStep(m_ptmin/m_fieldScale,m_maxdImpact,radiusPixelStart,radiusPixelEnd)/3.f : 0.f; 
    constexpr float radiusSctStart = 295.; ; 
    constexpr float radiusSctEnd = 560.; 
    const float binSizePhi_SSS = m_sct ? azimuthalStep(m_ptmin/m_fieldScale,m_maxdImpactSSS,radiusSctStart,radiusSctEnd)/3.f : 0.f; 
    m_inverseBinSizePhi = 1.f/std::max(binSizePhi_PPP, binSizePhi_SSS); 
  }
  else {
    float ptm = 400.; 
    if (m_ptmin/m_fieldScale < ptm) ptm = m_ptmin/m_fieldScale;
    m_inverseBinSizePhi = ptm /60.f;
  }

  if      (m_inverseBinSizePhi > inverseSizePhiMax) m_inverseBinSizePhi = inverseSizePhiMax;
  else if (m_inverseBinSizePhi < inverseSizePhiMin) m_inverseBinSizePhi = inverseSizePhiMin;

  m_maxPhiBin = static_cast<int>(twoPi*m_inverseBinSizePhi);
  if (m_maxPhiBin >=nPhiBinsMax) m_maxPhiBin = nPhiBinsMax-1;

  const int   nPhiBinsVertexMax  = arraySizePhiV;
  const float inverseBinSizePhiVertexMax = static_cast<float>(nPhiBinsVertexMax)/twoPi;
  m_inverseBinSizePhiVertex = m_ptmin/m_fieldScale/120.f;
  if (m_inverseBinSizePhiVertex > inverseBinSizePhiVertexMax) m_inverseBinSizePhiVertex = inverseBinSizePhiVertexMax;
  m_maxBinPhiVertex = static_cast<int>(twoPi*m_inverseBinSizePhiVertex);
  if (m_maxBinPhiVertex>=nPhiBinsVertexMax) m_maxBinPhiVertex = nPhiBinsVertexMax-1;
 

 
  for (int phiBin=0; phiBin<=m_maxPhiBin; ++phiBin) {
 
    int phiBelow = phiBin-1;
    if (phiBelow<0) phiBelow=m_maxPhiBin; 
 
    int phiAbove = phiBin+1;
    if (phiAbove>m_maxPhiBin) phiAbove=0; 
    
    for (int z=0; z<arraySizeZ; ++z) {

 
      int twoDbinSamePhi        = phiBin * arraySizeZ+z;
      int twoDbinLowerPhi        = phiBelow*arraySizeZ+z;
      int twoDbinHigherPhi        = phiAbove*  arraySizeZ+z;
 
      m_nNeighbourCellsBottom [twoDbinSamePhi]    = 3; 
      m_nNeighbourCellsTop [twoDbinSamePhi]    = 3;
 
      m_neighbourCellsBottom[twoDbinSamePhi][0] = twoDbinSamePhi; 
      m_neighbourCellsTop[twoDbinSamePhi][0] = twoDbinSamePhi;
 
      m_neighbourCellsBottom[twoDbinSamePhi][1] = twoDbinLowerPhi; 
      m_neighbourCellsTop[twoDbinSamePhi][1] = twoDbinLowerPhi;
 
      m_neighbourCellsBottom[twoDbinSamePhi][2] = twoDbinHigherPhi; 
      m_neighbourCellsTop[twoDbinSamePhi][2] = twoDbinHigherPhi;

      if (z==5) {
        m_nNeighbourCellsTop [twoDbinSamePhi]    = 9;
        // in the central z region, we include the two neighbouring 
        // z slices for the top neighbour search   
 
        m_neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi+1;
        m_neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
        m_neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
        m_neighbourCellsTop[twoDbinSamePhi][6] = twoDbinSamePhi-1;
        m_neighbourCellsTop[twoDbinSamePhi][7] = twoDbinLowerPhi-1;
        m_neighbourCellsTop[twoDbinSamePhi][8] = twoDbinHigherPhi-1;
      } 
      // z > 5: positive z values, |z| > 250mm 
      else if (z> 5) {
        // for the bottom SP search in positive non-central z, we include the 
        // neighbouring Z region on the left (towards the IP) in the bottom 
        // neighbour search 
        m_nNeighbourCellsBottom [twoDbinSamePhi]    = 6;
        m_neighbourCellsBottom[twoDbinSamePhi][3] = twoDbinSamePhi-1;
        m_neighbourCellsBottom[twoDbinSamePhi][4] = twoDbinLowerPhi-1;
        m_neighbourCellsBottom[twoDbinSamePhi][5] = twoDbinHigherPhi-1;
 
        if (z<10) {
          m_nNeighbourCellsTop [twoDbinSamePhi]    = 6;
          m_neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi+1;
          m_neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
          m_neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
        }
      } 
      // z < 5: negative z values, |z| > 250mm 
      else {
        m_nNeighbourCellsBottom [twoDbinSamePhi]    = 6;
        m_neighbourCellsBottom[twoDbinSamePhi][3] = twoDbinSamePhi+1;
        m_neighbourCellsBottom[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
        m_neighbourCellsBottom[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
 
        if (z>0) {
          // if there is a z region on the left (away from the IP), we include it in the top
          // neighbour search 
          m_nNeighbourCellsTop [twoDbinSamePhi]    = 6;
          m_neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi-1;
          m_neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi-1;
          m_neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi-1;
        }
      }

      if (z==3) {
        m_nNeighbourCellsBottom[twoDbinSamePhi]      = 9;
        m_neighbourCellsBottom[twoDbinSamePhi][6] = twoDbinSamePhi+2;
        m_neighbourCellsBottom[twoDbinSamePhi][7] = twoDbinLowerPhi+2;
        m_neighbourCellsBottom[twoDbinSamePhi][8] = twoDbinHigherPhi+2;
      } else if (z==7) {
        m_nNeighbourCellsBottom[twoDbinSamePhi]      = 9;
        m_neighbourCellsBottom[twoDbinSamePhi][6] = twoDbinSamePhi-2;
        m_neighbourCellsBottom[twoDbinSamePhi][7] = twoDbinLowerPhi-2;
        m_neighbourCellsBottom[twoDbinSamePhi][8] = twoDbinHigherPhi-2;
      }
    }
  }

  for (int phiBin=0; phiBin<=m_maxBinPhiVertex; ++phiBin) {
 
    int phiBinBelow = phiBin-1;
    if (phiBinBelow<0) phiBinBelow=m_maxBinPhiVertex;
 
    int phiBinTop = phiBin+1;
    if (phiBinTop>m_maxBinPhiVertex) phiBinTop=0;
    
    for (int zbin=0; zbin<arraySizeZV; ++zbin) {
      
      int twoDbinSamePhi  = phiBin *arraySizeZV+zbin;
      int twoDbinLowerPhi  = phiBinBelow*arraySizeZV+zbin;
      int twoDbinHigherPhi  = phiBinTop*arraySizeZV+zbin;

      m_nNeighboursVertexPhiZ[twoDbinSamePhi]    = 3;
      m_neighboursVertexPhiZ[twoDbinSamePhi][0] = twoDbinSamePhi;
      m_neighboursVertexPhiZ[twoDbinSamePhi][1] = twoDbinLowerPhi;
      m_neighboursVertexPhiZ[twoDbinSamePhi][2] = twoDbinHigherPhi;

      if (zbin>1) {
        m_nNeighboursVertexPhiZ[twoDbinSamePhi]    = 6;
        m_neighboursVertexPhiZ[twoDbinSamePhi][3] = twoDbinSamePhi-1;
        m_neighboursVertexPhiZ[twoDbinSamePhi][4] = twoDbinLowerPhi-1;
        m_neighboursVertexPhiZ[twoDbinSamePhi][5] = twoDbinHigherPhi-1;
      } 
      else if (zbin<1) {
        m_nNeighboursVertexPhiZ[twoDbinSamePhi]    = 6;
        m_neighboursVertexPhiZ[twoDbinSamePhi][3] = twoDbinSamePhi+1;
        m_neighboursVertexPhiZ[twoDbinSamePhi][4] = twoDbinLowerPhi+1;
        m_neighboursVertexPhiZ[twoDbinSamePhi][5] = twoDbinHigherPhi+1;
      }
    }
  }
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::buildBeamFrameWork(EventData& data) const
{ 
  SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey };
 
  const Amg::Vector3D &bsCentre = beamSpotHandle->beamPos();
  double     tx = std::tan(beamSpotHandle->beamTilt(0));
  double     ty = std::tan(beamSpotHandle->beamTilt(1));
 
  double phi   = std::atan2(ty,tx);
  double theta   = std::acos(1./std::sqrt(1.+tx*tx+ty*ty));
  double sinTheta = std::sin(theta);
  double cosTheta = std::cos(theta);
  double sinPhi = std::sin(phi);
  double cosPhi = std::cos(phi);
  
  data.xbeam[0] = static_cast<float>(bsCentre.x());
  data.xbeam[1] = static_cast<float>(cosTheta*cosPhi*cosPhi+sinPhi*sinPhi);
  data.xbeam[2] = static_cast<float>(cosTheta*sinPhi*cosPhi-sinPhi*cosPhi);
  data.xbeam[3] =-static_cast<float>(sinTheta*cosPhi               );
 
  data.ybeam[0] = static_cast<float>(bsCentre.y());
  data.ybeam[1] = static_cast<float>(cosTheta*cosPhi*sinPhi-sinPhi*cosPhi);
  data.ybeam[2] = static_cast<float>(cosTheta*sinPhi*sinPhi+cosPhi*cosPhi);
  data.ybeam[3] =-static_cast<float>(sinTheta*sinPhi               );
  
  data.zbeam[0] = static_cast<float>(bsCentre.z());
  data.zbeam[1] = static_cast<float>(sinTheta*cosPhi);
  data.zbeam[2] = static_cast<float>(sinTheta*sinPhi);
  data.zbeam[3] = static_cast<float>(cosTheta);
}

 
void  InDet::SiSpacePointsSeedMaker_ATLxk::convertToBeamFrameWork
(EventData& data, const Trk::SpacePoint*const& sp, std::array<float,3> & r) 
{
  r[0] = static_cast<float>(sp->globalPosition().x())-data.xbeam[0];
  r[1] = static_cast<float>(sp->globalPosition().y())-data.ybeam[0];
  r[2] = static_cast<float>(sp->globalPosition().z())-data.zbeam[0];
}
   
 
void InDet::SiSpacePointsSeedMaker_ATLxk::fillLists(EventData& data) const
{
  constexpr float twoPi = 2.*M_PI;
 
  int  firstRadialBin =     0;
  bool ibl = false;
 
  
    const std::map<float, int> ztoBin{
      {-2500.    ,0},
      {-1400.    ,1},
      {-925.     ,2},
      {-450.     ,3},
      {-250      ,4},
      { 250      ,5},
      { 450      ,6},
      { 925      ,7},
      { 1400     ,8},
      { 2500     ,9},
      { 100000   ,10},  
    };
 
  for (int radialBin=data.r_first; radialBin<m_nBinsR;  ++radialBin) {

    if (!data.r_map[radialBin]) continue;
    if (firstRadialBin == 0) firstRadialBin = radialBin;

    if (data.iteration) {
      if (!data.r_Sorted[radialBin].front()->spacepoint->clusterList().second) {
        if (radialBin * m_binSizeR < m_radiusCutIBL) ibl = true;
      } 


      else if (ibl) {
        break;
      } 
      else if (radialBin > 175) {
        break;
      }
    }
    // loop over the space points in the r-bin and sort them into the 2d phi-z binning 
    for (InDet::SiSpacePointForSeed* SP : data.r_Sorted[radialBin]) {
      
      float Phi = SP->phi();
      if (Phi<0.) Phi+=twoPi; // phi is defined in [0..2pi] for the binning 
      int phiBin = static_cast<int>(Phi*m_inverseBinSizePhi);
      if (phiBin < 0) {
        phiBin = m_maxPhiBin;
      } else if (phiBin > m_maxPhiBin) {
        phiBin = 0;
      }
 
      float Z = SP->z();
      int zBin{0};
      auto bound = ztoBin.lower_bound(Z); 
      if (bound == ztoBin.end()){
        --bound;  
      }
      zBin=bound->second;

      int twoDbin = phiBin*arraySizeZ+zBin;
      ++data.nsaz;
      // push our space point into the 2D binned array
      data.rfz_Sorted[twoDbin].push_back(SP);
      if (!data.rfz_map[twoDbin]++) data.rfz_index[data.nrfz++] = twoDbin;
    }
  }
  if (!m_sct && firstRadialBin && static_cast<float>(firstRadialBin)*m_binSizeR < m_radiusCutIBL) {
    data.maxScore = m_seedScoreThresholdPPPConfirmationSeed;
  }
 
  data.state = 0;
}
 
 
 
float InDet::SiSpacePointsSeedMaker_ATLxk::azimuthalStep(const float pTmin,const float maxd0,const float Rmin,const float Rmax) 
{
  // Tell clang to optimize assuming that FP exceptions can trap.
  // Otherwise, it can vectorize the division, which can lead to
  // spurious division-by-zero traps from unused vector lanes.
  CXXUTILS_TRAPPING_FP;

  float Rm = pTmin/.6f;

  float worstCaseD0 = maxd0;
  if (maxd0 > Rmin) worstCaseD0 = Rmin; 
 
  float sI = std::abs(std::asin(worstCaseD0/Rmin) - std::asin(worstCaseD0/Rmax)); 
  float sF = std::abs(std::asin(std::min(1.f,Rmax/(2.f*Rm))) -
                      std::asin(std::min(1.f,Rmin/(2.f*Rm))));
  return sI+sF; 
}
 
 
                                                                                          
// Erase space point information
 
void InDet::SiSpacePointsSeedMaker_ATLxk::erase(EventData& data) 
{
  for (int i=0; i<data.nrfz;  ++i) {
    int n = data.rfz_index[i];
    data.rfz_map[n] = 0;
    data.rfz_Sorted[n].clear();
  }
  
  for (int i=0; i<data.nrfzv; ++i) {
    int n = data.rfzv_index[i];
    data.rfzv_map[n] = 0;
    data.rfzv_Sorted[n].clear();
  }
  data.state = 0;
  data.nsaz  = 0;
  data.nsazv = 0;
  data.nrfz  = 0;
  data.nrfzv = 0;
}

// Test is space point used
 
 

// 2 space points seeds production
 
void InDet::SiSpacePointsSeedMaker_ATLxk::production2Sp(EventData& data) const
{
  if (data.nsazv<2) return;
 
  std::vector<InDet::SiSpacePointForSeed*>::iterator r0,r0e,r,re;
  int nseed = 0;
 
  // Loop thorugh all azimuthal regions
  //
  for (int f=data.fvNmin; f<=m_maxBinPhiVertex; ++f) {
 
    // For each azimuthal region loop through Z regions
    //
    int z = 0;
    if (!data.endlist) z = data.zMin;
    for (; z<arraySizeZV; ++z) {
      
      int a = f*arraySizeZV+z;
      if (!data.rfzv_map[a]) continue;
      r0  = data.rfzv_Sorted[a].begin();
      r0e = data.rfzv_Sorted[a].end  ();
 
      if (!data.endlist) {
        r0 = data.rMin;
        data.endlist = true;
      }
 
      // Loop through trigger space points
      //
      for (; r0!=r0e; ++r0) {
 
        float X  = (*r0)->x();
        float Y  = (*r0)->y();
        float R  = (*r0)->radius();
        if (R<m_r2minv) continue;
        if (R>m_r2maxv) break;
        float Z  = (*r0)->z();
        float ax = X/R;
        float ay = Y/R;
 
        // Bottom links production
        //
        int numberBottomCells = m_nNeighboursVertexPhiZ[a];
        for (int i=0; i<numberBottomCells; ++i) {
    
          int an = m_neighboursVertexPhiZ[a][i];
          if (!data.rfzv_map[an]) continue;
 
          r  = data.rfzv_Sorted[an].begin();
          re = data.rfzv_Sorted[an].end();
    
          for (; r!=re; ++r) {
      
            float Rb =(*r)->radius();
            if (Rb<m_r1minv) continue;
            if (Rb>m_r1maxv) break;
            float dR = R-Rb;
            if (dR<m_drminv) break;
            if (dR>m_drmax) continue;
            float dZ = Z-(*r)->z();
            float Tz = dZ/dR;
            if (Tz<data.dzdrmin || Tz>data.dzdrmax) continue;
            float Zo = Z-R*Tz;
 
            // Comparison with vertices Z coordinates
            //
            if (!isZCompatible(data, Zo, Rb, Tz)) continue;
 
            // Momentum cut
            //
            float dx =(*r)->x()-X;
            float dy =(*r)->y()-Y;
            float x  = dx*ax+dy*ay;
            float y  =-dx*ay+dy*ax;
            float xy = x*x+y*y;
            if (xy == 0.) continue;
            float r2 = 1.f/xy;
            float Ut = x*r2;
            float Vt = y*r2;
            float UR = Ut*R+1.f;
            if (UR == 0.) continue;
            float A  = Vt*R/UR;
            float B  = Vt-A*Ut;
            if (std::abs(B*data.K) > m_ipt*std::sqrt(1.f+A*A)) continue;
            ++nseed;
            newSeed(data, (*r), (*r0), Zo);
          }
        }
        if (nseed < m_maxsize) continue;
        data.endlist=false;
        data.rMin = (++r0);
        data.fvNmin=f;
        data.zMin=z;
        return;
      }
    }
  }
  data.endlist = true;
}

// Production 3 space points seeds 
 
void InDet::SiSpacePointsSeedMaker_ATLxk::production3Sp(EventData& data) const
{

  if (data.nsaz<3) return;

 
 
  const std::array<int,arraySizeZ> zBinIndex {5,6,7,8,9,10,4,3,2,1,0};

  std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_topCands;
  std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_endTopCands;
  std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_bottomCands;
  std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator,arraySizeNeighbourBins> iter_endBottomCands;
  
  int    nseed =    0;
  data.endlist = true;

  for (int phiBin=data.fNmin; phiBin<=m_maxPhiBin; ++phiBin) {
    
    int z = 0;
    if (!data.endlist) z = data.zMin;

    for (; z<arraySizeZ; ++z) {
      
      int phiZbin  = phiBin *arraySizeZ+zBinIndex[z];

      if (!data.rfz_map[phiZbin]) continue;

      int numberBottomCells = 0;
      int numberTopCells = 0;

      for (int neighbourCellNumber=0; neighbourCellNumber<m_nNeighbourCellsBottom[phiZbin]; ++neighbourCellNumber) {
  
        int theNeighbourCell =  m_neighbourCellsBottom[phiZbin][neighbourCellNumber];
        if (!data.rfz_map[theNeighbourCell]) continue;
        iter_bottomCands [numberBottomCells] = data.rfz_Sorted[theNeighbourCell].begin();
        iter_endBottomCands[numberBottomCells++] = data.rfz_Sorted[theNeighbourCell].end();
      } 

      for (int neighbourCellNumber=0; neighbourCellNumber<m_nNeighbourCellsTop[phiZbin]; ++neighbourCellNumber) {
 
        int theNeighbourCell =  m_neighbourCellsTop[phiZbin][neighbourCellNumber];
        if (!data.rfz_map[theNeighbourCell]) continue;
        iter_topCands [numberTopCells] = data.rfz_Sorted[theNeighbourCell].begin();
        iter_endTopCands[numberTopCells++] = data.rfz_Sorted[theNeighbourCell].end();
      } 

      if (!data.trigger) production3Sp       (data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed,zBinIndex[z]);
      else               production3SpTrigger(data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed);
    }
 
    if (nseed>=m_maxsize) {
      data.endlist=false;
      data.fNmin  =  phiBin+1;
      return;
    }
  }
  data.endlist = true;
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::production3Sp
(EventData& data,
 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_bottomCands ,
 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_endBottomCands,
 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_topCands ,
 std::array <std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & iter_endTopCands,
 const int numberBottomCells, const int numberTopCells, int& nseed, const int zbin) const
{ 

  std::vector<InDet::SiSpacePointForSeed*>::iterator iter_centralSP=iter_bottomCands[0];
  std::vector<InDet::SiSpacePointForSeed*>::iterator iter_otherSP;  


  bool isStrip = ((*iter_centralSP)->spacepoint->clusterList().second); 

  bool isBarrelRegion = (zbin >=4 && zbin <= 6); 
  bool isTransitionRegion = (zbin == 3 || zbin == 7);   // region in z = 450 - 925mm 
  float rmin =  40.;    
  float rmax = 140.;    
 
  if (isBarrelRegion){
    rmax = 100;   
  }
  if(isStrip) {
    rmin = 285.;    
    rmax = 450.;    
    if (isBarrelRegion){
      rmin = 335;  
      rmax = 450.; 
    }
    else if (isTransitionRegion){
      rmin = 335;  
      rmax = 550;  
    }
  }

  for(; iter_centralSP!=iter_endBottomCands[0]; ++iter_centralSP) {
    if((*iter_centralSP)->radius() > rmin) break;
  }

  iter_topCands[0] = iter_centralSP; 
  ++iter_topCands[0];

  const float ipt2K = data.ipt2K;
  const float ipt2C = data.ipt2C;
  const float COFK  = data.COFK;
  const float maxd0cut = m_maxdImpact;
  const float maxd0cutstrips = m_maxdImpactDecays;
  const float zmin  = data.zminU;
  const float zmax  = data.zmaxU;
  const float dzdrmax = data.dzdrmax;
  const float dzdrmin = data.dzdrmin;
  data.CmSp.clear();

  size_t SPcapacity = data.SP.size(); 

  for (; iter_centralSP!=iter_endBottomCands[0]; ++iter_centralSP) {
 
    const float& R  = (*iter_centralSP)->radius();
    if(R > rmax) break; 

    const float&        X    = (*iter_centralSP)->x();
    const float&        Y    = (*iter_centralSP)->y();
    const float&        Z    = (*iter_centralSP)->z();

    double absZ = std::abs(Z); 
    if (isStrip && absZ > 2650. ) continue;
    if (!isStrip && absZ > 600.) continue;
    if (isStrip && isTransitionRegion && absZ < 750. && R > 450.) continue; 

    size_t                 Nb   = 0;

    for (int cell=0; cell<numberBottomCells; ++cell) {
      for (iter_otherSP=iter_bottomCands[cell]; iter_otherSP!=iter_endBottomCands[cell]; ++iter_otherSP) {
        
        const float& Rb =(*iter_otherSP)->radius();
        float dR = R-Rb;

        if (dR > m_drmax) {
          iter_bottomCands[cell]=iter_otherSP;
          continue;
        }
        if (dR < m_drmin || (data.iteration && (*iter_otherSP)->spacepoint->clusterList().second)) break;

        const float dZdR = (Z-(*iter_otherSP)->z())/dR;
        const float absdZdR = std::abs(dZdR);
        if (absdZdR < dzdrmin or absdZdR > dzdrmax) continue;
  
        const float z0 = Z-R*dZdR;
          if(z0 > zmax || z0 < zmin) continue;
        data.SP[Nb] = (*iter_otherSP);
        if(m_writeNtuple) data.SP[Nb]->setDZDR(dZdR);
        if (++Nb==SPcapacity){
            data.resizeSPCont();
            SPcapacity=data.SP.size();
        }
      } 
    } 

    if (!Nb) continue;

    size_t Nt = Nb;
    

    for (int cell=0; cell<numberTopCells; ++cell) {
      for (iter_otherSP=iter_topCands[cell];iter_otherSP!=iter_endTopCands[cell]; ++iter_otherSP) {
  
        float Rt =(*iter_otherSP)->radius();
        float dR = Rt-R;

        if (dR<m_drmin) {
          iter_topCands[cell]=iter_otherSP;
          continue;
        }
        if (dR>m_drmax) break;

        float dZdR = ((*iter_otherSP)->z()-Z)/dR;
        float absdZdR = std::abs(dZdR);
          if (absdZdR < dzdrmin or absdZdR > dzdrmax) continue;

        float z0 = Z-R*dZdR;
          if(z0 > zmax || z0 < zmin) continue;
        data.SP[Nt] = (*iter_otherSP);
        if (m_writeNtuple)  data.SP[Nt]->setDZDR(dZdR);
        if (++Nt==SPcapacity) {
          data.resizeSPCont(); 
          SPcapacity=data.SP.size();
        }
      } 
    } 
    
    if (!(Nt-Nb)) continue;

    float covr0 = (*iter_centralSP)->covr ();
    float covz0 = (*iter_centralSP)->covz ();

    float ax    = X/R;
    float ay    = Y/R;

    for (size_t i=0; i<Nt; ++i) {
 
      InDet::SiSpacePointForSeed* sp = data.SP[i];

      float dx  = sp->x()-X;
      float dy  = sp->y()-Y;
      float dz  = sp->z()-Z;
      float x   = dx*ax+dy*ay;
      float y   = dy*ax-dx*ay;

      float r2  = 1.f/(x*x+y*y);
      float dr  = std::sqrt(r2);
      float tz  = dz*dr;

      if (i < Nb) tz = -tz;

      data.Tz[i]   = tz;      
      data.Zo[i]   = Z-R*tz;  
      data.R [i]   = dr;      
      data.U [i]   = x*r2;    
      data.V [i]   = y*r2;    
      data.Er[i]   = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;    
    }
 
    data.nOneSeeds = 0;
    data.mapOneSeeds_Pro.clear();

    for (size_t b=0; b<Nb; ++b) {
      
      float  Zob  = data.Zo[b];       
      float  Tzb  = data.Tz[b];       
      float  Erb  = data.Er[b];       
      float  Vb   = data.V [b];       
      float  Ub   = data.U [b];       
      float  Tzb2 = (1.f+Tzb*Tzb);     
      float sTzb2 = std::sqrt(Tzb2);       
      float sigmaSquaredScatteringPtDependent  = Tzb2*COFK;        
      float sigmaSquaredScatteringMinPt  = Tzb2*ipt2C;       
      float d0max  = maxd0cut;
      if (data.SP[b]->spacepoint->clusterList().second) d0max = maxd0cutstrips;

      for (size_t t=Nb;  t<Nt; ++t) {

        
        float meanOneOverTanTheta = (Tzb+data.Tz[t])/2.f; 
        float theta = 0.;
        if(m_writeNtuple){
          theta = std::atan(1.f/meanOneOverTanTheta);
        }
        float sigmaSquaredSpacePointErrors = Erb+data.Er[t]   
                        + 2.f * covz0 * data.R[t]*data.R[b]    
                        + 2.f * covr0 * data.R[t]*data.R[b] * meanOneOverTanTheta * meanOneOverTanTheta; // mixed term with r-uncertainy on central SP
        float remainingSquaredDelta  = (Tzb-data.Tz[t])*(Tzb-data.Tz[t]) - sigmaSquaredSpacePointErrors; 

        if (remainingSquaredDelta  - sigmaSquaredScatteringMinPt > 0 ) continue;    
 
        float deltaU  = data.U[t]-Ub;
        if (deltaU == 0.) continue;       
        float A   = (data.V[t]-Vb)/deltaU;  
        float B   = Vb-A*Ub;            
        float onePlusAsquare  = 1.f+A*A;
        float BSquare  = B*B;

        if (BSquare  > ipt2K*onePlusAsquare || remainingSquaredDelta*onePlusAsquare > BSquare*sigmaSquaredScatteringPtDependent) continue;
        float d0 = 0;         
        if(std::abs(B) < 1e-10) d0  = std::abs((A-B*R)*R);  
        else{
          float x0 = -A/(2.f*B);
          float rTrack = std::sqrt(onePlusAsquare/BSquare)*.5f;
          d0 = std::abs(-rTrack + std::sqrt(rTrack*rTrack +2.f*x0*R +R*R));
        }

        if (d0 <= d0max) {
          float dr = data.R[b];
          if (data.R[t] < data.R[b]) dr = data.R[t];
          data.SP[t]->setScorePenalty(std::abs((Tzb-data.Tz[t])/(dr*sTzb2)));
          data.SP[t]->setParam(d0);
 
          if(m_writeNtuple){
            data.SP[t]->setEta(-std::log(std::tan(0.5f*theta)));
            data.SP[t]->setPt(std::sqrt(onePlusAsquare/BSquare)/(1000.f*data.K)); 
          }
          data.CmSp.emplace_back(B/std::sqrt(onePlusAsquare), data.SP[t]);
 
        }
   
      }   
      if (!data.CmSp.empty()) {
        newOneSeedWithCurvaturesComparison(data, data.SP[b], (*iter_centralSP), Zob);
      }
    } 
    fillSeeds(data);
    nseed += data.fillOneSeeds;
   
  } 
}

// Production 3 space points seeds in ROI
 
void InDet::SiSpacePointsSeedMaker_ATLxk::production3SpTrigger
(EventData& data,
 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rb ,
 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rbe,
 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rt ,
 std::array<std::vector<InDet::SiSpacePointForSeed*>::iterator, arraySizeNeighbourBins> & rte,
 const int numberBottomCells, const int numberTopCells, int& nseed) const
{
  constexpr float twoPi = 2.*M_PI;
 
  std::vector<InDet::SiSpacePointForSeed*>::iterator r0=rb[0],r;
 
  float rmin =  40.;
  float rmax = 140.;
  if((*r0)->spacepoint->clusterList().second) {
    rmin = 280.;
    rmax = 540.;
  }
 
  for(; r0!=rbe[0]; ++r0) {if((*r0)->radius() > rmin) break;}
  rt[0] = r0; ++rt[0];
 
  float ipt2K = data.ipt2K;
  float ipt2C = data.ipt2C;
  float COFK  = data.COFK;
  float maxd0cut = m_maxdImpact;
  float maxd0cutstrips = m_maxdImpactSSS;
 
  data.CmSp.clear();
 
  // Loop through all trigger space points
  //
  for (; r0!=rbe[0]; ++r0) {
 
    data.nOneSeeds = 0;
    data.mapOneSeeds_Pro.clear();
  
    float R  = (*r0)->radius();
    if(R>rmax) break;
 
    const Trk::Surface* sur0 = (*r0)->sur();
    float               X    = (*r0)->x();
    float               Y    = (*r0)->y();
    float               Z    = (*r0)->z();
    int                 Nb   = 0;
 
    // Bottom links production
    //
    for (int i=0; i<numberBottomCells; ++i) {
 
      for (r=rb[i]; r!=rbe[i]; ++r) {
  
        float Rb =(*r)->radius();
 
        float dR = R-Rb;
        if (dR > m_drmax) {
          rb[i]=r;
          continue;
        }
        if ((*r)->sur()==sur0) continue;
 
        if (dR < m_drmin || (data.iteration && (*r)->spacepoint->clusterList().second)) break;
 
        // Comparison with  bottom and top Z 
        //
        float Tz = (Z-(*r)->z())/dR;
        float Zo = Z-R*Tz;
        if (Zo < data.zminB || Zo > data.zmaxB) continue;
        float Zu = Z+(550.f-R)*Tz;
        if (Zu < data.zminU || Zu > data.zmaxU) continue;
        data.SP[Nb] = (*r);
        if (++Nb==m_maxsizeSP) goto breakb;
      }
    }
  breakb:
    if (!Nb || Nb==m_maxsizeSP) continue;
    int Nt = Nb;
    
    // Top   links production
    //
    for (int i=0; i<numberTopCells; ++i) {
      
      for (r=rt[i]; r!=rte[i]; ++r) {
  
        float Rt =(*r)->radius();
        float dR = Rt-R;
  
        if (dR<m_drmin) {
          rt[i]=r;
          continue;
        }
        if (dR>m_drmax) break;
        if ((*r)->sur()==sur0) continue;
 
        // Comparison with  bottom and top Z 
        //
        float Tz = ((*r)->z()-Z)/dR;
        float Zo = Z-R*Tz;
        if (Zo < data.zminB || Zo > data.zmaxB) continue;
        float Zu = Z+(550.f-R)*Tz;
        if (Zu < data.zminU || Zu > data.zmaxU) continue;
        data.SP[Nt] = (*r);
        if (++Nt==m_maxsizeSP) goto breakt;
      }
    }
    
  breakt:
    if (!(Nt-Nb)) continue;
    float covr0 = (*r0)->covr ();
    float covz0 = (*r0)->covz ();
 
    float ax   = X/R;
    float ay   = Y/R;
    
    for (int i=0; i<Nt; ++i) {
 
      InDet::SiSpacePointForSeed* sp = data.SP[i];
 
      float dx  = sp->x()-X;
      float dy  = sp->y()-Y;
      float dz  = sp->z()-Z;
      float x   = dx*ax+dy*ay;
      float y   = dy*ax-dx*ay;
      float r2  = 1.f/(x*x+y*y);
      float dr  = std::sqrt(r2);
      float tz  = dz*dr;
      if (i < Nb) tz = -tz;
 
      data.Tz[i]   = tz;
      data.Zo[i] = Z-R*tz;
      data.R [i]   = dr;
      data.U [i]   = x*r2;
      data.V [i]   = y*r2;
      data.Er[i]   = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;
    }
    covr0      *= .5f;
    covz0      *= 2.f;
   
    // Three space points comparison
    //
    for (int b=0; b<Nb; ++b) {
    
      float  Zob  = data.Zo[b];
      float  Tzb  = data.Tz[b];
      float  Rb2r = data.R [b]*covr0;
      float  Rb2z = data.R [b]*covz0;
      float  Erb  = data.Er[b];
      float  Vb   = data.V [b];
      float  Ub   = data.U [b];
      float  Tzb2 = (1.f+Tzb*Tzb);
      float  CSA  = Tzb2*COFK;
      float ICSA  = Tzb2*ipt2C;
      float d0max  = maxd0cut;
      if (data.SP[b]->spacepoint->clusterList().second) d0max = maxd0cutstrips;
      
      for (int t=Nb;  t<Nt; ++t) {
  
        float dT  = ((Tzb-data.Tz[t])*(Tzb-data.Tz[t])-data.R[t]*Rb2z-(Erb+data.Er[t]))-(data.R[t]*Rb2r)*((Tzb+data.Tz[t])*(Tzb+data.Tz[t]));
        if ( dT > ICSA) continue;
 
        float deltaU  = data.U[t]-Ub;
        if (deltaU == 0.) continue;
        float A   = (data.V[t]-Vb)/deltaU;
        float onePlusAsquare  = 1.f+A*A;
        float B   = Vb-A*Ub;
        float BSquare  = B*B;
        if (BSquare  > ipt2K*onePlusAsquare || dT*onePlusAsquare > BSquare*CSA) continue;
 
        float Im  = std::abs((A-B*R)*R);
        if (Im > d0max) continue;
 
        // Azimuthal angle test
        //
        float y  = 1.;
        float x  = 2.f*B*R-A;
        float df = std::abs(std::atan2(ay*y-ax*x,ax*y+ay*x)-data.ftrig);
        if (df > M_PI) df = twoPi-df;
        if (df > data.ftrigW) continue;
        data.CmSp.emplace_back(B/std::sqrt(onePlusAsquare), data.SP[t]);
        data.SP[t]->setParam(Im);
      }
      if (!data.CmSp.empty()) {
        newOneSeedWithCurvaturesComparison(data, data.SP[b], (*r0), Zob);
      }
    }
    fillSeeds(data);
    nseed += data.fillOneSeeds;
  }
}

// New 3 space points pro seeds 
 
void InDet::SiSpacePointsSeedMaker_ATLxk::newOneSeed
(EventData& data,
 InDet::SiSpacePointForSeed*& p1, InDet::SiSpacePointForSeed*& p2,
 InDet::SiSpacePointForSeed*& p3, float z, float seedCandidateQuality) const
{
  float worstQualityInMap = std::numeric_limits<float>::min();
  InDet::SiSpacePointsProSeed* worstSeedSoFar = nullptr;
  if (!data.mapOneSeeds_Pro.empty()) {
    std::multimap<float,InDet::SiSpacePointsProSeed*>::reverse_iterator l = data.mapOneSeeds_Pro.rbegin();
    worstQualityInMap = (*l).first;
    worstSeedSoFar = (*l).second;
  }
  if (data.nOneSeeds < data.maxSeedsPerSP
      || (data.keepAllConfirmedSeeds  && worstQualityInMap <= seedCandidateQuality  && isConfirmedSeed(p1,p3,seedCandidateQuality) && data.nOneSeeds < data.seedPerSpCapacity)
      || (data.keepAllConfirmedSeeds  && worstQualityInMap >  seedCandidateQuality  && isConfirmedSeed(worstSeedSoFar->spacepoint0(),worstSeedSoFar->spacepoint2(),worstQualityInMap) && data.nOneSeeds < data.seedPerSpCapacity)
    ){
    data.OneSeeds_Pro[data.nOneSeeds].set(p1,p2,p3,z);
    data.mapOneSeeds_Pro.insert(std::make_pair(seedCandidateQuality, &data.OneSeeds_Pro[data.nOneSeeds]));
    ++data.nOneSeeds;
  } 
  else if (worstQualityInMap > seedCandidateQuality){
      worstSeedSoFar->set(p1,p2,p3,z);
      std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator 
        i = data.mapOneSeeds_Pro.insert(std::make_pair(seedCandidateQuality,worstSeedSoFar));
      for (++i; i!=data.mapOneSeeds_Pro.end(); ++i) {
        if ((*i).second==worstSeedSoFar) {
          data.mapOneSeeds_Pro.erase(i);
          return;
        }
      }
  }
}

// New 3 space points pro seeds production
namespace {
   inline
   float computeEta( float r, float z) {
      if (r <= 10e-9) return 0;
      auto asinh = [] (double x) { return std::log (std::sqrt (x*x+1) + x); };
      return asinh (z / r);
   }
}
 
void InDet::SiSpacePointsSeedMaker_ATLxk::newOneSeedWithCurvaturesComparison
(EventData& data, SiSpacePointForSeed*& SPb, SiSpacePointForSeed*& SP0, float Zob) const
{
  constexpr float curvatureInterval = .00003;
 
  bool  bottomSPisPixel = !SPb->spacepoint->clusterList().second;
  float bottomSPQuality   = SPb->quality();
  float centralSPQuality   = SP0->quality();

  if(data.CmSp.size() > 2) std::sort(data.CmSp.begin(), data.CmSp.end(), comCurvature());
      
  float bottomR=SPb->radius();
  float bottomZ=SPb->z();
 
  std::vector<std::pair<float,InDet::SiSpacePointForSeed*>>::iterator it_otherSP;
  std::vector<std::pair<float,InDet::SiSpacePointForSeed*>>::iterator it_commonTopSP = data.CmSp.begin(), ie = data.CmSp.end();
  std::vector<std::pair<float,InDet::SiSpacePointForSeed*>>::iterator it_startInnerLoop=it_commonTopSP;

  for (; it_commonTopSP!=ie; ++it_commonTopSP) {

    float seedIP   = (*it_commonTopSP).second->param();
    float seedQuality    = seedIP + (*it_commonTopSP).second->scorePenalty();
    float originalSeedQuality   = seedQuality;
 
    if(m_maxdImpact > 50){      //This only applies to LRT
      // Tell clang to optimize assuming that FP exceptions can trap.
      // Otherwise, it can vectorize the division, which can lead to
      // spurious division-by-zero traps from unused vector lanes.
      CXXUTILS_TRAPPING_FP;
 
      float topR=(*it_commonTopSP).second->radius();
      float topZ=(*it_commonTopSP).second->z();
 
      float Zot = std::abs(topR - bottomR) > 10e-9 ?
        bottomZ - (bottomR - originalSeedQuality) * ((topZ - bottomZ) / (topR - bottomR)) : bottomZ;
 
      float eta1 = computeEta(topR - bottomR, topZ - bottomZ);
      float eta0 = computeEta(seedIP, Zot);
 
      float deltaEta=std::abs(eta1-eta0); //For LLP daughters, the direction of the track is correlated with the direction of the LLP (which is correlated with the direction of the point of closest approach
      //calculate weighted average of d0 and deltaEta, normalized by their maximum values
      float f=std::min(0.5f,originalSeedQuality/200.f);  //0.5 and 200 are parameters chosen from a grid scan to optimize efficiency
      seedQuality*=(1.f-f)/300.f;
      seedQuality+=f*deltaEta/2.5f;
    }
 
    bool                topSPisPixel = !(*it_commonTopSP).second->spacepoint->clusterList().second;
    
    const Trk::Surface* surfaceTopSP  = (*it_commonTopSP).second->sur   ();
    float               radiusTopSP   = (*it_commonTopSP).second->radius();
    float               minCurvature  =(*it_commonTopSP).first-curvatureInterval;
    float               maxCurvature  =(*it_commonTopSP).first+curvatureInterval;
      

    if      (!bottomSPisPixel) seedQuality+=m_seedScoreBonusSSS;
    else if ( topSPisPixel) seedQuality+=m_seedScoreBonusPPP;
 
 
 
    for (it_otherSP=it_startInnerLoop;  it_otherSP!=ie; ++it_otherSP) {
      if (       it_otherSP == it_commonTopSP           ) continue;
      if ( (*it_otherSP).first < minCurvature       ) {
        it_startInnerLoop=it_otherSP; 
        ++it_startInnerLoop; 
        continue;
      }
      if ( (*it_otherSP).first > maxCurvature       ) break;
      if ( (*it_otherSP).second->sur()==surfaceTopSP) continue;
      float radiusOtherSP = (*it_otherSP).second->radius();
      if (std::abs(radiusOtherSP-radiusTopSP) < m_drmin) continue;
      // if we have a confirmation seed, we improve the score of the seed. 
      seedQuality += m_seedScoreBonusConfirmationSeed;
      // apply confirmation bonus only once
      break; 
    }

    if (seedQuality > data.maxScore) continue;
    
    if (bottomSPisPixel!=topSPisPixel) {
      if (seedQuality > 0. || 
        (seedQuality > bottomSPQuality && seedQuality > centralSPQuality && seedQuality > (*it_commonTopSP).second->quality()) 
      ) continue;
    }
    if (!isConfirmedSeed(SPb,it_commonTopSP->second,seedQuality)){
      double                maxdImpact = m_maxdImpact    - (m_dImpactCutSlopeUnconfirmedPPP * (*it_commonTopSP).second->scorePenalty()); 
      if (!bottomSPisPixel) maxdImpact = m_maxdImpactSSS - (m_dImpactCutSlopeUnconfirmedSSS * (*it_commonTopSP).second->scorePenalty()); 
      if (seedIP > maxdImpact) continue; 
    }
    newOneSeed(data, SPb, SP0, (*it_commonTopSP).second, Zob, seedQuality);
  } 
  data.CmSp.clear();
}

 
void InDet::SiSpacePointsSeedMaker_ATLxk::fillSeeds(EventData& data) const
{
 
  data.fillOneSeeds = 0;
 
  std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator it_firstSeedCandidate = data.mapOneSeeds_Pro.begin();
  std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator it_seedCandidate  = data.mapOneSeeds_Pro.begin();
  std::multimap<float,InDet::SiSpacePointsProSeed*>::iterator it_endSeedCandidates = data.mapOneSeeds_Pro.end();
  
  if (it_seedCandidate==it_endSeedCandidates) return;
 
  SiSpacePointsProSeed* theSeed{nullptr};

  for (; it_seedCandidate!=it_endSeedCandidates; ++it_seedCandidate) {

    float quality = (*it_seedCandidate).first;
    theSeed       = (*it_seedCandidate).second;
    if (it_seedCandidate!=it_firstSeedCandidate && theSeed->spacepoint0()->radius() < m_radiusCutIBL && quality > m_seedScoreThresholdPPPConfirmationSeed) continue;
    if (!theSeed->setQuality(quality)) continue;
    
    if (data.i_seede_Pro!=data.l_seeds_Pro.end()) {
      theSeed  = &(*data.i_seede_Pro++);
      *theSeed = *(*it_seedCandidate).second;
    } else {
      data.l_seeds_Pro.emplace_back(*(*it_seedCandidate).second);
      theSeed = &(data.l_seeds_Pro.back());
      data.i_seede_Pro = data.l_seeds_Pro.end();
    }
    
   ++data.fillOneSeeds;
  } 
}
 
const InDet::SiSpacePointsSeed* InDet::SiSpacePointsSeedMaker_ATLxk::next(const EventContext&, EventData& data) const
{
  if (not data.initialized) initializeEventData(data);
  
  if (data.nspoint==3) {
    do {
      if (data.i_seed_Pro==data.i_seede_Pro) {
        findNext(data);
        //cppcheck-suppress identicalInnerCondition 
        if (data.i_seed_Pro==data.i_seede_Pro) return nullptr;
      }
    } while (!(*data.i_seed_Pro++).set3(data.seedOutput));

    return &data.seedOutput;
  } else {
    if (data.i_seed_Pro==data.i_seede_Pro) {
      findNext(data);
      //cppcheck-suppress identicalInnerCondition
      if (data.i_seed_Pro==data.i_seede_Pro) return nullptr;
    } 
    (*data.i_seed_Pro++).set2(data.seedOutput);
    return &data.seedOutput;
  }
  return nullptr;
}

void InDet::SiSpacePointsSeedMaker_ATLxk::writeNtuple(const SiSpacePointsSeed* seed, const Trk::Track* track, int seedType, long eventNumber) const{
  
  if(m_writeNtuple) {
    std::lock_guard<std::mutex> lock(m_mutex);
 
    if(track != nullptr) {
      m_trackPt = (track->trackParameters()->front()->pT())/1000.f;
      m_trackEta = std::abs(track->trackParameters()->front()->eta());
    }
    else {
      m_trackPt = -1.;
      m_trackEta = -1.; 
    }
    m_d0           =   seed->d0();
    m_z0           =   seed->zVertex();
    m_eta          =   seed->eta();
    m_x1           =   seed->x1();
    m_x2           =   seed->x2();
    m_x3           =   seed->x3();
    m_y1           =   seed->y1();
    m_y2           =   seed->y2();
    m_y3           =   seed->y3();      
    m_z1           =   seed->z1();
    m_z2           =   seed->z2();
    m_z3           =   seed->z3();
    m_r1           =   seed->r1();
    m_r2           =   seed->r2();
    m_r3           =   seed->r3();
    m_type         =   seedType;
    m_dzdr_b       =   seed->dzdr_b();
    m_dzdr_t       =   seed->dzdr_t();
    m_pt           =   seed->pt();
    m_givesTrack   =   !(track == nullptr);
    m_eventNumber  =   eventNumber;
 
    // Ok: protected by mutex.
    TTree* outputTree ATLAS_THREAD_SAFE = m_outputTree;
    outputTree->Fill();
 
  }
 
}
 
bool InDet::SiSpacePointsSeedMaker_ATLxk::isZCompatible  
(EventData& data, const float Zv, const float R, const float T) const
{
  if (Zv < data.zminU or Zv > data.zmaxU) return false;
  if (not data.isvertex) return true;
  if (data.l_vertex.empty()) return false;
 
  float dZmin = std::numeric_limits<float>::max();
  for (const float& v: data.l_vertex) {
    float dZ = std::abs(v-Zv);
    if (dZ >= dZmin) break;
    dZmin = dZ;
  }
 
  //return dZmin < (m_dzver+m_dzdrver*R)*sqrt(1.+T*T);
  //(Minor) speed-up: Avoid calculation of sqrt, compare squares
  return dZmin*dZmin < (m_dzver+m_dzdrver*R)*(m_dzver+m_dzdrver*R)*(1.f+T*T);
}

 
InDet::SiSpacePointForSeed* InDet::SiSpacePointsSeedMaker_ATLxk::newSpacePoint
(EventData& data, const Trk::SpacePoint*const& sp) const
{
  InDet::SiSpacePointForSeed* sps = nullptr;

  std::array<float,3> r{0,0,0};
  convertToBeamFrameWork(data, sp, r);

  if (data.checketa) {
    float z = (std::abs(r[2])+m_zmax); 
    float x = r[0]*data.dzdrmin;
    float y = r[1]*data.dzdrmin;
    if ((z*z )<(x*x+y*y)) return sps;
  }
  if (data.i_spforseed!=data.l_spforseed.end()) {
    sps = &(*data.i_spforseed++);
    sps->set(sp,r);
  } else {
    data.l_spforseed.emplace_back(sp, r);
    sps = &(data.l_spforseed.back());
    data.i_spforseed = data.l_spforseed.end();
  }
      
  return sps;
}

// New 2 space points seeds 
 
void InDet::SiSpacePointsSeedMaker_ATLxk::newSeed
(EventData& data, InDet::SiSpacePointForSeed*& p1, InDet::SiSpacePointForSeed*& p2, float z) 
{
  InDet::SiSpacePointForSeed* p3 = nullptr;
 
  if (data.i_seede_Pro!=data.l_seeds_Pro.end()) {
    SiSpacePointsProSeed* s = &(*data.i_seede_Pro++);
    s->set(p1, p2, p3, z);
  } else {
    data.l_seeds_Pro.emplace_back(p1, p2, p3, z);
    data.i_seede_Pro = data.l_seeds_Pro.end();
  }
}
 
void InDet::SiSpacePointsSeedMaker_ATLxk::initializeEventData(EventData& data) const {
  int seedArrayPerSPSize = (m_maxOneSizePPP>m_maxOneSizeSSS ? m_maxOneSizePPP : m_maxOneSizeSSS); 
  if (m_alwaysKeepConfirmedStripSeeds || m_alwaysKeepConfirmedPixelSeeds)  seedArrayPerSPSize = 50; 
  data.initialize(EventData::ToolType::ATLxk,
                  m_maxsizeSP,
                  seedArrayPerSPSize,
                  0, 
                  m_nBinsR,
                  0, 
                  arraySizePhiZ,
                  arraySizePhiZV,
                  m_checketa);
}
 
bool InDet::SiSpacePointsSeedMaker_ATLxk::isConfirmedSeed(const InDet::SiSpacePointForSeed* bottomSP, const InDet::SiSpacePointForSeed* topSP, float quality) const{

    if (bottomSP->spacepoint->clusterList().second){
      return (quality < m_seedScoreThresholdSSSConfirmationSeed); 
    }
    else if (!topSP->spacepoint->clusterList().second){
      return (quality < m_seedScoreThresholdPPPConfirmationSeed); 
    }
    else return (quality < 0.); 
}
 
bool InDet::SiSpacePointsSeedMaker_ATLxk::getWriteNtupleBoolProperty() const{
    return m_writeNtuple;
}