UnbExtHitFill.cpp

Go to the documentation of this file.

#include "Containers/UnbExtHit.h"

namespace NAIA {

void UnbExtHitBaseData::Fill(AMSEventR *evPtr) {

  auto fillHit = [this](ExtHit type, TrRecHitR *hit, int imult, float beta) {
    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);
    m_hits[type].Charge[TrTrack::ChargeRecoType::HL] = hit->GetQH(2, beta);
    m_hits[type].Charge[TrTrack::ChargeRecoType::YJ] = hit->GetQYJ(2, beta);
  };

  // need at least a good beta and charge measurement
  if (tofBase->Beta.find(Tof::BetaType::BetaH) == end(tofBase->Beta) ||
      tofBase->Charge.find(Tof::ChargeType::Total) == end(tofBase->Charge))
    return;

  float _beta = tofBase->Beta.at(Tof::BetaType::BetaH);

  // Search unbiased hits on external layers
  TrRecHitR *hitL1 = nullptr;
  TrRecHitR *hitL9 = nullptr;

  if (floor(tofBase->Charge.at(Tof::ChargeType::Total) + 0.5) > 1) {
    // Case Helium or higher we use the highest charge hit
    Float_t prev_charge1 = 0;
    Float_t prev_charge9 = 0;
    for (int ihit = 0; ihit < evPtr->NTrRecHit(); ihit++) {
      if (evPtr->pTrRecHit(ihit)->OnlyX() || evPtr->pTrRecHit(ihit)->OnlyY())
        continue;
      if (evPtr->pTrRecHit(ihit)->GetLayerJ() == 1) {
        if (evPtr->pTrRecHit(ihit)->GetXCluster()->GetQ() > prev_charge1) {
          hitL1 = evPtr->pTrRecHit(ihit);
          prev_charge1 = evPtr->pTrRecHit(ihit)->GetXCluster()->GetQ();
        }
      } else if (evPtr->pTrRecHit(ihit)->GetLayerJ() == 9) {
        if (evPtr->pTrRecHit(ihit)->GetXCluster()->GetQ() > prev_charge9) {
          hitL9 = evPtr->pTrRecHit(ihit);
          prev_charge9 = evPtr->pTrRecHit(ihit)->GetXCluster()->GetQ();
        }
      }
    }
    // Solve multiplicity with tof standalone
    Int_t mult1 = -1, mult9 = -1;
    if (hitL1 != nullptr) {
      TVector3 betaInt = tofBase->InterpolateAtZ(hitL1->GetCoord().z());
      AMSPoint tofh_point = AMSPoint(betaInt[0], betaInt[1], betaInt[2]);
      hitL1->HitPointDist(tofh_point, mult1);
      // cout << mult1 << " " << hitL1->GetResolvedMultiplicity() << endl;
      fillHit(ExtHit::L1, hitL1, mult1, _beta);
    }
    if (hitL9 != nullptr) {
      TVector3 betaInt = tofBase->InterpolateAtZ(hitL9->GetCoord().z());
      AMSPoint tofh_point = AMSPoint(betaInt[0], betaInt[1], betaInt[2]);
      hitL9->HitPointDist(tofh_point, mult9);
      // cout << mult9 << " " << hitL9->GetResolvedMultiplicity() << endl;
      fillHit(ExtHit::L9, hitL9, mult9, _beta);
    }

  } else if (tofBase->Charge.at(Tof::ChargeType::Total) > 0) {
    TVector3 betaL1Int = tofBase->InterpolateAtZ(158.920);
    TVector3 betaL9Int = tofBase->InterpolateAtZ(-135.882);

    AMSPoint tofh_point_L1 = AMSPoint(betaL1Int[0], betaL1Int[1], betaL1Int[2]);
    AMSPoint tofh_point_L9 = AMSPoint(betaL9Int[0], betaL9Int[1], betaL9Int[2]);

    Float_t prev_dist1 = 99999;
    Float_t prev_dist9 = 99999;
    Int_t mult1 = -1, mult9 = -1;
    for (int ihit = 0; ihit < evPtr->NTrRecHit(); ihit++) {
      if (evPtr->pTrRecHit(ihit)->OnlyX() || evPtr->pTrRecHit(ihit)->OnlyY())
        continue;
      if (evPtr->pTrRecHit(ihit)->GetLayerJ() == 1) {
        if (evPtr->pTrRecHit(ihit)->HitPointDist(tofh_point_L1, mult1).norm() < prev_dist1) {
          hitL1 = evPtr->pTrRecHit(ihit);
          prev_dist1 = hitL1->HitPointDist(tofh_point_L1, mult1).norm();
        }
      } else if (evPtr->pTrRecHit(ihit)->GetLayerJ() == 9) {
        if (evPtr->pTrRecHit(ihit)->HitPointDist(tofh_point_L9, mult9).norm() < prev_dist9) {
          hitL9 = evPtr->pTrRecHit(ihit);
          prev_dist9 = hitL9->HitPointDist(tofh_point_L9, mult9).norm();
        }
      }
    }

    if (hitL1 != nullptr)
      fillHit(ExtHit::L1, hitL1, mult1, _beta);
    if (hitL9 != nullptr)
      fillHit(ExtHit::L9, hitL9, mult9, _beta);
  }
}
} // namespace NAIA