MCTruthFill.cpp

Go to the documentation of this file.

#include "Containers/MCTruth.h"
 
#include "TDatabasePDG.h"
 
#include <algorithm>
#include <optional>
 
namespace NAIA {
 
constexpr float amu = 0.93146;
 
MCParticle GetMCParticle(const MCEventgR &mcp) {
  // FIXME: It looks like G3 vs G4 is only determined by the absolute value of PartInfo.
  //        The sign is mainly used to distinguish between primary (positive) and secondary (negative) particles
  //        Except for the primary generation point. In that case PartInfo has a large value (and it could also be
  //        negative... KUDOS TO WHOEVER MADE THIS BRILLIANT AND IDIOTIC DECISION)
  int pdg{0};
  switch (std::abs(mcp.PartInfo)) {
  case 1:
    pdg = TDatabasePDG::Instance()->ConvertGeant3ToPdg(mcp.Particle);
    break;
  case 2:
    pdg = mcp.Particle;
    break;
  default:
    // FIXME: Check if the primary is always G3-coded
    pdg = TDatabasePDG::Instance()->ConvertGeant3ToPdg(mcp.Particle);
    break;
  }
 
  return {pdg,
          static_cast<int>(mcp.Charge),
          static_cast<unsigned int>(TMath::Nint(mcp.Mass / amu)),
          mcp.Mass,
          {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};
}
 
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;
  }
}
 
std::optional<MCTruthPlusData::TrMCHit> GetTrMCHit(TrTrackR *pTrTrack, unsigned int jLayer, float tof_beta,
                                                   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");
 
  if (!pTrTrack->TestHitLayerJ(jLayer))
    return {};
 
  AMSPlaneM plm = pTrTrack->GetHitCooLJN(jLayer);
  TVector3 mglob = plm.getMGlobal();
  AMSPoint point = AMSPoint(mglob.x(), mglob.y(), mglob.z());
 
  // Get all MCClusters within the match radius
  std::vector<TrMCClusterR> mcHits{};
  std::copy_if(begin(event->TrMCCluster()), end(event->TrMCCluster()), std::back_inserter(mcHits),
               [&point, &match_radius](TrMCClusterR &hit) { return hit.GetXgl().dist(point) < match_radius; });
 
  // Fail if no hits within match_radius
  if (mcHits.empty())
    return {};
 
  decltype(mcHits)::iterator mcclusterIt;
 
  if (pTrTrack->GetInnerQ(tof_beta) > 1.5) {
    // Get the cluster with the maximum edep within match_radius
    mcclusterIt = std::max_element(begin(mcHits), end(mcHits),
                                   [](TrMCClusterR &cl1, TrMCClusterR &cl2) { return cl1.GetEdep() < cl2.GetEdep(); });
  } else {
    // Get the cluster closest to the track
    mcclusterIt = std::min_element(begin(mcHits), end(mcHits), [&point](TrMCClusterR &cl1, TrMCClusterR &cl2) {
      return cl1.GetXgl().dist(point) < cl2.GetXgl().dist(point);
    });
  }
 
  MCTruthPlusData::TrMCHit mcHit;
  mcHit.pID = TDatabasePDG::Instance()->ConvertGeant3ToPdg(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;
  }
};
 
bool 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)");
  }
 
  constexpr float z_tolerance = 2.0f;
 
  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) < z_tolerance; });
    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]};
 
      Primary.Position.resize(position + 1);
      Primary.Position[position] = TVector3{mcp.Coo[0], mcp.Coo[1], mcp.Coo[2]};
    }
  }
 
  return true;
}
 
bool 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];
  for (MCEventgR &mcp : evPtr->MCEventg()) {
    // don't store the primary
    if (mcp.parentID == 0)
      continue;
 
    if (mcp.Momentum / firstPrim.Momentum > 0.005) {
      Secondaries.push_back(GetMCParticle(mcp));
    }
  }
 
  return true;
}
 
bool MCTruthPlusData::FillTrMCHit(TrTrackR *pTrTrack, float tof_beta) {
  for (unsigned int ij = 0; ij < 9; ij++) {
    if (auto mcHit = GetTrMCHit(pTrTrack, ij + 1, tof_beta); mcHit)
      TrackMCHits[ij] = mcHit.value();
  }
  return true;
}
} // namespace NAIA