CaTripletConstructorSW.cxx

Go to the documentation of this file.

/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergey Gorbunov [committer] */

 
#include "CaTripletConstructorSW.h"
 
#include "AlgoFairloggerCompat.h"
#include "CaDebugger.h"
#include "CaDefines.h"
#include "CaFramework.h"
#include "CaGridArea.h"
#include "KfTrackKalmanFilter.h"
#include "KfTrackParam.h"
 
#include <algorithm>
#include <iostream>
 
 
using namespace cbm::algo::ca;
 
TripletConstructorSW::TripletConstructorSW(const ca::Framework& framework, const ca::Parameters<fvec>& pars,
                                           WindowData& wData, const fscal mass, const ca::TrackingMode& mode)
  : fFramework(framework)
  , fParameters(pars)
  , fSetup(fParameters.GetActiveSetup())
  , fField(fSetup.GetField())
  , fSwMaps(fParameters.GetSearchWindows())
  , frWData(wData)
  , fDefaultMass(mass)
  , fTrackingMode(mode)
{
  // FIXME: SZh 24.08.2022
  //        This approach is suitable only for a case, when all the stations inside a magnetic field come right before
  //        all the stations outside of the field!
  fNfieldStations =
    std::lower_bound(fSetup.GetActiveLayers().cbegin(), fSetup.GetActiveLayers().cbegin() + fSetup.GetNofLayers(),
                     0,  // we are looking for the first zero element
                     [](const auto& s, int edge) { return s.IsInField() > edge; })
    - fSetup.GetActiveLayers().cbegin();
 
  fIsTargetField = !frWData.TargB().IsZero();
 
  if constexpr (fDebugDublets) {
    utils::Debugger::Instance().AddNtuple("doubletSW",
                                          "event:iter:sta:p:q:vx:vy:vz:x0:y0:z0:x1:y1:z1:ax:ay:az:adx:ady:dx:dy:dz");
  }
 
  if constexpr (fDebugTriplets) {
    //utils::Debugger::Instance().AddNtuple(
    //"triplets", "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf:chi2time:ndfTime");
    utils::Debugger::Instance().AddNtuple(
      "tripletSW", "event:iter:sta:jumpL:jumpM:p:q:vx:vy:vz:x0:y0:z0:x1:y1:z1:x2:y2:z2:ax:ay:az:adx:ady:dx:dy:dz");
  }
 
  //utils::Debugger::Instance().AddNtuple(
  //"fit", "good:event:iter:sta:jumpL:jumpM:p:q:vx:vy:vz:x:y:z:chi2:chi2kf:chi2kfS:chi2kfT");
}
 
 
bool TripletConstructorSW::InitStations(int istal, int istam, int istar)
{
  fIstaL = istal;
  fIstaM = istam;
  fIstaR = istar;
 
  if (fIstaM >= fSetup.GetNofLayers()) {
    return false;
  }
  fStaL = &fSetup.GetActiveLayer(fIstaL);
  fStaM = &fSetup.GetActiveLayer(fIstaM);
  fStaR = &fSetup.GetActiveLayer(fIstaR);
 
  // FIXME: The following blocks do nothing => remove
  {  // two stations for approximating the field between the target and the left hit
    const int sta1 = (fNfieldStations <= 1) ? 1 : std::clamp(fIstaL, 1, fNfieldStations - 1);
    const int sta0 = sta1 / 2;  // station between fIstaL and the target
 
    assert(0 <= sta0 && sta0 < sta1 && sta1 < fSetup.GetNofLayers());
  }
 
  {  // three stations for approximating the field between the left and the right hit
 
    int sta0 = fIstaL;
    int sta1 = fIstaM;
    int sta2 = fIstaM + 1;
    if (sta2 >= fSetup.GetNofLayers()) {
      sta2 = fIstaM;
      sta1 = fIstaM - 1;
      sta0 = (sta1 <= 0) ? 0 : std::clamp(sta0, 0, sta1 - 1);
    }
    if (fSetup.GetNofLayers() >= 3) {
      assert(0 <= sta0 && sta0 < sta1 && sta1 < sta2 && sta2 < fSetup.GetNofLayers());
    }
  }
  return true;
}
 
 
void TripletConstructorSW::CreateTripletsForHit(Vector<ca::Triplet>& tripletsOut, int istal, int istam, int istar,
                                                ca::HitIndex_t iHitL)
{
  tripletsOut.clear();
 
  if (!InitStations(istal, istam, istar) || (fIstaR >= fSetup.GetNofLayers())) {
    return;
  }
 
  fIhitL           = iHitL;
  const auto& hitl = frWData.Hit(fIhitL);
 
  {  
    const bool matchMc = fParameters.GetConfig().DevMatchDoubletsViaMc();
    //std::cout << "MC for hit " << fIhitL << ", " << hitl.Id() << " " << ca::Framework::GetMcTrackIdForCaHit(hitl.Id())
    //        << std::endl;
    McMatch mc = matchMc ? fFramework.GetMcMatchForCaHit(hitl.Id()) : McMatch();
    fDoubletData.clear();
    if (mc.IsEmpty() && matchMc) {
      return;
    }
 
    const DoubletSearchWindowMap& sa = fSwMaps.DoubletSw(frWData.CurrentIterationIndex(), fIstaM, istam - istal - 1);
 
    SearchWindowMap::SearchWindow area2D =
      sa.GetSearchWindow(hitl.X(), hitl.Y(), hitl.Z(), hitl.RangeX(), hitl.RangeY());
 
    CollectHits(fDoubletData, area2D, hitl.T(), hitl.RangeT(), fIstaM,
                sqrt(frWData.CurrentIteration()->GetDoubletChi2Cut()) / 3.5,
                fParameters.GetConfig().GetMaxDoubletsPerSinglet(), mc);
 
    if constexpr (fDebugDublets) {
      if (!mc.IsEmpty() && mc.GetNofLinks() == 1) {
        int trueHitId = FindClosestHitWithMc(area2D, fIstaM, mc);
        if (trueHitId >= 0) {
          const ca::Hit& trueHit = frWData.Hit(trueHitId);
          const McTrack& mctr    = fFramework.GetMcData()->GetTrack(mc.GetLink(0));
          float x                = 0.5 * (area2D.xMax + area2D.xMin);
          float y                = 0.5 * (area2D.yMax + area2D.yMin);
          float z                = fSetup.GetActiveLayer(fIstaM).GetZref()[0];
          utils::Debugger::Instance().FillNtuple(
            "doubletSW", mctr.GetEventId(), frWData.CurrentIterationIndex(), fIstaM, mctr.GetP(), mctr.GetCharge(),
            mctr.GetStartX(), mctr.GetStartY(), mctr.GetStartZ(), hitl.X(), hitl.Y(), hitl.Z(), trueHit.X(),
            trueHit.Y(), trueHit.Z(), x, y, z, 0.5 * (area2D.xMax - area2D.xMin), 0.5 * (area2D.yMax - area2D.yMin),
            trueHit.X() - x, trueHit.Y() - y, trueHit.Z() - z);
        }
      }
    }
  }
 
  SuppressDoubletClones();
 
  FindTripletHits();
  FitTriplets(tripletsOut);
 
  //tripletsOut.clear();
  //FitTriplets(tripletsOut);
}
 
 
void TripletConstructorSW::SuppressDoubletClones()
{
  // ---- Add the middle hits to parameters estimation ----
  Vector<ca::HitIndex_t>& hitsM = fDoubletData;
 
  const ca::Hit& hitl = frWData.Hit(fIhitL);
 
  auto it2 = hitsM.begin();
  for (auto it = hitsM.begin(); it != hitsM.end(); it++) {
 
    const ca::HitIndex_t indM = *it;
    const ca::Hit& hitm       = frWData.Hit(indM);
 
    if (frWData.IsHitSuppressed(indM)) {
      continue;
    }
 
    // check if there are other hits close to the doublet on the same station
    if (ca::kMcbm != fTrackingMode) {
      // TODO: SG: adjust cuts, put them to parameters
      fscal dz = hitm.Z() - hitl.Z();
      fscal tx = (hitm.X() - hitl.X()) / dz;
      fscal ty = (hitm.Y() - hitl.Y()) / dz;
 
      for (auto itClone = it + 1; itClone != hitsM.end(); itClone++) {
 
        const int indClone      = *itClone;
        const ca::Hit& hitClone = frWData.Hit(indClone);
 
        const fscal dzc = hitClone.Z() - hitm.Z();
 
        if (fStaM->IsTimeMeasured()) {
          const fscal dt = hitm.T() - hitClone.T();
          if (!(fabs(dt) <= hitm.RangeT() + hitClone.RangeT())) {
            continue;
          }
        }
 
        const fscal dx = hitm.X() + tx * dzc - hitClone.X();
        if (!(fabs(dx) <= hitm.RangeX() + hitClone.RangeX())) {
          continue;
        }
 
        const fscal dy = hitm.Y() + ty * dzc - hitClone.Y();
        if (!(fabs(dy) <= hitm.RangeY() + hitClone.RangeY())) {
          continue;
        }
 
        if (fParameters.GetConfig().DevSuppressOverlapHitsViaMc()) {
          auto mc = fFramework.GetMcMatchForCaHit(hitl.Id());
          if ((mc != fFramework.GetMcMatchForCaHit(hitm.Id()))
              || (mc != fFramework.GetMcMatchForCaHit(hitClone.Id()))) {
            continue;
          }
        }
 
        frWData.SuppressHit(indClone);
      }
    }
 
    *it2 = indM;
    it2++;
  }  // it
  hitsM.shrink(std::distance(hitsM.begin(), it2));
}
 
 
void TripletConstructorSW::FindTripletHits()
{
  // ---- Find the triplets (right hit). Reformat data into portions of triplets. ----
 
  Vector<ca::HitIndex_t>& hitsM = std::get<0>(fTripletData);
  Vector<ca::HitIndex_t>& hitsR = std::get<1>(fTripletData);
  Vector<fscal>& msM            = std::get<2>(fTripletData);
  {
    const int maxTriplets = fDoubletData.size() * fParameters.GetConfig().GetMaxTripletsPerDoublet();
    hitsM.clear();
    hitsR.clear();
    msM.clear();
    hitsM.reserve(maxTriplets);
    hitsR.reserve(maxTriplets);
    msM.reserve(maxTriplets);
  }
 
  const TripletSearchWindowMap& swMap =
    fSwMaps.TripletSw(frWData.CurrentIterationIndex(), fIstaR, fIstaR - fIstaL - 2, fIstaR - fIstaM - 1);
 
 
  const double maxSlope       = frWData.CurrentIteration()->GetMaxSlope();
  const double tripletChi2Cut = frWData.CurrentIteration()->GetTripletChi2Cut();
 
  Vector<ca::HitIndex_t> collectedHits;
  collectedHits.reserve(fParameters.GetConfig().GetMaxTripletsPerDoublet());
 
  const auto& sta  = fSetup.GetActiveLayer(fIstaR);
  const auto& grid = frWData.Grid(fIstaR);
 
  for (unsigned int i2 = 0; i2 < fDoubletData.size(); i2++) {
 
    const ca::Hit& h0 = frWData.Hit(fIhitL);
    const ca::Hit& h1 = frWData.Hit(fDoubletData[i2]);
 
    // TODO: better estimation of the extrapolated tx for the low-p tracks
    fscal dz = h1.Z() - h0.Z();
    fscal tx = (h1.X() - h0.X()) / dz;
    fscal ty = (h1.Y() - h0.Y()) / dz;
    if ((fabs(tx) > maxSlope) || (fabs(ty) > maxSlope)) {
      continue;
    }
 
    McMatch mc;
 
    if (fParameters.GetConfig().DevMatchTripletsViaMc()) {
      mc       = fFramework.GetMcMatchForCaHit(h0.Id());
      auto mc1 = fFramework.GetMcMatchForCaHit(h1.Id());
      if (mc.IsEmpty() || mc != mc1) {
        continue;
      }
    }
 
    auto [area2D, ms] = swMap.GetSearchWindowAndMs(h0.X(), h0.Y(), h0.Z(), h0.RangeX(), h0.RangeY(), h1.X(), h1.Y(),
                                                   h1.Z(), h1.RangeX(), h1.RangeY());
 
    {  // add correction to the hit differences in Z
      fscal dz1 = grid.GetMinZ() - sta.GetZref()[0];
      fscal dz2 = grid.GetMaxZ() - sta.GetZref()[0];
      fscal dx1 = tx * dz1;
      fscal dx2 = tx * dz2;
      fscal dy1 = ty * dz1;
      fscal dy2 = ty * dz2;
      assert(std::min(dx1, dx2) <= 0.f && std::max(dx1, dx2) >= 0.f);
      assert(std::min(dy1, dy2) <= 0.f && std::max(dy1, dy2) >= 0.f);
      area2D.xMin += std::min(dx1, dx2);
      area2D.xMax += std::max(dx1, dx2);
      area2D.yMin += std::min(dy1, dy2);
      area2D.yMax += std::max(dy1, dy2);
    }
 
    collectedHits.clear();
    CollectHits(collectedHits, area2D, h1.T(), h1.RangeT(), fIstaR, sqrt(tripletChi2Cut) / 3.5,
                fParameters.GetConfig().GetMaxTripletsPerDoublet(), mc);
 
    if constexpr (fDebugTriplets) {  // debug
      auto mc1 = fFramework.GetMcMatchForCaHit(h0.Id());
      auto mc2 = fFramework.GetMcMatchForCaHit(h1.Id());
      if (!mc1.IsEmpty() && mc1 == mc2 && mc1.GetNofLinks() == 1 && mc2.GetNofLinks() == 1) {
 
        int trueHitId = FindClosestHitWithMc(area2D, fIstaR, mc2);
        if (trueHitId >= 0) {
          const ca::Hit& trueHit = frWData.Hit(trueHitId);
          const McTrack& mctr    = fFramework.GetMcData()->GetTrack(mc2.GetLink(0));
 
          float x = 0.5 * (area2D.xMax + area2D.xMin);
          float y = 0.5 * (area2D.yMax + area2D.yMin);
          float z = fSetup.GetActiveLayer(fIstaR).GetZref()[0];
          //const auto& staR = fParameters.GetStation(fIstaR);
 
          utils::Debugger::Instance().FillNtuple(
            "tripletSW", mctr.GetEventId(), frWData.CurrentIterationIndex(), fIstaR, fIstaR - fIstaM - 1,
            fIstaR - fIstaL - 2, mctr.GetP(), mctr.GetCharge(), mctr.GetStartX(), mctr.GetStartY(), mctr.GetStartZ(),
            h0.X(), h0.Y(), h0.Z(), h1.X(), h1.Y(), h1.Z(), trueHit.X(), trueHit.Y(), trueHit.Z(), x, y, z,
            0.5 * (area2D.xMax - area2D.xMin), 0.5 * (area2D.yMax - area2D.yMin), trueHit.X() - x, trueHit.Y() - y,
            trueHit.Z() - z);
        }
      }
    }
 
    if (static_cast<int>(collectedHits.size()) >= fParameters.GetConfig().GetMaxTripletsPerDoublet()) {
      // reject already created triplets for this doublet
 
      // TODO: collect this cases and make a common warning at the end of the data processing
      //      LOG(debug) << "Ca: GetMaxTripletPerDoublets==" << fParameters.GetMaxTripletPerDoublets()
      //                 << " reached in FindTripletHits()";
 
      collectedHits.clear();
    }
 
    for (ca::HitIndex_t& irh : collectedHits) {
      hitsM.push_back(fDoubletData[i2]);
      hitsR.push_back(irh);
      msM.push_back(ms);
    }
  }  // i2
}
 
 
void TripletConstructorSW::FitTriplets(Vector<ca::Triplet>& tripletsOut)
{
 
  Vector<ca::HitIndex_t>& hitsM = std::get<0>(fTripletData);
  Vector<ca::HitIndex_t>& hitsR = std::get<1>(fTripletData);
  Vector<fscal>& msM            = std::get<2>(fTripletData);
  assert(hitsM.size() == hitsR.size());
 
  //tripletsOut.clear();
  //tripletsOut.reserve(hitsM.size());
 
  const bool isMomentumFitted = ((fStaL->IsInField()) || (fStaM->IsInField()) || (fStaR->IsInField()));
 
  if constexpr (fDebugTriplets) {
    //cbm::ca::tools::Debugger::Instance().AddNtuple(
    //  "triplets", "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf:chi2time:ndfTime");
  }
 
  const ca::HitIndex_t ihitl = fIhitL;
  const ca::Hit& hitL        = frWData.Hit(ihitl);
  auto fldValueL             = fField.GetFieldValue(fIstaL, hitL.X(), hitL.Y());
  const auto& target         = fSetup.GetTarget();
  const auto fldValueTarget  = fField.GetPrimVertexField().Get(target.GetX(), target.GetY(), target.GetZ());
 
  fscal c_light = fscal(1.e-5 * kf::defs::SpeedOfLight<double>);
 
  for (size_t i3 = 0; i3 < hitsM.size(); ++i3) {
 
    // prepare data
 
    const ca::HitIndex_t ihitm = hitsM[i3];
    const ca::HitIndex_t ihitr = hitsR[i3];
 
    const ca::Hit& hitM = frWData.Hit(ihitm);
    const ca::Hit& hitR = frWData.Hit(ihitr);
 
    // find the field along the track
 
    auto fldValueM =
      fField.GetFieldValue(fIstaM, hitM.X(), hitM.Y());  // TODO: SZh 24.08.2022: check if fIstaM is in field region
 
    auto fld =
      kf::FieldRegion<fscal>::CopyConvert(kf::FieldRegion<fvec>(fldValueL,                                         //
                                                                fldValueM,                                         //
                                                                fField.GetFieldValue(fIstaR, hitR.X(), hitR.Y()))  //
      );
 
    //fscal hT = target.GetZ()[0] - hitM.Z();
    fscal hL = hitL.Z() - hitM.Z();
    fscal hR = hitR.Z() - hitM.Z();
 
    fscal JxL4 = 0.;
    fscal JyL4 = 0.;
    fscal msX = 0., msY = 0.;
 
    {  // left hit
      fscal tx = (hitL.X() - hitM.X()) / hL;
      fscal ty = (hitL.Y() - hitM.Y()) / hL;
      fscal xx = tx * tx;
      fscal yy = ty * ty;
      fscal xy = tx * ty;
      fscal cL = c_light * sqrt(1.f + xx + yy);
 
      auto [Sx, Sy, Sz] = fld.GetDoubleIntegrals(hitM.X(), hitM.Y(), hitM.Z(), hitL.X(), hitL.Y(), hitL.Z());
 
      JxL4 = cL * (Sx * xy - Sy * (xx + 1.f) + Sz * ty);
      JyL4 = cL * (Sx * (yy + 1.f) - Sy * xy - Sz * tx);
 
      fscal a = msM[i3];
      msX     = (1.f + xx) * a;
      msY     = (1.f + yy) * a;
    }
 
    fscal JxR4 = 0.;
    fscal JyR4 = 0.;
    {  // right hit
      fscal tx = (hitR.X() - hitM.X()) / hR;
      fscal ty = (hitR.Y() - hitM.Y()) / hR;
      fscal xx = tx * tx;
      fscal yy = ty * ty;
      fscal xy = tx * ty;
      fscal cL = c_light * sqrt(1.f + xx + yy);
 
      auto [Sx, Sy, Sz] = fld.GetDoubleIntegrals(hitM.X(), hitM.Y(), hitM.Z(), hitR.X(), hitR.Y(), hitR.Z());
 
      JxR4 = cL * (Sx * xy - Sy * (xx + 1.f) + Sz * ty);
      JyR4 = cL * (Sx * (yy + 1.f) - Sy * xy - Sz * tx);
    }
 
    fscal qp, Cqp;
    {  // fit q/p
      fscal d = hR * JxL4 - hL * JxR4;
      qp      = (hL * (hitM.X() - hitR.X()) + hR * (hitL.X() - hitM.X())) / d;
      Cqp = (hL * hL * hitR.dX2() + (hL - hR) * (hL - hR) * hitM.dX2() + hR * hR * hitL.dX2() + hL * hL * hR * hR * msX)
            / d / d;
    }
 
    fscal x = hitL.X();
    fscal y = hitL.Y();
    fscal z = hitL.Z();
 
    fscal dz01 = 1. / (hitM.Z() - hitL.Z());
 
    fscal tx    = (hitM.X() - hitL.X()) * dz01;
    fscal Ctx   = 1.;
    fscal chi2x = 0.;
 
    fscal ty    = (hitM.Y() - hitL.Y()) * dz01;
    fscal Cty   = 1.;
    fscal chi2y = 0.;
 
 
    {  // get chi2 from y fit
      fscal yL = hitL.Y() - JyL4 * qp;
      fscal lS = hitL.dY2() + JyL4 * JyL4 * Cqp;
 
      fscal yM = hitM.Y();
      fscal mS = hitM.dY2();
 
      fscal yR = hitR.Y() - JyR4 * qp;
      fscal rS = hitR.dY2() + JyR4 * JyR4 * Cqp;
 
      fscal zeta = (hL * (yM - yR) + hR * (yL - yM));
      fscal c    = hL * hL * rS + (hL - hR) * (hL - hR) * mS + hR * hR * lS + hR * hR * hL * hL * msY;
 
      chi2y = zeta * zeta / c;
    }
 
    fscal chi2Tgt = 0.;
 
    {  // chi2 to the target
 
      auto fldTarget = kf::FieldRegion<fscal>::CopyConvert(kf::FieldRegion<fvec>(fldValueTarget, fldValueL, fldValueM));
 
      fscal xTarget = target.GetX()[0];
      fscal yTarget = target.GetY()[0];
      fscal zTarget = target.GetZ()[0];
 
      fscal h  = (zTarget - z);
      fscal xx = tx * tx;
      fscal yy = ty * ty;
      fscal xy = tx * ty;
 
      fscal cL = c_light * sqrt(fscal(1.) + xx + yy);
 
      fscal Sx, Sy, Sz;
      std::tie(Sx, Sy, Sz) = fldTarget.GetDoubleIntegrals(x, y, z, xTarget, yTarget, zTarget);
 
      fscal JxTgt4 = cL * (Sx * xy - Sy * (xx + fscal(1.)) + Sz * ty);
      fscal JyTgt4 = cL * (Sx * (yy + fscal(1.)) - Sy * xy - Sz * tx);
 
      const kf::MeasurementXy<fvec>& m = frWData.TargetMeasurement();
 
      fscal zeta0 = x + h * tx + JxTgt4 * qp - xTarget;
      fscal zeta1 = y + h * ty + JyTgt4 * qp - yTarget;
 
      fscal S00 = m.Dx2()[0] + JxTgt4 * JxTgt4 * Cqp;
      fscal S10 = m.Dxy()[0] + JxTgt4 * JyTgt4 * Cqp;
      fscal S11 = m.Dy2()[0] + JyTgt4 * JyTgt4 * Cqp;
 
      fscal si = fscal(1.) / (S00 * S11 - S10 * S10);
 
      fscal S00tmp = S00;
      S00          = si * S11;
      S10          = -si * S10;
      S11          = si * S00tmp;
 
      chi2Tgt = zeta0 * zeta0 * S00 + fscal(2.) * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
    }
 
 
    // TODO: SG: normalize chi2, separate cuts on time and space
 
    fscal chi2 = chi2x + chi2y;
 
    if (chi2 > frWData.CurrentIteration()->GetTripletFinalChi2Cut()) {
      continue;
    }
 
    if (!std::isfinite(chi2) || chi2 < 0) {
      continue;
    }
 
    chi2 += chi2Tgt;
 
    tripletsOut.emplace_back(ihitl, ihitm, ihitr, fIstaL, fIstaM, fIstaR, 0, 0, 0, chi2, qp, Cqp, tx, Ctx, ty, Cty,
                             isMomentumFitted);
  }  //i3
}  // FitTriplets
 
 
void TripletConstructorSW::CollectHits(Vector<ca::HitIndex_t>& collectedHits,
                                       const SearchWindowMap::SearchWindow& area2D, fscal time, fscal timeRange,
                                       const int iSta, const fscal areaExtension, const int maxNhits, const McMatch& mc)
{
  collectedHits.clear();
  collectedHits.reserve(maxNhits);
 
  const auto& sta  = fSetup.GetActiveLayer(iSta);
  const auto& grid = frWData.Grid(iSta);
 
  //const fvec maxDZ = frWData.CurrentIteration()->GetMaxDZ();
 
  fscal areaX  = 0.5 * (area2D.xMin + area2D.xMax);
  fscal areaDx = 0.5 * (area2D.xMax - area2D.xMin);
  fscal areaY  = 0.5 * (area2D.yMin + area2D.yMax);
  fscal areaDy = 0.5 * (area2D.yMax - area2D.yMin);
  areaDx       = areaExtension * areaDx;
  areaDy       = areaExtension * areaDy;
 
  fscal dxExtended = areaDx + grid.GetMaxRangeX();  // + maxDZ * kf::utils::fabs(T.Tx()))[0];
  fscal dyExtended = areaDy + grid.GetMaxRangeY();  // + maxDZ * kf::utils::fabs(T.Ty()))[0];
  ca::GridArea area(grid, areaX, areaY, dxExtended, dyExtended);
 
  if constexpr (fDebugCollectHits) {
    LOG(info) << "search area: x " << areaX << " +-" << areaDx << ", y " << areaY << " +-" << areaDy;
  }
  if (fParameters.GetConfig().DevIgnoreHitSearchAreas()) {
    area.DoLoopOverEntireGrid();
  }
 
  // loop over station hits (index incremented in condition)
  for (ca::HitIndex_t ih = 0; area.GetNextObjectId(ih) && ((int) collectedHits.size() < maxNhits);) {
    if (frWData.IsHitSuppressed(ih)) {
      continue;
    }
 
    const ca::Hit& hit = frWData.Hit(ih);
    if constexpr (fDebugCollectHits) {
      LOG(info) << "hit in the search area: " << hit.ToString();
      LOG(info) << "   check the hit.. ";
    }
 
 
    if (!mc.IsEmpty()) {  // match via MC
      if (mc != fFramework.GetMcMatchForCaHit(hit.Id())) {
        if constexpr (fDebugCollectHits) {
          LOG(info) << "   hit mc does not match";
        }
        continue;
      }
    }
 
    {
      float dx = hit.X() - areaX;
      float dy = hit.Y() - areaY;
      if (fabs(dx) > areaDx + hit.RangeX() || fabs(dy) > areaDy + hit.RangeY()) {
        continue;
      }
    }
 
    // check time-boundaries
 
    //TODO: move hardcoded cuts to parameters
    if (sta.IsTimeMeasured()) {
      if (fabs(time - hit.T()) > 1.4 * (timeRange + hit.RangeT())) {
        if constexpr (fDebugCollectHits) {
          LOG(info) << "   hit time does not match";
        }
        continue;
      }
    }
 
    collectedHits.push_back(ih);
 
  }  // loop over the hits in the area
}
 
 
int TripletConstructorSW::FindClosestHitWithMc(const SearchWindowMap::SearchWindow& area2D, const int iSta,
                                               const McMatch& mc)
{
 
  if (mc.IsEmpty()) {
    return -1;
  }
 
  const auto& grid = frWData.Grid(iSta);
 
  fscal areaX = 0.5 * (area2D.xMin + area2D.xMax);
  fscal areaY = 0.5 * (area2D.yMin + area2D.yMax);
 
  ca::GridArea area(grid, areaX, areaY, 10000., 10000.);
 
  area.DoLoopOverEntireGrid();
 
  int iret         = -1;
  double bestDist2 = 1.e10;
  // loop over station hits (index incremented in condition)
  for (ca::HitIndex_t ih = 0; area.GetNextObjectId(ih);) {
    if (frWData.IsHitSuppressed(ih)) {
      continue;
    }
 
    const ca::Hit& hit = frWData.Hit(ih);
 
    if (mc == fFramework.GetMcMatchForCaHit(hit.Id())) {
      double dx    = areaX - hit.X();
      double dy    = areaY - hit.Y();
      double dist2 = dx * dx + dy * dy;
      if (dist2 < bestDist2) {
        bestDist2 = dist2;
        iret      = ih;
      }
    }
  }
  return iret;
}