TrdKFill.cpp

Go to the documentation of this file.

#include "Containers/TrdK.h"

#include <iostream>

namespace NAIA {
using namespace TrdK;

static constexpr float threshold = 15;
static constexpr float TRDCenter = 115;

void TrdKBaseData::ResizeMembers() {
  NHits.resize(TrdK::numQualTypes);
  Amps.resize(TrdK::numQualTypes);
  Charge.resize(TrdK::numChargeTypes);
}

void TrdKBaseData::Fill(TrdKCluster *trdkPtr, TrTrackR *pTrTrack) {
  AlignStatus = trdkPtr->IsReadAlignmentOK;
  CalibStatus = trdkPtr->IsReadCalibOK;
  if (!AlignStatus || !CalibStatus)
    return;

  if (trdkPtr->NHits() < 1)
    return;

  double _tmpLR[3], _tmpL[3];
  float trd_pathl;

  int isValid;
  int nhits;
  float amps;
  if (!_isStandalone) {
    isValid = trdkPtr->GetLikelihoodRatio_TrTrack(threshold, _tmpLR, _tmpL, nhits, trd_pathl, amps);
  } else {
    isValid = trdkPtr->GetLikelihoodRatio_TRDRefit(threshold, _tmpLR, _tmpL, nhits, trd_pathl, amps);
  }

  ResizeMembers();

  NHits[QualType::OnTrack] = nhits;
  Amps[QualType::OnTrack] = amps;

  AMSPoint trd_point;
  AMSDir trd_dir;
  if (isValid) {
    Likelihood.resize(TrdK::numLikelihoodTypes-1);
    for (auto ltype : likelihoodTypes) {
      // we deal with deuterium later...
      if (ltype == LikelihoodType::Deuterium)
        continue;

      Likelihood[ltype] = _tmpL[ltype];
    }

    LikelihoodRatio.resize(TrdK::numLikelihoodRTypes);
    for (auto ltype : likelihoodRTypes) {
      LikelihoodRatio[ltype] = _tmpLR[ltype];
    }
  }

  ChargeStatus = trdkPtr->CalculateTRDCharge(0, _beta);
  Charge[ChargeType::Total] = trdkPtr->GetTRDCharge();
  Charge[ChargeType::Upper] = trdkPtr->GetTRDChargeUpper();
  Charge[ChargeType::Lower] = trdkPtr->GetTRDChargeLower();

  trdkPtr->GetOffTrackHit_TrTrack(NHits[QualType::OffTrack], Amps[QualType::OffTrack]);

  trdkPtr->FitTRDTrack();
  trdkPtr->GetTRDRefittedTrack(trd_point, trd_dir);

  if (!_isStandalone) {
    int status = trdkPtr->GetTrTrackExtrapolation(trd_point, trd_dir);
    if (status >= 0) {
      for (unsigned int i = 0; i < trdkPtr->NHits(); i++) {
        TrdKHit *hit = trdkPtr->GetHit(i);
        if (!hit)
          continue;
        float tube_length = hit->Tube_Track_3DLength(&trd_point, &trd_dir);
        if (tube_length <= 0.02)
          continue;
        float amplitude = hit->TRDHit_Amp;
        if (amplitude <= 0)
          continue;
        if (amplitude > 32767)
          amplitude = 32767;
        unsigned int ilay = hit->TRDHit_Layer;
        if (amplitude <= Edep[ilay])
          continue; // MaxADC
        Edep[ilay] = amplitude;
        Pathlength[ilay] = tube_length;
      }
    }
  }

  Cleanliness = trdkPtr->GetCleanliness();
  IPChi2 = trdkPtr->GetIPChi2();

  if (isValid && !_isStandalone && pTrTrack) {
    // here do deuterium...
    auto fit_ID_maxspan =
        std::max(trkBase->_fit_ID.at(TrTrack::Fit::Choutko).at(trkBase->GetBestSpan(TrTrack::Fit::Choutko)), 0);
    // trick for getting the deuteron likelihood
    pTrTrack->Interpolate(TRDCenter, trd_point, trd_dir, fit_ID_maxspan);
    float rigidity = pTrTrack->GetRigidity(fit_ID_maxspan);
    trdkPtr->SetTrTrack(&trd_point, &trd_dir, 2 * rigidity);

    // this is the reason why YOU DON'T USE REFERENCES OR POINTERS FOR FUNCTION RETURNS!
    int hitDummy;
    float pathlDummy, ampDummy;
    trdkPtr->GetLikelihoodRatio_TrTrack(threshold, _tmpLR, _tmpL, hitDummy, pathlDummy, ampDummy, -1, 0);

    if (Likelihood.size() == TrdK::numLikelihoodTypes - 1) {
      Likelihood.resize(TrdK::numLikelihoodTypes);
      Likelihood[LikelihoodType::Deuterium] = _tmpL[1];
    }
  }

  m_pointAndDir.emplace_back(trd_point.x(), trd_point.y(), trd_point.z());
  m_pointAndDir.emplace_back(trd_dir.x(), trd_dir.y(), trd_dir.z());
}

} // namespace NAIA