#include "HEP/Detectors/DetectorChamber.h"
#include <cmath>

namespace uLib {

MuonEvent DetectorChamber::ProjectMuonEvent(const MuonEvent &muon) const {
  MuonEvent projectedMuon = muon;

  // Transform the local projection plane to world coordinates
  HLine3f worldPlane = this->GetWorldProjectionPlane();
  HPoint3f P = worldPlane.origin;
  HVector3f N = worldPlane.direction;

  HPoint3f X_in = muon.LineIn().origin;
  HPoint3f X_out = muon.LineOut().origin;

  // Calculate squared distances to the plane normal point for comparison
  // Actually, we should probably follow the user's description literally:
  // "closest ... with the projection plane ( so the colsest direction point with the point of the normal defining the plane )"
  // This could mean point-to-plane or point-to-point. 
  // Given "closest with the projection plane", point-to-plane is more natural.
  // However, "closest direction point with the point of the normal" strongly suggests point-to-point distance.
  // Let's use distance to the plane for the first part and keep the logic consistent.
  
  float dist_in = std::abs((X_in - P).dot(N));
  float dist_out = std::abs((X_out - P).dot(N));

  const HLine3f &chosenLine = (dist_in <= dist_out) ? muon.LineIn() : muon.LineOut();
  HPoint3f X_chosen = chosenLine.origin;

  // Project X_chosen into the plane defined by P and normal N
  // X_proj = X_chosen - ((X_chosen - P) . N / (N . N)) * N
  float dot = (X_chosen - P).dot(N);
  float n_sq = N.dot(N);

  HPoint3f X_proj = X_chosen;
  if (n_sq > 0) {
    X_proj = X_chosen - (dot / n_sq) * N;
  }

  // Define the projected line with projected origin and original direction
  HLine3f projectedLine;
  projectedLine.origin = X_proj;
  projectedLine.direction = chosenLine.direction;

  // Set both input and output lines of the projected muon to the same projected line
  projectedMuon.LineIn() = projectedLine;
  projectedMuon.LineOut() = projectedLine;

  return projectedMuon;
}

} // namespace uLib