140 lines
5.8 KiB
C++
140 lines
5.8 KiB
C++
/*//////////////////////////////////////////////////////////////////////////////
|
|
// 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);
|
|
|
|
// Test Case 4: Transformed DetectorChamber
|
|
DetectorChamber chamber2;
|
|
chamber2.SetPosition(Vector3f(0, 0, 100)); // Move chamber to z=100
|
|
// chamber2.GetProjectionPlane has default origin (0,0,0) and direction (0,0,1)
|
|
// In world coordinates, this plane is at z = 100 + 0 = 100.
|
|
|
|
muon.LineIn().origin = HPoint3f(50, 60, 50); // dist to world plane (z=100) is 50
|
|
muon.LineOut().origin = HPoint3f(50, 60, 200); // dist to world plane (z=100) is 100
|
|
|
|
projected = chamber2.ProjectMuonEvent(muon);
|
|
expectedPos = HPoint3f(50, 60, 100);
|
|
|
|
std::cout << "\nTest Case 4: Transformed DetectorChamber (active world matrix)" << std::endl;
|
|
std::cout << "Projected Position: " << projected.LineIn().origin.transpose() << std::endl;
|
|
std::cout << "Expected Position: " << expectedPos.transpose() << std::endl;
|
|
|
|
bool posOk4 = (projected.LineIn().origin - expectedPos).norm() < 1e-5;
|
|
TEST1(posOk4);
|
|
|
|
END_TESTING;
|
|
}
|