MCTruthFill.cpp

Go to the documentation of this file.

#include "Containers/MCTruth.h"

#include "TDatabasePDG.h"

#include <algorithm>

namespace NAIA {

static float amu = 0.93146;

MCParticle GetMCParticle(const MCEventgR &mcp) {
  MCParticle particle{TDatabasePDG::Instance()->ConvertGeant3ToPdg(mcp.Particle),
                      static_cast<int>(mcp.Charge),
                      TMath::Nint(mcp.Mass / amu),
                      TVector3{mcp.Coo[0], mcp.Coo[1], mcp.Coo[2]},
                      {mcp.Momentum * TVector3{mcp.Dir[0], mcp.Dir[1], mcp.Dir[2]}},
              mcp.trkID,
              mcp.parentID,
                      mcp.Nskip >> 24,
                      mcp.Nskip & 0xFFFFFF};

  return particle;
}

unsigned int TkIDToLayerJ(int tkID) {
  tkID = abs(tkID);

  unsigned int layerOld = tkID / 100;
  switch (layerOld) {
  case 8:
    return 1;
    break;
  case 9:
    return 9;
  default:
    return layerOld + 1;
  }
}

MCTruthPlusData::TrMCHit GetTrMCHit(TrTrackR *pTrTrack, unsigned int jLayer, double match_radius = 0.1) {
  AMSEventR *event = AMSEventR::Head();
  if (!event)
    throw std::runtime_error("Invalid pointer for AMSEventR::Head");
  if (event->nMCEventgC() == 0)
    throw std::runtime_error("No MCEventg information in this event");

  MCTruthPlusData::TrMCHit mcHit;

  TrRecHitR *hit = pTrTrack->GetHitLJ(jLayer);
  if (!hit)
    return mcHit;

  AMSPoint point = hit->GetCoord();

  auto mcclusterIt =
      std::min_element(begin(event->TrMCCluster()), end(event->TrMCCluster()),
                       [&point](auto &cl1, auto &cl2) { return cl1.GetXgl().dist(point) < cl2.GetXgl().dist(point); });

  if (mcclusterIt == end(event->TrMCCluster()))
    return mcHit;

  // Fail if minimal distance is not in match_radius
  auto rmin = mcclusterIt->GetXgl().dist(point);
  if (rmin > match_radius)
    return mcHit;

  mcHit.pID = mcclusterIt->GetPart();
  for (unsigned int ic = 0; ic < 3; ic++) {
    mcHit.Coo[ic] = mcclusterIt->GetXgl()[ic];
    mcHit.Mom[ic] = mcclusterIt->GetMom()[ic];
  }

  return mcHit;
}

auto particleSorter = [](const MCEventgR &a, const MCEventgR &b) -> bool {
  if (a.parentID < b.parentID)
    return true;
  else if (a.parentID > b.parentID)
    return false;
  else {
    if (a.Coo[2] > b.Coo[2])
      return true;
    else if (a.Coo[2] < b.Coo[2])
      return false;
    else
      return a.Momentum > b.Momentum;
  }
};

void MCTruthBaseData::Fill(AMSEventR *evPtr) {

  // sort all MC particles by interaction point/momentum
  std::sort(begin(evPtr->MCEventg()), end(evPtr->MCEventg()), particleSorter);

  auto primIt = std::find_if(begin(evPtr->MCEventg()), end(evPtr->MCEventg()),
                             [](const MCEventgR &mcp) { return mcp.Coo[2] > 160.0; });

  if (primIt == end(evPtr->MCEventg())) {
    throw std::runtime_error("No generation point found for primary particle (with CooZ > 160 cm)");
  }

  Primary = GetMCParticle(*primIt);
  for (const auto &mcp : evPtr->MCEventg()) {
    if (mcp.parentID != 0)
      break;

    auto hIt = std::find_if(std::next(begin(MCTruth::MCHeightsZ)), end(MCTruth::MCHeightsZ),
                            [&mcp](float height) { return fabs(mcp.Coo[2] - height) < 3.5e-1; });
    if (hIt != end(MCTruth::MCHeightsZ)) {
      auto position = std::distance(begin(MCTruth::MCHeightsZ), hIt);
      Primary.Momentum.resize(position + 1);
      Primary.Momentum[position] = mcp.Momentum * TVector3{mcp.Dir[0], mcp.Dir[1], mcp.Dir[2]};
    }
  }
}

void MCTruthPlusData::Fill(AMSEventR *evPtr) {

  // this should already be sorted, but you never know...
  if (!std::is_sorted(begin(evPtr->MCEventg()), end(evPtr->MCEventg()), particleSorter)) {
    // sort all MC particles by interaction point/momentum
    std::sort(begin(evPtr->MCEventg()), end(evPtr->MCEventg()), particleSorter);
  }

  const MCEventgR &firstPrim = evPtr->MCEventg()[0];

  const auto lastPrimIt = std::find_if(rbegin(evPtr->MCEventg()), rend(evPtr->MCEventg()),
                                       [](const MCEventgR &mcp) { return mcp.parentID == 0; });

  if (lastPrimIt == rend(evPtr->MCEventg()))
    return;

  for (MCEventgR &mcp : evPtr->MCEventg()) {
    if (mcp.Momentum / firstPrim.Momentum > 0.01 && mcp.Coo[2] <= lastPrimIt->Coo[2]) {
      Secondaries.push_back(GetMCParticle(mcp));
    }
  }
}

void MCTruthPlusData::FillTrMCHit(TrTrackR *pTrTrack) {
  for (unsigned int ij = 0; ij < 9; ij++) {
    auto mcHit = GetTrMCHit(pTrTrack, ij + 1);
    if (mcHit.pID > 0)
      TrackMCHits[ij] = mcHit;
  }
}
} // namespace NAIA