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 <cmath>
#include <iostream>
 
namespace NAIA {
using namespace TrTrack;
 
static std::array<unsigned long int, numSpans> patterns = {3, 5, 6, 7, 0, 0};
std::array<int, numFits> fitAlgoCode{1, 6, 7, 6, 11, 17};
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:
    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] += std::lround(9 * powersOfTen[il - 1]);
    }
 
    if (il > 2) {
      patterns[Span::LowerHalfInner] += std::lround(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 = 1;
    else
      pars.refit = 3;
  } else {
    pars.refit = 3; // MC
  }
  return pars;
}
 
namespace {
// these never change so let's make them static and cache them
AMSPlane pl_TOI_top(TVector3(0, 0, 195), TVector3(1, 0, 0), TVector3(0, 1, 0));
AMSPlane pl_TOI_bot(TVector3(0, 0, -195), TVector3(1, 0, 0), TVector3(0, 1, 0));
AMSPlane pl_RICH(TVector3(0, 0, -70), TVector3(1, 0, 0), TVector3(0, 1, 0));
 
// FIXME: this is horrible. I know.
AMSPlane pl_TrCenter(TVector3(0, 0, 0), TVector3(1, 0, 0), TVector3(0, 1, 0));
std::array<AMSPlane, numFitPositionHeights> pl_TrPlanes{
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[0]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[1]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[2]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[3]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[4]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[5]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[6]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[7]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[8]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[9]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[10]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[11]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[12]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[13]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[14]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[15]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
    AMSPlane(TVector3(0, 0, fitPositionHeightZ[16]), TVector3(1, 0, 0), TVector3(0, 1, 0)),
};
} // namespace
 
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;
 
  // more accurate HighZ charge
  if (tmp_inner_charge > 1.5) {
    auto tmp_inner_charge_YJ = trackPtr->GetQYJ_all(2, _beta, temp_fit_ID).InnerQ;
    if (tmp_inner_charge_YJ > 8.5)
      tmp_inner_charge = tmp_inner_charge_YJ;
  }
 
  FitCharge = tmp_inner_charge;
 
  auto fitPars = InitTrkReFit(tmp_inner_charge, _isMC, FitCharge);
  _tmpCharge = fitPars.charge;
 
  ComputePatterns(trackPtr);
 
  // Rigidity
  for (auto fit : fitTypes) {
    if (type == TrTrackFillType::Standalone) {
      switch (fit) {
      case Fit::KalmanElectron:
      case Fit::ChoutkoNoMS:
      case Fit::GBLNoMS:
        // for standalone we only consider Choutko, Kalman and GBL
        continue;
      default:
        break;
      }
    }
 
    // Kalman with electron mass hyp. to be done in dedicated fill method
    if (fit == Fit::KalmanElectron)
      continue;
 
    for (auto span : spanTypes) {
      // For NoMS fits we only care about InnerOnly
      if ((fit == Fit::ChoutkoNoMS || fit == Fit::GBLNoMS) && span != Span::InnerOnly) {
        continue;
      }
 
      if (type == TrTrackFillType::Standalone) {
        switch (span) {
        case Span::LowerHalfInner:
        case Span::UpperHalfInner:
        case Span::InnerL9:
          // for standalone we only consider InnerOnly, InnerL1 and FullSpan
          continue;
        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);
        if (_fit_ID[fit][span] > 0) {
          double rigidity{0};
          trackPtr->Interpolate(pl_TrCenter, 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]));
        }
        Logging::UnmuteOutput();
      }
    }
  }
  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;
  }
 
  // Charges
  auto stdCharge = trackPtr->GetQ_all(_beta, base_fit_ID);
  auto stdChargeInner = trackPtr->GetInnerQ_all(_beta, base_fit_ID);
 
  auto hlCharge = trackPtr->GetQH_all(2, _beta, base_fit_ID);
  auto hlChargeInner = trackPtr->GetInnerQH_all(2, _beta, base_fit_ID);
 
  auto yj_opt = TrChargeYJ::DefaultCorOpt | TrChargeYJ::kPow;
 
  // FIXME: As soon as Yi Jia fixes charge reconstruction for Z=1 and Z=2, we can remove this option
  if (_tmpCharge < 2.5) {
    yj_opt |= TrChargeYJ::kSingle;
  }
 
  auto yjCharge = trackPtr->GetQYJ_all(2, _beta, base_fit_ID);
 
  //  - Inner tracker
  Charge[ChargeRecoType::STD] = stdCharge.Mean;
  Charge[ChargeRecoType::HL] = hlCharge.Mean;
  // Charge[ChargeRecoType::YJ] does not support "global" charge;
  ChargeRMS[ChargeRecoType::STD] = stdCharge.RMS;
  ChargeRMS[ChargeRecoType::HL] = hlCharge.RMS;
  // ChargeRMS[ChargeRecoType::YJ] does not support "global" charge;
  InnerCharge[ChargeRecoType::STD] = stdChargeInner.Mean;
  InnerCharge[ChargeRecoType::HL] = hlChargeInner.Mean;
  InnerCharge[ChargeRecoType::YJ] = yjCharge.InnerQ;
  InnerChargeRMS[ChargeRecoType::STD] = stdChargeInner.RMS;
  InnerChargeRMS[ChargeRecoType::HL] = hlChargeInner.RMS;
  InnerChargeRMS[ChargeRecoType::YJ] = yjCharge.InnerQRMS;
 
  // geometry and position
  for (int ijl = 0; ijl < 9; ijl++) {
 
    if (!trackPtr->GetHitLJ(ijl + 1)) {
      continue;
    }
 
    if (trackPtr->GetHitLJ(ijl + 1)->OnlyY()) {
      m_trackPattern |= (static_cast<uint32_t>(HitClusterAssociation::Y) << (2 * ijl));
    } else {
      // Only-X case does not exist, right?
      m_trackPattern |= (static_cast<uint32_t>(HitClusterAssociation::XY) << (2 * ijl));
    }
 
    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) {
      if (type == TrTrackFillType::Standalone) {
        switch (fit) {
        case Fit::KalmanElectron:
        case Fit::ChoutkoNoMS:
        case Fit::GBLNoMS:
          // for standalone we only consider Choutko, Kalman and GBL
          continue;
        default:
          break;
        }
      }
 
      for (auto span : spanTypes) {
        // For NoMS fits we only care about InnerOnly
        if ((fit == Fit::ChoutkoNoMS || fit == Fit::GBLNoMS) && span != Span::InnerOnly) {
          continue;
        }
 
        if (type == TrTrackFillType::Standalone) {
          switch (span) {
          case Span::LowerHalfInner:
          case Span::UpperHalfInner:
          case Span::InnerL9:
            // 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
          trackPtr->Interpolate(pl_TrPlanes[fitPos], rigidity, fit_ID);
          TVector3 pglob_fitpos = pl_TrPlanes[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;
 
    Logging::MuteOutput();
    for (auto span : spanTypes) {
      // For NoMS fits we only care about InnerOnly
      if ((fit == Fit::ChoutkoNoMS || fit == Fit::GBLNoMS) && span != Span::InnerOnly) {
        continue;
      }
 
      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;
 
        AMSPlane &pl_TOI = _beta < 0 ? pl_TOI_bot : pl_TOI_top;
        double rigidity;
        if (fit == Fit::Kalman) {
          // TOI rigidity for both downward and upward going particles.
          trackPtr->Interpolate(pl_TOI, rigidity, _fit_ID[fit][span], fitMass, _tmpCharge, 1, (_beta > 0) ? 1 : -1);
          // NOTE: To keep the same convention for all rigidities measured in AMS we invert the result for upgoing
          // particles. We need to do this because Kalman keeps track of the propagation direction so an upgoing proton
          // is reconstructed with a positive TOI rigidity while all the "other" rigidities are reconstructed negative.
          RigidityTOI[fit][span] = static_cast<float>(rigidity * (_beta > 0 ? 1.0f : -1.0f));
          // Rigidity at RICH
          trackPtr->Interpolate(pl_RICH, rigidity, _fit_ID[fit][span], fitMass, _tmpCharge, 1);
          RigidityRICH[fit][span] = static_cast<float>(rigidity);
        }
        const TrTrackPar &trackPar = trackPtr->gTrTrackPar(_fit_ID[fit][span]);
        InvRigErr[fit][span] = static_cast<float>(trackPar.ErrRinv);
 
        // Theta and Phi
        trackPtr->Interpolate(pl_TOI, rigidity, _fit_ID[fit][span], fitMass, _tmpCharge, 1, (_beta > 0) ? 1 : -1);
        Theta[fit][span] = static_cast<float>(pl_TOI.getD().Theta());
        Phi[fit][span] = static_cast<float>(pl_TOI.getD().Phi());
      }
    }
    Logging::UnmuteOutput();
 
    // using default fit and inner tracker charge to restore coordinate
    trackPtr->RecalcHitCoordinates(0, false, false, _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);
 
        Logging::MuteOutput();
        int pfit_ID = trackPtr->iTrTrackPar(fitAlgoCode[fit], partialPattern, 30 + refit, fitMass, _tmpCharge);
        Logging::UnmuteOutput();
        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;
        trackPtr->Interpolate(pl_TrPlanes[ijl - 1], rigidity, pfit_ID);
        TVector3 pglob_fitpos = pl_TrPlanes[ijl - 1].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)) {
    const 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;
 
        if (!trackPtr->GetHitLJ(ijl + 1)->OnlyY()) {
          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
          trackPtr->Interpolate(pl_TrPlanes[fitPos], rigidity, fit_ID);
          TVector3 pglob_fitpos = pl_TrPlanes[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();
            if (!trackPtr->GetHitLJ(fitPos + 1)->OnlyY()) {
              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;
    trackPtr->Interpolate(pl_TrPlanes[jLayer - 1], rigidity, _fit_ID[Fit::Kalman][Span::InnerOnly]);
    TVector3 pglob = pl_TrPlanes[jLayer - 1].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;
      trackPtr->Interpolate(pl_TrCenter, 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) {
      AMSPlane &pl_TOI = _beta < 0 ? pl_TOI_bot : pl_TOI_top;
      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;
      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
      Theta[fit][span] = static_cast<float>(pl_TOI.getD().Theta());
      Phi[fit][span] = static_cast<float>(pl_TOI.getD().Phi());
 
      const TrTrackPar &trackPar = trackPtr->gTrTrackPar(fit_ID);
      InvRigErr[fit][span] = static_cast<float>(trackPar.ErrRinv);
    }
  }
 
  return true;
}
} // namespace NAIA