TofFill.cpp

Go to the documentation of this file.

#include "Containers/Tof.h"

namespace NAIA {
using namespace Tof;

bool TofBaseData::Fill(BetaR *betaPtr, BetaHR *betahPtr, bool standalone) {
  if (!betahPtr)
    return false;

  if (betaPtr && !standalone) {
    Beta[BetaType::Evgeni] = betaPtr->Beta;
  }

  auto betah = _isMC ? betahPtr->GetMCBeta() : betahPtr->GetBeta();
  Beta[BetaType::BetaH] = betah;

  int nlayers = 0;
  float qrms = 0;
  if (standalone || betahPtr->IsGoodQPathL(0) || betahPtr->IsGoodQPathL(1))
    Charge[ChargeType::Upper] =
        betahPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon,
                       1000 * betahPtr->IsGoodQPathL(0) + 100 * betahPtr->IsGoodQPathL(1), betah);
  if (standalone || betahPtr->IsGoodQPathL(2) || betahPtr->IsGoodQPathL(3))
    Charge[ChargeType::Lower] = betahPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon,
                                               10 * betahPtr->IsGoodQPathL(2) + 1 * betahPtr->IsGoodQPathL(3), betah);
  if (((betahPtr->IsGoodQPathL(0) || betahPtr->IsGoodQPathL(1)) &&
       (betahPtr->IsGoodQPathL(2) || betahPtr->IsGoodQPathL(3))) ||
      standalone)
    Charge[ChargeType::Total] = betahPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon,
                                               1000 * betahPtr->IsGoodQPathL(0) + 100 * betahPtr->IsGoodQPathL(1) +
                                                   10 * betahPtr->IsGoodQPathL(2) + 1 * betahPtr->IsGoodQPathL(3),
                                               betah);

  ChargeNoPL[ChargeType::Upper] = betahPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon, 1100, betah);
  ChargeNoPL[ChargeType::Lower] = betahPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon, 11, betah);
  ChargeNoPL[ChargeType::Total] = betahPtr->GetQ(nlayers, qrms, 2, TofClusterHR::DefaultQOptIon, 1111, betah);

  double tofh_time;
  AMSPoint tofh_point;
  AMSDir tofh_dir;
  betahPtr->TInterpolate(tofRefZ, tofh_point, tofh_dir, tofh_time, false);
  m_pointAndDir.emplace_back(tofh_point.x(), tofh_point.y(), tofh_point.z());
  m_pointAndDir.emplace_back(tofh_dir.x(), tofh_dir.y(), tofh_dir.z());

  return true;
}

void TofPlusData::ResizeMembers() { LayerGoodPathl.resize(4); }

bool TofPlusData::Fill(BetaR *betaPtr, BetaHR *betahPtr, bool standalone) {
  auto *evPtr = AMSEventR::Head();

  if (!betahPtr)
    return false;

  if (betaPtr && !standalone) {
    BetaPattern[BetaType::Evgeni] = static_cast<short>(betaPtr->Pattern);
  }

  short tofhbetapatt = 0;
  for (int i = 0; i < 4; i++) {
    if (!(betahPtr->TestExistHL(i))) {
      tofhbetapatt = static_cast<short>(tofhbetapatt | (1 << i));
    }
  }
  BetaPattern[BetaType::BetaH] = tofhbetapatt;

  Chi2Coo = betahPtr->GetNormChi2C();
  Chi2Time = betahPtr->GetNormChi2T();

  NTrkClusters = static_cast<short>(betahPtr->GetSumHit());
  NBetaClusters = static_cast<short>(betahPtr->GetUseHit());

  // this is fucking bad. 4 lines to get one variable. (FOR FUCK'S SAKE STOP USING REFERENCES FOR OUTPUT!!!)
  int nhit;
  float qrms;
  betahPtr->GetQ(nhit, qrms);
  NChargeClusters = static_cast<short>(nhit);

  // same here... gbatch has the mother of all bad APIs
  ChargeLikelihoodEst = betahPtr->GetZ(nhit, ChargeLikelihoodProb);

  ResizeMembers();
  for (int ip = 0; ip < 4; ip++) {
    LayerGoodPathl[ip] = betahPtr->IsGoodQPathL(ip);

    float layer_charge = betahPtr->GetQL(ip);
    if (layer_charge > 0.0f) {
      LayerCharge[ip] = layer_charge;
    }
    if (betahPtr->GetClusterHL(ip)) {
      LayerEdep[ip] = betahPtr->GetClusterHL(ip)->GetEdep();
    }

    // taken from dbar code
    if (betahPtr->TestExistHL(ip)) {
      float ltime = betahPtr->GetTime(ip); // ns
      for (auto &cluster : evPtr->TofClusterH()) {
        if (cluster.Layer != ip || (cluster.NBetaHUsed() > 0) || !cluster.IsGoodTime())
          continue;

        float dt = cluster.Time - ltime; // later cluster has positive dt
        float tcut = ip < 2 ? 10 : 4;    // 10ns for top layers, 4 ns for bottom layers

        int itm = -1;
        if (std::fabs(dt) <= tcut)
          itm = 0; // around-time hits ("in time")
        else if (dt > tcut)
          itm = 1; // later hits ("off time")

        switch (itm) {
        case 0:
          NClusters[ip][BetaClusterType::OnTime]++;
          AverageDTime[ip][BetaClusterType::OnTime] += dt;
          AverageEdep[ip][BetaClusterType::OnTime] += cluster.GetEdep();
          break;
        case 1:
          NClusters[ip][BetaClusterType::OffTime]++;
          AverageDTime[ip][BetaClusterType::OffTime] += dt;
          AverageEdep[ip][BetaClusterType::OffTime] += cluster.GetEdep();
          break;
        default:
          break;
        }
      }
    }
  }
  return true;
}

} // namespace NAIA