MCParticle.hpp

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <array>
#include <cstdint>
#include <vector>
 
namespace SHiP {
 
struct MCParticle {
  std::int32_t pdg_code{0};
  std::array<double, 3> vertex{0, 0, 0};    
  std::array<double, 3> momentum{0, 0, 0};  
  double energy{0};                         
  double time{0};                           
  std::int32_t mother_id{-1};               
  std::vector<std::int32_t> mothers{};      
  std::int32_t status{1};                   
 
  bool operator==(MCParticle const&) const = default;
};
 
inline bool mothers_are_consistent(std::vector<MCParticle> const& particles) {
  auto const size = static_cast<std::int32_t>(particles.size());
  for (std::int32_t i = 0; i < size; ++i) {
    auto const& mothers = particles[i].mothers;
    if (particles[i].mother_id < -1 || particles[i].mother_id >= size) {
      return false;
    }
    if (mothers.empty()) {
      continue;
    }
    if (mothers.front() != particles[i].mother_id) {
      return false;
    }
    for (auto m = mothers.begin(); m != mothers.end(); ++m) {
      if (*m < 0 || *m >= size || *m == i) {
        return false;
      }
      // Mother lists hold a handful of entries at most, so scanning the
      // prefix already seen beats allocating a set.
      if (std::find(mothers.begin(), m, *m) != m) {
        return false;
      }
    }
  }
  return true;
}
 
inline bool mothers_are_populated(std::vector<MCParticle> const& particles) {
  for (auto const& p : particles) {
    if (p.mother_id >= 0 && p.mothers.empty()) {
      return false;
    }
  }
  return mothers_are_consistent(particles);
}
 
}  // namespace SHiP