/*//////////////////////////////////////////////////////////////////////////////
// CMT Cosmic Muon Tomography project //////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

  Copyright (c) 2014, Universita' degli Studi di Padova, INFN sez. di Padova
  All rights reserved

  Authors: Andrea Rigoni Garola < andrea.rigoni@pd.infn.it >

  ------------------------------------------------------------------
  This library is free software;  you  can  redistribute  it  and/or
  modify it  under the  terms  of  the  GNU  Lesser  General  Public
  License as published  by  the  Free  Software  Foundation;  either
  version 3.0 of the License, or (at your option) any later version.

  This library is  distributed in  the hope that it will  be useful,
  but  WITHOUT ANY WARRANTY;  without  even  the implied warranty of
  MERCHANTABILITY  or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.

  You should have received a copy of  the GNU Lesser General  Public
  License along with this library.

//////////////////////////////////////////////////////////////////////////////*/

#include "HEP/Detectors/DetectorChamber.h"
#include "testing-prototype.h"
#include <iostream>

using namespace uLib;

int main() {
    BEGIN_TESTING(DetectorChamber Projection);

    DetectorChamber chamber;

    // Define a horizontal plane at z = 100
    HLine3f plane;
    plane.origin = HPoint3f(0, 0, 100);
    plane.direction = HVector3f(0, 0, 1); // Normal to the plane
    chamber.SetProjectionPlane(plane);

    // Create a muon with two segments
    MuonEvent muon;
    
    // Segment 1 (muon.LineIn()): origin is at (10, 20, 50), direction along Z
    muon.LineIn().origin = HPoint3f(10, 20, 50);
    muon.LineIn().direction = HVector3f(0, 0, 1);
    
    // Segment 2 (muon.LineOut()): origin is at (10, 20, 200), direction along Z
    muon.LineOut().origin = HPoint3f(10, 20, 200);
    muon.LineOut().direction = HVector3f(0, 0, 1);
    
    // distance_in = |50 - 100| = 50
    // distance_out = |200 - 100| = 100
    // LineIn is closer to the plane (z=100)
    
    MuonEvent projected = chamber.ProjectMuonEvent(muon);

    // Expected:
    // chosenLine = LineIn (it is closer)
    // X_chosen = (10, 20, 50)
    // X_proj = (10, 20, 50) - (( (10, 20, 50) - (0, 0, 100) ) . (0, 0, 1)) * (0, 0, 1)
    // X_proj = (10, 20, 50) - (-50) * (0, 0, 1) = (10, 20, 100)
    
    HPoint3f expectedPos(10, 20, 100);
    
    std::cout << "Test Case 1: LineIn is closer" << std::endl;
    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;

    // Check if the projected position is correct
    // norm() includes the 4th component, which for HVector3f (diff of points) should be 0.
    bool posOk1 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
    TEST1(posOk1);

    // Check if LineIn and LineOut are the same
    bool linesMatch1 = (projected.LineIn().origin == projected.LineOut().origin) &&
                       (projected.LineIn().direction == projected.LineOut().direction);
    TEST1(linesMatch1);
    
    // Test Case 2: LineOut is closer
    muon.LineIn().origin = HPoint3f(30, 40, 0);      // dist = 100
    muon.LineOut().origin = HPoint3f(30, 40, 110);   // dist = 10
    
    projected = chamber.ProjectMuonEvent(muon);
    expectedPos = HPoint3f(30, 40, 100); // projection of (30,40,110) onto z=100
    
    std::cout << "\nTest Case 2: LineOut is closer" << std::endl;
    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;

    bool posOk2 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
    TEST1(posOk2);
    
    // Test Case 3: Oblique plane
    // Plane through (0,0,0) with normal (1,1,1) (normalized)
    plane.origin = HPoint3f(0, 0, 0);
    plane.direction = HVector3f(1, 1, 1).normalized();
    chamber.SetProjectionPlane(plane);
    
    muon.LineIn().origin = HPoint3f(1, 1, 1); // dist = (1,1,1) . (1,1,1).norm()
    muon.LineIn().direction = HVector3f(0, 0, 1);
    muon.LineOut().origin = HPoint3f(10, 10, 10); // dist = 10 * sqrt(3)
    
    projected = chamber.ProjectMuonEvent(muon);
    // X_chosen = (1,1,1)
    // X_proj = (1,1,1) - ((1,1,1) . N) * N   where N = (1,1,1)/sqrt(3)
    // X_proj = (1,1,1) - (sqrt(3)) * (1,1,1)/sqrt(3) = (1,1,1) - (1,1,1) = (0,0,0)
    expectedPos = HPoint3f(0, 0, 0);
    
    std::cout << "\nTest Case 3: Oblique plane" << std::endl;
    std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
    std::cout << "Expected Position:  " << expectedPos.transpose() << std::endl;
    
    bool posOk3 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
    TEST1(posOk3);

    END_TESTING;
}