UnbExtHitFill.cpp

Go to the documentation of this file.

#include "Containers/UnbExtHit.h"

namespace NAIA {

bool UnbExtHitBaseData::Fill(AMSEventR *evPtr, bool use_trd) {

  auto fillHit = [this](ExtHit type, TrRecHitR *hit, int imult, float beta, const TVector3 &interp_direction) {
    float tof_dxdz = interp_direction.x() / interp_direction.z();
    float tof_dydz = interp_direction.y() / interp_direction.z();

    m_hits[type].ChargeStatus = hit->GetQStatus();
    m_hits[type].HitPosX = hit->GetCoord(imult).x();
    m_hits[type].HitPosY = hit->GetCoord(imult).y();
    m_hits[type].Charge[TrTrack::ChargeRecoType::STD] = hit->GetQ(2, beta, 0, 0, imult, tof_dxdz, tof_dydz);
    m_hits[type].Charge[TrTrack::ChargeRecoType::HL] = hit->GetQH(2, beta, 0, imult, tof_dxdz, tof_dydz);
    m_hits[type].Charge[TrTrack::ChargeRecoType::YJ] = hit->GetQYJ(2, beta, 0, imult, tof_dxdz, tof_dydz);
  };

  constexpr auto tof_beta_type = Tof::BetaType::BetaH;
  constexpr auto tof_charge_type = Tof::ChargeType::Upper;

  // need at least a good beta and charge measurement
  if (!ContainsKeys(tofBase->Beta, tof_beta_type) || !ContainsKeys(tofBase->Charge, tof_charge_type))
    return false;

  float tof_beta = tofBase->Beta.at(tof_beta_type);
  const auto &interp_direction = use_trd ? trdBase->GetLineDirection() : tofBase->GetLineDirection();

  // Search unbiased hits on external layers
  TrRecHitR *hitL1 = nullptr;
  TrRecHitR *hitL9 = nullptr;
  int multL1{0};
  int multL9{0};

  if (floor(tofBase->Charge.at(tof_charge_type) + 0.5) > 1) {
    // Case Helium or higher we use the highest charge hit
    float prev_charge1 = 0;
    float prev_charge9 = 0;
    for (TrRecHitR &hit : evPtr->TrRecHit()) {
      if (hit.OnlyY() || (hit.GetLayerJ() != 1 && hit.GetLayerJ() != 9))
        continue;

      // Solve multiplicity with tof standalone
      int imult{0};
      TVector3 vInterp =
          use_trd ? trdBase->InterpolateAtZ(hit.GetCoord().z()) : tofBase->InterpolateAtZ(hit.GetCoord().z());
      AMSPoint interp_point = AMSPoint(vInterp[0], vInterp[1], vInterp[2]);
      hit.HitPointDist(interp_point, imult); // again, how can I know imult is an output variable?

      float interp_dxdz = interp_direction.x() / interp_direction.z();
      float interp_dydz = interp_direction.y() / interp_direction.z();

      float hit_charge = hit.GetQYJ(2, tof_beta, 0, imult, interp_dxdz, interp_dydz);

      if (hit.GetLayerJ() == 1) {
        if (hit_charge > prev_charge1) {
          hitL1 = &hit;
          multL1 = imult;
          prev_charge1 = hit_charge;
        }
      } else if (hit.GetLayerJ() == 9) {
        if (hit_charge > prev_charge9) {
          hitL9 = &hit;
          multL9 = imult;
          prev_charge9 = hit_charge;
        }
      }
    }

    if (hitL1 != nullptr) {
      fillHit(ExtHit::L1, hitL1, multL1, tof_beta, interp_direction);
    }
    if (hitL9 != nullptr) {
      fillHit(ExtHit::L9, hitL9, multL9, tof_beta, interp_direction);
    }

  } else if (tofBase->Charge.at(tof_charge_type) > 0) {
    static constexpr auto l1Z = TrTrack::fitPositionHeightZ[TrTrack::FitPositionHeight::Layer1];
    static constexpr auto l9Z = TrTrack::fitPositionHeightZ[TrTrack::FitPositionHeight::Layer9];

    TVector3 vL1Int = use_trd ? trdBase->InterpolateAtZ(l1Z) : tofBase->InterpolateAtZ(l1Z);
    TVector3 vL9Int = use_trd ? trdBase->InterpolateAtZ(l9Z) : tofBase->InterpolateAtZ(l9Z);

    AMSPoint interp_point_L1 = AMSPoint(vL1Int[0], vL1Int[1], vL1Int[2]);
    AMSPoint interp_point_L9 = AMSPoint(vL9Int[0], vL9Int[1], vL9Int[2]);

    float prev_dist1 = std::numeric_limits<float>::max();
    float prev_dist9 = std::numeric_limits<float>::max();
    int mult1 = -1, mult9 = -1;

    for (TrRecHitR &hit : evPtr->TrRecHit()) {
      if (hit.OnlyY() || (hit.GetLayerJ() != 1 && hit.GetLayerJ() != 9))
        continue;

      if (hit.GetLayerJ() == 1) {
        if (hit.HitPointDist(interp_point_L1, mult1).norm() < prev_dist1) {
          hitL1 = &hit;
          prev_dist1 = hitL1->HitPointDist(interp_point_L1, mult1).norm();
        }
      } else if (hit.GetLayerJ() == 9) {
        if (hit.HitPointDist(interp_point_L9, mult9).norm() < prev_dist9) {
          hitL9 = &hit;
          prev_dist9 = hitL9->HitPointDist(interp_point_L9, mult9).norm();
        }
      }
    }

    if (hitL1 != nullptr)
      fillHit(ExtHit::L1, hitL1, mult1, tof_beta, interp_direction);
    if (hitL9 != nullptr)
      fillHit(ExtHit::L9, hitL9, mult9, tof_beta, interp_direction);
  }

  return true;
}
} // namespace NAIA