TrTrackFill.cpp

Go to the documentation of this file.

#include "Containers/TrTrack.h"
#include "Utility/Logging/Logging.h"

#include "NtpMaker/NtpTools.h"

#include "bcorr.h"

#include <iostream>

namespace NAIA {
using namespace TrTrack;

static std::array<unsigned long int, numSpans> patterns = {3, 5, 6, 7, 0, 0, 3};
std::array<int, numFits> fitAlgoCode{1, 6, 7};
static bool TrackHasExtHits(TrTrackR *trackPtr, Span span) {
  switch (span) {
  case Span::InnerOnly:
    return true;
  case Span::InnerL1:
    return trackPtr->TestHitLayerJHasXY(1);
  case Span::InnerL9:
    return trackPtr->TestHitLayerJHasXY(9);
  case Span::FullSpan:
    return (trackPtr->TestHitLayerJHasXY(1) && trackPtr->TestHitLayerJHasXY(9));
  case Span::UpperHalfInner:
  case Span::LowerHalfInner:
  case Span::InnerNoMS:
    return true;
  default:
    throw std::runtime_error("Span not supported");
  }
}

static constexpr std::array<double, 9> powersOfTen = {1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8};

void ComputePatterns(TrTrackR *trackPtr) {
  patterns[Span::UpperHalfInner] = patterns[Span::LowerHalfInner] = 0;

  for (int il = 2; il < 9; il++) {
    if (!trackPtr->TestHitLayerJ(il))
      continue;

    if (il < 7) {
      patterns[Span::UpperHalfInner] += static_cast<unsigned long int>(9 * powersOfTen[il - 1]);
    }

    if (il > 2) {
      patterns[Span::LowerHalfInner] += static_cast<unsigned long int>(9 * powersOfTen[il - 1]);
    }
  }
}

bool TrTrackBaseData::_FitIDExists(Fit fit, Span span) const {
  auto itr = _fit_ID.find(fit);
  if (itr == end(_fit_ID)) {
    return false;
  }

  auto fit_ID_itr = itr->second.find(span);
  if (fit_ID_itr == end(itr->second)) {
    return false;
  }

  return (fit_ID_itr->second > 0);
}

Span TrTrackBaseData::_GetBestSpan(Fit fit) const {
  Span bestspan = _FitIDExists(fit, Span::FullSpan)
                      ? Span::FullSpan
                      : (_FitIDExists(fit, Span::InnerL1) > 0
                             ? Span::InnerL1
                             : (_FitIDExists(fit, Span::InnerL9) > 0 ? Span::InnerL9 : Span::InnerOnly));
  return bestspan;
}

struct ReFitParameters {
  int refit;
  float mass;
  int charge;
};

ReFitParameters InitTrkReFit(float inner_charge, bool _isMC, float production_charge) {
  ReFitParameters pars{};

  pars.charge = static_cast<int>(inner_charge + 0.5f);
  if (pars.charge < 1)
    pars.charge = 1;
  pars.mass = TrFit::Mproton; // Proton
  if (pars.charge >= 2)
    pars.mass = TrFit::Mhelium / 2.0f * pars.charge; // Helium+Ion

  if (!_isMC) { // Data
    if (pars.charge == static_cast<int>(production_charge + 0.5f))
      pars.refit = 0;
    else
      pars.refit = 3;
  } else {
    pars.refit = 3; // MC
  }
  return pars;
}

bool TrTrackBaseData::Fill(TrTrackR *trackPtr, TrTrackFillType type) {

  // some helper default variables
  int temp_fit_ID = 0;
  if (type == TrTrackFillType::Standalone) {
    _beta = 1.0f;
    temp_fit_ID = -1;
  }

  // New Charge-refit Logic from Qi
  float tmp_inner_charge = trackPtr->GetQ_all(_beta, temp_fit_ID).Mean;
  if (tmp_inner_charge > 1.5) { // more accurate HighZ charge
    if (trackPtr->GetQYJ_all(2, _beta, temp_fit_ID).InnerQ > 8.5)
      tmp_inner_charge = trackPtr->GetQYJ_all(2, _beta, temp_fit_ID).InnerQ;
  }

  FitCharge = tmp_inner_charge;

  auto fitPars = InitTrkReFit(tmp_inner_charge, _isMC, FitCharge);
  _tmpCharge = fitPars.charge;

  ComputePatterns(trackPtr);

  // Rigidity
  for (auto fit : fitTypes) {
    // Kalman with electron mass hyp. to be done in dedicated fill method
    if (fit == Fit::KalmanElectron)
      continue;

    for (auto span : spanTypes) {
      if (type == TrTrackFillType::Standalone) {
        switch (span) {
        case Span::LowerHalfInner:
        case Span::UpperHalfInner:
        case Span::InnerL9:
        case Span::InnerNoMS:
          // for standalone we only consider InnerOnly, InnerL1 and FullSpan
          continue;
          break;
        default:
          break;
        }
      }

      if (TrackHasExtHits(trackPtr, span)) {
        Logging::MuteOutput(); // this can get quite verbose...
        _fit_ID[fit][span] = trackPtr->iTrTrackPar(fitAlgoCode[fit], static_cast<int>(patterns[span]),
                                                   30 + fitPars.refit, fitPars.mass, _tmpCharge);
        Logging::UnmuteOutput();
        if (_fit_ID[fit][span] > 0) {
          double rigidity;
          AMSPlane pl(TVector3(0, 0, 0), TVector3(1, 0, 0), TVector3(0, 1, 0));
          trackPtr->Interpolate(pl, rigidity, _fit_ID[fit][span]);
          Rigidity[fit][span] = static_cast<float>(rigidity);
          TrChiSq[fit][span][Side::X] =
              static_cast<float>(trackPtr->GetChisqX(_fit_ID[fit][span]) / trackPtr->GetNdofX(_fit_ID[fit][span]));
          TrChiSq[fit][span][Side::Y] =
              static_cast<float>(trackPtr->GetChisqY(_fit_ID[fit][span]) / trackPtr->GetNdofY(_fit_ID[fit][span]));
        }
      }
    }
  }
  trackPtr->RecalcHitCoordinates(0, false, false, _tmpCharge); // using default fit and inner tracker charge to restore
                                                               // coordinate

  int base_fit_ID = FitIDExists(Fit::Choutko, Span::InnerOnly) ? _fit_ID[Fit::Choutko][Span::InnerOnly] : -1;
  if (type == TrTrackFillType::Standalone) {
    base_fit_ID = -1;
  }

  // Charge
  //  - Inner tracker
  Charge[ChargeRecoType::STD] = trackPtr->GetQ_all(_beta, base_fit_ID).Mean;
  Charge[ChargeRecoType::HL] = trackPtr->GetQH_all(2, _beta, base_fit_ID).Mean;
  // Charge[ChargeRecoType::YJ] does not support "global" charge;
  ChargeRMS[ChargeRecoType::STD] = trackPtr->GetQ_all(_beta, base_fit_ID).RMS;
  ChargeRMS[ChargeRecoType::HL] = trackPtr->GetQH_all(2, _beta, base_fit_ID).RMS;
  // ChargeRMS[ChargeRecoType::YJ] does not support "global" charge;
  InnerCharge[ChargeRecoType::STD] = trackPtr->GetInnerQ_all(_beta, base_fit_ID).Mean;
  InnerCharge[ChargeRecoType::HL] = trackPtr->GetInnerQH_all(2, _beta, base_fit_ID).Mean;
  InnerCharge[ChargeRecoType::YJ] = trackPtr->GetQYJ_all(2, _beta, base_fit_ID).InnerQ;
  InnerChargeRMS[ChargeRecoType::STD] = trackPtr->GetInnerQ_all(_beta, base_fit_ID).RMS;
  InnerChargeRMS[ChargeRecoType::HL] = trackPtr->GetInnerQH_all(2, _beta, base_fit_ID).RMS;
  InnerChargeRMS[ChargeRecoType::YJ] = trackPtr->GetQYJ_all(2, _beta, base_fit_ID).InnerQRMS;

  // geometry and position
  for (int ijl = 0; ijl < 9; ijl++) {
    if (!trackPtr->GetHitLJ(ijl + 1))
      continue;
    AMSPlaneM plm = trackPtr->GetHitCooLJN(ijl + 1);
    TVector3 mglob = plm.getMGlobal();
    if (trackPtr->GetHitLJ(ijl + 1)->GetXCluster())
      TrTrackHitPos[ijl][0] = static_cast<float>(mglob.x());
    if (trackPtr->GetHitLJ(ijl + 1)->GetYCluster())
      TrTrackHitPos[ijl][1] = static_cast<float>(mglob.y());
    if (trackPtr->GetHitLJ(ijl + 1))
      TrTrackHitPos[ijl][2] = static_cast<float>(mglob.z());
  }

  if (type == TrTrackFillType::SecondTrack) {
    return true;
  }

  //  - single layers
  for (int ijl = 0; ijl < 9; ijl++) {
    Fit fit = Fit::Choutko;
    Span span = Span::InnerOnly;

    if (trackPtr->GetHitLJ(ijl + 1)) {
      LayerChargeStatus[ijl] = trackPtr->GetHitLJ(ijl + 1)->GetQStatus();

      LayerChargeXY[ijl][ChargeRecoType::STD] = trackPtr->GetHitLJ(ijl + 1)->GetQ(2, _beta, Rigidity[fit][span]);
      LayerChargeXY[ijl][ChargeRecoType::HL] = trackPtr->GetLayerJQH(ijl + 1, 2, _beta, _fit_ID[fit][span]);
      LayerChargeXY[ijl][ChargeRecoType::YJ] = trackPtr->GetLayerQYJ(ijl + 1, 2, _beta, _fit_ID[fit][span]);
    }
  }

  for (const auto fitPos : fitPositionHeights) {
    // Define a Plane from the Hits Coordinates
    AMSPlaneM plm_res;
    TVector3 mglob;
    if (fitPos <= TrTrack::FitPositionHeight::Layer9) { // Just for the Layers
      plm_res = trackPtr->GetHitCooLJN(fitPos + 1);
      mglob = plm_res.getMGlobal();
    }
    for (auto fit : fitTypes) {
      for (auto span : spanTypes) {
        if (type == TrTrackFillType::Standalone) {
          switch (span) {
          case Span::LowerHalfInner:
          case Span::UpperHalfInner:
          case Span::InnerL9:
          case Span::InnerNoMS:
            // for standalone we only consider InnerOnly, InnerL1 and FullSpan
            continue;
            break;
          default:
            break;
          }
        }

        if (FitIDExists(fit, span)) {
          int fit_ID = _fit_ID.at(fit).at(span);
          if (fit_ID < 0)
            continue;

          double rigidity;

          // Define a Plane from a given Nominal HeightZ to store the FitPositions
          AMSPlane plm_fitpos(TVector3(0, 0, fitPositionHeightZ[fitPos]), TVector3(1, 0, 0), TVector3(0, 1, 0));
          trackPtr->Interpolate(plm_fitpos, rigidity, fit_ID);
          TVector3 pglob_fitpos = plm_fitpos.getP();
          TrTrackFitPos[fitPos][fit][span][Side::X] = static_cast<float>(pglob_fitpos.x());
          TrTrackFitPos[fitPos][fit][span][Side::Y] = static_cast<float>(pglob_fitpos.y());
        }
      }
    }
  }

  return true;
}

bool TrTrackPlusData::Fill(TrTrackR *trackPtr) {
  const auto &_beta = trkBase->_beta;
  const auto &_tmpCharge = trkBase->_tmpCharge;

  if (trkBase->InnerCharge.empty())
    return false;

  auto *pev = AMSEventR::Head();

  int refit = ((trkBase->_isMC) || (_tmpCharge != trackPtr->GetAdvancedFitCharge())) ? 3 : 1;
  float fitMass = (_tmpCharge >= 2) ? TrFit::Mhelium / 2 * _tmpCharge : TrFit::Mproton;

  auto getPartialPattern = [](TrTrackR *track, unsigned int jLayer) {
    unsigned int pattern = 0;
    for (unsigned int ijl = 1; ijl < 10; ijl++) {
      if (track->TestHitLayerJ(ijl) && ijl != jLayer)
        pattern += 9 * powersOfTen[ijl - 1];
    }
    return pattern;
  };

  for (unsigned int ijl = 2; ijl < 9; ijl++) {
    TrackFeetDistance[ijl - 1] = pev->GetTkFeetDist(ijl, trackPtr);
  }

  // Rigidity
  for (auto fit : fitTypes) {

    // Kalman with electron mass hyp. to be done in dedicated fill method
    if (fit == Fit::KalmanElectron)
      continue;

    for (auto span : spanTypes) {
      if (trkBase->FitIDExists(fit, span)) {
        _fit_ID[fit][span] =
            trackPtr->iTrTrackPar(fitAlgoCode[fit], static_cast<int>(patterns[span]), 30 + refit, fitMass, _tmpCharge);
        if (_fit_ID[fit][span] < 0)
          continue;

        double rigidity;
        if (fit == Fit::Kalman) {
          // TOI rigidity for both downward and upward going particles.
          float zpl = (_beta > 0) ? 195 : -195;
          AMSPlane pl_TOI(TVector3(0, 0, zpl), TVector3(1, 0, 0), TVector3(0, 1, 0));
          trackPtr->Interpolate(pl_TOI, rigidity, _fit_ID[fit][span], fitMass, _tmpCharge, 1, (_beta > 0) ? 1 : -1);
          RigidityTOI[fit][span] = static_cast<float>(rigidity);
          // Rigidity at RICH
          AMSPlane pl_RICH(TVector3(0, 0, -70), TVector3(1, 0, 0), TVector3(0, 1, 0));
          trackPtr->Interpolate(pl_RICH, rigidity, _fit_ID[fit][span], fitMass, _tmpCharge, 1);
          RigidityRICH[fit][span] = static_cast<float>(rigidity);
        }
        // Theta and Phi
        AMSPlane pl(TVector3(0, 0, 0), TVector3(1, 0, 0), TVector3(0, 1, 0));
        trackPtr->Interpolate(pl, rigidity, _fit_ID[fit][span]);
        Theta[fit][span] = static_cast<float>(pl.getD().Theta());
        Phi[fit][span] = static_cast<float>(pl.getD().Phi());

        const TrTrackPar &trackPar = trackPtr->gTrTrackPar(_fit_ID[fit][span]);
        InvRigErr[fit][span] = static_cast<float>(trackPar.ErrRinv);
      }
    }
    // using default fit and inner tracker charge to restore coordinate
    trackPtr->RecalcHitCoordinates(0, 0, 0, _tmpCharge);

    // doesn't make sense to try if inner only fit failed...
    if (trkBase->FitIDExists(fit, Span::InnerOnly)) {
      for (unsigned int ijl = 1; ijl < 10; ijl++) {
        if (!trackPtr->TestHitLayerJ(ijl))
          continue;
        auto partialPattern = getPartialPattern(trackPtr, ijl);

        int pfit_ID = trackPtr->iTrTrackPar(fitAlgoCode[fit], partialPattern, 30 + refit, fitMass, _tmpCharge);
        if (pfit_ID < 0)
          continue;

        const TrTrackPar &partialTrackPar = trackPtr->gTrTrackPar(pfit_ID);
        PartialRigidity[ijl - 1][fit] = partialTrackPar.Rigidity;
        PartialInvRigErr[ijl - 1][fit] = partialTrackPar.ErrRinv;
        PartialTrChiSq[ijl - 1][fit][Side::X] = static_cast<float>(partialTrackPar.ChisqX / partialTrackPar.NdofX);
        PartialTrChiSq[ijl - 1][fit][Side::Y] = static_cast<float>(partialTrackPar.ChisqY / partialTrackPar.NdofY);

        // Define a Plane from a given Nominal HeightZ to store the Partial FitPositions
        double rigidity;
        AMSPlane plm_fitpos(TVector3(0, 0, fitPositionHeightZ[ijl - 1]), TVector3(1, 0, 0), TVector3(0, 1, 0));
        trackPtr->Interpolate(plm_fitpos, rigidity, pfit_ID);
        TVector3 pglob_fitpos = plm_fitpos.getP();
        PartialTrTrackFitPos[ijl - 1][fit][Side::X] = static_cast<float>(pglob_fitpos.x());
        PartialTrTrackFitPos[ijl - 1][fit][Side::Y] = static_cast<float>(pglob_fitpos.y());

        // Define a Plane from the Hits Coordinates to calculate the Partial Residuals
        // by using the real Height of the event
        AMSPlaneM plm_res = trackPtr->GetHitCooLJN(ijl + 1);
        TVector3 mglob = plm_res.getMGlobal();
        trackPtr->Interpolate(plm_res, rigidity, pfit_ID);
        TVector3 pglob_res = plm_res.getP();
        PartialTrTrackResidual[ijl - 1][fit][Side::X] = static_cast<float>(mglob.x() - pglob_res.x());
        PartialTrTrackResidual[ijl - 1][fit][Side::Y] = static_cast<float>(mglob.y() - pglob_res.y());
      }
    }
  }
  trackPtr->RecalcHitCoordinates(0, 0, 0, _tmpCharge); // using default fit and inner tracker charge to restore
                                                       // coordinate

  auto dir_fit = TrTrack::Fit::Choutko;
  auto dir_span = trkBase->GetBestSpan(dir_fit);
  if (trkBase->FitIDExists(dir_fit, dir_span)) {
    AMSDir trk_dir{Theta[dir_fit][dir_span], Phi[dir_fit][dir_span]};
    // Once again, I come to you to complain about the worst API ever designed: gbatch
    // 5 LOC when 2 would have been enough
    double tmp;
    AMSEventR::Head()->GetStoermerCutoff(tmp, -1, trk_dir);
    DirectionalStoermerCutoff[0] = tmp;
    AMSEventR::Head()->GetStoermerCutoff(tmp, 1, trk_dir);
    DirectionalStoermerCutoff[1] = tmp;
  }

  //  - single layers
  constexpr Fit base_fit = Fit::Choutko;
  constexpr Span base_span = Span::InnerOnly;
  if (trkBase->FitIDExists(base_fit, base_span)) {
    int base_fit_ID = _fit_ID.at(base_fit).at(base_span);

    for (int ijl = 0; ijl < 9; ijl++) {
      if (trackPtr->GetHitLJ(ijl + 1)) {
        float rigidity = (base_fit_ID > 0) ? trkBase->Rigidity.at(base_fit).at(base_span) : 0;

        LayerCharge[ijl][ChargeRecoType::STD][Side::X] = trackPtr->GetHitLJ(ijl + 1)->GetQ(0, _beta, rigidity);
        LayerCharge[ijl][ChargeRecoType::HL][Side::X] = trackPtr->GetLayerJQH(ijl + 1, 0, _beta, base_fit_ID);
        LayerCharge[ijl][ChargeRecoType::YJ][Side::X] = trackPtr->GetLayerQYJ(ijl + 1, 0, _beta, base_fit_ID);

        if (trackPtr->GetHitLJ(ijl + 1)->GetXCluster())
          LayerEdep[ijl][Side::X] = trackPtr->GetHitLJ(ijl + 1)->GetXCluster()->GetEdep();
        if (trackPtr->GetHitLJ(ijl + 1)->GetYCluster())
          LayerEdep[ijl][Side::Y] = trackPtr->GetHitLJ(ijl + 1)->GetYCluster()->GetEdep();

        LayerCharge[ijl][ChargeRecoType::STD][Side::Y] = trackPtr->GetHitLJ(ijl + 1)->GetQ(1, _beta, rigidity);
        LayerCharge[ijl][ChargeRecoType::HL][Side::Y] = trackPtr->GetLayerJQH(ijl + 1, 1, _beta, base_fit_ID);
        LayerCharge[ijl][ChargeRecoType::YJ][Side::Y] = trackPtr->GetLayerQYJ(ijl + 1, 1, _beta, base_fit_ID);
      }
    }
  }

  for (const auto fitPos : fitPositionHeights) {
    // Define a Plane from the Hits Coordinates
    AMSPlaneM plm_res;
    TVector3 mglob;
    if (fitPos <= TrTrack::FitPositionHeight::Layer9) { // Just for the Layers
      plm_res = trackPtr->GetHitCooLJN(fitPos + 1);
      mglob = plm_res.getMGlobal();
    }
    for (auto fit : fitTypes) {
      for (auto span : spanTypes) {
        if (trkBase->FitIDExists(fit, span)) {
          int fit_ID = _fit_ID.at(fit).at(span);
          if (fit_ID < 0)
            continue;
          double rigidity;

          // Define a Plane from a given Nominal HeightZ to store the FitPositions
          AMSPlane plm_fitpos(TVector3(0, 0, fitPositionHeightZ[fitPos]), TVector3(1, 0, 0), TVector3(0, 1, 0));
          trackPtr->Interpolate(plm_fitpos, rigidity, fit_ID);
          TVector3 pglob_fitpos = plm_fitpos.getP();

          // From the Plane defined by the Hits Cooordinates, calculate the Residuals using the real Height associated
          // with the event
          if (fitPos <= TrTrack::FitPositionHeight::Layer9) { // Just for the Layers
            if (!trackPtr->TestHitLayerJ(fitPos + 1))
              continue;
            trackPtr->Interpolate(plm_res, rigidity, fit_ID);
            TVector3 pglob_res = plm_res.getP();
            TrTrackResidual[fitPos][fit][span][Side::X] = static_cast<float>(mglob.x() - pglob_res.x());
            TrTrackResidual[fitPos][fit][span][Side::Y] = static_cast<float>(mglob.y() - pglob_res.y());
          }
        }
      }
    }
  }

  static constexpr std::array<float, 5> window = {0.1, 1., 2., 5., 10.};

  // from dbar code
  for (int jLayer = 1; jLayer <= 9; jLayer++) {

    ClusterSignalRatio[jLayer - 1] = pev->GetTrackerRawSignalRatio(jLayer, 10);

    // identify the clusters used in trtrack
    TrClusterR *trackClusterPtrX = trackPtr->GetHitLJ(jLayer) ? trackPtr->GetHitLJ(jLayer)->GetXCluster() : nullptr;
    TrClusterR *trackClusterPtrY = trackPtr->GetHitLJ(jLayer) ? trackPtr->GetHitLJ(jLayer)->GetYCluster() : nullptr;

    // Get track interpolation to jLayer
    double rigidity;
    AMSPlane plm(TVector3(0, 0, fitPositionHeightZ[jLayer - 1]), TVector3(1, 0, 0), TVector3(0, 1, 0));
    trackPtr->Interpolate(plm, rigidity, _fit_ID[Fit::Kalman][Span::InnerOnly]);
    TVector3 pglob = plm.getP();
    AMSPoint fitcoo = AMSPoint(pglob.x(), pglob.y(), pglob.z());

    // TODO: Check if zero threshold is correct
    static constexpr float noiseThreshold = 0;

    for (int icl = 0; icl < pev->nTrCluster(); icl++) {
      auto *cluster = pev->pTrCluster(icl);

      if (!cluster || cluster->GetLayerJ() != jLayer)
        continue;

      float edep = 1000 * cluster->GetEdep();
      if (cluster == trackClusterPtrX) {
        NClusters[jLayer - 1][DistanceFromTrack::OnTrack][Side::X]++;
        ClustersEdep[jLayer - 1][DistanceFromTrack::OnTrack][Side::X] += edep;
        continue;
      }
      if (cluster == trackClusterPtrY) {
        NClusters[jLayer - 1][DistanceFromTrack::OnTrack][Side::Y]++;
        ClustersEdep[jLayer - 1][DistanceFromTrack::OnTrack][Side::Y] += edep;
        continue;
      }

      if (edep < noiseThreshold)
        continue; // apply a threshold to account for noise

      // Deal with multiplicities: Choose closest to track
      float clusterDistanceToTrack = std::numeric_limits<float>::max();
      for (int im = 0; im < cluster->GetMultiplicity(); im++) {
        AMSPlaneM pl_mult = cluster->GetGCoordN(im);
        TVector3 mglob_mult = pl_mult.getMGlobal();

        clusterDistanceToTrack =
            std::min(clusterDistanceToTrack,
                     fabs(static_cast<float>(fitcoo[cluster->GetSide()] - mglob_mult[cluster->GetSide()])));
      }

      Side side = cluster->GetSide() == 0 ? Side::X : Side::Y;

      NClusters[jLayer - 1][DistanceFromTrack::AllLayer][side]++;
      ClustersEdep[jLayer - 1][DistanceFromTrack::AllLayer][side] += edep;
      if (edep > MaxClusterEdep[jLayer - 1][side]) {
        MaxClusterEdep[jLayer - 1][side] = edep;
        MaxClusterDistance[jLayer - 1][side] = clusterDistanceToTrack;
      }

      // FIXME: I don't really like using an int instead of the right enum
      for (unsigned int iw = 0; iw < window.size(); iw++) {
        if (std::fabs(clusterDistanceToTrack) < window[iw]) {
          NClusters[jLayer - 1][static_cast<DistanceFromTrack>(iw + 1)][side]++;
          ClustersEdep[jLayer - 1][static_cast<DistanceFromTrack>(iw + 1)][side] += edep;
        }
      }
    }
  }
  return true;
}

bool TrTrackBaseData::FillElectronVars(TrTrackR *trackPtr, float zEcalCOG, bool refitKalman) {
  // Rigidity
  auto fit = Fit::KalmanElectron;
  auto span = _GetBestSpan(Fit::Kalman);

  if (TrackHasExtHits(trackPtr, span) && refitKalman) {
    Logging::MuteOutput(); // this can get quite verbose...
    _fit_ID[fit][span] = trackPtr->iTrTrackPar(fitAlgoCode[fit], static_cast<int>(patterns[span]), 32,
                                               static_cast<float>(TrFit::Melectron), 1);
    Logging::UnmuteOutput();
    if (_fit_ID[fit][span] > 0) {
      double rigidity;
      AMSPlane pl(TVector3(0, 0, 0), TVector3(1, 0, 0), TVector3(0, 1, 0));
      trackPtr->Interpolate(pl, rigidity, _fit_ID[fit][span], TrFit::Melectron, 1);
      Rigidity[fit][span] = static_cast<float>(rigidity);
      TrChiSq[fit][span][Side::X] =
          static_cast<float>(trackPtr->GetChisqX(_fit_ID[fit][span]) / trackPtr->GetNdofX(_fit_ID[fit][span]));
      TrChiSq[fit][span][Side::Y] =
          static_cast<float>(trackPtr->GetChisqY(_fit_ID[fit][span]) / trackPtr->GetNdofY(_fit_ID[fit][span]));
    }
  } else {
    return false;
  }

  if (zEcalCOG < std::numeric_limits<float>::max()) {
    for (auto fit : fitTypes) {
      for (auto span : spanTypes) {
        if (FitIDExists(fit, span)) {
          int fit_ID = _fit_ID.at(fit).at(span);
          double rigidity;
          AMSPlane plm(TVector3(0, 0, zEcalCOG), TVector3(1, 0, 0), TVector3(0, 1, 0));
          trackPtr->Interpolate(plm, rigidity, fit_ID);
          TVector3 pglob = plm.getP();
          TrTrackFitPos[FitPositionHeight::EcalCOG][fit][span][Side::X] = static_cast<float>(pglob.x());
          TrTrackFitPos[FitPositionHeight::EcalCOG][fit][span][Side::Y] = static_cast<float>(pglob.y());
        }
      }
    }
  }

  return true;
}

bool TrTrackPlusData::FillElectronVars(TrTrackR *trackPtr, bool refitKalman) {
  const auto &base_fit_ID = trkBase->_fit_ID;
  const auto &_beta = trkBase->_beta;

  {
    auto fit = Fit::KalmanElectron;
    auto span = trkBase->GetBestSpan(Fit::Kalman);
    if (trkBase->FitIDExists(fit, span) && refitKalman) {
      int fit_ID = base_fit_ID.at(fit).at(span);
      // TOI rigidity for both downward and upward going particles.
      double rigidity;
      float zpl = (_beta > 0) ? 195 : -195;
      AMSPlane pl_TOI(TVector3(0, 0, zpl), TVector3(1, 0, 0), TVector3(0, 1, 0));
      trackPtr->Interpolate(pl_TOI, rigidity, fit_ID, TrFit::Melectron, 1, 1, (_beta > 0) ? 1 : -1);
      RigidityTOI[fit][span] = static_cast<float>(rigidity);
      // Rigidity at RICH
      AMSPlane pl_RICH(TVector3(0, 0, -70), TVector3(1, 0, 0), TVector3(0, 1, 0));
      trackPtr->Interpolate(pl_RICH, rigidity, fit_ID, TrFit::Melectron, 1, 1);
      RigidityRICH[fit][span] = static_cast<float>(rigidity);
      // Theta and Phi
      AMSPlane pl(TVector3(0, 0, 0), TVector3(1, 0, 0), TVector3(0, 1, 0));
      trackPtr->Interpolate(pl, rigidity, fit_ID, TrFit::Melectron, 1);
      Theta[fit][span] = static_cast<float>(pl.getD().Theta());
      Phi[fit][span] = static_cast<float>(pl.getD().Phi());

      const TrTrackPar &trackPar = trackPtr->gTrTrackPar(fit_ID);
      InvRigErr[fit][span] = static_cast<float>(trackPar.ErrRinv);
    }
  }

  return true;
}
} // namespace NAIA