CaSearchWindowMapContainerFactory.h

Go to the documentation of this file.

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

 
#pragma once
 
#include "CaConfig.h"
#include "CaDefs.h"
#include "CaSearchWindowMapContainer.h"
#include "KfFieldRegion.h"
#include "KfSetup.h"
#include "KfTrackKalmanFilter.h"
 
#include <vector>
 
#include <fmt/format.h>
 
namespace cbm::algo::ca
{
  class SearchWindowMapContainerFactory {
   public:
    template<typename F>
    static SearchWindowMapContainer Create(const kf::Setup<F>& actSetup,  //
                                           const kf::FieldFn_t& fieldFn,  //
                                           const ca::Config& config)
    {
      using kf::utils::simd::Cast;
      SearchWindowMapContainer res;
 
      constexpr int MaxNofTripletGaps = SearchWindowMapContainer::kMaxNofTripletGaps;
 
      // FIXME: Define these parameters in ca::Config::SwFactory
      constexpr int BinsX{10};
      constexpr int BinsY{20};
      constexpr double MaxExtrapolationStep{10.};  // [cm]
 
      auto field    = kf::FieldRegion<double>(actSetup.GetField().GetFieldType(), fieldFn);
      auto target   = kf::Target<double>(actSetup.GetTarget());
      auto stations = std::vector<kf::ActiveLayer<double>>{};
      stations.reserve(actSetup.GetNofLayers());
      for (const auto& layer : actSetup.GetActiveLayers()) {  // Note: copy layers in double precision
        stations.emplace_back(layer);
      }
 
      {
        // Validate input
        // FIXME: maybe move to another function (Validate(config, actSetup, ....))
        // 1. Check, if first station is far enough from the target
        double zFstSta = Cast<F, double>(actSetup.GetActiveLayer(0).GetZref());
        if (zFstSta - target.GetZ() < 1.e-2) {
          throw std::runtime_error(fmt::format("The first station (z={}) is too close to or before the target (z={})",
                                               zFstSta, target.GetZ()));
        }
      }
 
      // Doublet search window maps
      for (int iIter = 0; iIter < config.GetNofIterations(); ++iIter) {
        const auto& iter = config.GetIteration(iIter);
        for (int iStaM = 0; iStaM < actSetup.GetNofLayers(); ++iStaM) {
          const auto& staM = stations[iStaM];  // middle station
          for (int iGap = 0; iGap < MaxNofTripletGaps; ++iGap) {
            auto& swMap = res.DoubletSw(iIter, iStaM, iGap);
            swMap.Init(target.GetX(), target.GetY(), target.GetZ(),  //
                       staM.GetZref(),                               //
                       staM.GetXmin(), staM.GetXmax(), BinsX,        //
                       staM.GetYmin(), staM.GetYmax(), BinsY);
 
            // left station
            int iStaL = iStaM - iGap - 1;
            if (iStaL < 0) continue;
            if (iGap > iter.GetMaxStationGap()) continue;
            const auto& staL = stations[iStaL];
 
            // 0.3 -> 0.333 ?
            double stepXm = 0.3 * (staM.GetXmax() - staM.GetXmin()) / BinsX;
            double stepYm = 0.3 * (staM.GetYmax() - staM.GetYmin()) / BinsY;
            for (double xm = staM.GetXmin(); xm <= staM.GetXmax(); xm += stepXm) {
              for (double ym = staM.GetYmin(); ym <= staM.GetYmax(); ym += stepYm) {
                // Extrapolate a track from the target area through the left station
                const double dzm = staM.GetZref() - target.GetZ();
 
                kf::TrackKalmanFilter<double, kf::FilterSettings<kf::LinearizationQp>> fit;
                fit.SetMaxExtrapolationStep(MaxExtrapolationStep);
                fit.SetParticleMass(iter.GetElectronFlag() ? constants::phys::ElectronMass : constants::phys::MuonMass);
                fit.Linearization().qp = 0.;  // zero linearization with qp = 0
 
                // Initial track parameters (at target.GetZ())
                auto& T  = fit.Tr();
                T.X()    = target.GetX();
                T.Y()    = target.GetY();
                T.Z()    = target.GetZ();
                T.Tx()   = (xm - target.GetX()) / dzm;
                T.Ty()   = (ym - target.GetY()) / dzm;
                T.Qp()   = 0.;
                T.Time() = 0.;
                T.Vi()   = constants::phys::SpeedOfLightInv;
 
                const double slopeVar = iter.GetMaxSlopePV() * iter.GetMaxSlopePV() / 9.;
                const double maxQp    = iter.GetMaxQp();
                const double qpVar    = maxQp * maxQp / 9.;
 
                T.ResetErrors(0., 0., slopeVar, slopeVar, qpVar, 0., 1.e+10);
                T.InitVelocityRange(1. / iter.GetMaxQp());
                T.C00() = iter.GetTargetPosSigmaX() * iter.GetTargetPosSigmaX();
                T.C10() = 0.;
                T.C11() = iter.GetTargetPosSigmaY() * iter.GetTargetPosSigmaY();
 
                // Extrapolate to the left station (to staL.GetZref())
                fit.ExtrapolateLineInField(staL.GetZref(), field);
 
                // Add the left hit to the covariance matrix (NOTE: initial hit errors are ignored)
                fit.FilterXY(kf::MeasurementXy<double>(T.X(), T.Y(), 1.e-8, 1.e-8, 0., 1., 1.));
 
                // Add material effects
                fit.Linearization().qp = maxQp;
                fit.MultipleScattering(actSetup.GetMaterial(iStaL).GetThicknessX0(T.GetX(), T.GetY()));
                fit.Linearization().qp = 0.;
 
                // Extrapolate to the middle station (to staM.GetZref())
                fit.ExtrapolateLineInField(staM.GetZref(), field);
 
                double dxM = 3.5 * std::sqrt(T.C00());
                double dyM = 3.5 * std::sqrt(T.C11());
                assert(std::isfinite(dxM));
                assert(std::isfinite(dyM));
                swMap.GetMap().Update(xm, ym, dxM, dyM, 0.);
              }  // loop over ym
            }    // loop over xm
          }      // loop over iGap
        }        // loop over iStaM
      }          // loop over iIter
 
 
      // Triplet search window maps
      for (int iIter = 0; iIter < config.GetNofIterations(); ++iIter) {
        const auto& iter = config.GetIteration(iIter);
        for (int iStaR = 0; iStaR < actSetup.GetNofLayers(); ++iStaR) {
          const auto& staR = stations[iStaR];
          for (int iGapL = 0; iGapL < MaxNofTripletGaps; ++iGapL) {
            for (int iGapM = 0; iGapM < MaxNofTripletGaps; ++iGapM) {
              auto& swMap = res.TripletSw(iIter, iStaR, iGapL, iGapM);
              swMap.Init(target.GetX(), target.GetY(), target.GetZ(),  //
                         staR.GetZref(),                               //
                         staR.GetXmin(), staR.GetXmax(), BinsX,        //
                         staR.GetYmin(), staR.GetYmax(), BinsY);
 
              int iStaL = iStaR - iGapL - 2;
              int iStaM = iStaR - iGapM - 1;
              if (iStaL < 0 || iStaL >= iStaM) continue;
              if (iGapL > iter.GetMaxStationGap() || iGapM > iter.GetMaxStationGap()) continue;
 
              const auto& staL = stations[iStaL];
              const auto& staM = stations[iStaM];
 
              // 0.3 -> 0.333 ?
              double stepXr = 0.3 * (staR.GetXmax() - staR.GetXmin()) / BinsX;
              double stepYr = 0.3 * (staR.GetYmax() - staR.GetYmin()) / BinsY;
 
              for (double xr = staR.GetXmin(); xr <= staR.GetXmax(); xr += stepXr) {
                for (double yr = staR.GetYmin(); yr <= staR.GetYmax(); yr += stepYr) {
                  // track from the target area
                  const double dzr = staR.GetZref() - target.GetZ();
 
                  kf::TrackKalmanFilter<double, kf::FilterSettings<kf::LinearizationQp>> fit;
                  fit.SetMaxExtrapolationStep(MaxExtrapolationStep);
                  fit.SetParticleMass(iter.GetElectronFlag() ? constants::phys::ElectronMass
                                                             : constants::phys::MuonMass);
                  fit.Linearization().qp = 0.;
 
                  // Initial track parameters (z = target.GetZ())
                  auto& T  = fit.Tr();
                  T.X()    = target.GetX();
                  T.Y()    = target.GetY();
                  T.Z()    = target.GetZ();
                  T.Tx()   = (xr - target.GetX()) / dzr;
                  T.Ty()   = (yr - target.GetY()) / dzr;
                  T.Qp()   = 0.;
                  T.Time() = 0.;
                  T.Vi()   = constants::phys::SpeedOfLightInv;
 
                  const double slopeVar = iter.GetMaxSlopePV() * iter.GetMaxSlopePV() / 9.;
                  const double maxQp    = iter.GetMaxQp();
                  const double qpVar    = maxQp * maxQp / 9.;
 
                  T.ResetErrors(0., 0., slopeVar, slopeVar, qpVar, 0., 1.e+10);
                  T.InitVelocityRange(1. / maxQp);
                  T.C00() = iter.GetTargetPosSigmaX() * iter.GetTargetPosSigmaX();
                  T.C10() = 0.;
                  T.C11() = iter.GetTargetPosSigmaY() * iter.GetTargetPosSigmaY();
 
                  // Extrapolate to the left station (z = staL.GetZref())
                  fit.ExtrapolateLineInField(staL.GetZref(), field);
 
                  // Add the left hit XY-measurement without errors
                  fit.FilterXY(kf::MeasurementXy<double>(T.X(), T.Y(), 1.e-8, 1.e-8, 0., 1., 1.));
 
                  // Add material effects
                  fit.Linearization().qp = maxQp;
                  fit.MultipleScattering(actSetup.GetMaterial(iStaL).GetThicknessX0(T.GetX(), T.GetY()));
                  fit.Linearization().qp = 0.;
                  T.Qp() = 0.;
 
                  // Extrapolate to the middle station (z = staM.GetZref())
                  fit.ExtrapolateLineInField(staM.GetZref(), field);
 
                  // Add the middle hit without errors to the covariance matrix
                  fit.FilterXY(kf::MeasurementXy<double>(T.X(), T.Y(), 1.e-8, 1.e-8, 0., 1., 1.));
 
                  // Add material effects
                  fit.Linearization().qp = maxQp;
                  double radThick = actSetup.GetMaterial(iStaM).GetThicknessX0(T.GetX(), T.GetY());
                  {
                    double dz = staR.GetZref() - staM.GetZref();
                    double rX = T.GetX() + dz * T.Tx();
                    double rY = T.GetY() + dz * T.Ty();
                    radThick += actSetup.GetMaterial(iStaR).GetThicknessX0(rX, rY);
                  }
                  fit.MultipleScattering(radThick);
 
                  double ms;
                  {
                    double tx = T.Tx();
                    double ty = T.Ty();
                    double t  = sqrt(1. + tx * tx + ty * ty);
                    double lg = 0.0136 * (1. + 0.038 * log(radThick * t));
                    if (lg < 0.) lg = 0.;
                    double qp   = maxQp / 7.;  // ??? correction is needed some reason. Investigate.
                    double s0   = lg * qp * t;
                    double mass = iter.GetElectronFlag() ? constants::phys::ElectronMass : constants::phys::MuonMass;
                    ms          = (1. + mass * mass * qp * qp) * s0 * s0 * t * radThick;
                  }
                  fit.Linearization().qp = 0.;
                  T.Qp() = 0.;
 
                  // Extrapolate to the right station (z = staR.GetZref())
                  fit.ExtrapolateLineInField(staR.GetZref(), field);
 
                  double dxR = 3.5 * std::sqrt(T.C00());
                  double dyR = 3.5 * std::sqrt(T.C11());
                  assert(std::isfinite(dxR));
                  assert(std::isfinite(dyR));
                  swMap.GetMap().Update(xr, yr, dxR, dyR, ms);
                }  // loop over yr
              }    // loop over xr
            }      // loop over iGapM
          }        // loop over iGapL
        }          // loop over iStaR
      }            // loop over iIter
 
      return res;
    }
  };
}  // namespace cbm::algo::ca